hydro.h 28.1 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
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2016 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/>.
 *
 ******************************************************************************/
#ifndef SWIFT_PRESSURE_ENTROPY_HYDRO_H
#define SWIFT_PRESSURE_ENTROPY_HYDRO_H

/**
23
24
25
 * @file PressureEntropy/hydro.h
 * @brief Pressure-Entropy implementation of SPH (Non-neighbour loop
 * equations)
26
27
28
29
 *
 * The thermal variable is the entropy (S) and the entropy is smoothed over
 * contact discontinuities to prevent spurious surface tension.
 *
30
31
 * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
 * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
32
33
34
 */

#include "adiabatic_index.h"
35
#include "approx_math.h"
36
#include "cosmology.h"
37
#include "dimension.h"
Matthieu Schaller's avatar
Matthieu Schaller committed
38
#include "entropy_floor.h"
39
40
#include "equation_of_state.h"
#include "hydro_properties.h"
41
#include "hydro_space.h"
42
#include "kernel_hydro.h"
43
#include "minmax.h"
44

45
46
#include "./hydro_parameters.h"

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

58
  return gas_internal_energy_from_entropy(p->rho_bar, xp->entropy_full);
59
60
61
}

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

  return gas_internal_energy_from_entropy(p->rho_bar * cosmo->a3_inv,
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
                                          xp->entropy_full);
}
/**
 * @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_bar, p->entropy);
}

/**
 * @brief Returns the physical internal energy 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_internal_energy(const struct part *restrict p,
                                           const struct cosmology *cosmo) {

100
101
  return gas_internal_energy_from_entropy(p->rho_bar * cosmo->a3_inv,
                                          p->entropy);
102
103
104
105
106
107
108
109
}

/**
 * @brief Returns the comoving pressure of a particle
 *
 * @param p The particle of interest
 */
__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
110
    const struct part *restrict p) {
111

112
  return gas_pressure_from_entropy(p->rho_bar, p->entropy);
113
114
115
}

/**
116
117
118
119
120
121
122
123
124
125
126
 * @brief Returns the physical pressure of a particle
 *
 * @param p The particle of interest
 */
__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_bar * cosmo->a3_inv, p->entropy);
}

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

136
  return xp->entropy_full;
137
138
139
}

/**
140
141
 * @brief Returns the physical entropy of a particl at the last
 * time the particle was kicked.
142
143
 *
 * @param p The particle of interest.
144
 * @param xp The extended data of the particle of interest.
145
146
147
 * @param cosmo The cosmological model.
 */
__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
148
149
150
151
152
153
154
155
    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;
}

156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
/**
 * @brief Sets the physical internal energy of a particle
 *
 * @param p The particle of interest.
 * @param xp The extended particle data.
 * @param cosmo Cosmology data structure
 * @param u The physical entropy
 */
__attribute__((always_inline)) INLINE static void
hydro_set_physical_internal_energy(struct part *p, struct xpart *xp,
                                   const struct cosmology *cosmo,
                                   const float u) {
  error("To be implemented");
}

/**
 * @brief Sets the drifted physical internal energy of a particle
 *
 * @param p The particle of interest.
 * @param xp The extended particle data.
 * @param cosmo Cosmology data structure
 * @param u The physical entropy
 */
__attribute__((always_inline)) INLINE static void
hydro_set_drifted_physical_internal_energy(struct part *p,
                                           const struct cosmology *cosmo,
                                           const float u) {
  error("To be implemented");
}

186
187
188
189
190
191
192
193
194
195
196
197
/**
 * @brief Update the value of the viscosity alpha for the scheme.
 *
 * @param p the particle of interest
 * @param alpha the new value for the viscosity coefficient.
 */
__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
    struct part *restrict p, float alpha) {
  /* This scheme has fixed alpha */
}

/**
198
 * @brief Update the value of the diffusive coefficients to the
199
200
201
202
203
 *        feedback reset value for the scheme.
 *
 * @param p the particle of interest
 */
__attribute__((always_inline)) INLINE static void
204
hydro_diffusive_feedback_reset(struct part *restrict p) {
205
206
207
  /* This scheme has fixed alpha */
}

208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
/**
 * @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) {
230
231
232
233
234
235
236
237

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

/**
 * @brief Returns the comoving sound speed of a particle
238
239
240
 *
 * @param p The particle of interest
 */
241
242
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_soundspeed(const struct part *restrict p) {
243
244
245
246

  return p->force.soundspeed;
}

247
/**
248
 * @brief Returns the physical sound speed of a particle
249
250
 *
 * @param p The particle of interest
251
 * @param cosmo The cosmological model.
252
 */
253
254
255
256
257
258
259
260
261
262
263
264
265
__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
 *
 * @param p The particle of interest
 */
__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
266
267
268
269
270
    const struct part *restrict p) {

  return p->rho;
}

271
272
273
274
275
276
277
278
279
280
281
/**
 * @brief Returns the physical density of a particle
 *
 * @param p The particle of interest
 */
__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;
}

282
283
284
285
286
287
288
289
290
291
292
/**
 * @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;
}

293
294
295
296
297
298
299
300
301
302
303
304
/**
 * @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;
}

305
306
307
308
309
/**
 * @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.
310
311
 * @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.
312
313
314
 * @param v (return) The velocities at the current time.
 */
__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
315
316
317
318
319
320
321
322
323
    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;
324
325
}

326
/**
327
 * @brief Returns the time derivative of internal energy of a particle
328
 *
329
 * We assume a constant density.
330
 *
331
 * @param p The particle of interest
332
 */
333
334
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_internal_energy_dt(const struct part *restrict p) {
335

336
  return gas_internal_energy_from_entropy(p->rho_bar, p->entropy_dt);
337
338
}

339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
/**
 * @brief Returns the time derivative of physical internal energy of a particle
 *
 * We assume a constant density.
 *
 * @param p The particle of interest.
 * @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_bar * cosmo->a3_inv,
                                          p->entropy_dt);
}

355
/**
356
 * @brief Returns the time derivative of internal energy of a particle
357
 *
358
 * We assume a constant density.
359
 *
360
361
 * @param p The particle of interest.
 * @param du_dt The new time derivative of the internal energy.
362
 */
363
364
__attribute__((always_inline)) INLINE static void
hydro_set_comoving_internal_energy_dt(struct part *restrict p, float du_dt) {
365

366
  p->entropy_dt = gas_entropy_from_internal_energy(p->rho_bar, du_dt);
367
368
}

369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
/**
 * @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_bar * cosmo->a3_inv, du_dt);
}
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
/**
 * @brief Sets the physical entropy of a particle
 *
 * @param p The particle of interest.
 * @param xp The extended particle data.
 * @param cosmo Cosmology data structure
 * @param entropy The physical entropy
 */
__attribute__((always_inline)) INLINE static void hydro_set_physical_entropy(
    struct part *p, struct xpart *xp, const struct cosmology *cosmo,
    const float entropy) {

  /* Note there is no conversion from physical to comoving entropy */
  xp->entropy_full = entropy;
}
400

401
402
403
/**
 * @brief Computes the hydro time-step of a given particle
 *
404
405
406
407
 * @param p Pointer to the particle data.
 * @param xp Pointer to the extended particle data.
 * @param hydro_properties The constants used in the scheme.
 * @param cosmo The current cosmological model.
408
409
410
 */
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
    const struct part *restrict p, const struct xpart *restrict xp,
411
412
    const struct hydro_props *restrict hydro_properties,
    const struct cosmology *restrict cosmo) {
413

414
  const float CFL = hydro_properties->CFL_condition;
415
416

  /* CFL condition */
417
418
  const float dt_cfl = 2.f * kernel_gamma * CFL * cosmo->a * p->h /
                       (cosmo->a_factor_sound_speed * p->force.v_sig);
419
420
421
422

  return dt_cfl;
}

423
424
425
426
427
428
429
430
431
432
/**
 * @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) {}

433
434
435
436
437
438
439
/**
 * @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
440
 * @param hs #hydro_space containing hydro specific space information.
441
442
 */
__attribute__((always_inline)) INLINE static void hydro_init_part(
443
    struct part *restrict p, const struct hydro_space *hs) {
444

445
  p->rho = 0.f;
446
  p->rho_bar = 0.f;
447
448
  p->density.wcount = 0.f;
  p->density.wcount_dh = 0.f;
449
450
  p->density.rho_dh = 0.f;
  p->density.pressure_dh = 0.f;
451

452
453
454
455
  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;
456
457
458
459
460
461
462
463
464
}

/**
 * @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
465
 * @param cosmo The current cosmological model.
466
467
 */
__attribute__((always_inline)) INLINE static void hydro_end_density(
468
    struct part *restrict p, const struct cosmology *cosmo) {
469
470
471

  /* Some smoothing length multiples. */
  const float h = p->h;
472
473
474
  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) */
475
476
477

  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
478
  p->rho_bar += p->mass * p->entropy_one_over_gamma * kernel_root;
479
480
  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
  p->density.pressure_dh -=
481
      hydro_dimension * p->mass * p->entropy_one_over_gamma * kernel_root;
482
  p->density.wcount += kernel_root;
483
  p->density.wcount_dh -= hydro_dimension * kernel_root;
484
485
486

  /* Finish the calculation by inserting the missing h-factors */
  p->rho *= h_inv_dim;
487
  p->rho_bar *= h_inv_dim;
488
489
  p->density.rho_dh *= h_inv_dim_plus_one;
  p->density.pressure_dh *= h_inv_dim_plus_one;
490
491
  p->density.wcount *= h_inv_dim;
  p->density.wcount_dh *= h_inv_dim_plus_one;
492

493
  const float rho_inv = 1.f / p->rho;
494
  const float a_inv2 = cosmo->a2_inv;
495
  const float entropy_minus_one_over_gamma = 1.f / p->entropy_one_over_gamma;
496

497
498
  /* Final operation on the weighted density */
  p->rho_bar *= entropy_minus_one_over_gamma;
499
500

  /* Finish calculation of the velocity curl components */
501
502
503
  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;
504
505

  /* Finish calculation of the velocity divergence */
506
  p->density.div_v *= h_inv_dim_plus_one * rho_inv;
507
508
}

509
510
511
512
513
/**
 * @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
514
 * @param cosmo The current cosmological model.
515
516
 */
__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
517
518
    struct part *restrict p, struct xpart *restrict xp,
    const struct cosmology *cosmo) {
519
520
521
522
523
524
525
526
527

  /* 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->rho_bar = p->mass * kernel_root * h_inv_dim;
528
  p->density.wcount = kernel_root * h_inv_dim;
529
530
531
532
533
534
535
536
537
  p->density.rho_dh = 0.f;
  p->density.wcount_dh = 0.f;
  p->density.pressure_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;
}

538
539
540
/**
 * @brief Prepare a particle for the force calculation.
 *
541
542
543
544
545
546
 * 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.
547
548
549
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
550
 * @param cosmo The current cosmological model.
Matthieu Schaller's avatar
Matthieu Schaller committed
551
 * @param hydro_props Hydrodynamic properties.
552
553
 * @param dt_alpha The time-step used to evolve non-cosmological quantities such
 *                 as the artificial viscosity.
554
555
 */
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
556
    struct part *restrict p, struct xpart *restrict xp,
Josh Borrow's avatar
Josh Borrow committed
557
558
    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
    const float dt_alpha) {
559

560
  const float fac_mu = cosmo->a_factor_mu;
561
562
563
564
565
566
567
568
569

  /* Compute the norm of the curl */
  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]);

  /* Compute the norm of div v */
  const float abs_div_v = fabsf(p->density.div_v);

570
  /* Compute the pressure */
571
  const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
572

573
  /* Compute the sound speed from the pressure*/
574
575
  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);

576
577
  /* Compute the Balsara switch */
  const float balsara =
Josh Borrow's avatar
Josh Borrow committed
578
579
      hydro_props->viscosity.alpha * abs_div_v /
      (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu / p->h);
580

Matthieu Schaller's avatar
Matthieu Schaller committed
581
582
583
584
585
  /* Divide the pressure by the density squared to get the SPH term */
  const float rho_bar_inv = 1.f / p->rho_bar;
  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;

  /* Compute "grad h" term (note we use rho here and not rho_bar !)*/
586
587
588
589
590
  float rho_dh = p->density.rho_dh;
  /* Ignore changing-kernel effects when h ~= h_max */
  if (p->h > 0.9999f * hydro_props->h_max) {
    rho_dh = 0.f;
  }
591
  const float rho_inv = 1.f / p->rho;
592
  const float rho_dh_term =
593
      1.f / (1.f + hydro_dimension_inv * p->h * rho_dh * rho_inv);
594
  const float pressure_dh_term =
595
      p->density.pressure_dh * rho_inv * p->h * hydro_dimension_inv;
596

597
  const float grad_h_term = rho_dh_term * pressure_dh_term;
598

599
600
  /* Update variables. */
  p->force.soundspeed = soundspeed;
601
  p->force.P_over_rho2 = P_over_rho2;
602
603
  p->force.balsara = balsara;
  p->force.f = grad_h_term;
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
}

/**
 * @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 *restrict p) {

  /* 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. */
  p->entropy_dt = 0.0f;
  p->force.h_dt = 0.0f;

  /* Reset maximal signal velocity */
627
  p->force.v_sig = p->force.soundspeed;
628
629
}

630
631
632
633
634
635
/**
 * @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.
Loic Hausammann's avatar
Loic Hausammann committed
636
 * @param cosmo The cosmological model.
637
638
 */
__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
Loic Hausammann's avatar
Loic Hausammann committed
639
640
    struct part *restrict p, const struct xpart *restrict xp,
    const struct cosmology *cosmo) {
641
642
643
644
645
646
647
648

  /* 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->entropy = xp->entropy_full;
Matthieu Schaller's avatar
Matthieu Schaller committed
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663

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

  /* 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 */
  const float rho_inv = 1.f / p->rho;
  const float P_over_rho2 = pressure * rho_inv * rho_inv;

  /* Update variables */
  p->force.soundspeed = soundspeed;
  p->force.P_over_rho2 = P_over_rho2;
  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
664
665
}

666
667
668
669
670
/**
 * @brief Predict additional particle fields forward in time when drifting
 *
 * @param p The particle
 * @param xp The extended data of the particle
671
672
 * @param dt_drift The drift time-step for positions.
 * @param dt_therm The drift time-step for thermal quantities.
Matthieu Schaller's avatar
Matthieu Schaller committed
673
674
675
 * @param cosmo The cosmological model.
 * @param hydro_props The properties of the hydro scheme.
 * @param floor_props The properties of the entropy floor.
676
677
 */
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
678
    struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
Matthieu Schaller's avatar
Matthieu Schaller committed
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
    float dt_therm, const struct cosmology *cosmo,
    const struct hydro_props *hydro_props,
    const struct entropy_floor_properties *floor_props) {

  /* Predict the entropy */
  p->entropy += p->entropy_dt * dt_therm;

  /* Check against entropy floor */
  const float floor_A = entropy_floor(p, cosmo, floor_props);

  /* Check against absolute minimum */
  const float min_u = hydro_props->minimal_internal_energy;
  const float min_A =
      gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, min_u);

  p->entropy = max(p->entropy, floor_A);
  p->entropy = max(p->entropy, min_A);
696

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

699
  /* Predict smoothing length */
700
  const float w1 = p->force.h_dt * h_inv * dt_drift;
701
  if (fabsf(w1) < 0.2f)
702
    p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
703
  else
704
705
706
707
708
709
710
711
712
713
714
715
716
    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) */
    p->rho_bar *= approx_expf(w2);
  } else {
    p->rho *= expf(w2);
    p->rho_bar *= expf(w2);
  }

  /* Compute the pressure */
717
  const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
718

719
720
721
722
  /* Compute the new sound speed */
  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);

  /* Divide the pressure by the density squared to get the SPH term */
723
  const float rho_bar_inv = 1.f / p->rho_bar;
724
  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
725

726
  /* Update the variables */
727
  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
728
  p->force.soundspeed = soundspeed;
729
  p->force.P_over_rho2 = P_over_rho2;
730
731
732
733
734
735
736
737
}

/**
 * @brief Finishes the force calculation.
 *
 * Multiplies the forces and accelerationsby the appropiate constants
 *
 * @param p The particle to act upon
738
 * @param cosmo The current cosmological model.
739
740
 */
__attribute__((always_inline)) INLINE static void hydro_end_force(
741
    struct part *restrict p, const struct cosmology *cosmo) {
742
743
744

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

745
746
  p->entropy_dt = 0.5f * cosmo->a2_inv *
                  gas_entropy_from_internal_energy(p->rho_bar, p->entropy_dt);
747
748
749
750
751
752
753
}

/**
 * @brief Kick the additional variables
 *
 * @param p The particle to act upon
 * @param xp The particle extended data to act upon
754
 * @param dt_therm The time-step for this kick (for thermodynamic quantities)
755
756
 * @param cosmo The cosmological model.
 * @param hydro_props The constants used in the scheme
Matthieu Schaller's avatar
Matthieu Schaller committed
757
 * @param floor_props The properties of the entropy floor.
758
759
 */
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
760
    struct part *restrict p, struct xpart *restrict xp, float dt_therm,
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
761
    float dt_grav, float dt_hydro, float dt_kick_corr,
Matthieu Schaller's avatar
Matthieu Schaller committed
762
763
    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
    const struct entropy_floor_properties *floor_props) {
764

Matthieu Schaller's avatar
Matthieu Schaller committed
765
766
  /* Integrate the entropy forward in time */
  const float delta_entropy = p->entropy_dt * dt_therm;
767

Matthieu Schaller's avatar
Matthieu Schaller committed
768
769
770
  /* Do not decrease the entropy by more than a factor of 2 */
  xp->entropy_full =
      max(xp->entropy_full + delta_entropy, 0.5f * xp->entropy_full);
771

Matthieu Schaller's avatar
Matthieu Schaller committed
772
773
  /* Check against entropy floor */
  const float floor_A = entropy_floor(p, cosmo, floor_props);
774

Matthieu Schaller's avatar
Matthieu Schaller committed
775
776
777
778
  /* Check against absolute minimum */
  const float min_u = hydro_props->minimal_internal_energy;
  const float min_A =
      gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, min_u);
779

Matthieu Schaller's avatar
Matthieu Schaller committed
780
781
782
783
784
785
786
  /* Take highest of both limits */
  const float entropy_min = max(min_A, floor_A);

  if (xp->entropy_full < entropy_min) {
    xp->entropy_full = entropy_min;
    p->entropy_dt = 0.f;
  }
787
788
789
790
791
792
793
794
}

/**
 *  @brief Converts hydro quantity of a particle at the start of a run
 *
 * Requires the density to be known
 *
 * @param p The particle to act upon
795
796
 * @param xp The extended data.
 * @param cosmo The cosmological model.
797
798
 */
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
799
    struct part *restrict p, struct xpart *restrict xp,
800
    const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
801
802

  /* We read u in the entropy field. We now get S from u */
803
804
  xp->entropy_full =
      gas_entropy_from_internal_energy(p->rho_bar * cosmo->a3_inv, p->entropy);
805
  p->entropy = xp->entropy_full;
806
807
808
  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);

  /* Compute the pressure */
809
  const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
810

811
812
813
814
  /* Compute the sound speed */
  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);

  /* Divide the pressure by the density squared to get the SPH term */
815
  const float rho_bar_inv = 1.f / p->rho_bar;
816
817
  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;

818
  p->force.soundspeed = soundspeed;
819
  p->force.P_over_rho2 = P_over_rho2;
820
821
}

822
823
824
825
826
827
828
829
830
831
832
833
834
/**
 * @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;
835
  p->wakeup = time_bin_not_awake;
836
837
838
839
840
  p->rho_bar = 0.f;
  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
  xp->v_full[0] = p->v[0];
  xp->v_full[1] = p->v[1];
  xp->v_full[2] = p->v[2];
841
842
843
  xp->a_grav[0] = 0.f;
  xp->a_grav[1] = 0.f;
  xp->a_grav[2] = 0.f;
844
845

  hydro_reset_acceleration(p);
846
  hydro_init_part(p, NULL);
847
848
}

849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
/**
 * @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;
}

866
867
868
869
870
871
872
873
874
875
/**
 * @brief Operations performed when a particle gets removed from the
 * simulation volume.
 *
 * @param p The particle.
 * @param xp The extended particle data.
 */
__attribute__((always_inline)) INLINE static void hydro_remove_part(
    const struct part *p, const struct xpart *xp) {}

876
#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */