hydro.h 14.2 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 *
 * 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/>.
 *
 ******************************************************************************/
19
20
21
22
23
24
25
26
27
28
29
30
31
#ifndef SWIFT_MINIMAL_HYDRO_H
#define SWIFT_MINIMAL_HYDRO_H

/**
 * @file Minimal/hydro.h
 * @brief Minimal conservative implementation of SPH (Non-neighbour loop
 * equations)
 *
 * The thermal variable is the internal energy (u). Simple constant
 * viscosity term without switches is implemented. No thermal conduction
 * term is implemented.
 *
 * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
32
33
 * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
 * Physics, 2012, Volume 231, Issue 3, pp. 759-794.
34
 */
35

36
#include "adiabatic_index.h"
37
#include "approx_math.h"
38
#include "dimension.h"
39
#include "equation_of_state.h"
40
41
#include "hydro_properties.h"
#include "kernel_hydro.h"
42
#include "minmax.h"
43
44
45
46
47
48
49
50
51
52
53
54
55
56

/**
 * @brief Returns the internal energy of a particle
 *
 * For implementations where the main thermodynamic variable
 * is not internal energy, this function computes the internal
 * energy from the thermodynamic variable.
 *
 * @param p The particle of interest
 * @param dt Time since the last kick
 */
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
    const struct part *restrict p, float dt) {

57
  return p->u + p->u_dt * dt;
58
59
60
61
62
63
64
65
66
67
68
}

/**
 * @brief Returns the pressure of a particle
 *
 * @param p The particle of interest
 * @param dt Time since the last kick
 */
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
    const struct part *restrict p, float dt) {

69
70
71
  const float u = p->u + p->u_dt * dt;

  return gas_pressure_from_internal_energy(p->rho, u);
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
}

/**
 * @brief Returns the entropy of a particle
 *
 * For implementations where the main thermodynamic variable
 * is not entropy, this function computes the entropy from
 * the thermodynamic variable.
 *
 * @param p The particle of interest
 * @param dt Time since the last kick
 */
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
    const struct part *restrict p, float dt) {

87
88
89
  const float u = p->u + p->u_dt * dt;

  return gas_entropy_from_internal_energy(p->rho, u);
90
91
92
93
94
95
96
97
98
99
100
}

/**
 * @brief Returns the sound speed of a particle
 *
 * @param p The particle of interest
 * @param dt Time since the last kick
 */
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
    const struct part *restrict p, float dt) {

101
  return p->force.soundspeed;
102
}
103

104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
/**
 * @brief Returns the density of a particle
 *
 * @param p The particle of interest
 */
__attribute__((always_inline)) INLINE static float hydro_get_density(
    const struct part *restrict p) {

  return p->rho;
}

/**
 * @brief Returns the mass of a particle
 *
 * @param p The particle of interest
 */
__attribute__((always_inline)) INLINE static float hydro_get_mass(
    const struct part *restrict p) {

  return p->mass;
}

126
127
128
129
/**
 * @brief Modifies the thermal state of a particle to the imposed internal
 * energy
 *
130
131
132
 * This overwrites the current state of the particle but does *not* change its
 * time-derivatives. Internal energy, pressure, sound-speed and signal velocity
 * will be updated.
133
134
135
136
137
138
139
140
 *
 * @param p The particle
 * @param u The new internal energy
 */
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
    struct part *restrict p, float u) {

  p->u = u;
141

142
  /* Compute the new pressure */
143
  const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
144

145
146
  /* Compute the new sound speed */
  const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
147

148
149
150
151
152
  /* Update the signal velocity */
  const float v_sig_old = p->force.v_sig;
  const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
  const float v_sig = max(v_sig_old, v_sig_new);

153
  p->force.soundspeed = soundspeed;
154
  p->force.pressure = pressure;
155
  p->force.v_sig = v_sig;
156
157
158
159
160
}

/**
 * @brief Modifies the thermal state of a particle to the imposed entropy
 *
161
162
163
 * This overwrites the current state of the particle but does *not* change its
 * time-derivatives. Internal energy, pressure, sound-speed and signal velocity
 * will be updated.
164
165
166
167
168
169
170
171
 *
 * @param p The particle
 * @param S The new entropy
 */
__attribute__((always_inline)) INLINE static void hydro_set_entropy(
    struct part *restrict p, float S) {

  p->u = gas_internal_energy_from_entropy(p->rho, S);
172
173
174

  /* Compute the pressure */
  const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
175

176
177
  /* Compute the new sound speed */
  const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
178

179
180
181
182
183
  /* Update the signal velocity */
  const float v_sig_old = p->force.v_sig;
  const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
  const float v_sig = max(v_sig_old, v_sig_new);

184
  p->force.soundspeed = soundspeed;
185
  p->force.pressure = pressure;
186
  p->force.v_sig = v_sig;
187
188
}

189
190
191
/**
 * @brief Computes the hydro time-step of a given particle
 *
192
193
194
 * This function returns the time-step of a particle given its hydro-dynamical
 * state. A typical time-step calculation would be the use of the CFL condition.
 *
195
196
 * @param p Pointer to the particle data
 * @param xp Pointer to the extended particle data
197
 * @param hydro_properties The SPH parameters
198
199
200
 *
 */
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
201
202
    const struct part *restrict p, const struct xpart *restrict xp,
    const struct hydro_props *restrict hydro_properties) {
203
204

  const float CFL_condition = hydro_properties->CFL_condition;
205
206

  /* CFL condition */
207
208
  const float dt_cfl =
      2.f * kernel_gamma * CFL_condition * p->h / p->force.v_sig;
209

210
  return dt_cfl;
211
212
}

213
214
215
216
/**
 * @brief Initialises the particles for the first time
 *
 * This function is called only once just after the ICs have been
217
218
 * read in to do some conversions or assignments between the particle
 * and extended particle fields.
219
220
221
222
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
 */
223
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
224
    struct part *restrict p, struct xpart *restrict xp) {
225

226
227
228
229
230
  p->ti_begin = 0;
  p->ti_end = 0;
  xp->v_full[0] = p->v[0];
  xp->v_full[1] = p->v[1];
  xp->v_full[2] = p->v[2];
231
232
}

233
234
235
236
/**
 * @brief Prepares a particle for the density calculation.
 *
 * Zeroes all the relevant arrays in preparation for the sums taking place in
237
238
 * the various density loop over neighbours. Typically, all fields of the
 * density sub-structure of a particle get zeroed in here.
239
240
241
 *
 * @param p The particle to act upon
 */
242
__attribute__((always_inline)) INLINE static void hydro_init_part(
243
    struct part *restrict p) {
244

245
246
247
  p->density.wcount = 0.f;
  p->density.wcount_dh = 0.f;
  p->rho = 0.f;
248
  p->density.rho_dh = 0.f;
249
250
251
252
253
254
255
}

/**
 * @brief Finishes the density calculation.
 *
 * Multiplies the density and number of neighbours by the appropiate constants
 * and add the self-contribution term.
256
257
258
 * Additional quantities such as velocity gradients will also get the final
 *terms
 * added to them here.
259
260
261
 *
 * @param p The particle to act upon
 */
262
__attribute__((always_inline)) INLINE static void hydro_end_density(
263
    struct part *restrict p) {
264
265
266

  /* Some smoothing length multiples. */
  const float h = p->h;
267
268
269
  const float h_inv = 1.0f / h;                       /* 1/h */
  const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
  const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
270

271
272
  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
273
  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
274
275
276
  p->density.wcount += kernel_root;

  /* Finish the calculation by inserting the missing h-factors */
277
  p->rho *= h_inv_dim;
278
  p->density.rho_dh *= h_inv_dim_plus_one;
279
  p->density.wcount *= kernel_norm;
280
  p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
281
282
283
284
285
}

/**
 * @brief Prepare a particle for the force calculation.
 *
286
287
288
289
290
291
 * This function is called in the ghost task to convert some quantities coming
 * from the density loop over neighbours into quantities ready to be used in the
 * force loop over neighbours. Quantities are typically read from the density
 * sub-structure and written to the force sub-structure.
 * Examples of calculations done here include the calculation of viscosity term
 * constants, thermal conduction terms, hydro conversions, etc.
292
293
294
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
295
296
 * @param ti_current The current time (on the timeline)
 * @param timeBase The minimal time-step size
297
 */
298
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
299
300
    struct part *restrict p, struct xpart *restrict xp, int ti_current,
    double timeBase) {
301

302
  /* Compute the pressure */
303
304
  const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
  const float pressure = hydro_get_pressure(p, half_dt);
305
306

  /* Compute the sound speed */
307
  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
308

309
  /* Compute the "grad h" term */
310
  const float rho_inv = 1.f / p->rho;
311
312
  const float grad_h_term =
      1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
313

314
315
  /* Update variables. */
  p->force.f = grad_h_term;
316
  p->force.pressure = pressure;
317
  p->force.soundspeed = soundspeed;
318
319
320
321
322
323
}

/**
 * @brief Reset acceleration fields of a particle
 *
 * Resets all hydro acceleration and time derivative fields in preparation
324
 * for the sums taking  place in the various force tasks.
325
326
327
 *
 * @param p The particle to act upon
 */
328
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
329
    struct part *restrict p) {
330
331
332
333
334
335
336

  /* Reset the acceleration. */
  p->a_hydro[0] = 0.0f;
  p->a_hydro[1] = 0.0f;
  p->a_hydro[2] = 0.0f;

  /* Reset the time derivatives. */
337
  p->u_dt = 0.0f;
338
  p->force.h_dt = 0.0f;
339
340
341
342
343
344
  p->force.v_sig = 0.0f;
}

/**
 * @brief Predict additional particle fields forward in time when drifting
 *
345
346
347
 * Additional hydrodynamic quantites are drifted forward in time here. These
 * include thermal quantities (thermal energy or total energy or entropy, ...).
 *
348
349
350
351
352
353
 * @param p The particle.
 * @param xp The extended data of the particle.
 * @param dt The drift time-step.
 * @param t0 The time at the start of the drift (on the timeline).
 * @param t1 The time at the end of the drift (on the timeline).
 * @param timeBase The minimal time-step size.
354
355
 */
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
    struct part *restrict p, const struct xpart *restrict xp, float dt, int t0,
    int t1, double timeBase) {

  const float h_inv = 1.f / p->h;

  /* Predict smoothing length */
  const float w1 = p->force.h_dt * h_inv * dt;
  if (fabsf(w1) < 0.2f)
    p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
  else
    p->h *= expf(w1);

  /* Predict density */
  const float w2 = -hydro_dimension * w1;
  if (fabsf(w2) < 0.2f)
    p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
  else
    p->rho *= expf(w2);
374

375
376
  /* Drift the pressure */
  const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
377
378
379
  const float pressure = hydro_get_pressure(p, dt_entr);

  /* Compute the new sound speed */
380
  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
381
382
383

  p->force.pressure = pressure;
  p->force.soundspeed = soundspeed;
384
385
386
387
388
}

/**
 * @brief Finishes the force calculation.
 *
389
390
391
 * Multiplies the force and accelerations by the appropiate constants
 * and add the self-contribution term. In most cases, there is nothing
 * to do here.
392
393
394
 *
 * @param p The particle to act upon
 */
395
__attribute__((always_inline)) INLINE static void hydro_end_force(
396
397
    struct part *restrict p) {

398
  p->force.h_dt *= p->h * hydro_dimension_inv;
399
}
400
401
402
403

/**
 * @brief Kick the additional variables
 *
404
405
406
 * Additional hydrodynamic quantites are kicked forward in time here. These
 * include thermal quantities (thermal energy or total energy or entropy, ...).
 *
407
 * @param p The particle to act upon
408
 * @param xp The particle extended data to act upon
409
 * @param dt The time-step for this kick
410
 * @param half_dt The half time-step for this kick
411
 */
412
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
413
414
    struct part *restrict p, struct xpart *restrict xp, float dt,
    float half_dt) {
415

416
417
418
419
420
421
  /* Do not decrease the energy by more than a factor of 2*/
  const float u_change = p->u_dt * dt;
  if (u_change > -0.5f * p->u)
    p->u += u_change;
  else
    p->u *= 0.5f;
422

423
424
  /* Do not 'overcool' when timestep increases */
  if (p->u + p->u_dt * half_dt < 0.5f * p->u) p->u_dt = -0.5f * p->u / half_dt;
425
426
427

  /* Compute the pressure */
  const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
428

429
430
  /* Compute the sound speed */
  const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
431

432
  p->force.pressure = pressure;
433
  p->force.soundspeed = soundspeed;
434
}
435
436

/**
437
 * @brief Converts hydro quantity of a particle at the start of a run
438
 *
439
440
441
442
 * This function is called once at the end of the engine_init_particle()
 * routine (at the start of a calculation) after the densities of
 * particles have been computed.
 * This can be used to convert internal energy into entropy for instance.
443
444
445
 *
 * @param p The particle to act upon
 */
446
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
447
448
449
450
    struct part *restrict p) {

  /* Compute the pressure */
  const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
451
452
453
454

  /* Compute the sound speed */
  const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);

455
  p->force.pressure = pressure;
456
  p->force.soundspeed = soundspeed;
457
}
458
459

#endif /* SWIFT_MINIMAL_HYDRO_H */