hydro.h 19.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
 * @brief Returns the comoving internal energy of a particle
46
47
48
 *
 * @param p The particle of interest
 */
49
50
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_internal_energy(const struct part *restrict p) {
51

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

55
/**
56
 * @brief Returns the physical internal energy of a particle
57
58
59
60
61
62
63
64
65
66
67
 *
 * @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 gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
}

68
/**
69
 * @brief Returns the comoving pressure of a particle
70
71
72
 *
 * @param p The particle of interest
 */
73
__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
74
    const struct part *restrict p) {
75

76
  return gas_pressure_from_entropy(p->rho, p->entropy);
77
78
}

79
80
81
82
83
84
85
86
87
88
89
/**
 * @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 * cosmo->a3_inv, p->entropy);
}

90
91
92
93
94
95
96
97
98
99
100
/**
 * @brief Returns the comoving entropy of a particle
 *
 * @param p The particle of interest.
 */
__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
    const struct part *restrict p) {

  return p->entropy;
}

101
/**
102
 * @brief Returns the physical entropy of a particle
103
 *
104
105
 * @param p The particle of interest.
 * @param cosmo The cosmological model.
106
 */
107
__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
108
    const struct part *restrict p, const struct cosmology *cosmo) {
109

110
111
  /* Note: no cosmological conversion required here with our choice of
   * coordinates. */
112
  return p->entropy;
113
114
115
}

/**
116
 * @brief Returns the comoving sound speed of a particle
117
118
119
 *
 * @param p The particle of interest
 */
120
121
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_soundspeed(const struct part *restrict p) {
122
123
124

  return p->force.soundspeed;
}
125

126
127
128
129
130
131
132
133
134
135
136
137
138
/**
 * @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;
}

139
/**
140
 * @brief Returns the comoving density of a particle
141
142
143
 *
 * @param p The particle of interest
 */
144
__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
145
146
147
148
149
    const struct part *restrict p) {

  return p->rho;
}

150
151
152
153
154
155
156
157
158
159
160
161
/**
 * @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;
}

162
163
164
165
166
167
168
169
170
171
172
/**
 * @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;
}

173
174
175
176
177
178
179
180
181
182
183
184
/**
 * @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;
}

185
186
187
188
189
/**
 * @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.
190
191
 * @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.
192
193
194
 * @param v (return) The velocities at the current time.
 */
__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
195
196
197
198
199
200
201
202
203
    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;
204
205
}

206
/**
207
 * @brief Returns the time derivative of internal energy of a particle
208
 *
209
 * We assume a constant density.
210
 *
211
 * @param p The particle of interest
212
 */
213
214
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
    const struct part *restrict p) {
215

216
  return gas_internal_energy_from_entropy(p->rho, p->entropy_dt);
217
218
219
}

/**
220
 * @brief Returns the time derivative of internal energy of a particle
221
 *
222
 * We assume a constant density.
223
 *
224
225
 * @param p The particle of interest.
 * @param du_dt The new time derivative of the internal energy.
226
 */
227
228
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
    struct part *restrict p, float du_dt) {
229

230
  p->entropy_dt = gas_entropy_from_internal_energy(p->rho, du_dt);
231
232
}

233
234
235
236
237
/**
 * @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
238
 * @param hydro_properties The constants used in the scheme
239
 * @param cosmo The cosmological model.
240
241
 */
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
242
    const struct part *restrict p, const struct xpart *restrict xp,
243
244
    const struct hydro_props *restrict hydro_properties,
    const struct cosmology *restrict cosmo) {
245

246
  const float CFL = hydro_properties->CFL_condition;
247

248
  /* CFL condition */
249
  const float dt_cfl = 2.f * kernel_gamma * CFL * cosmo->a * p->h /
250
                       (cosmo->a_factor_sound_speed * p->force.v_sig);
251

252
  return dt_cfl;
253
254
}

255
256
257
258
259
260
261
262
263
264
/**
 * @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) {}

265
266
267
/**
 * @brief Prepares a particle for the density calculation.
 *
268
 * Zeroes all the relevant arrays in preparation for the sums taking place in
269
270
271
 * the variaous density tasks
 *
 * @param p The particle to act upon
272
 * @param hs #hydro_space containing hydro specific space information.
273
 */
274
__attribute__((always_inline)) INLINE static void hydro_init_part(
275
    struct part *restrict p, const struct hydro_space *hs) {
276

277
#ifdef DEBUG_INTERACTIONS_SPH
278
  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS; ++i) p->ids_ngbs_density[i] = -1;
279
280
  p->num_ngb_density = 0;
#endif
281

282
  p->rho = 0.f;
283
284
  p->density.wcount = 0.f;
  p->density.wcount_dh = 0.f;
285
  p->density.rho_dh = 0.f;
286
  p->density.div_v = 0.f;
287
288
289
  p->density.rot_v[0] = 0.f;
  p->density.rot_v[1] = 0.f;
  p->density.rot_v[2] = 0.f;
290
291
292
}

/**
293
 * @brief Finishes the density calculation.
294
 *
295
 * Multiplies the density and number of neighbours by the appropiate constants
296
297
298
 * and add the self-contribution term.
 *
 * @param p The particle to act upon
299
 * @param cosmo The current cosmological model.
300
 */
301
__attribute__((always_inline)) INLINE static void hydro_end_density(
302
    struct part *restrict p, const struct cosmology *cosmo) {
303
304
305

  /* Some smoothing length multiples. */
  const float h = p->h;
306
307
308
  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) */
309

310
311
  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
312
  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
313
  p->density.wcount += kernel_root;
314
  p->density.wcount_dh -= hydro_dimension * kernel_root;
315

316
  /* Finish the calculation by inserting the missing h-factors */
317
  p->rho *= h_inv_dim;
318
  p->density.rho_dh *= h_inv_dim_plus_one;
319
320
  p->density.wcount *= h_inv_dim;
  p->density.wcount_dh *= h_inv_dim_plus_one;
321

322
  const float rho_inv = 1.f / p->rho;
323
  const float a_inv2 = cosmo->a2_inv;
Matthieu Schaller's avatar
Matthieu Schaller committed
324

325
326
327
328
  /* 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;
329

330
331
  /* Finish calculation of the (physical) velocity divergence */
  p->density.div_v *= h_inv_dim_plus_one * a_inv2 * rho_inv;
332
333
}

334
335
336
337
338
/**
 * @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
339
 * @param cosmo The current cosmological model.
340
341
 */
__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
342
343
    struct part *restrict p, struct xpart *restrict xp,
    const struct cosmology *cosmo) {
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360

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

361
362
363
364
365
366
367
/**
 * @brief Prepare a particle for the force calculation.
 *
 * Computes viscosity term, conduction term and smoothing length gradient terms.
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
368
 * @param cosmo The current cosmological model.
369
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
370
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
371
372
    struct part *restrict p, struct xpart *restrict xp,
    const struct cosmology *cosmo) {
373

374
  const float fac_mu = cosmo->a_factor_mu;
375

376
377
378
  /* Inverse of the physical density */
  const float rho_inv = 1.f / p->rho;

379
  /* Compute the norm of the curl */
380
381
382
  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]);
383

Matthieu Schaller's avatar
Matthieu Schaller committed
384
385
386
  /* Compute the norm of div v */
  const float abs_div_v = fabsf(p->density.div_v);

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

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

393
  /* Divide the pressure by the density squared to get the SPH term */
394
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
395

396
  /* Compute the Balsara switch */
Matthieu Schaller's avatar
Matthieu Schaller committed
397
  const float balsara =
398
      abs_div_v / (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
Matthieu Schaller's avatar
Matthieu Schaller committed
399

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

404
  /* Update variables. */
405
  p->force.f = omega_inv;
Matthieu Schaller's avatar
Matthieu Schaller committed
406
  p->force.P_over_rho2 = P_over_rho2;
407
408
  p->force.soundspeed = soundspeed;
  p->force.balsara = balsara;
409
410
411
412
413
414
415
416
417
418
}

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

422
#ifdef DEBUG_INTERACTIONS_SPH
423
  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS; ++i) p->ids_ngbs_force[i] = -1;
424
425
  p->num_ngb_force = 0;
#endif
426

427
  /* Reset the acceleration. */
Matthieu Schaller's avatar
Matthieu Schaller committed
428
429
430
  p->a_hydro[0] = 0.0f;
  p->a_hydro[1] = 0.0f;
  p->a_hydro[2] = 0.0f;
431

432
  /* Reset the time derivatives. */
433
  p->entropy_dt = 0.0f;
434
  p->force.h_dt = 0.0f;
435

436
  /* Reset maximal signal velocity */
437
  p->force.v_sig = p->force.soundspeed;
438
439
}

440
441
442
443
444
445
446
/**
 * @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.
 */
447
448
449
450
451
452
453
__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];
454
455
456

  /* Re-set the entropy */
  p->entropy = xp->entropy_full;
457
458
}

459
460
461
462
463
/**
 * @brief Predict additional particle fields forward in time when drifting
 *
 * @param p The particle
 * @param xp The extended data of the particle
464
465
 * @param dt_drift The drift time-step for positions.
 * @param dt_therm The drift time-step for thermal quantities.
466
 */
467
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
468
469
    struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
    float dt_therm) {
470
471
472
473

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

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

487
  /* Predict the entropy */
488
  p->entropy += p->entropy_dt * dt_therm;
Matthieu Schaller's avatar
Matthieu Schaller committed
489

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

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

496
497
  /* Divide the pressure by the density squared to get the SPH term */
  const float rho_inv = 1.f / p->rho;
498
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
Matthieu Schaller's avatar
Matthieu Schaller committed
499

500
501
  /* Update variables */
  p->force.soundspeed = soundspeed;
502
  p->force.P_over_rho2 = P_over_rho2;
503
504
505
}

/**
506
 * @brief Finishes the force calculation.
507
 *
508
 * Multiplies the forces and accelerationsby the appropiate constants
509
510
 *
 * @param p The particle to act upon
511
 * @param cosmo The current cosmological model.
512
 */
513
__attribute__((always_inline)) INLINE static void hydro_end_force(
514
    struct part *restrict p, const struct cosmology *cosmo) {
515

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

518
519
  p->entropy_dt = 0.5f * cosmo->a2_inv *
                  gas_entropy_from_internal_energy(p->rho, p->entropy_dt);
520
521
}

Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
522
523
524
525
/**
 * @brief Kick the additional variables
 *
 * @param p The particle to act upon
526
 * @param xp The particle extended data to act upon
527
 * @param dt_therm The time-step for this kick (for thermodynamic quantities)
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
528
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
529
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
530
    struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
531

532
  /* Do not decrease the entropy by more than a factor of 2 */
533
  if (dt_therm > 0. && p->entropy_dt * dt_therm < -0.5f * xp->entropy_full) {
Matthieu Schaller's avatar
Matthieu Schaller committed
534
535
536
    /* message("Warning! Limiting entropy_dt. Possible cooling error.\n
     * entropy_full = %g \n entropy_dt * dt =%g \n", */
    /* 	     xp->entropy_full,p->entropy_dt * dt); */
537
    p->entropy_dt = -0.5f * xp->entropy_full / dt_therm;
Matthieu Schaller's avatar
Matthieu Schaller committed
538
  }
539
  xp->entropy_full += p->entropy_dt * dt_therm;
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
540

541
  /* Compute the pressure */
542
  const float pressure = gas_pressure_from_entropy(p->rho, xp->entropy_full);
543

544
545
546
547
  /* 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 */
548
  const float rho_inv = 1.f / p->rho;
549
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
550
551

  p->force.soundspeed = soundspeed;
552
  p->force.P_over_rho2 = P_over_rho2;
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
553
554
}

555
/**
556
 *  @brief Converts hydro quantity of a particle at the start of a run
557
558
559
 *
 * Requires the density to be known
 *
560
561
562
 * @param p The particle to act upon.
 * @param xp The extended data.
 * @param cosmo The cosmological model.
563
 */
564
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
565
566
    struct part *restrict p, struct xpart *restrict xp,
    const struct cosmology *cosmo) {
567

568
569
570
571
  /* We read u in the entropy field. We now get (comoving) S from (physical) u
   * and (physical) rho. Note that comoving S == physical S */
  xp->entropy_full =
      gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, p->entropy);
572
  p->entropy = xp->entropy_full;
Matthieu Schaller's avatar
Matthieu Schaller committed
573

574
575
576
  /* Compute the pressure */
  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);

577
578
579
580
  /* 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 */
581
  const float rho_inv = 1.f / p->rho;
582
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
583
584

  p->force.soundspeed = soundspeed;
585
  p->force.P_over_rho2 = P_over_rho2;
586
}
587

Matthieu Schaller's avatar
Matthieu Schaller committed
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
/**
 * @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;
  xp->v_full[0] = p->v[0];
  xp->v_full[1] = p->v[1];
  xp->v_full[2] = p->v[2];
604
605
606
  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
607
608
609
  xp->entropy_full = p->entropy;

  hydro_reset_acceleration(p);
610
  hydro_init_part(p, NULL);
Matthieu Schaller's avatar
Matthieu Schaller committed
611
612
}

613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
/**
 * @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;
}

630
#endif /* SWIFT_GADGET2_HYDRO_H */