hydro.h 19.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 "cosmology.h"
39
#include "dimension.h"
40
#include "equation_of_state.h"
41
#include "hydro_properties.h"
42
#include "hydro_space.h"
43
#include "kernel_hydro.h"
44
#include "minmax.h"
45
46

/**
47
 * @brief Returns the comoving internal energy of a particle
48
49
50
51
52
53
54
 *
 * 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
 */
55
56
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_internal_energy(const struct part *restrict p) {
57

58
  return p->u;
59
60
61
}

/**
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
 * @brief Returns the physical 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 and converts it to
 * physical coordinates.
 *
 * @param p The particle of interest.
 * @param cosmo The cosmological model.
 */
__attribute__((always_inline)) INLINE static float
hydro_get_physical_internal_energy(const struct part *restrict p,
                                   const struct cosmology *cosmo) {

  return p->u * cosmo->a_factor_internal_energy;
}

/**
 * @brief Returns the comoving pressure of a particle
 *
 * Computes the pressure based on the particle's properties.
83
84
85
 *
 * @param p The particle of interest
 */
86
__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
87
    const struct part *restrict p) {
88

89
  return gas_pressure_from_internal_energy(p->rho, p->u);
90
91
92
}

/**
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
 * @brief Returns the physical pressure of a particle
 *
 * Computes the pressure based on the particle's properties and
 * convert it to physical coordinates.
 *
 * @param p The particle of interest
 * @param cosmo The cosmological model.
 */
__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
    const struct part *restrict p, const struct cosmology *cosmo) {

  return cosmo->a_factor_pressure *
         gas_pressure_from_internal_energy(p->rho, p->u);
}

/**
 * @brief Returns the comoving entropy of a particle
110
111
112
113
114
115
116
 *
 * 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
 */
117
__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
118
    const struct part *restrict p) {
119

120
  return gas_entropy_from_internal_energy(p->rho, p->u);
121
122
123
}

/**
124
125
126
127
128
129
 * @brief Returns the physical entropy of a particle
 *
 * For implementations where the main thermodynamic variable
 * is not entropy, this function computes the entropy from
 * the thermodynamic variable and converts it to
 * physical coordinates.
130
131
 *
 * @param p The particle of interest
132
 * @param cosmo The cosmological model.
133
 */
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
    const struct part *restrict p, const struct cosmology *cosmo) {

  /* Note: no cosmological conversion required here with our choice of
   * coordinates. */
  return gas_entropy_from_internal_energy(p->rho, p->u);
}

/**
 * @brief Returns the comoving sound speed of a particle
 *
 * @param p The particle of interest
 */
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_soundspeed(const struct part *restrict p) {
149

150
  return p->force.soundspeed;
151
}
152

153
/**
154
155
156
157
158
159
160
161
162
163
164
165
166
167
 * @brief Returns the physical sound speed of a particle
 *
 * @param p The particle of interest
 * @param cosmo The cosmological model.
 */
__attribute__((always_inline)) INLINE static float
hydro_get_physical_soundspeed(const struct part *restrict p,
                              const struct cosmology *cosmo) {

  return cosmo->a_factor_sound_speed * p->force.soundspeed;
}

/**
 * @brief Returns the comoving density of a particle
168
169
170
 *
 * @param p The particle of interest
 */
171
__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
172
173
174
175
176
    const struct part *restrict p) {

  return p->rho;
}

177
178
179
180
181
182
183
184
185
186
187
188
/**
 * @brief Returns the comoving density of a particle.
 *
 * @param p The particle of interest
 * @param cosmo The cosmological model.
 */
__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
    const struct part *restrict p, const struct cosmology *cosmo) {

  return cosmo->a3_inv * p->rho;
}

189
190
191
192
193
194
195
196
197
198
199
/**
 * @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;
}

200
201
202
203
204
205
206
207
208
209
210
211
/**
 * @brief Sets the mass of a particle
 *
 * @param p The particle of interest
 * @param m The mass to set.
 */
__attribute__((always_inline)) INLINE static void hydro_set_mass(
    struct part *restrict p, float m) {

  p->mass = m;
}

212
213
214
215
216
/**
 * @brief Returns the velocities drifted to the current time of a particle.
 *
 * @param p The particle of interest
 * @param xp The extended data of the particle.
217
218
 * @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
 * @param dt_kick_grav The time (for gravity accelerations) since the last kick.
219
220
221
 * @param v (return) The velocities at the current time.
 */
__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
222
223
224
225
226
227
228
229
230
    const struct part *restrict p, const struct xpart *xp, float dt_kick_hydro,
    float dt_kick_grav, float v[3]) {

  v[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro +
         xp->a_grav[0] * dt_kick_grav;
  v[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro +
         xp->a_grav[1] * dt_kick_grav;
  v[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro +
         xp->a_grav[2] * dt_kick_grav;
231
232
}

233
/**
234
 * @brief Returns the time derivative of internal energy of a particle
235
 *
236
 * We assume a constant density.
237
 *
238
 * @param p The particle of interest
239
 */
240
241
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
    const struct part *restrict p) {
242

243
  return p->u_dt;
244
245
246
}

/**
247
 * @brief Returns the time derivative of internal energy of a particle
248
 *
249
 * We assume a constant density.
250
 *
251
252
 * @param p The particle of interest.
 * @param du_dt The new time derivative of the internal energy.
253
 */
254
255
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
    struct part *restrict p, float du_dt) {
256

257
  p->u_dt = du_dt;
258
}
259
260
261
/**
 * @brief Computes the hydro time-step of a given particle
 *
262
263
264
 * 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.
 *
265
266
 * @param p Pointer to the particle data
 * @param xp Pointer to the extended particle data
267
 * @param hydro_properties The SPH parameters
268
 * @param cosmo The cosmological model.
269
270
 */
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
271
    const struct part *restrict p, const struct xpart *restrict xp,
272
273
    const struct hydro_props *restrict hydro_properties,
    const struct cosmology *restrict cosmo) {
274
275

  const float CFL_condition = hydro_properties->CFL_condition;
276
277

  /* CFL condition */
278
279
  const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
                       (cosmo->a_factor_sound_speed * p->force.v_sig);
280

281
  return dt_cfl;
282
283
}

284
285
286
287
288
289
290
291
292
293
/**
 * @brief Does some extra hydro operations once the actual physical time step
 * for the particle is known.
 *
 * @param p The particle to act upon.
 * @param dt Physical time step of the particle during the next step.
 */
__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
    struct part *p, float dt) {}

294
295
296
297
/**
 * @brief Prepares a particle for the density calculation.
 *
 * Zeroes all the relevant arrays in preparation for the sums taking place in
298
299
 * the various density loop over neighbours. Typically, all fields of the
 * density sub-structure of a particle get zeroed in here.
300
301
 *
 * @param p The particle to act upon
302
 * @param hs #hydro_space containing hydro specific space information.
303
 */
304
__attribute__((always_inline)) INLINE static void hydro_init_part(
305
    struct part *restrict p, const struct hydro_space *hs) {
306

307
308
309
  p->density.wcount = 0.f;
  p->density.wcount_dh = 0.f;
  p->rho = 0.f;
310
  p->density.rho_dh = 0.f;
311
312
313
314
315
316
317
}

/**
 * @brief Finishes the density calculation.
 *
 * Multiplies the density and number of neighbours by the appropiate constants
 * and add the self-contribution term.
318
 * Additional quantities such as velocity gradients will also get the final
319
320
321
 * terms added to them here.
 *
 * Also adds/multiplies the cosmological terms if need be.
322
323
 *
 * @param p The particle to act upon
324
 * @param cosmo The cosmological model.
325
 */
326
__attribute__((always_inline)) INLINE static void hydro_end_density(
327
    struct part *restrict p, const struct cosmology *cosmo) {
328
329
330

  /* Some smoothing length multiples. */
  const float h = p->h;
331
332
333
  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) */
334

335
336
  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
337
  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
338
  p->density.wcount += kernel_root;
339
  p->density.wcount_dh -= hydro_dimension * kernel_root;
340
341

  /* Finish the calculation by inserting the missing h-factors */
342
  p->rho *= h_inv_dim;
343
  p->density.rho_dh *= h_inv_dim_plus_one;
344
  p->density.wcount *= h_inv_dim;
345
346
347
348
349
350
  p->density.wcount_dh *= h_inv_dim_plus_one;
}

/**
 * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
 *
351
352
353
354
 * In the desperate case where a particle has no neighbours (likely because
 * of the h_max ceiling), set the particle fields to something sensible to avoid
 * NaNs in the next calculations.
 *
355
356
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
357
 * @param cosmo The cosmological model.
358
359
 */
__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
360
361
    struct part *restrict p, struct xpart *restrict xp,
    const struct cosmology *cosmo) {
362
363
364
365
366
367
368
369
370
371
372

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

  /* Re-set problematic values */
  p->rho = p->mass * kernel_root * h_inv_dim;
  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
  p->density.rho_dh = 0.f;
  p->density.wcount_dh = 0.f;
373
374
375
376
377
}

/**
 * @brief Prepare a particle for the force calculation.
 *
378
379
380
381
382
383
 * 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.
384
385
386
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
387
 * @param cosmo The current cosmological model.
388
 */
389
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
390
391
    struct part *restrict p, struct xpart *restrict xp,
    const struct cosmology *cosmo) {
392

393
  /* Compute the pressure */
394
  const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
395
396

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

399
  /* Compute the "grad h" term */
400
  const float rho_inv = 1.f / p->rho;
401
402
  const float grad_h_term =
      1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
403

404
405
  /* Update variables. */
  p->force.f = grad_h_term;
406
  p->force.pressure = pressure;
407
  p->force.soundspeed = soundspeed;
408
409
410
411
412
413
}

/**
 * @brief Reset acceleration fields of a particle
 *
 * Resets all hydro acceleration and time derivative fields in preparation
414
 * for the sums taking  place in the various force tasks.
415
416
417
 *
 * @param p The particle to act upon
 */
418
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
419
    struct part *restrict p) {
420
421
422
423
424
425
426

  /* 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. */
427
  p->u_dt = 0.0f;
428
  p->force.h_dt = 0.0f;
429
430
431
  p->force.v_sig = 0.0f;
}

432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
/**
 * @brief Sets the values to be predicted in the drifts to their values at a
 * kick time
 *
 * @param p The particle.
 * @param xp The extended data of this particle.
 */
__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
    struct part *restrict p, const struct xpart *restrict xp) {

  /* Re-set the predicted velocities */
  p->v[0] = xp->v_full[0];
  p->v[1] = xp->v_full[1];
  p->v[2] = xp->v_full[2];

  /* Re-set the entropy */
  p->u = xp->u_full;
}

451
452
453
/**
 * @brief Predict additional particle fields forward in time when drifting
 *
454
455
456
 * Additional hydrodynamic quantites are drifted forward in time here. These
 * include thermal quantities (thermal energy or total energy or entropy, ...).
 *
457
458
459
 * Note the different time-step sizes used for the different quantities as they
 * include cosmological factors.
 *
460
461
 * @param p The particle.
 * @param xp The extended data of the particle.
462
463
 * @param dt_drift The drift time-step for positions.
 * @param dt_therm The drift time-step for thermal quantities.
464
465
 */
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
466
467
    struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
    float dt_therm) {
468
469
470
471

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

  /* Predict smoothing length */
472
  const float w1 = p->force.h_dt * h_inv * dt_drift;
473
474
475
476
477
478
479
480
481
482
483
  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);
484

485
  /* Predict the internal energy */
486
  p->u += p->u_dt * dt_therm;
487
488
489

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

  /* Compute the new sound speed */
492
  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
493
494
495

  p->force.pressure = pressure;
  p->force.soundspeed = soundspeed;
496
497
498
499
500
}

/**
 * @brief Finishes the force calculation.
 *
501
 * Multiplies the force and accelerations by the appropiate constants
502
 * and add the self-contribution term. In most cases, there is little
503
 * to do here.
504
 *
505
506
 * Cosmological terms are also added/multiplied here.
 *
507
 * @param p The particle to act upon
508
 * @param cosmo The current cosmological model.
509
 */
510
__attribute__((always_inline)) INLINE static void hydro_end_force(
511
    struct part *restrict p, const struct cosmology *cosmo) {
512

513
  p->force.h_dt *= p->h * hydro_dimension_inv;
514
}
515
516
517
518

/**
 * @brief Kick the additional variables
 *
519
520
521
 * Additional hydrodynamic quantites are kicked forward in time here. These
 * include thermal quantities (thermal energy or total energy or entropy, ...).
 *
522
523
524
 * @param p The particle to act upon.
 * @param xp The particle extended data to act upon.
 * @param dt_therm The time-step for this kick (for thermodynamic quantities).
525
 */
526
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
527
    struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
528

529
  /* Do not decrease the energy by more than a factor of 2*/
530
531
  if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * xp->u_full) {
    p->u_dt = -0.5f * xp->u_full / dt_therm;
532
  }
533
  xp->u_full += p->u_dt * dt_therm;
534
535

  /* Compute the pressure */
536
  const float pressure = gas_pressure_from_internal_energy(p->rho, xp->u_full);
537

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

541
  p->force.pressure = pressure;
542
  p->force.soundspeed = soundspeed;
543
}
544
545

/**
546
 * @brief Converts hydro quantity of a particle at the start of a run
547
 *
548
549
550
551
 * 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.
552
553
 *
 * @param p The particle to act upon
554
 * @param xp The extended particle to act upon
555
 * @param cosmo The cosmological model.
556
 */
557
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
558
559
    struct part *restrict p, struct xpart *restrict xp,
    const struct cosmology *cosmo) {
560
561
562

  /* Compute the pressure */
  const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
563
564
565
566

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

567
  p->force.pressure = pressure;
568
  p->force.soundspeed = soundspeed;
569
}
570

571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
/**
 * @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 or assignments between the particle
 * and extended particle fields.
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
 */
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
    struct part *restrict p, struct xpart *restrict xp) {

  p->time_bin = 0;
  xp->v_full[0] = p->v[0];
  xp->v_full[1] = p->v[1];
  xp->v_full[2] = p->v[2];
588
589
590
  xp->a_grav[0] = 0.f;
  xp->a_grav[1] = 0.f;
  xp->a_grav[2] = 0.f;
591
592
593
  xp->u_full = p->u;

  hydro_reset_acceleration(p);
594
  hydro_init_part(p, NULL);
595
596
}

597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
/**
 * @brief Overwrite the initial internal energy of a particle.
 *
 * Note that in the cases where the thermodynamic variable is not
 * internal energy but gets converted later, we must overwrite that
 * field. The conversion to the actual variable happens later after
 * the initial fake time-step.
 *
 * @param p The #part to write to.
 * @param u_init The new initial internal energy.
 */
__attribute__((always_inline)) INLINE static void
hydro_set_init_internal_energy(struct part *p, float u_init) {

  p->u = u_init;
}

614
#endif /* SWIFT_MINIMAL_HYDRO_H */