hydro.h 7.04 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
/*******************************************************************************
 * This file is part of SWIFT.
 * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

/**
 * @brief Computes the hydro time-step of a given particle
 *
 * @param p Pointer to the particle data
 * @param xp Pointer to the extended particle data
 *
 */
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
    struct part* p, struct xpart* xp) {

  /* CFL condition */
31
  const float dt_cfl = const_cfl * p->h / p->v_sig;
32

33
34
35
36
  /* /\* Limit change in h *\/ */
  /* float dt_h_change = (p->h_dt != 0.0f) */
  /*                         ? fabsf(const_ln_max_h_change * p->h / p->h_dt) */
  /*                         : FLT_MAX; */
37

38
39
40
41
  /* /\* Limit change in u *\/ */
  /* float dt_u_change = (p->force.u_dt != 0.0f) */
  /*                         ? fabsf(const_max_u_change * p->u / p->force.u_dt) */
  /*                         : FLT_MAX; */
42

43
44
  //  return fminf(dt_cfl, fminf(dt_h_change, dt_u_change));
  return dt_cfl; 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
}



/**
 * @brief Prepares a particle for the density calculation.
 *
 * Zeroes all the relevant arrays in preparation for the sums taking place in 
 * the variaous density tasks
 *
 * @param p The particle to act upon
 */
__attribute__((always_inline)) INLINE static void hydro_init_part(struct part* p) {
  p->density.wcount = 0.f;
  p->density.wcount_dh = 0.f;
  p->rho = 0.f;
  p->rho_dh = 0.f;
62
63
64
65
66
  p->div_v = 0.f;
  p->curl_v = 0.f;
  p->rot_v[0] = 0.f;
  p->rot_v[1] = 0.f;
  p->rot_v[2] = 0.f;
67
68
69
70
71
72
73
74
75
}

/**
 * @brief Finishes the density calculation. 
 *
 * Multiplies the density and number of neighbours by the appropiate constants 
 * and add the self-contribution term.
 *
 * @param p The particle to act upon
76
 * @param time The current time
77
 */
78
__attribute__((always_inline)) INLINE static void hydro_end_density(struct part* p, float time) {
79
80
81
82
83
84
85

  /* Some smoothing length multiples. */
  const float h = p->h;
  const float ih = 1.0f / h;
  const float ih2 = ih * ih;
  const float ih4 = ih2 * ih2;
  
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
  p->rho_dh -= 3.0f * p->mass * kernel_root * kernel_igamma;
  p->density.wcount += kernel_root;
  
  /* Finish the calculation by inserting the missing h-factors */
  p->rho *= ih * ih2;
  p->rho_dh  *= ih4;
  p->density.wcount *= (4.0f / 3.0 * M_PI * kernel_gamma3);
  p->density.wcount_dh *= ih * (4.0f / 3.0 * M_PI * kernel_gamma3);
  
  /* Compute the derivative term */
  p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh / p->rho);

  /* Finish calculation of the velocity curl */
101
102
103
104
105
  p->rot_v[0] *= ih4;
  p->rot_v[1] *= ih4;
  p->rot_v[2] *= ih4;

  /* And compute the norm of the curl */
106
107
108
109
110
  p->curl_v = sqrtf(p->rot_v[0] * p->rot_v[0] +
		    p->rot_v[1] * p->rot_v[1] +
		    p->rot_v[2] * p->rot_v[2]) / p->rho;
  
  /* Finish calculation of the velocity divergence */
111
  p->div_v *= ih4 / p->rho;
112
113
114
  
  /* Compute the pressure */
  const float dt = time - 0.5f*(p->t_begin + p->t_end);
115
  p->pressure = (p->entropy + p->entropy_dt * dt) * powf(p->rho, const_hydro_gamma);
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

}


/**
 * @brief Prepare a particle for the force calculation.
 *
 * Computes viscosity term, conduction term and smoothing length gradient terms.
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
 */
__attribute__((always_inline)) INLINE static void hydro_prepare_force(struct part* p,
								      struct xpart* xp) {

  /* Some smoothing length multiples. */
132
133
134
135
  /* const float h = p->h; */
  /* const float ih = 1.0f / h; */
  /* const float ih2 = ih * ih; */
  /* const float ih4 = ih2 * ih2; */
136
137

  
138
139
140
141
142
143
144
145
146
147
148
  /* /\* Pre-compute some stuff for the balsara switch. *\/ */
  /* const float normDiv_v = fabs(p->density.div_v / p->rho * ih4); */
  /* const float normCurl_v = sqrtf(p->density.curl_v[0] * p->density.curl_v[0] + */
  /* 				 p->density.curl_v[1] * p->density.curl_v[1] + */
  /* 				 p->density.curl_v[2] * p->density.curl_v[2]) / */
  /*   p->rho * ih4; */

  /* /\* Compute this particle's sound speed. *\/ */
  /* const float u = p->u; */
  /* const float fc = p->force.c =  */
  /*   sqrtf(const_hydro_gamma * (const_hydro_gamma - 1.0f) * u); */
149
  
150
151
152
  /* /\* Compute the P/Omega/rho2. *\/ */
  /* xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho; */
  /* p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (p->rho * xp->omega); */
153
  
154
155
156
  /* /\* Balsara switch *\/ */
  /* p->force.balsara = */
  /*   normDiv_v / (normDiv_v + normCurl_v + 0.0001f * fc * ih); */
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
}



/**
 * @brief Reset acceleration fields of a particle
 *
 * Resets all hydro acceleration and time derivative fields in preparation
 * for the sums taking place in the variaous force tasks
 *
 * @param p The particle to act upon
 */
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(struct part* p) {

  /* Reset the acceleration. */
  p->a[0] = 0.0f;
  p->a[1] = 0.0f;
  p->a[2] = 0.0f;
  
  /* Reset the time derivatives. */
177
178
179
180
  p->entropy_dt = 0.0f;
  
  /* Reset minimal signal velocity */
  p->v_sig = FLT_MAX;
181
182
183
184
185
186
187
188
}


/**
 * @brief Predict additional particle fields forward in time when drifting
 *
 * @param p The particle
 * @param xp The extended data of the particle
189
190
 * @param t0 The time at the start of the drift
 * @param t1 The time at the end of the drift
191
192
193
 */
__attribute__((always_inline)) INLINE static void hydro_predict_extra(struct part* p,
								      struct xpart* xp,
194
195
196
197
198
199
200
201
202
203
204
								      float t0,
								      float t1) {

  const float dt = t1 - t0;
  
  p->rho *= expf(-p->div_v * dt);
  p->h *= expf( 0.33333333f * p->div_v * dt);

  const float dt_entr = t1 - 0.5f*(p->t_begin + p->t_end);
  p->pressure = (p->entropy + p->entropy_dt * dt_entr) * powf(p->rho, const_hydro_gamma);
  
205
206
207
208
209
210
211
212
213
214
215
216
217
}



/**
 * @brief Finishes the force calculation. 
 *
 * Multiplies the forces and accelerationsby the appropiate constants 
 *
 * @param p The particle to act upon
 */
__attribute__((always_inline)) INLINE static void hydro_end_force(struct part* p) {

218
219
  p->entropy_dt *= (const_hydro_gamma - 1.f) * powf(p->rho, -(const_hydro_gamma - 1.f));
  
220
221
222
223
224
225
226
227
228
229
230
231
232
}

/**
 * @brief Converts hydro quantity of a particle 
 *
 * Requires the density to be known
 *
 * @param p The particle to act upon
 */
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(struct part* p) {

  p->entropy = (const_hydro_gamma - 1.f) * p->entropy * pow(p->rho, -(const_hydro_gamma - 1.f));

233
}