hydro.h 16.9 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_DEFAULT_HYDRO_H
#define SWIFT_DEFAULT_HYDRO_H
21

22
#include "adiabatic_index.h"
23
#include "approx_math.h"
24
25
#include "cosmology.h"
#include "dimension.h"
26
#include "equation_of_state.h"
27
#include "hydro_space.h"
28
#include "kernel_hydro.h"
29
#include "minmax.h"
30

31
32
#include <float.h>

33
/**
34
 * @brief Returns the comoving internal energy of a particle
35
36
37
 *
 * @param p The particle of interest
 */
38
39
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_internal_energy(const struct part *restrict p) {
40
41
42
43
44

  return p->u;
}

/**
45
46
47
48
49
50
51
52
53
54
55
56
57
58
 * @brief Returns the physical internal energy of a particle
 *
 * @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 cosmo->a_factor_internal_energy * p->u;
}

/**
 * @brief Returns the comoving pressure of a particle
59
60
61
 *
 * @param p The particle of interest
 */
62
__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
63
    const struct part *restrict p) {
64
65
66
67
68

  return gas_pressure_from_internal_energy(p->rho, p->u);
}

/**
69
70
71
72
73
74
75
76
77
78
79
80
81
82
 * @brief Returns the physical pressure of a particle
 *
 * @param p The particle of interest
 * @param cosmo The cosmological model.
 */
__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
    const struct part *restrict p, const struct cosmology *cosmo) {

  return cosmo->a_factor_pressure *
         gas_pressure_from_internal_energy(p->rho, p->u);
}

/**
 * @brief Returns the comoving entropy of a particle
83
84
85
 *
 * @param p The particle of interest
 */
86
__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
87
    const struct part *restrict p) {
88
89
90
91
92

  return gas_entropy_from_internal_energy(p->rho, p->u);
}

/**
93
 * @brief Returns the physical entropy of a particle
94
95
 *
 * @param p The particle of interest
96
 * @param cosmo The cosmological model.
97
 */
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
    const struct part *restrict p, const struct cosmology *cosmo) {

  /* Note: no cosmological conversion required here with our choice of
   * coordinates. */
  return gas_entropy_from_internal_energy(p->rho, p->u);
}

/**
 * @brief Returns the comoving sound speed of a particle
 *
 * @param p The particle of interest
 */
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_soundspeed(const struct part *restrict p) {
113
114
115
116

  return p->force.soundspeed;
}

117
/**
118
 * @brief Returns the physical sound speed of a particle
119
120
 *
 * @param p The particle of interest
121
 * @param cosmo The cosmological model.
122
 */
123
124
125
126
127
128
129
130
131
132
133
134
135
__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(
136
137
138
139
140
    const struct part *restrict p) {

  return p->rho;
}

141
142
143
144
145
146
147
148
149
150
151
/**
 * @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;
}

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

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

175
176
177
178
179
/**
 * @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.
180
181
 * @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.
182
183
184
 * @param v (return) The velocities at the current time.
 */
__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
185
186
187
188
189
190
191
192
193
    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;
194
195
}

196
/**
197
 * @brief Returns the time derivative of internal energy of a particle
198
 *
199
 * We assume a constant density.
200
 *
201
 * @param p The particle of interest
202
 */
203
204
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
    const struct part *restrict p) {
205

206
  return p->force.u_dt;
207
208
209
}

/**
210
 * @brief Returns the time derivative of internal energy of a particle
211
 *
212
 * We assume a constant density.
213
 *
214
215
 * @param p The particle of interest.
 * @param du_dt The new time derivative of the internal energy.
216
 */
217
218
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
    struct part *restrict p, float du_dt) {
219

220
  p->force.u_dt = du_dt;
221
222
}

223
224
225
226
227
/**
 * @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
228
229
 * @param hydro_properties The constants used in the scheme
 * @param cosmo The cosmological model.
230
231
 */
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
232
    const struct part *restrict p, const struct xpart *restrict xp,
233
234
    const struct hydro_props *restrict hydro_properties,
    const struct cosmology *restrict cosmo) {
235

236
  const float CFL = hydro_properties->CFL_condition;
237
238

  /* CFL condition */
239
240
  const float dt_cfl = 2.f * kernel_gamma * CFL * cosmo->a * p->h /
                       (cosmo->a_factor_sound_speed * p->force.v_sig);
241
242

  /* Limit change in u */
Matthieu Schaller's avatar
Matthieu Schaller committed
243
244
245
  const float dt_u_change =
      (p->force.u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->force.u_dt)
                              : FLT_MAX;
246

247
  return min(dt_cfl, dt_u_change);
248
249
}

250
251
252
253
254
255
256
257
258
259
/**
 * @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) {}

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

272
273
274
275
276
  p->density.wcount = 0.f;
  p->density.wcount_dh = 0.f;
  p->rho = 0.f;
  p->rho_dh = 0.f;
  p->density.div_v = 0.f;
277
278
279
  p->density.rot_v[0] = 0.f;
  p->density.rot_v[1] = 0.f;
  p->density.rot_v[2] = 0.f;
280
281
282
}

/**
283
 * @brief Finishes the density calculation.
284
 *
285
 * Multiplies the density and number of neighbours by the appropiate constants
286
287
288
 * and add the self-contribution term.
 *
 * @param p The particle to act upon
289
 * @param cosmo The current cosmological model.
290
 */
291
__attribute__((always_inline)) INLINE static void hydro_end_density(
292
    struct part *restrict p, const struct cosmology *cosmo) {
293
294
295

  /* Some smoothing length multiples. */
  const float h = p->h;
296
297
298
  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) */
299

300
301
  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
302
  p->rho_dh -= hydro_dimension * p->mass * kernel_root;
303
304
305
  p->density.wcount += kernel_root;

  /* Finish the calculation by inserting the missing h-factors */
306
307
  p->rho *= h_inv_dim;
  p->rho_dh *= h_inv_dim_plus_one;
308
  p->density.wcount *= kernel_norm;
309
  p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
310
311

  const float irho = 1.f / p->rho;
312
  const float a_inv2 = cosmo->a2_inv;
313

314
  /* Finish calculation of the velocity curl components */
315
316
317
  p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * irho;
  p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * irho;
  p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * irho;
318
319

  /* Finish calculation of the velocity divergence */
320
  p->density.div_v *= h_inv_dim_plus_one * a_inv2 * irho;
321
322
}

323
324
325
326
327
/**
 * @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
328
 * @param cosmo The current cosmological model.
329
330
 */
__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
331
332
    struct part *restrict p, struct xpart *restrict xp,
    const struct cosmology *cosmo) {
333

334
335
336
337
338
  /* 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 */

339
  /* Re-set problematic values */
340
341
  p->rho = p->mass * kernel_root * h_inv_dim;
  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
342
343
344
345
346
347
348
349
  p->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;
}

350
351
352
353
354
355
356
/**
 * @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
357
 * @param cosmo The current cosmological model.
358
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
359
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
360
361
362
363
    struct part *restrict p, struct xpart *restrict xp,
    const struct cosmology *cosmo) {

  const float fac_mu = cosmo->a_factor_mu;
364
365
366

  /* Some smoothing length multiples. */
  const float h = p->h;
367
  const float h_inv = 1.0f / h;
368
369

  /* Pre-compute some stuff for the balsara switch. */
370
  const float normDiv_v = fabs(p->density.div_v);
371
  const float normRot_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
Matthieu Schaller's avatar
Matthieu Schaller committed
372
                                p->density.rot_v[1] * p->density.rot_v[1] +
373
                                p->density.rot_v[2] * p->density.rot_v[2]);
374
375
376

  /* Compute this particle's sound speed. */
  const float u = p->u;
Matthieu Schaller's avatar
Matthieu Schaller committed
377
378
  const float fc = p->force.soundspeed =
      sqrtf(hydro_gamma * hydro_gamma_minus_one * u);
379

380
381
382
  /* Compute the derivative term */
  p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh / p->rho);

383
  /* Compute the P/Omega/rho2. */
384
  xp->omega = 1.0f + hydro_dimension_inv * h * p->rho_dh / p->rho;
385
  p->force.P_over_rho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
386
387

  /* Balsara switch */
388
389
  p->force.balsara =
      normDiv_v / (normDiv_v + normRot_v + 0.0001f * fac_mu * fc * h_inv);
390
391

  /* Viscosity parameter decay time */
392
393
  /* const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
   */
394
395

  /* Viscosity source term */
396
  /* const float S = max(-normDiv_v, 0.f); */
397
398

  /* Compute the particle's viscosity parameter time derivative */
399
400
  /* const float alpha_dot = (const_viscosity_alpha_min - p->alpha) / tau + */
  /*                         (const_viscosity_alpha_max - p->alpha) * S; */
401
402

  /* Update particle's viscosity paramter */
403
  /* p->alpha += alpha_dot * (p->ti_end - p->ti_begin) * timeBase; */  // MATTHIEU
404
405
406
407
408
409
410
411
412
413
}

/**
 * @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
 */
414
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
415
    struct part *restrict p) {
416
417

  /* Reset the acceleration. */
Matthieu Schaller's avatar
Matthieu Schaller committed
418
419
420
  p->a_hydro[0] = 0.0f;
  p->a_hydro[1] = 0.0f;
  p->a_hydro[2] = 0.0f;
421

422
423
  /* Reset the time derivatives. */
  p->force.u_dt = 0.0f;
424
  p->force.h_dt = 0.0f;
425
426
427
  p->force.v_sig = 0.0f;
}

428
429
430
431
432
433
434
435
436
437
438
439
440
441
/**
 * @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];
442
443

  p->u = xp->u_full;
444
445
}

446
447
448
449
450
/**
 * @brief Predict additional particle fields forward in time when drifting
 *
 * @param p The particle
 * @param xp The extended data of the particle
451
452
 * @param dt_drift The drift time-step for positions.
 * @param dt_therm The drift time-step for thermal quantities.
453
 */
454
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
455
456
    struct part *restrict p, struct xpart *restrict xp, float dt_drift,
    float dt_therm) {
457
458
  float u, w;

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

  /* Predict smoothing length */
462
  const float w1 = p->force.h_dt * h_inv * dt_drift;
463
464
465
466
467
468
469
470
471
472
473
  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);
Matthieu Schaller's avatar
Matthieu Schaller committed
474

475
  /* Predict internal energy */
476
  w = p->force.u_dt / p->u * dt_therm;
477
478
  if (fabsf(w) < 0.2f)
    u = p->u *= approx_expf(w);
479
480
  else
    u = p->u *= expf(w);
481

482
  /* Predict gradient term */
483
  p->force.P_over_rho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
484
485
486
}

/**
487
 * @brief Finishes the force calculation.
488
 *
489
 * Multiplies the forces and accelerationsby the appropiate constants
490
491
 *
 * @param p The particle to act upon
492
 * @param cosmo The current cosmological model.
493
 */
494
__attribute__((always_inline)) INLINE static void hydro_end_force(
495
    struct part *restrict p, const struct cosmology *cosmo) {
496
  p->force.h_dt *= p->h * hydro_dimension_inv;
497
}
498
499
500
501
502

/**
 * @brief Kick the additional variables
 *
 * @param p The particle to act upon
503
 * @param xp The particle extended data to act upon
504
 * @param dt_therm The time-step for this kick (for thermodynamic quantities)
505
 */
506
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
507
    struct part *restrict p, struct xpart *restrict xp, float dt_therm) {}
508
509

/**
510
 *  @brief Converts hydro quantity of a particle at the start of a run
511
512
513
514
 *
 * Requires the density to be known
 *
 * @param p The particle to act upon
515
516
 * @param xp The extended data.
 * @param cosmo The cosmological model.
517
 */
518
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
519
520
    struct part *restrict p, struct xpart *restrict xp,
    const struct cosmology *cosmo) {}
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537

/**
 * @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];
538
539
540
  xp->a_grav[0] = 0.f;
  xp->a_grav[1] = 0.f;
  xp->a_grav[2] = 0.f;
541
542
543
  xp->u_full = p->u;

  hydro_reset_acceleration(p);
544
  hydro_init_part(p, NULL);
545
}
546

547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
/**
 * @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->u = u_init;
}

564
#endif /* SWIFT_DEFAULT_HYDRO_H */