hydro.h 8.17 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/*******************************************************************************
 * 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/>.
 *
 ******************************************************************************/

20
#include "adiabatic_index.h"
21
22
#include "approx_math.h"

23
24
25
26
27
28
29
30
/**
 * @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(
31
32
    const struct part *restrict p, const struct xpart *restrict xp,
    const struct hydro_props *restrict hydro_properties) {
33
34

  const float CFL_condition = hydro_properties->CFL_condition;
35
36

  /* CFL condition */
37
38
  const float dt_cfl =
      2.f * kernel_gamma * CFL_condition * p->h / p->force.v_sig;
39
40

  /* Limit change in u */
Matthieu Schaller's avatar
Matthieu Schaller committed
41
42
43
  const float dt_u_change =
      (p->force.u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->force.u_dt)
                              : FLT_MAX;
44

45
  return fminf(dt_cfl, dt_u_change);
46
47
}

48
49
50
51
52
53
54
55
56
/**
 * @brief Initialises the particles for the first time
 *
 * This function is called only once just after the ICs have been
 * read in to do some conversions.
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
 */
57
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
58
    struct part *restrict p, struct xpart *restrict xp) {}
59

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

/**
81
 * @brief Finishes the density calculation.
82
 *
83
 * Multiplies the density and number of neighbours by the appropiate constants
84
85
86
 * and add the self-contribution term.
 *
 * @param p The particle to act upon
Matthieu Schaller's avatar
Matthieu Schaller committed
87
 * @param time The current time
88
 */
89
__attribute__((always_inline)) INLINE static void hydro_end_density(
90
    struct part *restrict p, float time) {
91
92
93
94
95
96

  /* 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;
97

98
99
  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
100
  p->rho_dh -= 3.0f * p->mass * kernel_root;
101
102
103
104
105
  p->density.wcount += kernel_root;

  /* Finish the calculation by inserting the missing h-factors */
  p->rho *= ih * ih2;
  p->rho_dh *= ih4;
106
107
  p->density.wcount *= kernel_norm;
  p->density.wcount_dh *= ih * kernel_gamma * kernel_norm;
108
109
110
111
112

  const float irho = 1.f / p->rho;

  /* Compute the derivative term */
  p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
113
114
115
116
117
118
119
120

  /* Finish calculation of the velocity curl components */
  p->density.rot_v[0] *= ih4 * irho;
  p->density.rot_v[1] *= ih4 * irho;
  p->density.rot_v[2] *= ih4 * irho;

  /* Finish calculation of the velocity divergence */
  p->density.div_v *= ih4 * irho;
121
122
123
124
125
126
127
128
129
}

/**
 * @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
130
 * @param time The current time
131
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
132
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
133
134
    struct part *restrict p, struct xpart *restrict xp, int ti_current,
    double timeBase) {
135
136
137
138
139
140

  /* Some smoothing length multiples. */
  const float h = p->h;
  const float ih = 1.0f / h;

  /* Pre-compute some stuff for the balsara switch. */
141
  const float normDiv_v = fabs(p->density.div_v);
142
  const float normRot_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
Matthieu Schaller's avatar
Matthieu Schaller committed
143
                                p->density.rot_v[1] * p->density.rot_v[1] +
144
                                p->density.rot_v[2] * p->density.rot_v[2]);
145
146
147

  /* Compute this particle's sound speed. */
  const float u = p->u;
Matthieu Schaller's avatar
Matthieu Schaller committed
148
149
  const float fc = p->force.soundspeed =
      sqrtf(hydro_gamma * hydro_gamma_minus_one * u);
150
151
152

  /* Compute the P/Omega/rho2. */
  xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho;
153
  p->force.P_over_rho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
154
155

  /* Balsara switch */
156
  p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * ih);
157
158

  /* Viscosity parameter decay time */
159
  const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
160
161
162
163
164
165

  /* Viscosity source term */
  const float S = fmaxf(-normDiv_v, 0.f);

  /* Compute the particle's viscosity parameter time derivative */
  const float alpha_dot = (const_viscosity_alpha_min - p->alpha) / tau +
166
                          (const_viscosity_alpha_max - p->alpha) * S;
167
168

  /* Update particle's viscosity paramter */
169
  p->alpha += alpha_dot * (p->ti_end - p->ti_begin) * timeBase;
170
171
172
173
174
175
176
177
178
179
}

/**
 * @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
 */
180
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
181
    struct part *restrict p) {
182
183

  /* Reset the acceleration. */
Matthieu Schaller's avatar
Matthieu Schaller committed
184
185
186
  p->a_hydro[0] = 0.0f;
  p->a_hydro[1] = 0.0f;
  p->a_hydro[2] = 0.0f;
187

188
189
  /* Reset the time derivatives. */
  p->force.u_dt = 0.0f;
190
  p->force.h_dt = 0.0f;
191
192
193
194
195
196
197
198
  p->force.v_sig = 0.0f;
}

/**
 * @brief Predict additional particle fields forward in time when drifting
 *
 * @param p The particle
 * @param xp The extended data of the particle
Matthieu Schaller's avatar
Matthieu Schaller committed
199
200
 * @param t0 The time at the start of the drift
 * @param t1 The time at the end of the drift
201
 * @param timeBase The minimal time-step size
202
 */
203
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
204
205
    struct part *restrict p, struct xpart *restrict xp, int t0, int t1,
    double timeBase) {
206
207
  float u, w;

208
  const float dt = (t1 - t0) * timeBase;
Matthieu Schaller's avatar
Matthieu Schaller committed
209

210
211
  /* Predict internal energy */
  w = p->force.u_dt / p->u * dt;
212
213
  if (fabsf(w) < 0.2f)
    u = p->u *= approx_expf(w);
214
215
  else
    u = p->u *= expf(w);
216

217
  /* Predict gradient term */
218
  p->force.P_over_rho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
219
220
221
}

/**
222
 * @brief Finishes the force calculation.
223
 *
224
 * Multiplies the forces and accelerationsby the appropiate constants
225
226
227
 *
 * @param p The particle to act upon
 */
228
__attribute__((always_inline)) INLINE static void hydro_end_force(
229
230
231
    struct part *restrict p) {
  p->force.h_dt *= p->h * 0.333333333f;
}
232
233
234
235
236

/**
 * @brief Kick the additional variables
 *
 * @param p The particle to act upon
237
238
239
 * @param xp The particle extended data to act upon
 * @param dt The time-step for this kick
 * @param half_dt The half time-step for this kick
240
 */
241
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
242
243
    struct part *restrict p, struct xpart *restrict xp, float dt,
    float half_dt) {}
244
245

/**
246
 *  @brief Converts hydro quantity of a particle at the start of a run
247
248
249
250
251
 *
 * Requires the density to be known
 *
 * @param p The particle to act upon
 */
252
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
253
    struct part *restrict p) {}
254
255
256
257
258

/**
 * @brief Returns the internal energy of a particle
 *
 * @param p The particle of interest
259
 * @param dt Time since the last kick
260
 */
261
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
262
    const struct part *restrict p, float dt) {
263
264
265

  return p->u;
}