hydro.h 6.65 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
23
24
25
26
27
28
29
30
/*******************************************************************************
 * 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/>.
 *
 ******************************************************************************/

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

  /* CFL condition */
Matthieu Schaller's avatar
Matthieu Schaller committed
31
  const float dt_cfl = 2.f * const_cfl * kernel_gamma * p->h / p->force.v_sig;
32
33

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

38
  return fminf(dt_cfl, dt_u_change);
39
40
41
42
43
}

/**
 * @brief Prepares a particle for the density calculation.
 *
44
 * Zeroes all the relevant arrays in preparation for the sums taking place in
45
46
47
48
 * the variaous density tasks
 *
 * @param p The particle to act upon
 */
49
50
__attribute__((always_inline))
    INLINE static void hydro_init_part(struct part* p) {
51
52
53
54
55
56
57
58
59
60
61
  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;
  p->density.curl_v[0] = 0.f;
  p->density.curl_v[1] = 0.f;
  p->density.curl_v[2] = 0.f;
}

/**
62
 * @brief Finishes the density calculation.
63
 *
64
 * Multiplies the density and number of neighbours by the appropiate constants
65
66
67
 * and add the self-contribution term.
 *
 * @param p The particle to act upon
Matthieu Schaller's avatar
Matthieu Schaller committed
68
 * @param time The current time
69
 */
70
__attribute__((always_inline))
Matthieu Schaller's avatar
Matthieu Schaller committed
71
    INLINE static void hydro_end_density(struct part* p, float time) {
72
73
74
75
76
77

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

79
80
81
  /* Final operation on the density. */
  p->rho = ih * ih2 * (p->rho + p->mass * kernel_root);
  p->rho_dh = (p->rho_dh - 3.0f * p->mass * kernel_root) * ih4;
82
83
84
85
  p->density.wcount =
      (p->density.wcount + kernel_root) * (4.0f / 3.0 * M_PI * kernel_gamma3);
  p->density.wcount_dh =
      p->density.wcount_dh * ih * (4.0f / 3.0 * M_PI * kernel_gamma3);
86
87
88
89
90
91
92
93
94
}

/**
 * @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
95
 * @param time The current time
96
 */
97
__attribute__((always_inline))
98
INLINE static void hydro_prepare_force(struct part* p, struct xpart* xp, float time) {
99
100
101
102
103
104
105
106
107
108

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

  /* Pre-compute some stuff for the balsara switch. */
  const float normDiv_v = fabs(p->density.div_v / p->rho * ih4);
  const float normCurl_v = sqrtf(p->density.curl_v[0] * p->density.curl_v[0] +
109
110
111
                                 p->density.curl_v[1] * p->density.curl_v[1] +
                                 p->density.curl_v[2] * p->density.curl_v[2]) /
                           p->rho * ih4;
112
113
114

  /* Compute this particle's sound speed. */
  const float u = p->u;
115
116
  const float fc = p->force.c =
      sqrtf(const_hydro_gamma * (const_hydro_gamma - 1.0f) * u);
117
118
119
120
121
122

  /* Compute the P/Omega/rho2. */
  xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho;
  p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (p->rho * xp->omega);

  /* Balsara switch */
123
  p->force.balsara = normDiv_v / (normDiv_v + normCurl_v + 0.0001f * fc * ih);
124
125
126
127
128
129
130
131
132

  /* Viscosity parameter decay time */
  const float tau = h / (2.f * const_viscosity_length * p->force.c);

  /* Viscosity source term */
  const float S = fmaxf(-normDiv_v, 0.f);

  /* Compute the particle's viscosity parameter time derivative */
  const float alpha_dot = (const_viscosity_alpha_min - p->alpha) / tau +
133
                          (const_viscosity_alpha_max - p->alpha) * S;
134
135
136
137
138
139
140
141
142
143
144
145
146

  /* Update particle's viscosity paramter */
  p->alpha += alpha_dot * (p->t_end - p->t_begin);
}

/**
 * @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
 */
147
148
__attribute__((always_inline))
    INLINE static void hydro_reset_acceleration(struct part* p) {
149
150

  /* Reset the acceleration. */
Matthieu Schaller's avatar
Matthieu Schaller committed
151
152
153
  p->a_hydro[0] = 0.0f;
  p->a_hydro[1] = 0.0f;
  p->a_hydro[2] = 0.0f;
154

155
156
  /* Reset the time derivatives. */
  p->force.u_dt = 0.0f;
157
  p->h_dt = 0.0f;
158
159
160
161
162
163
164
165
  p->force.v_sig = 0.0f;
}

/**
 * @brief Predict additional particle fields forward in time when drifting
 *
 * @param p The particle
 * @param xp The extended data of the particle
Matthieu Schaller's avatar
Matthieu Schaller committed
166
167
 * @param t0 The time at the start of the drift
 * @param t1 The time at the end of the drift
168
 */
169
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
170
    struct part* p, struct xpart* xp, float t0, float t1) {
171
172
  float u, w;

173
  const float dt = t1 - t0;
Matthieu Schaller's avatar
Matthieu Schaller committed
174

175
176
177
178
  /* Predict internal energy */
  w = p->force.u_dt / p->u * dt;
  if (fabsf(w) < 0.01f) /* 1st order expansion of exp(w) */
    u = p->u *=
179
        1.0f + w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f * w)));
180
181
  else
    u = p->u *= expf(w);
182

183
184
185
186
187
  /* Predict gradient term */
  p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (p->rho * xp->omega);
}

/**
188
 * @brief Finishes the force calculation.
189
 *
190
 * Multiplies the forces and accelerationsby the appropiate constants
191
192
193
 *
 * @param p The particle to act upon
 */
194
__attribute__((always_inline))
195
    INLINE static void hydro_end_force(struct part* p) {}
196

197
198


199
200
201
202
203
204
/**
 * @brief Kick the additional variables
 *
 * @param p The particle to act upon
 */
__attribute__((always_inline))
205
   INLINE static void hydro_kick_extra(struct part* p, float dt) {}
206
207
208
209
210
211
212
213
214

/**
 * @brief Converts hydro quantity of a particle
 *
 * Requires the density to be known
 *
 * @param p The particle to act upon
 */
__attribute__((always_inline))
Matthieu Schaller's avatar
Matthieu Schaller committed
215
    INLINE static void hydro_convert_quantities(struct part* p) {}