hydro_iact.h 13.3 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 *
 * 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/>.
 *
 ******************************************************************************/

#include "adiabatic_index.h"
#include "hydro_gradients.h"
#include "riemann.h"
23
#include "voronoi_algorithm.h"
24
25

/**
26
 * @brief Calculate the Voronoi cell by interacting particle pi and pj
27
 *
28
 * This method wraps around voronoi_cell_interact().
29
30
31
32
33
34
35
36
37
38
39
 *
 * @param r2 Squared distance between particle i and particle j.
 * @param dx Distance vector between the particles (dx = pi->x - pj->x).
 * @param hi Smoothing length of particle i.
 * @param hj Smoothing length of particle j.
 * @param pi Particle i.
 * @param pj Particle j.
 */
__attribute__((always_inline)) INLINE static void runner_iact_density(
    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {

40
41
42
43
44
45
46
  float mindx[3];

  voronoi_cell_interact(&pi->cell, dx, pj->id);
  mindx[0] = -dx[0];
  mindx[1] = -dx[1];
  mindx[2] = -dx[2];
  voronoi_cell_interact(&pj->cell, mindx, pi->id);
47
48
49
}

/**
50
 * @brief Calculate the Voronoi cell by interacting particle pi with pj
51
 *
52
 * This method wraps around voronoi_cell_interact().
53
54
55
56
57
58
59
60
61
62
63
 *
 * @param r2 Squared distance between particle i and particle j.
 * @param dx Distance vector between the particles (dx = pi->x - pj->x).
 * @param hi Smoothing length of particle i.
 * @param hj Smoothing length of particle j.
 * @param pi Particle i.
 * @param pj Particle j.
 */
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {

64
  voronoi_cell_interact(&pi->cell, dx, pj->id);
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
}

/**
 * @brief Calculate the gradient interaction between particle i and particle j
 *
 * This method wraps around hydro_gradients_collect, which can be an empty
 * method, in which case no gradients are used.
 *
 * @param r2 Squared distance between particle i and particle j.
 * @param dx Distance vector between the particles (dx = pi->x - pj->x).
 * @param hi Smoothing length of particle i.
 * @param hj Smoothing length of particle j.
 * @param pi Particle i.
 * @param pj Particle j.
 */
__attribute__((always_inline)) INLINE static void runner_iact_gradient(
    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {

  hydro_gradients_collect(r2, dx, hi, hj, pi, pj);
}

/**
 * @brief Calculate the gradient interaction between particle i and particle j:
 * non-symmetric version
 *
 * This method wraps around hydro_gradients_nonsym_collect, which can be an
 * empty method, in which case no gradients are used.
 *
 * @param r2 Squared distance between particle i and particle j.
 * @param dx Distance vector between the particles (dx = pi->x - pj->x).
 * @param hi Smoothing length of particle i.
 * @param hj Smoothing length of particle j.
 * @param pi Particle i.
 * @param pj Particle j.
 */
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {

  hydro_gradients_nonsym_collect(r2, dx, hi, hj, pi, pj);
}

/**
 * @brief Common part of the flux calculation between particle i and j
 *
 * Since the only difference between the symmetric and non-symmetric version
 * of the flux calculation  is in the update of the conserved variables at the
 * very end (which is not done for particle j if mode is 0 and particle j is
 * active), both runner_iact_force and runner_iact_nonsym_force call this
 * method, with an appropriate mode.
 *
115
116
117
118
119
 * This method retrieves the oriented surface area and face midpoint for the
 * Voronoi face between pi and pj (if it exists). It uses the midpoint position
 * to reconstruct the primitive quantities (if gradients are used) at the face
 * and then uses the face quantities to estimate a flux through the face using
 * a Riemann solver.
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
 *
 * This method also calculates the maximal velocity used to calculate the time
 * step.
 *
 * @param r2 Squared distance between particle i and particle j.
 * @param dx Distance vector between the particles (dx = pi->x - pj->x).
 * @param hi Smoothing length of particle i.
 * @param hj Smoothing length of particle j.
 * @param pi Particle i.
 * @param pj Particle j.
 */
__attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
    int mode) {

  float r = sqrtf(r2);
136
  int k;
137
  float A;
138
  float xij_i[3];
139
140
141
142
143
144
  float vmax, dvdotdx;
  float vi[3], vj[3], vij[3];
  float Wi[5], Wj[5];
  float dti, dtj, mindt;
  float n_unit[3];

145
  A = voronoi_get_face(&pi->cell, pj->id, xij_i);
146
  if (A == 0.0f) {
147
148
149
150
    /* this neighbour does not share a face with the cell, return */
    return;
  }

151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
  /* Initialize local variables */
  for (k = 0; k < 3; k++) {
    vi[k] = pi->force.v_full[k]; /* particle velocities */
    vj[k] = pj->force.v_full[k];
  }
  Wi[0] = pi->primitives.rho;
  Wi[1] = pi->primitives.v[0];
  Wi[2] = pi->primitives.v[1];
  Wi[3] = pi->primitives.v[2];
  Wi[4] = pi->primitives.P;
  Wj[0] = pj->primitives.rho;
  Wj[1] = pj->primitives.v[0];
  Wj[2] = pj->primitives.v[1];
  Wj[3] = pj->primitives.v[2];
  Wj[4] = pj->primitives.P;

  dti = pi->force.dt;
  dtj = pj->force.dt;

  /* calculate the maximal signal velocity */
171
172
173
  vmax = 0.0f;
  if (Wi[0] > 0.) {
    vmax += gas_soundspeed_from_pressure(Wi[0], Wi[4]);
174
  }
175
176
177
178
179

  if (Wj[0] > 0.) {
    vmax += gas_soundspeed_from_pressure(Wj[0], Wj[4]);
  }

180
181
182
183
184
  dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
            (Wi[3] - Wj[3]) * dx[2];
  if (dvdotdx > 0.) {
    vmax -= dvdotdx / r;
  }
185

186
187
188
189
190
191
192
193
194
195
  pi->timestepvars.vmax = fmaxf(pi->timestepvars.vmax, vmax);
  if (mode == 1) {
    pj->timestepvars.vmax = fmaxf(pj->timestepvars.vmax, vmax);
  }

  /* The flux will be exchanged using the smallest time step of the two
   * particles */
  mindt = fminf(dti, dtj);

  /* compute the normal vector of the interface */
196
197
198
  for (k = 0; k < 3; ++k) {
    n_unit[k] = -dx[k] / r;
  }
199
200

  /* Compute interface velocity */
201
202
203
204
205
206
207
  float fac = (vi[0] - vj[0]) * (xij_i[0] + 0.5f * dx[0]) +
              (vi[1] - vj[1]) * (xij_i[1] + 0.5f * dx[1]) +
              (vi[2] - vj[2]) * (xij_i[2] + 0.5f * dx[2]);
  fac /= r;
  vij[0] = 0.5f * (vi[0] + vj[0]) - fac * dx[0];
  vij[1] = 0.5f * (vi[1] + vj[1]) - fac * dx[1];
  vij[2] = 0.5f * (vi[2] + vj[2]) - fac * dx[2];
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
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
258
259
260
261
262
263
264

  /* Boost the primitive variables to the frame of reference of the interface */
  /* Note that velocities are indices 1-3 in W */
  Wi[1] -= vij[0];
  Wi[2] -= vij[1];
  Wi[3] -= vij[2];
  Wj[1] -= vij[0];
  Wj[2] -= vij[1];
  Wj[3] -= vij[2];

  hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj, mindt);

  /* we don't need to rotate, we can use the unit vector in the Riemann problem
   * itself (see GIZMO) */

  if (Wi[0] < 0.0f || Wj[0] < 0.0f || Wi[4] < 0.0f || Wj[4] < 0.0f) {
    printf("mindt: %g\n", mindt);
    printf("WL: %g %g %g %g %g\n", pi->primitives.rho, pi->primitives.v[0],
           pi->primitives.v[1], pi->primitives.v[2], pi->primitives.P);
#ifdef USE_GRADIENTS
    printf("dWL: %g %g %g %g %g\n", dWi[0], dWi[1], dWi[2], dWi[3], dWi[4]);
#endif
    printf("gradWL[0]: %g %g %g\n", pi->primitives.gradients.rho[0],
           pi->primitives.gradients.rho[1], pi->primitives.gradients.rho[2]);
    printf("gradWL[1]: %g %g %g\n", pi->primitives.gradients.v[0][0],
           pi->primitives.gradients.v[0][1], pi->primitives.gradients.v[0][2]);
    printf("gradWL[2]: %g %g %g\n", pi->primitives.gradients.v[1][0],
           pi->primitives.gradients.v[1][1], pi->primitives.gradients.v[1][2]);
    printf("gradWL[3]: %g %g %g\n", pi->primitives.gradients.v[2][0],
           pi->primitives.gradients.v[2][1], pi->primitives.gradients.v[2][2]);
    printf("gradWL[4]: %g %g %g\n", pi->primitives.gradients.P[0],
           pi->primitives.gradients.P[1], pi->primitives.gradients.P[2]);
    printf("WL': %g %g %g %g %g\n", Wi[0], Wi[1], Wi[2], Wi[3], Wi[4]);
    printf("WR: %g %g %g %g %g\n", pj->primitives.rho, pj->primitives.v[0],
           pj->primitives.v[1], pj->primitives.v[2], pj->primitives.P);
#ifdef USE_GRADIENTS
    printf("dWR: %g %g %g %g %g\n", dWj[0], dWj[1], dWj[2], dWj[3], dWj[4]);
#endif
    printf("gradWR[0]: %g %g %g\n", pj->primitives.gradients.rho[0],
           pj->primitives.gradients.rho[1], pj->primitives.gradients.rho[2]);
    printf("gradWR[1]: %g %g %g\n", pj->primitives.gradients.v[0][0],
           pj->primitives.gradients.v[0][1], pj->primitives.gradients.v[0][2]);
    printf("gradWR[2]: %g %g %g\n", pj->primitives.gradients.v[1][0],
           pj->primitives.gradients.v[1][1], pj->primitives.gradients.v[1][2]);
    printf("gradWR[3]: %g %g %g\n", pj->primitives.gradients.v[2][0],
           pj->primitives.gradients.v[2][1], pj->primitives.gradients.v[2][2]);
    printf("gradWR[4]: %g %g %g\n", pj->primitives.gradients.P[0],
           pj->primitives.gradients.P[1], pj->primitives.gradients.P[2]);
    printf("WR': %g %g %g %g %g\n", Wj[0], Wj[1], Wj[2], Wj[3], Wj[4]);
    error("Negative density or pressure!\n");
  }

  float totflux[5];
  riemann_solve_for_flux(Wi, Wj, n_unit, vij, totflux);

  /* Update conserved variables */
  /* eqn. (16) */
265
266
267
268
269
  pi->conserved.flux.mass -= mindt * A * totflux[0];
  pi->conserved.flux.momentum[0] -= mindt * A * totflux[1];
  pi->conserved.flux.momentum[1] -= mindt * A * totflux[2];
  pi->conserved.flux.momentum[2] -= mindt * A * totflux[3];
  pi->conserved.flux.energy -= mindt * A * totflux[4];
270

271
#ifndef SHADOWFAX_TOTAL_ENERGY
272
273
274
  float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
                       pi->primitives.v[1] * pi->primitives.v[1] +
                       pi->primitives.v[2] * pi->primitives.v[2]);
275
276
277
278
279
  pi->conserved.flux.energy += mindt * A * totflux[1] * pi->primitives.v[0];
  pi->conserved.flux.energy += mindt * A * totflux[2] * pi->primitives.v[1];
  pi->conserved.flux.energy += mindt * A * totflux[3] * pi->primitives.v[2];
  pi->conserved.flux.energy -= mindt * A * totflux[0] * ekin;
#endif
280
281
282
283
284
285
286
287
288
289
290
291
292

  /* here is how it works:
     Mode will only be 1 if both particles are ACTIVE and they are in the same
     cell. In this case, this method IS the flux calculation for particle j, and
     we HAVE TO UPDATE it.
     Mode 0 can mean several things: it can mean that particle j is INACTIVE, in
     which case we NEED TO UPDATE it, since otherwise the flux is lost from the
     system and the conserved variable is not conserved.
     It can also mean that particle j sits in another cell and is ACTIVE. In
     this case, the flux exchange for particle j is done TWICE and we SHOULD NOT
     UPDATE particle j.
     ==> we update particle j if (MODE IS 1) OR (j IS INACTIVE)
  */
293
  if (mode == 1 || pj->force.active == 0) {
294
295
296
297
298
    pj->conserved.flux.mass += mindt * A * totflux[0];
    pj->conserved.flux.momentum[0] += mindt * A * totflux[1];
    pj->conserved.flux.momentum[1] += mindt * A * totflux[2];
    pj->conserved.flux.momentum[2] += mindt * A * totflux[3];
    pj->conserved.flux.energy += mindt * A * totflux[4];
299

300
#ifndef SHADOWFAX_TOTAL_ENERGY
301
302
303
    ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
                   pj->primitives.v[1] * pj->primitives.v[1] +
                   pj->primitives.v[2] * pj->primitives.v[2]);
304
305
306
307
308
    pj->conserved.flux.energy -= mindt * A * totflux[1] * pj->primitives.v[0];
    pj->conserved.flux.energy -= mindt * A * totflux[2] * pj->primitives.v[1];
    pj->conserved.flux.energy -= mindt * A * totflux[3] * pj->primitives.v[2];
    pj->conserved.flux.energy += mindt * A * totflux[0] * ekin;
#endif
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
  }
}

/**
 * @brief Flux calculation between particle i and particle j
 *
 * This method calls runner_iact_fluxes_common with mode 1.
 *
 * @param r2 Squared distance between particle i and particle j.
 * @param dx Distance vector between the particles (dx = pi->x - pj->x).
 * @param hi Smoothing length of particle i.
 * @param hj Smoothing length of particle j.
 * @param pi Particle i.
 * @param pj Particle j.
 */
__attribute__((always_inline)) INLINE static void runner_iact_force(
    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {

  runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1);
}

/**
 * @brief Flux calculation between particle i and particle j: non-symmetric
 * version
 *
 * This method calls runner_iact_fluxes_common with mode 0.
 *
 * @param r2 Squared distance between particle i and particle j.
 * @param dx Distance vector between the particles (dx = pi->x - pj->x).
 * @param hi Smoothing length of particle i.
 * @param hj Smoothing length of particle j.
 * @param pi Particle i.
 * @param pj Particle j.
 */
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {

  runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0);
}