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

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

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

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

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

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

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

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

  return p->force.soundspeed;
}
86

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

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

119
  return gas_internal_energy_from_entropy(p->rho, p->entropy_dt);
120
121
122
}

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

133
  p->entropy_dt = gas_entropy_from_internal_energy(p->rho, du_dt);
134
135
}

136
137
138
139
140
141
142
143
/**
 * @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(
144
145
    const struct part *restrict p, const struct xpart *restrict xp,
    const struct hydro_props *restrict hydro_properties) {
146

147
  const float CFL_condition = hydro_properties->CFL_condition;
148

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

153
  return dt_cfl;
154
155
}

156
157
158
159
160
161
162
163
164
165
/**
 * @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) {}

166
167
168
/**
 * @brief Prepares a particle for the density calculation.
 *
169
 * Zeroes all the relevant arrays in preparation for the sums taking place in
170
171
172
 * the variaous density tasks
 *
 * @param p The particle to act upon
173
 * @param hs #hydro_space containing hydro specific space information.
174
 */
175
__attribute__((always_inline)) INLINE static void hydro_init_part(
176
    struct part *restrict p, const struct hydro_space *hs) {
177

178
#ifdef DEBUG_INTERACTIONS_SPH
179
  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS; ++i) p->ids_ngbs_density[i] = -1;
180
181
182
  p->num_ngb_density = 0;
#endif
  
183
  p->rho = 0.f;
184
185
  p->density.wcount = 0.f;
  p->density.wcount_dh = 0.f;
186
  p->density.rho_dh = 0.f;
187
  p->density.div_v = 0.f;
188
189
190
  p->density.rot_v[0] = 0.f;
  p->density.rot_v[1] = 0.f;
  p->density.rot_v[2] = 0.f;
191
192
193
}

/**
194
 * @brief Finishes the density calculation.
195
 *
196
 * Multiplies the density and number of neighbours by the appropiate constants
197
198
199
200
 * and add the self-contribution term.
 *
 * @param p The particle to act upon
 */
201
__attribute__((always_inline)) INLINE static void hydro_end_density(
202
    struct part *restrict p) {
203
204
205

  /* Some smoothing length multiples. */
  const float h = p->h;
206
207
208
  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) */
209

210
211
  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
212
  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
213
  p->density.wcount += kernel_root;
214
  p->density.wcount_dh -= hydro_dimension * kernel_root;
215

216
  /* Finish the calculation by inserting the missing h-factors */
217
  p->rho *= h_inv_dim;
218
  p->density.rho_dh *= h_inv_dim_plus_one;
219
220
  p->density.wcount *= h_inv_dim;
  p->density.wcount_dh *= h_inv_dim_plus_one;
221

222
  const float rho_inv = 1.f / p->rho;
Matthieu Schaller's avatar
Matthieu Schaller committed
223

224
  /* Finish calculation of the velocity curl components */
225
226
227
  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;
228

229
  /* Finish calculation of the velocity divergence */
230
  p->density.div_v *= h_inv_dim_plus_one * rho_inv;
231
232
}

233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
/**
 * @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
 */
__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
    struct part *restrict p, struct xpart *restrict xp) {

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

258
259
260
261
262
263
264
/**
 * @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
265
266
 * @param ti_current The current time (on the timeline)
 * @param timeBase The minimal time-step size
267
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
268
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
Matthieu Schaller's avatar
Matthieu Schaller committed
269
    struct part *restrict p, struct xpart *restrict xp) {
270

271
272
  const float fac_mu = 1.f; /* Will change with cosmological integration */

273
274
275
  /* Inverse of the physical density */
  const float rho_inv = 1.f / p->rho;

276
  /* Compute the norm of the curl */
277
278
279
  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]);
280

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

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

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

290
  /* Divide the pressure by the density squared to get the SPH term */
291
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
292

293
  /* Compute the Balsara switch */
Matthieu Schaller's avatar
Matthieu Schaller committed
294
295
296
  const float balsara =
      abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);

297
  /* Compute the "grad h" term */
298
  const float omega_inv =
299
300
      1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);

301
  /* Update variables. */
302
  p->force.f = omega_inv;
Matthieu Schaller's avatar
Matthieu Schaller committed
303
  p->force.P_over_rho2 = P_over_rho2;
304
305
  p->force.soundspeed = soundspeed;
  p->force.balsara = balsara;
306
307
308
309
310
311
312
313
314
315
}

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

319
#ifdef DEBUG_INTERACTIONS_SPH
320
  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS; ++i) p->ids_ngbs_force[i] = -1;
321
322
323
  p->num_ngb_force = 0;
#endif
  
324
  /* Reset the acceleration. */
Matthieu Schaller's avatar
Matthieu Schaller committed
325
326
327
  p->a_hydro[0] = 0.0f;
  p->a_hydro[1] = 0.0f;
  p->a_hydro[2] = 0.0f;
328

329
  /* Reset the time derivatives. */
330
  p->entropy_dt = 0.0f;
331
  p->force.h_dt = 0.0f;
332

333
  /* Reset maximal signal velocity */
334
  p->force.v_sig = p->force.soundspeed;
335
336
}

337
338
339
340
341
342
343
/**
 * @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.
 */
344
345
346
347
348
349
350
__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];
351
352
353

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

356
357
358
359
360
/**
 * @brief Predict additional particle fields forward in time when drifting
 *
 * @param p The particle
 * @param xp The extended data of the particle
361
 * @param dt The drift time-step.
362
363
 * @param t0 The time at the start of the drift
 * @param t1 The time at the end of the drift
364
 * @param timeBase The minimal time-step size
365
 */
366
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
Matthieu Schaller's avatar
Matthieu Schaller committed
367
    struct part *restrict p, const struct xpart *restrict xp, float dt) {
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383

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

  /* Predict smoothing length */
  const float w1 = p->force.h_dt * h_inv * dt;
  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);
384

385
386
  /* Predict the entropy */
  p->entropy += p->entropy_dt * dt;
Matthieu Schaller's avatar
Matthieu Schaller committed
387

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

391
392
  /* Compute the new sound speed */
  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
393

394
395
  /* Divide the pressure by the density squared to get the SPH term */
  const float rho_inv = 1.f / p->rho;
396
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
Matthieu Schaller's avatar
Matthieu Schaller committed
397

398
399
  /* Update variables */
  p->force.soundspeed = soundspeed;
400
  p->force.P_over_rho2 = P_over_rho2;
401
402
403
}

/**
404
 * @brief Finishes the force calculation.
405
 *
406
 * Multiplies the forces and accelerationsby the appropiate constants
407
408
409
 *
 * @param p The particle to act upon
 */
410
__attribute__((always_inline)) INLINE static void hydro_end_force(
411
    struct part *restrict p) {
412

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

415
416
  p->entropy_dt =
      0.5f * gas_entropy_from_internal_energy(p->rho, p->entropy_dt);
417
418
}

Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
419
420
421
422
/**
 * @brief Kick the additional variables
 *
 * @param p The particle to act upon
423
424
 * @param xp The particle extended data to act upon
 * @param dt The time-step for this kick
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
425
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
426
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
427
    struct part *restrict p, struct xpart *restrict xp, float dt) {
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
428

429
  /* Do not decrease the entropy by more than a factor of 2 */
430
  if (dt > 0. && p->entropy_dt * dt < -0.5f * xp->entropy_full) {
Matthieu Schaller's avatar
Matthieu Schaller committed
431
432
433
434
435
    /* 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); */
    p->entropy_dt = -0.5f * xp->entropy_full / dt;
  }
436
  xp->entropy_full += p->entropy_dt * dt;
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
437

438
  /* Compute the pressure */
439
  const float pressure = gas_pressure_from_entropy(p->rho, xp->entropy_full);
440

441
442
443
444
  /* 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 */
445
  const float rho_inv = 1.f / p->rho;
446
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
447
448

  p->force.soundspeed = soundspeed;
449
  p->force.P_over_rho2 = P_over_rho2;
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
450
451
}

452
/**
453
 *  @brief Converts hydro quantity of a particle at the start of a run
454
455
456
457
458
 *
 * Requires the density to be known
 *
 * @param p The particle to act upon
 */
459
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
460
    struct part *restrict p, struct xpart *restrict xp) {
461

462
  /* We read u in the entropy field. We now get S from u */
463
464
  xp->entropy_full = gas_entropy_from_internal_energy(p->rho, p->entropy);
  p->entropy = xp->entropy_full;
Matthieu Schaller's avatar
Matthieu Schaller committed
465

466
467
468
  /* Compute the pressure */
  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);

469
470
471
472
  /* 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 */
473
  const float rho_inv = 1.f / p->rho;
474
  const float P_over_rho2 = pressure * rho_inv * rho_inv;
475
476

  p->force.soundspeed = soundspeed;
477
  p->force.P_over_rho2 = P_over_rho2;
478
}
479

Matthieu Schaller's avatar
Matthieu Schaller committed
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
/**
 * @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];
  xp->entropy_full = p->entropy;

  hydro_reset_acceleration(p);
499
  hydro_init_part(p, NULL);
Matthieu Schaller's avatar
Matthieu Schaller committed
500
501
}

502
#endif /* SWIFT_GADGET2_HYDRO_H */