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

22
23
24
25
26
27
28
29
30
31
32
33
/**
 * @file Gadget2/hydro.h
 * @brief SPH interaction functions following the Gadget-2 version of SPH.
 *
 * The interactions computed here are the ones presented in the Gadget-2 paper
 * Springel, V., MNRAS, Volume 364, Issue 4, pp. 1105-1134.
 * We use the same numerical coefficients as the Gadget-2 code. When used with
 * the Spline-3 kernel, the results should be equivalent to the ones obtained
 * with Gadget-2 up to the rounding errors and interactions missed by the
 * Gadget-2 tree-code neighbours search.
 */

34
#include "adiabatic_index.h"
35
#include "approx_math.h"
36
#include "cosmology.h"
37
#include "dimension.h"
38
#include "equation_of_state.h"
39
#include "hydro_properties.h"
40
#include "hydro_space.h"
41
#include "kernel_hydro.h"
42
#include "minmax.h"
43
44

/**
45
46
 * @brief Returns the comoving internal energy of a particle at the last
 * time the particle was kicked.
47
48
 *
 * @param p The particle of interest
49
 * @param xp The extended data of the particle of interest.
50
 */
51
__attribute__((always_inline)) INLINE static float
52
53
hydro_get_comoving_internal_energy(const struct part *restrict p,
                                   const struct xpart *restrict xp) {
54

55
  return gas_internal_energy_from_entropy(p->rho, xp->entropy_full);
56
57
}

58
/**
59
60
 * @brief Returns the physical internal energy of a particle at the last
 * time the particle was kicked.
61
62
 *
 * @param p The particle of interest.
63
 * @param xp The extended data of the particle of interest.
64
65
66
67
 * @param cosmo The cosmological model.
 */
__attribute__((always_inline)) INLINE static float
hydro_get_physical_internal_energy(const struct part *restrict p,
68
                                   const struct xpart *restrict xp,
69
70
                                   const struct cosmology *cosmo) {

71
72
  return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv,
                                          xp->entropy_full);
73
74
}

75
/**
76
77
78
79
80
81
82
83
84
85
86
87
88
89
 * @brief Returns the comoving internal energy of a particle drifted to the
 * current time.
 *
 * @param p The particle of interest
 */
__attribute__((always_inline)) INLINE static float
hydro_get_drifted_comoving_internal_energy(const struct part *restrict p) {

  return gas_internal_energy_from_entropy(p->rho, p->entropy);
}

/**
 * @brief Returns the physical internal energy of a particle drifted to the
 * current time.
90
91
92
93
94
 *
 * @param p The particle of interest.
 * @param cosmo The cosmological model.
 */
__attribute__((always_inline)) INLINE static float
95
96
hydro_get_drifted_physical_internal_energy(const struct part *restrict p,
                                           const struct cosmology *cosmo) {
97

98
  return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
99
100
}

101
/**
102
 * @brief Returns the comoving pressure of a particle
103
104
105
 *
 * @param p The particle of interest
 */
106
__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
107
    const struct part *restrict p) {
108

109
  return gas_pressure_from_entropy(p->rho, p->entropy);
110
111
}

112
113
114
/**
 * @brief Returns the physical pressure of a particle
 *
115
116
 * @param p The particle of interest.
 * @param cosmo The cosmological model.
117
118
119
120
121
122
123
 */
__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
    const struct part *restrict p, const struct cosmology *cosmo) {

  return gas_pressure_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
}

124
/**
125
126
 * @brief Returns the comoving entropy of a particle at the last
 * time the particle was kicked.
127
128
 *
 * @param p The particle of interest.
129
 * @param xp The extended data of the particle of interest.
130
131
 */
__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
132
    const struct part *restrict p, const struct xpart *restrict xp) {
133

134
  return xp->entropy_full;
135
136
}

137
/**
138
139
 * @brief Returns the physical entropy of a particl at the last
 * time the particle was kicked.
140
 *
141
 * @param p The particle of interest.
142
 * @param xp The extended data of the particle of interest.
143
 * @param cosmo The cosmological model.
144
 */
145
__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
    const struct part *restrict p, const struct xpart *restrict xp,
    const struct cosmology *cosmo) {

  /* Note: no cosmological conversion required here with our choice of
   * coordinates. */
  return xp->entropy_full;
}

/**
 * @brief Returns the comoving entropy of a particle drifted to the
 * current time.
 *
 * @param p The particle of interest.
 */
__attribute__((always_inline)) INLINE static float
hydro_get_drifted_comoving_entropy(const struct part *restrict p) {

  return p->entropy;
}

/**
 * @brief Returns the physical entropy of a particle drifted to the
 * current time.
 *
 * @param p The particle of interest.
 * @param cosmo The cosmological model.
 */
__attribute__((always_inline)) INLINE static float
hydro_get_drifted_physical_entropy(const struct part *restrict p,
                                   const struct cosmology *cosmo) {
176

177
178
  /* Note: no cosmological conversion required here with our choice of
   * coordinates. */
179
  return p->entropy;
180
181
182
}

/**
183
 * @brief Returns the comoving sound speed of a particle
184
185
186
 *
 * @param p The particle of interest
 */
187
188
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_soundspeed(const struct part *restrict p) {
189
190
191

  return p->force.soundspeed;
}
192

193
194
195
196
197
198
199
200
201
202
203
204
205
/**
 * @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;
}

206
/**
207
 * @brief Returns the comoving density of a particle
208
209
210
 *
 * @param p The particle of interest
 */
211
__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
212
213
214
215
216
    const struct part *restrict p) {

  return p->rho;
}

217
218
219
220
221
222
223
224
225
226
227
228
/**
 * @brief Returns the physical 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 p->rho * cosmo->a3_inv;
}

229
230
231
232
233
234
235
236
237
238
239
/**
 * @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;
}

240
241
242
243
244
245
246
247
248
249
250
251
/**
 * @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;
}

252
253
254
255
256
/**
 * @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.
257
258
 * @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.
259
260
261
 * @param v (return) The velocities at the current time.
 */
__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
262
263
264
265
266
267
268
269
270
    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;
271
272
}

273
/**
274
 * @brief Returns the time derivative of co-moving internal energy of a particle
275
 *
276
 * We assume a constant density.
277
 *
278
 * @param p The particle of interest
279
 */
280
281
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_internal_energy_dt(const struct part *restrict p) {
282

283
  return gas_internal_energy_from_entropy(p->rho, p->entropy_dt);
284
285
286
}

/**
287
 * @brief Returns the time derivative of physical internal energy of a particle
288
 *
289
 * We assume a constant density.
290
 *
291
 * @param p The particle of interest.
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
 * @param cosmo The cosmological model.
 */
__attribute__((always_inline)) INLINE static float
hydro_get_physical_internal_energy_dt(const struct part *restrict p,
                                      const struct cosmology *cosmo) {

  return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv,
                                          p->entropy_dt);
}

/**
 * @brief Sets the time derivative of the co-moving internal energy of a
 * particle
 *
 * We assume a constant density for the conversion to entropy.
 *
 * @param p The particle of interest.
309
 * @param du_dt The new time derivative of the internal energy.
310
 */
311
312
__attribute__((always_inline)) INLINE static void
hydro_set_comoving_internal_energy_dt(struct part *restrict p, float du_dt) {
313

314
  p->entropy_dt = gas_entropy_from_internal_energy(p->rho, du_dt);
315
316
}

317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
/**
 * @brief Sets the time derivative of the physical internal energy of a particle
 *
 * We assume a constant density for the conversion to entropy.
 *
 * @param p The particle of interest.
 * @param cosmo Cosmology data structure
 * @param du_dt The time derivative of the internal energy.
 */
__attribute__((always_inline)) INLINE static void
hydro_set_physical_internal_energy_dt(struct part *restrict p,
                                      const struct cosmology *restrict cosmo,
                                      float du_dt) {
  p->entropy_dt =
      gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, du_dt);
}

334
335
336
337
338
/**
 * @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
339
 * @param hydro_properties The constants used in the scheme
340
 * @param cosmo The cosmological model.
341
342
 */
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
343
    const struct part *restrict p, const struct xpart *restrict xp,
344
345
    const struct hydro_props *restrict hydro_properties,
    const struct cosmology *restrict cosmo) {
346

347
  const float CFL = hydro_properties->CFL_condition;
348

349
  /* CFL condition */
350
  const float dt_cfl = 2.f * kernel_gamma * CFL * cosmo->a * p->h /
351
                       (cosmo->a_factor_sound_speed * p->force.v_sig);
352

353
  return dt_cfl;
354
355
}

356
357
358
359
360
361
362
363
364
365
/**
 * @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) {}

366
367
368
/**
 * @brief Prepares a particle for the density calculation.
 *
369
 * Zeroes all the relevant arrays in preparation for the sums taking place in
370
371
372
 * the variaous density tasks
 *
 * @param p The particle to act upon
373
 * @param hs #hydro_space containing hydro specific space information.
374
 */
375
__attribute__((always_inline)) INLINE static void hydro_init_part(
376
    struct part *restrict p, const struct hydro_space *hs) {
377

378
#ifdef DEBUG_INTERACTIONS_SPH
379
  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS; ++i) p->ids_ngbs_density[i] = -1;
380
381
  p->num_ngb_density = 0;
#endif
382

383
  p->rho = 0.f;
384
385
  p->density.wcount = 0.f;
  p->density.wcount_dh = 0.f;
386
  p->density.rho_dh = 0.f;
387
  p->density.div_v = 0.f;
388
389
390
  p->density.rot_v[0] = 0.f;
  p->density.rot_v[1] = 0.f;
  p->density.rot_v[2] = 0.f;
391
392
393
}

/**
394
 * @brief Finishes the density calculation.
395
 *
396
 * Multiplies the density and number of neighbours by the appropiate constants
397
398
399
 * and add the self-contribution term.
 *
 * @param p The particle to act upon
400
 * @param cosmo The current cosmological model.
401
 */
402
__attribute__((always_inline)) INLINE static void hydro_end_density(
403
    struct part *restrict p, const struct cosmology *cosmo) {
404
405
406

  /* Some smoothing length multiples. */
  const float h = p->h;
407
408
409
  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) */
410

411
412
  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
413
  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
414
  p->density.wcount += kernel_root;
415
  p->density.wcount_dh -= hydro_dimension * kernel_root;
416

417
  /* Finish the calculation by inserting the missing h-factors */
418
  p->rho *= h_inv_dim;
419
  p->density.rho_dh *= h_inv_dim_plus_one;
420
421
  p->density.wcount *= h_inv_dim;
  p->density.wcount_dh *= h_inv_dim_plus_one;
422

423
  const float rho_inv = 1.f / p->rho;
424
  const float a_inv2 = cosmo->a2_inv;
Matthieu Schaller's avatar
Matthieu Schaller committed
425

426
427
428
429
  /* Finish calculation of the (physical) velocity curl components */
  p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
  p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
  p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
430

431
432
  /* Finish calculation of the (physical) velocity divergence */
  p->density.div_v *= h_inv_dim_plus_one * a_inv2 * rho_inv;
433
434
}

435
436
437
438
439
/**
 * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
440
 * @param cosmo The current cosmological model.
441
442
 */
__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
443
444
    struct part *restrict p, struct xpart *restrict xp,
    const struct cosmology *cosmo) {
445
446
447
448
449
450
451
452

  /* 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;
453
  p->density.wcount = kernel_root * h_inv_dim;
454
455
456
457
458
459
460
461
  p->density.rho_dh = 0.f;
  p->density.wcount_dh = 0.f;
  p->density.div_v = 0.f;
  p->density.rot_v[0] = 0.f;
  p->density.rot_v[1] = 0.f;
  p->density.rot_v[2] = 0.f;
}

462
463
464
/**
 * @brief Prepare a particle for the force calculation.
 *
465
466
467
468
469
470
 * 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.
471
472
473
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
474
 * @param cosmo The current cosmological model.
Matthieu Schaller's avatar
Matthieu Schaller committed
475
 * @param hydro_props Hydrodynamic properties.
476
477
 * @param dt_alpha The time-step used to evolve non-cosmological quantities such
 *                 as the artificial viscosity.
478
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
479
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
480
    struct part *restrict p, struct xpart *restrict xp,
Josh Borrow's avatar
Josh Borrow committed
481
482
    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
    const float dt_alpha) {
483

484
  const float fac_Balsara_eps = cosmo->a_factor_Balsara_eps;
485

486
  /* Inverse of the co-moving density */
487
488
  const float rho_inv = 1.f / p->rho;

489
490
491
  /* Inverse of the smoothing length */
  const float h_inv = 1.f / p->h;

492
  /* Compute the norm of the curl */
493
494
495
  const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
                             p->density.rot_v[1] * p->density.rot_v[1] +
                             p->density.rot_v[2] * p->density.rot_v[2]);
496

497
498
499
  /* Compute the norm of div v including the Hubble flow term */
  const float div_physical_v = p->density.div_v + 3.f * cosmo->H;
  const float abs_div_physical_v = fabsf(div_physical_v);
Matthieu Schaller's avatar
Matthieu Schaller committed
500

501
  /* Compute the pressure */
502
  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
Matthieu Schaller's avatar
Matthieu Schaller committed
503

504
505
  /* Compute the sound speed */
  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
Matthieu Schaller's avatar
Matthieu Schaller committed
506

507
  /* Divide the pressure by the density squared to get the SPH term */
508
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
509

510
  /* Compute the Balsara switch */
Matthieu Schaller's avatar
Matthieu Schaller committed
511
512
  /* Pre-multiply in the AV factor; hydro_props are not passed to the iact
   * functions */
513
514
515
  const float balsara = hydro_props->viscosity.alpha * abs_div_physical_v /
                        (abs_div_physical_v + curl_v +
                         0.0001f * fac_Balsara_eps * soundspeed * h_inv);
Matthieu Schaller's avatar
Matthieu Schaller committed
516

517
  /* Compute the "grad h" term */
518
  const float omega_inv =
519
520
      1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);

521
  /* Update variables. */
522
  p->force.f = omega_inv;
Matthieu Schaller's avatar
Matthieu Schaller committed
523
  p->force.P_over_rho2 = P_over_rho2;
524
525
  p->force.soundspeed = soundspeed;
  p->force.balsara = balsara;
526
527
528
529
530
531
532
533
534
535
}

/**
 * @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
 */
536
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
537
    struct part *restrict p) {
538

539
#ifdef DEBUG_INTERACTIONS_SPH
540
  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS; ++i) p->ids_ngbs_force[i] = -1;
541
542
  p->num_ngb_force = 0;
#endif
543

544
  /* Reset the acceleration. */
Matthieu Schaller's avatar
Matthieu Schaller committed
545
546
547
  p->a_hydro[0] = 0.0f;
  p->a_hydro[1] = 0.0f;
  p->a_hydro[2] = 0.0f;
548

549
  /* Reset the time derivatives. */
550
  p->entropy_dt = 0.0f;
551
  p->force.h_dt = 0.0f;
552

553
  /* Reset maximal signal velocity */
554
  p->force.v_sig = p->force.soundspeed;
555
556
}

557
558
559
560
561
562
563
/**
 * @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.
 */
564
565
566
567
568
569
570
__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];
571
572
573

  /* Re-set the entropy */
  p->entropy = xp->entropy_full;
574
575
}

576
577
578
579
580
/**
 * @brief Predict additional particle fields forward in time when drifting
 *
 * @param p The particle
 * @param xp The extended data of the particle
581
582
 * @param dt_drift The drift time-step for positions.
 * @param dt_therm The drift time-step for thermal quantities.
583
 */
584
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
585
586
    struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
    float dt_therm) {
587
588
589
590

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

  /* Predict smoothing length */
591
  const float w1 = p->force.h_dt * h_inv * dt_drift;
592
593
594
595
596
597
598
599
600
601
602
  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);
603

604
605
    /* Predict the entropy */
#ifdef SWIFT_DEBUG_CHECKS
606
607
  if (p->entropy + p->entropy_dt * dt_therm <= 0)
    error(
608
        "Negative entropy for particle id %llu old entropy %.5e d_entropy %.5e "
609
610
        "entropy_dt %.5e dt therm %.5e",
        p->id, p->entropy, p->entropy_dt * dt_therm, p->entropy_dt, dt_therm);
611
#endif
612
  p->entropy += p->entropy_dt * dt_therm;
Matthieu Schaller's avatar
Matthieu Schaller committed
613

614
615
  /* Re-compute the pressure */
  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
616

617
618
  /* Compute the new sound speed */
  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
619

620
621
  /* Divide the pressure by the density squared to get the SPH term */
  const float rho_inv = 1.f / p->rho;
622
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
Matthieu Schaller's avatar
Matthieu Schaller committed
623

624
625
  /* Update variables */
  p->force.soundspeed = soundspeed;
626
  p->force.P_over_rho2 = P_over_rho2;
627
628
629
}

/**
630
 * @brief Finishes the force calculation.
631
 *
632
 * Multiplies the forces and accelerationsby the appropiate constants
633
634
 *
 * @param p The particle to act upon
635
 * @param cosmo The current cosmological model.
636
 */
637
__attribute__((always_inline)) INLINE static void hydro_end_force(
638
    struct part *restrict p, const struct cosmology *cosmo) {
639

640
  p->force.h_dt *= p->h * hydro_dimension_inv;
641

642
643
  p->entropy_dt =
      0.5f * gas_entropy_from_internal_energy(p->rho, p->entropy_dt);
644
645
}

Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
646
647
648
649
/**
 * @brief Kick the additional variables
 *
 * @param p The particle to act upon
650
 * @param xp The particle extended data to act upon
651
 * @param dt_therm The time-step for this kick (for thermodynamic quantities)
652
653
654
 * @param dt_grav The time-step for this kick (for gravity forces)
 * @param dt_hydro The time-step for this kick (for hydro forces)
 * @param dt_kick_corr The time-step for this kick (for correction of the kick)
655
656
 * @param cosmo The cosmological model.
 * @param hydro_props The constants used in the scheme
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
657
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
658
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
659
    struct part *restrict p, struct xpart *restrict xp, float dt_therm,
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
660
    float dt_grav, float dt_hydro, float dt_kick_corr,
661
    const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
662

663
  /* Do not decrease the entropy by more than a factor of 2 */
664
665
666
667
  if (dt_therm > 0. && p->entropy_dt * dt_therm < -0.5f * xp->entropy_full) {
    p->entropy_dt = -0.5f * xp->entropy_full / dt_therm;
  }
  xp->entropy_full += p->entropy_dt * dt_therm;
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
668

669
  /* Apply the minimal energy limit */
670
671
672
673
674
675
676
  const float physical_density = p->rho * cosmo->a3_inv;
  const float min_physical_energy = hydro_props->minimal_internal_energy;
  const float min_physical_entropy =
      gas_entropy_from_internal_energy(physical_density, min_physical_energy);
  const float min_comoving_entropy = min_physical_entropy; /* A' = A */
  if (xp->entropy_full < min_comoving_entropy) {
    xp->entropy_full = min_comoving_entropy;
677
678
679
    p->entropy_dt = 0.f;
  }

680
  /* Compute the pressure */
681
  const float pressure = gas_pressure_from_entropy(p->rho, xp->entropy_full);
682

683
684
685
686
  /* Compute the new sound speed */
  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);

  /* Divide the pressure by the density squared to get the SPH term */
687
  const float rho_inv = 1.f / p->rho;
688
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
689
690

  p->force.soundspeed = soundspeed;
691
  p->force.P_over_rho2 = P_over_rho2;
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
692
693
}

694
/**
695
 *  @brief Converts hydro quantity of a particle at the start of a run
696
697
698
 *
 * Requires the density to be known
 *
699
700
701
 * @param p The particle to act upon.
 * @param xp The extended data.
 * @param cosmo The cosmological model.
702
 * @param hydro_props The constants used in the scheme.
703
 */
704
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
705
    struct part *restrict p, struct xpart *restrict xp,
706
    const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
707

708
709
  /* We read u in the entropy field. We now get (comoving) A from (physical) u
   * and (physical) rho. Note that comoving A (A') == physical A */
710
711
  xp->entropy_full =
      gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, p->entropy);
712
  p->entropy = xp->entropy_full;
Matthieu Schaller's avatar
Matthieu Schaller committed
713

714
  /* Apply the minimal energy limit */
715
716
717
718
719
720
721
722
  const float physical_density = p->rho * cosmo->a3_inv;
  const float min_physical_energy = hydro_props->minimal_internal_energy;
  const float min_physical_entropy =
      gas_entropy_from_internal_energy(physical_density, min_physical_energy);
  const float min_comoving_entropy = min_physical_entropy; /* A' = A */
  if (xp->entropy_full < min_comoving_entropy) {
    xp->entropy_full = min_comoving_entropy;
    p->entropy = min_comoving_entropy;
723
724
725
    p->entropy_dt = 0.f;
  }

726
727
728
  /* Compute the pressure */
  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);

729
730
731
732
  /* Compute the sound speed */
  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);

  /* Divide the pressure by the density squared to get the SPH term */
733
  const float rho_inv = 1.f / p->rho;
734
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
735
736

  p->force.soundspeed = soundspeed;
737
  p->force.P_over_rho2 = P_over_rho2;
738
}
739

Matthieu Schaller's avatar
Matthieu Schaller committed
740
741
742
743
744
745
746
747
748
749
750
751
752
/**
 * @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
 */
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
    struct part *restrict p, struct xpart *restrict xp) {

  p->time_bin = 0;
753
  p->wakeup = time_bin_not_awake;
Matthieu Schaller's avatar
Matthieu Schaller committed
754
755
756
  xp->v_full[0] = p->v[0];
  xp->v_full[1] = p->v[1];
  xp->v_full[2] = p->v[2];
757
758
759
  xp->a_grav[0] = 0.f;
  xp->a_grav[1] = 0.f;
  xp->a_grav[2] = 0.f;
Matthieu Schaller's avatar
Matthieu Schaller committed
760
761
762
  xp->entropy_full = p->entropy;

  hydro_reset_acceleration(p);
763
  hydro_init_part(p, NULL);
Matthieu Schaller's avatar
Matthieu Schaller committed
764
765
}

766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
/**
 * @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->entropy = u_init;
}

783
#endif /* SWIFT_GADGET2_HYDRO_H */