hydro.h 14.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
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

/**
 * @brief Returns the internal energy of a particle
 *
 * @param p The particle of interest
 */
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
48
    const struct part *restrict p) {
49

50
  return gas_internal_energy_from_entropy(p->rho_bar, p->entropy);
51
52
53
54
55
56
57
58
}

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

61
  return gas_pressure_from_entropy(p->rho_bar, p->entropy);
62
63
64
65
66
67
68
69
}

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

72
  return p->entropy;
73
74
75
76
77
78
79
80
}

/**
 * @brief Returns the sound speed of a particle
 *
 * @param p The particle of interest
 */
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
81
    const struct part *restrict p) {
82
83
84
85

  return p->force.soundspeed;
}

86
/**
87
 * @brief Returns the physical density of a particle
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
 *
 * @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;
}

108
/**
109
 * @brief Returns the time derivative of internal energy of a particle
110
 *
111
 * We assume a constant density.
112
 *
113
 * @param p The particle of interest
114
 */
115
116
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
    const struct part *restrict p) {
117

118
  return gas_internal_energy_from_entropy(p->rho_bar, p->entropy_dt);
119
120
121
}

/**
122
 * @brief Returns the time derivative of internal energy of a particle
123
 *
124
 * We assume a constant density.
125
 *
126
127
 * @param p The particle of interest.
 * @param du_dt The new time derivative of the internal energy.
128
 */
129
130
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
    struct part *restrict p, float du_dt) {
131

132
  p->entropy_dt = gas_entropy_from_internal_energy(p->rho_bar, du_dt);
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
}

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

166
  p->rho = 0.f;
167
  p->rho_bar = 0.f;
168
169
  p->density.wcount = 0.f;
  p->density.wcount_dh = 0.f;
170
171
  p->density.rho_dh = 0.f;
  p->density.pressure_dh = 0.f;
172

173
174
175
176
  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;
177
178
179
180
181
182
183
184
185
186
187
}

/**
 * @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
 */
__attribute__((always_inline)) INLINE static void hydro_end_density(
188
    struct part *restrict p) {
189
190
191

  /* Some smoothing length multiples. */
  const float h = p->h;
192
193
194
  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) */
195
196
197

  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
198
  p->rho_bar += p->mass * p->entropy_one_over_gamma * kernel_root;
199
200
  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
  p->density.pressure_dh -=
201
      hydro_dimension * p->mass * p->entropy_one_over_gamma * kernel_root;
202
203
204
205
  p->density.wcount += kernel_root;

  /* Finish the calculation by inserting the missing h-factors */
  p->rho *= h_inv_dim;
206
  p->rho_bar *= h_inv_dim;
207
208
  p->density.rho_dh *= h_inv_dim_plus_one;
  p->density.pressure_dh *= h_inv_dim_plus_one;
209
210
211
  p->density.wcount *= kernel_norm;
  p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;

212
213
  const float rho_inv = 1.f / p->rho;
  const float entropy_minus_one_over_gamma = 1.f / p->entropy_one_over_gamma;
214

215
216
  /* Final operation on the weighted density */
  p->rho_bar *= entropy_minus_one_over_gamma;
217
218

  /* Finish calculation of the velocity curl components */
219
220
221
  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;
222
223

  /* Finish calculation of the velocity divergence */
224
  p->density.div_v *= h_inv_dim_plus_one * rho_inv;
225
226
227
228
229
230
231
232
233
234
235
}

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

238
239
240
241
242
243
244
245
246
247
  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);

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

251
  /* Compute the sound speed from the pressure*/
252
253
  const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);

254
255
256
257
  /* 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
258
259
260
261
262
  /* 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 !)*/
263
264
265
  const float rho_inv = 1.f / p->rho;
  const float rho_dh =
      1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
266
267
  const float pressure_dh =
      p->density.pressure_dh * rho_inv * p->h * hydro_dimension_inv;
268
269

  const float grad_h_term = rho_dh * pressure_dh;
270

271
272
  /* Update variables. */
  p->force.soundspeed = soundspeed;
273
  p->force.P_over_rho2 = P_over_rho2;
274
275
  p->force.balsara = balsara;
  p->force.f = grad_h_term;
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
}

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

302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
/**
 * @brief Sets the values to be predicted in the drifts to their values at a
 * kick time
 *
 * @param p The particle.
 * @param xp The extended data of this particle.
 */
__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
    struct part *restrict p, const struct xpart *restrict xp) {

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

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

321
322
323
324
325
/**
 * @brief Predict additional particle fields forward in time when drifting
 *
 * @param p The particle
 * @param xp The extended data of the particle
326
 * @param dt The drift time-step.
327
328
 */
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
329
    struct part *restrict p, const struct xpart *restrict xp, float dt) {
330

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

333
334
  /* Predict smoothing length */
  const float w1 = p->force.h_dt * h_inv * dt;
335
  if (fabsf(w1) < 0.2f)
336
    p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
337
  else
338
339
340
341
342
343
344
345
346
347
348
349
    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);
  }

350
351
  /* Predict the entropy */
  p->entropy += p->entropy_dt * dt;
352
353

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

356
357
358
359
  /* 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 */
360
  const float rho_bar_inv = 1.f / p->rho_bar;
361
  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
362

363
  /* Update the variables */
364
  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
365
  p->force.soundspeed = soundspeed;
366
  p->force.P_over_rho2 = P_over_rho2;
367
368
369
370
371
372
373
374
375
376
377
378
379
380
}

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

Matthieu Schaller's avatar
Matthieu Schaller committed
381
  p->entropy_dt =
382
      0.5f * gas_entropy_from_internal_energy(p->rho_bar, p->entropy_dt);
383
384
385
386
387
388
389
390
391
392
393
}

/**
 * @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(
394
    struct part *restrict p, struct xpart *restrict xp, float dt) {
395
396

  /* Do not decrease the entropy (temperature) by more than a factor of 2*/
397
398
399
400
  if (p->entropy_dt < -0.5f * xp->entropy_full / dt) {
    p->entropy_dt = -0.5f * xp->entropy_full / dt;
  }
  xp->entropy_full += p->entropy_dt * dt;
401
402

  /* Compute the pressure */
403
404
  const float pressure =
      gas_pressure_from_entropy(p->rho_bar, xp->entropy_full);
405

406
407
408
409
  /* 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 */
410
  const float rho_bar_inv = 1.f / p->rho_bar;
411
  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
412
413
414

  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
  p->force.soundspeed = soundspeed;
415
  p->force.P_over_rho2 = P_over_rho2;
416
417
418
419
420
421
422
423
424
425
}

/**
 *  @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(
426
    struct part *restrict p, struct xpart *restrict xp) {
427
428

  /* We read u in the entropy field. We now get S from u */
429
430
  xp->entropy_full = gas_entropy_from_internal_energy(p->rho_bar, p->entropy);
  p->entropy = xp->entropy_full;
431
432
433
  p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);

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

436
437
438
439
  /* 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 */
440
  const float rho_bar_inv = 1.f / p->rho_bar;
441
442
  const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;

443
  p->force.soundspeed = soundspeed;
444
  p->force.P_over_rho2 = P_over_rho2;
445
446
}

447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
/**
 * @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;
  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];

  hydro_reset_acceleration(p);
  hydro_init_part(p);
}

470
#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */