hydro.h 12.1 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 *
 * 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
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#ifndef SWIFT_MINIMAL_HYDRO_H
#define SWIFT_MINIMAL_HYDRO_H

/**
 * @file Minimal/hydro.h
 * @brief Minimal conservative implementation of SPH (Non-neighbour loop
 * equations)
 *
 * The thermal variable is the internal energy (u). Simple constant
 * viscosity term without switches is implemented. No thermal conduction
 * term is implemented.
 *
 * This corresponds to equations (43), (44), (45), (101), (103)  and (104) with
 * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
 * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
 * pp. 759-794.
 */
36

37
#include "adiabatic_index.h"
38
#include "approx_math.h"
39
#include "dimension.h"
40
#include "equation_of_state.h"
41
42
#include "hydro_properties.h"
#include "kernel_hydro.h"
43
44
45
46
47
48
49
50
51
52
53
54
55
56

/**
 * @brief Returns the internal energy of a particle
 *
 * For implementations where the main thermodynamic variable
 * is not internal energy, this function computes the internal
 * energy from the thermodynamic variable.
 *
 * @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) {

57
  return p->u + p->u_dt * dt;
58
59
60
61
62
63
64
65
66
67
68
}

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

69
70
71
  const float u = p->u + p->u_dt * dt;

  return gas_pressure_from_internal_energy(p->rho, u);
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
}

/**
 * @brief Returns the entropy of a particle
 *
 * For implementations where the main thermodynamic variable
 * is not entropy, this function computes the entropy from
 * the thermodynamic variable.
 *
 * @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) {

87
88
89
  const float u = p->u + p->u_dt * dt;

  return gas_entropy_from_internal_energy(p->rho, u);
90
91
92
93
94
95
96
97
98
99
100
}

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

101
102
103
  const float u = p->u + p->u_dt * dt;

  return gas_soundspeed_from_internal_energy(p->rho, u);
104
}
105

106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
/**
 * @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;
}

128
129
130
131
132
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
/**
 * @brief Modifies the thermal state of a particle to the imposed internal
 * energy
 *
 * This overrides the current state of the particle but does *not* change its
 * time-derivatives
 *
 * @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) {

  p->u = u;
}

/**
 * @brief Modifies the thermal state of a particle to the imposed entropy
 *
 * This overrides the current state of the particle but does *not* change its
 * time-derivatives
 *
 * @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->u = gas_internal_energy_from_entropy(p->rho, S);
}

159
160
161
/**
 * @brief Computes the hydro time-step of a given particle
 *
162
163
164
 * This function returns the time-step of a particle given its hydro-dynamical
 * state. A typical time-step calculation would be the use of the CFL condition.
 *
165
166
 * @param p Pointer to the particle data
 * @param xp Pointer to the extended particle data
167
 * @param hydro_properties The SPH parameters
168
169
170
 *
 */
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
171
172
    const struct part *restrict p, const struct xpart *restrict xp,
    const struct hydro_props *restrict hydro_properties) {
173
174

  const float CFL_condition = hydro_properties->CFL_condition;
175
176

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

180
  return dt_cfl;
181
182
}

183
184
185
186
/**
 * @brief Initialises the particles for the first time
 *
 * This function is called only once just after the ICs have been
187
188
 * read in to do some conversions or assignments between the particle
 * and extended particle fields.
189
190
191
192
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
 */
193
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
194
    struct part *restrict p, struct xpart *restrict xp) {
195

196
197
198
199
200
  p->ti_begin = 0;
  p->ti_end = 0;
  xp->v_full[0] = p->v[0];
  xp->v_full[1] = p->v[1];
  xp->v_full[2] = p->v[2];
201
202
}

203
204
205
206
/**
 * @brief Prepares a particle for the density calculation.
 *
 * Zeroes all the relevant arrays in preparation for the sums taking place in
207
208
 * the various density loop over neighbours. Typically, all fields of the
 * density sub-structure of a particle get zeroed in here.
209
210
211
 *
 * @param p The particle to act upon
 */
212
__attribute__((always_inline)) INLINE static void hydro_init_part(
213
    struct part *restrict p) {
214

215
216
217
218
219
220
221
222
223
224
225
  p->density.wcount = 0.f;
  p->density.wcount_dh = 0.f;
  p->rho = 0.f;
  p->rho_dh = 0.f;
}

/**
 * @brief Finishes the density calculation.
 *
 * Multiplies the density and number of neighbours by the appropiate constants
 * and add the self-contribution term.
226
227
228
 * Additional quantities such as velocity gradients will also get the final
 *terms
 * added to them here.
229
230
231
232
 *
 * @param p The particle to act upon
 * @param time The current time
 */
233
__attribute__((always_inline)) INLINE static void hydro_end_density(
234
    struct part *restrict p, float time) {
235
236
237

  /* Some smoothing length multiples. */
  const float h = p->h;
238
239
240
  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) */
241

242
243
  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
244
  p->rho_dh -= hydro_dimension * p->mass * kernel_root;
245
246
247
  p->density.wcount += kernel_root;

  /* Finish the calculation by inserting the missing h-factors */
248
249
  p->rho *= h_inv_dim;
  p->rho_dh *= h_inv_dim_plus_one;
250
  p->density.wcount *= kernel_norm;
251
  p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
252
253
254
255

  const float irho = 1.f / p->rho;

  /* Compute the derivative term */
256
  p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
257
258
259
260
261
}

/**
 * @brief Prepare a particle for the force calculation.
 *
262
263
264
265
266
267
 * This function is called in the ghost task to convert some quantities coming
 * from the density loop over neighbours into quantities ready to be used in the
 * force loop over neighbours. Quantities are typically read from the density
 * sub-structure and written to the force sub-structure.
 * Examples of calculations done here include the calculation of viscosity term
 * constants, thermal conduction terms, hydro conversions, etc.
268
269
270
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
271
272
 * @param ti_current The current time (on the timeline)
 * @param timeBase The minimal time-step size
273
 */
274
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
275
276
    struct part *restrict p, struct xpart *restrict xp, int ti_current,
    double timeBase) {
277

278
279
280
281
  const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
  const float pressure = hydro_get_pressure(p, half_dt);

  p->force.pressure = pressure;
282
283
284
285
286
287
}

/**
 * @brief Reset acceleration fields of a particle
 *
 * Resets all hydro acceleration and time derivative fields in preparation
288
 * for the sums taking  place in the various force tasks.
289
290
291
 *
 * @param p The particle to act upon
 */
292
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
293
    struct part *restrict p) {
294
295
296
297
298
299
300

  /* 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. */
301
  p->u_dt = 0.0f;
302
  p->force.h_dt = 0.0f;
303
304
305
306
307
308
  p->force.v_sig = 0.0f;
}

/**
 * @brief Predict additional particle fields forward in time when drifting
 *
309
310
311
 * Additional hydrodynamic quantites are drifted forward in time here. These
 * include thermal quantities (thermal energy or total energy or entropy, ...).
 *
312
313
314
315
316
317
 * @param p The particle.
 * @param xp The extended data of the particle.
 * @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).
 * @param timeBase The minimal time-step size.
318
319
 */
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
    struct part *restrict p, const struct xpart *restrict xp, float dt, int t0,
    int t1, double timeBase) {

  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);
338

339
340
341
  /* Drift the pressure */
  const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
  p->force.pressure = hydro_get_pressure(p, dt_entr);
342
343
344
345
346
}

/**
 * @brief Finishes the force calculation.
 *
347
348
349
 * Multiplies the force and accelerations by the appropiate constants
 * and add the self-contribution term. In most cases, there is nothing
 * to do here.
350
351
352
 *
 * @param p The particle to act upon
 */
353
__attribute__((always_inline)) INLINE static void hydro_end_force(
354
355
    struct part *restrict p) {

356
  p->force.h_dt *= p->h * hydro_dimension_inv;
357
}
358
359
360
361

/**
 * @brief Kick the additional variables
 *
362
363
364
 * Additional hydrodynamic quantites are kicked forward in time here. These
 * include thermal quantities (thermal energy or total energy or entropy, ...).
 *
365
 * @param p The particle to act upon
366
 * @param xp The particle extended data to act upon
367
 * @param dt The time-step for this kick
368
 * @param half_dt The half time-step for this kick
369
 */
370
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
371
372
    struct part *restrict p, struct xpart *restrict xp, float dt,
    float half_dt) {
373

374
375
376
377
378
379
  /* Do not decrease the energy by more than a factor of 2*/
  const float u_change = p->u_dt * dt;
  if (u_change > -0.5f * p->u)
    p->u += u_change;
  else
    p->u *= 0.5f;
380

381
382
  /* Do not 'overcool' when timestep increases */
  if (p->u + p->u_dt * half_dt < 0.5f * p->u) p->u_dt = -0.5f * p->u / half_dt;
383
}
384
385

/**
386
 * @brief Converts hydro quantity of a particle at the start of a run
387
 *
388
389
390
391
 * This function is called once at the end of the engine_init_particle()
 * routine (at the start of a calculation) after the densities of
 * particles have been computed.
 * This can be used to convert internal energy into entropy for instance.
392
393
394
 *
 * @param p The particle to act upon
 */
395
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
396
    struct part *restrict p) {}
397
398

#endif /* SWIFT_MINIMAL_HYDRO_H */