hydro.h 16.3 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
37
38
39
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "kernel_hydro.h"
40
#include "minmax.h"
41
42
43
44
45
46
47
48
49
50
51
52

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

  const float entropy = p->entropy + p->entropy_dt * dt;

53
  return gas_internal_energy_from_entropy(p->rho_bar, entropy);
54
55
56
57
58
59
60
61
62
63
64
65
66
}

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

  const float entropy = p->entropy + p->entropy_dt * dt;

67
  return gas_pressure_from_entropy(p->rho_bar, entropy);
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
}

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

  return p->entropy + p->entropy_dt * dt;
}

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

  return p->force.soundspeed;
}

94
/**
95
 * @brief Returns the physical density of a particle
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
 *
 * @param p The particle of interest
 */
__attribute__((always_inline)) INLINE static float hydro_get_density(
    const struct part *restrict p) {

  return p->rho;
}

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

  return p->mass;
}

116
117
118
119
/**
 * @brief Modifies the thermal state of a particle to the imposed internal
 * energy
 *
120
121
122
 * This overwrites the current state of the particle but does *not* change its
 * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
 * updated.
123
124
125
126
127
128
129
 *
 * @param p The particle
 * @param u The new internal energy
 */
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
    struct part *restrict p, float u) {

130
  p->entropy = gas_entropy_from_internal_energy(p->rho_bar, u);
131
132
133
  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);

  /* Compute the pressure */
134
  const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
135
136

  /* Compute the sound speed from the pressure*/
137
  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
138
139
140
141
142
143

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

144
  const float rho_bar_inv = 1.f / p->rho_bar;
145

146
147
  p->force.soundspeed = soundspeed;
  p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
148
  p->force.v_sig = v_sig;
149
150
151
152
153
}

/**
 * @brief Modifies the thermal state of a particle to the imposed entropy
 *
154
155
156
 * This overwrites the current state of the particle but does *not* change its
 * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
 * updated.
157
158
159
160
161
162
163
164
 *
 * @param p The particle
 * @param S The new entropy
 */
__attribute__((always_inline)) INLINE static void hydro_set_entropy(
    struct part *restrict p, float S) {

  p->entropy = S;
165
166
167
  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);

  /* Compute the pressure */
168
  const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
169
170

  /* Compute the sound speed from the pressure*/
171
  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
172
173
174
175
176
177

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

178
  const float rho_bar_inv = 1.f / p->rho_bar;
179

180
181
  p->force.soundspeed = soundspeed;
  p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
182
  p->force.v_sig = v_sig;
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
}

/**
 * @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
 *
 */
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
    const struct part *restrict p, const struct xpart *restrict xp,
    const struct hydro_props *restrict hydro_properties) {

  const float CFL_condition = hydro_properties->CFL_condition;

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

  return dt_cfl;
}

/**
 * @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->ti_begin = 0;
  p->ti_end = 0;
219
220
  p->rho_bar = 0.f;
  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
  xp->v_full[0] = p->v[0];
  xp->v_full[1] = p->v[1];
  xp->v_full[2] = p->v[2];
}

/**
 * @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
 */
__attribute__((always_inline)) INLINE static void hydro_init_part(
    struct part *restrict p) {
236

237
  p->rho = 0.f;
238
  p->rho_bar = 0.f;
239
240
  p->density.wcount = 0.f;
  p->density.wcount_dh = 0.f;
241
242
  p->density.rho_dh = 0.f;
  p->density.pressure_dh = 0.f;
243

244
245
246
247
  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;
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
}

/**
 * @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
 * @param ti_current The current time (on the integer timeline)
 */
__attribute__((always_inline)) INLINE static void hydro_end_density(
    struct part *restrict p, int ti_current) {

  /* Some smoothing length multiples. */
  const float h = p->h;
264
265
266
  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) */
267
268
269

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

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

284
285
  const float rho_inv = 1.f / p->rho;
  const float entropy_minus_one_over_gamma = 1.f / p->entropy_one_over_gamma;
286

287
288
  /* Final operation on the weighted density */
  p->rho_bar *= entropy_minus_one_over_gamma;
289
290

  /* Finish calculation of the velocity curl components */
291
292
293
  p->density.rot_v[0] *= h_inv_dim_plus_one * rho_inv;
  p->density.rot_v[1] *= h_inv_dim_plus_one * rho_inv;
  p->density.rot_v[2] *= h_inv_dim_plus_one * rho_inv;
294
295

  /* Finish calculation of the velocity divergence */
296
  p->density.div_v *= h_inv_dim_plus_one * rho_inv;
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
}

/**
 * @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
 * @param ti_current The current time (on the timeline)
 * @param timeBase The minimal time-step size
 */
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
    struct part *restrict p, struct xpart *restrict xp, int ti_current,
    double timeBase) {

313
314
315
316
317
318
319
320
321
322
  const float fac_mu = 1.f; /* Will change with cosmological integration */

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

323
324
  /* Compute the pressure */
  const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
Matthieu Schaller's avatar
Matthieu Schaller committed
325
326
  const float entropy = hydro_get_entropy(p, half_dt);
  const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
327

328
  /* Compute the sound speed from the pressure*/
329
330
  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);

331
332
333
334
  /* Compute the Balsara switch */
  const float balsara =
      abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);

Matthieu Schaller's avatar
Matthieu Schaller committed
335
336
337
338
339
  /* 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 !)*/
340
341
342
343
344
345
346
347
  const float rho_inv = 1.f / p->rho;
  const float entropy_minus_one_over_gamma = 1.f / p->entropy_one_over_gamma;
  const float rho_dh =
      1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
  const float pressure_dh = p->density.pressure_dh * rho_inv * p->h *
                            hydro_dimension_inv * entropy_minus_one_over_gamma;

  const float grad_h_term = rho_dh * pressure_dh;
348

349
350
  /* Update variables. */
  p->force.soundspeed = soundspeed;
351
  p->force.P_over_rho2 = P_over_rho2;
352
353
  p->force.balsara = balsara;
  p->force.f = grad_h_term;
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
}

/**
 * @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 */
  p->force.v_sig = 0.0f;
}

/**
 * @brief Predict additional particle fields forward in time when drifting
 *
 * @param p The particle
 * @param xp The extended data of the particle
385
386
387
 * @param dt The drift time-step.
 * @param t0 The time at the start of the drift (on the timeline).
 * @param t1 The time at the end of the drift (on the timeline).
388
389
390
 * @param timeBase The minimal time-step size
 */
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
391
392
    struct part *restrict p, const struct xpart *restrict xp, float dt, int t0,
    int t1, double timeBase) {
393

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

396
397
  /* Predict smoothing length */
  const float w1 = p->force.h_dt * h_inv * dt;
398
  if (fabsf(w1) < 0.2f)
399
    p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
400
  else
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
    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);
  }

  /* Drift the entropy */
  const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
  const float entropy = hydro_get_entropy(p, dt_entr);

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

420
421
422
423
  /* 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 */
424
  const float rho_bar_inv = 1.f / p->rho_bar;
425
  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
426

427
428
  /* Update the variables */
  p->entropy_one_over_gamma = pow_one_over_gamma(entropy);
429
  p->force.soundspeed = soundspeed;
430
  p->force.P_over_rho2 = P_over_rho2;
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
}

/**
 * @brief Finishes the force calculation.
 *
 * Multiplies the forces and accelerationsby the appropiate constants
 *
 * @param p The particle to act upon
 */
__attribute__((always_inline)) INLINE static void hydro_end_force(
    struct part *restrict p) {

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

  p->entropy_dt *=
446
      0.5f * hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho_bar);
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
}

/**
 * @brief Kick the additional variables
 *
 * @param p The particle to act upon
 * @param xp The particle extended data to act upon
 * @param dt The time-step for this kick
 * @param half_dt The half time-step for this kick
 */
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
    struct part *restrict p, struct xpart *restrict xp, float dt,
    float half_dt) {

  /* Do not decrease the entropy (temperature) by more than a factor of 2*/
  const float entropy_change = p->entropy_dt * dt;
  if (entropy_change > -0.5f * p->entropy)
    p->entropy += entropy_change;
  else
    p->entropy *= 0.5f;

  /* Do not 'overcool' when timestep increases */
  if (p->entropy + p->entropy_dt * half_dt < 0.5f * p->entropy)
    p->entropy_dt = -0.5f * p->entropy / half_dt;
471
472
473
474

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

475
476
477
478
  /* 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 */
479
  const float rho_bar_inv = 1.f / p->rho_bar;
480
  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
481
482
483

  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
  p->force.soundspeed = soundspeed;
484
  p->force.P_over_rho2 = P_over_rho2;
485
486
487
488
489
490
491
492
493
494
495
496
497
}

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

  /* We read u in the entropy field. We now get S from u */
498
  p->entropy = gas_entropy_from_internal_energy(p->rho_bar, p->entropy);
499
500
501
  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);

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

504
505
506
507
  /* 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 */
508
  const float rho_bar_inv = 1.f / p->rho_bar;
509
510
  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;

511
  p->force.soundspeed = soundspeed;
512
  p->force.P_over_rho2 = P_over_rho2;
513
514
515
}

#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */