hydro.h 7.82 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
19
 *
 * 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/>.
 *
 ******************************************************************************/

20
#include "approx_math.h"
21
#include "gamma.h"
22

23
24
25
/**
 * @brief Computes the hydro time-step of a given particle
 *
26
27
28
 * 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.
 *
29
30
 * @param p Pointer to the particle data
 * @param xp Pointer to the extended particle data
31
 * @param hydro_properties The SPH parameters
32
33
34
 *
 */
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
35
    const struct part* p, const struct xpart* xp,
36
37
38
    const struct hydro_props* hydro_properties) {

  const float CFL_condition = hydro_properties->CFL_condition;
39
40

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

44
  return dt_cfl;
45
46
}

47
48
49
50
/**
 * @brief Initialises the particles for the first time
 *
 * This function is called only once just after the ICs have been
51
52
 * read in to do some conversions or assignments between the particle
 * and extended particle fields.
53
54
55
56
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
 */
57
58
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
    struct part* p, struct xpart* xp) {
59
60
61
62

  xp->u_full = p->u;
}

63
64
65
66
/**
 * @brief Prepares a particle for the density calculation.
 *
 * Zeroes all the relevant arrays in preparation for the sums taking place in
67
68
 * the various density loop over neighbours. Typically, all fields of the
 * density sub-structure of a particle get zeroed in here.
69
70
71
 *
 * @param p The particle to act upon
 */
72
73
__attribute__((always_inline)) INLINE static void hydro_init_part(
    struct part* p) {
74
75
76
77
78
79
80
81
82
83
84
  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.
85
86
87
 * Additional quantities such as velocity gradients will also get the final
 *terms
 * added to them here.
88
89
90
91
 *
 * @param p The particle to act upon
 * @param time The current time
 */
92
93
__attribute__((always_inline)) INLINE static void hydro_end_density(
    struct part* p, float time) {
94
95
96
97
98
99
100

  /* Some smoothing length multiples. */
  const float h = p->h;
  const float ih = 1.0f / h;
  const float ih2 = ih * ih;
  const float ih4 = ih2 * ih2;

101
102
  /* Final operation on the density (add self-contribution). */
  p->rho += p->mass * kernel_root;
103
  p->rho_dh -= 3.0f * p->mass * kernel_root;
104
105
106
107
108
109
  p->density.wcount += kernel_root;

  /* Finish the calculation by inserting the missing h-factors */
  p->rho *= ih * ih2;
  p->rho_dh *= ih4;
  p->density.wcount *= (4.0f / 3.0f * M_PI * kernel_gamma3);
110
111
112
113
114
115
  p->density.wcount_dh *= ih * (4.0f / 3.0f * M_PI * kernel_gamma4);

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

  /* Compute the derivative term */
  p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
116
117
118
119
120
}

/**
 * @brief Prepare a particle for the force calculation.
 *
121
122
123
124
125
126
 * 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.
127
128
129
 *
 * @param p The particle to act upon
 * @param xp The extended particle data to act upon
130
131
 * @param ti_current The current time (on the timeline)
 * @param timeBase The minimal time-step size
132
 */
133
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
134
    struct part* p, struct xpart* xp, int ti_current, double timeBase) {
135

136
  p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
137
138
139
140
141
142
}

/**
 * @brief Reset acceleration fields of a particle
 *
 * Resets all hydro acceleration and time derivative fields in preparation
143
 * for the sums taking  place in the various force tasks.
144
145
146
 *
 * @param p The particle to act upon
 */
147
148
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
    struct part* p) {
149
150
151
152
153
154
155

  /* 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. */
156
  p->u_dt = 0.0f;
157
158
159
160
161
162
163
  p->h_dt = 0.0f;
  p->force.v_sig = 0.0f;
}

/**
 * @brief Predict additional particle fields forward in time when drifting
 *
164
165
166
 * Additional hydrodynamic quantites are drifted forward in time here. These
 * include thermal quantities (thermal energy or total energy or entropy, ...).
 *
167
168
 * @param p The particle
 * @param xp The extended data of the particle
169
170
171
 * @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
172
173
 */
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
174
    struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {
175

176
  p->u = xp->u_full;
177

178
  /* Need to recompute the pressure as well */
179
  p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
180
181
182
183
184
}

/**
 * @brief Finishes the force calculation.
 *
185
186
187
 * Multiplies the force and accelerations by the appropiate constants
 * and add the self-contribution term. In most cases, there is nothing
 * to do here.
188
189
190
 *
 * @param p The particle to act upon
 */
191
192
__attribute__((always_inline)) INLINE static void hydro_end_force(
    struct part* p) {}
193
194
195
196

/**
 * @brief Kick the additional variables
 *
197
198
199
 * Additional hydrodynamic quantites are kicked forward in time here. These
 * include thermal quantities (thermal energy or total energy or entropy, ...).
 *
200
 * @param p The particle to act upon
201
 * @param xp The particle extended data to act upon
202
 * @param dt The time-step for this kick
203
 * @param half_dt The half time-step for this kick
204
 */
205
206
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
    struct part* p, struct xpart* xp, float dt, float half_dt) {
207
208
209
210
211
212
213

  /* Kick in momentum space */
  xp->u_full += p->u_dt * dt;

  /* Get the predicted internal energy */
  p->u = xp->u_full - half_dt * p->u_dt;
}
214
215

/**
216
 * @brief Converts hydro quantity of a particle at the start of a run
217
 *
218
219
220
221
 * 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.
222
223
224
 *
 * @param p The particle to act upon
 */
225
226
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
    struct part* p) {}
227
228
229
230
231
232
233
234
235

/**
 * @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
236
 * @param dt Time since the last kick
237
 */
238
239
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
    const struct part* p, float dt) {
240
241
242

  return p->u;
}