hydro_iact.h 18.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/*******************************************************************************
 * This file is part of SWIFT.
 * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *                    Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
 *
 * 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/>.
 *
 ******************************************************************************/
21
22
#ifndef SWIFT_GIZMO_HYDRO_IACT_H
#define SWIFT_GIZMO_HYDRO_IACT_H
23

24
#include "adiabatic_index.h"
25
#include "hydro_gradients.h"
26
27
#include "riemann.h"

28
29
#define GIZMO_VOLUME_CORRECTION

30
31
32
33
34
35
36
37
38
39
/**
 * @brief Calculate the volume interaction between particle i and particle j
 *
 * The volume is in essence the same as the weighted number of neighbours in a
 * classical SPH density calculation.
 *
 * We also calculate the components of the matrix E, which is used for second
 * order accurate gradient calculations and for the calculation of the interface
 * surface areas.
 *
40
41
42
43
44
 * @param r2 Comoving squared distance between particle i and particle j.
 * @param dx Comoving distance vector between the particles (dx = pi->x -
 * pj->x).
 * @param hi Comoving smoothing-length of particle i.
 * @param hj Comoving smoothing-length of particle j.
45
46
 * @param pi Particle i.
 * @param pj Particle j.
47
48
 * @param a Current scale factor.
 * @param H Current Hubble parameter.
49
 */
50
__attribute__((always_inline)) INLINE static void runner_iact_density(
51
52
    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
    struct part *restrict pj, float a, float H) {
53
54
55
56
57
58
59
60
61
62
63
64
65

  float r = sqrtf(r2);
  float xi, xj;
  float h_inv;
  float wi, wj, wi_dx, wj_dx;
  int k, l;

  /* Compute density of pi. */
  h_inv = 1.0 / hi;
  xi = r * h_inv;
  kernel_deval(xi, &wi, &wi_dx);

  pi->density.wcount += wi;
66
  pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
67
68
69
70
71
72

  /* these are eqns. (1) and (2) in the summary */
  pi->geometry.volume += wi;
  for (k = 0; k < 3; k++)
    for (l = 0; l < 3; l++) pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;

73
74
75
76
  pi->geometry.centroid[0] -= dx[0] * wi;
  pi->geometry.centroid[1] -= dx[1] * wi;
  pi->geometry.centroid[2] -= dx[2] * wi;

77
78
79
80
81
82
  /* Compute density of pj. */
  h_inv = 1.0 / hj;
  xj = r * h_inv;
  kernel_deval(xj, &wj, &wj_dx);

  pj->density.wcount += wj;
83
  pj->density.wcount_dh -= (hydro_dimension * wj + xj * wj_dx);
84
85
86
87
88

  /* these are eqns. (1) and (2) in the summary */
  pj->geometry.volume += wj;
  for (k = 0; k < 3; k++)
    for (l = 0; l < 3; l++) pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
89
90
91
92

  pj->geometry.centroid[0] += dx[0] * wj;
  pj->geometry.centroid[1] += dx[1] * wj;
  pj->geometry.centroid[2] += dx[2] * wj;
93
94
}

95
96
97
98
99
100
101
102
103
104
105
/**
 * @brief Calculate the volume interaction between particle i and particle j:
 * non-symmetric version
 *
 * The volume is in essence the same as the weighted number of neighbours in a
 * classical SPH density calculation.
 *
 * We also calculate the components of the matrix E, which is used for second
 * order accurate gradient calculations and for the calculation of the interface
 * surface areas.
 *
106
107
108
109
110
 * @param r2 Comoving squared distance between particle i and particle j.
 * @param dx Comoving distance vector between the particles (dx = pi->x -
 * pj->x).
 * @param hi Comoving smoothing-length of particle i.
 * @param hj Comoving smoothing-length of particle j.
111
112
 * @param pi Particle i.
 * @param pj Particle j.
113
114
 * @param a Current scale factor.
 * @param H Current Hubble parameter.
115
 */
116
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
117
118
    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
    const struct part *restrict pj, float a, float H) {
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

  float r;
  float xi;
  float h_inv;
  float wi, wi_dx;
  int k, l;

  /* Get r and r inverse. */
  r = sqrtf(r2);

  h_inv = 1.0 / hi;
  xi = r * h_inv;
  kernel_deval(xi, &wi, &wi_dx);

  pi->density.wcount += wi;
134
  pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
135
136
137
138
139

  /* these are eqns. (1) and (2) in the summary */
  pi->geometry.volume += wi;
  for (k = 0; k < 3; k++)
    for (l = 0; l < 3; l++) pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
140
141
142
143

  pi->geometry.centroid[0] -= dx[0] * wi;
  pi->geometry.centroid[1] -= dx[1] * wi;
  pi->geometry.centroid[2] -= dx[2] * wi;
144
145
}

146
147
148
149
150
151
/**
 * @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.
 *
152
153
154
155
156
 * @param r2 Comoving squared distance between particle i and particle j.
 * @param dx Comoving distance vector between the particles (dx = pi->x -
 * pj->x).
 * @param hi Comoving smoothing-length of particle i.
 * @param hj Comoving smoothing-length of particle j.
157
158
 * @param pi Particle i.
 * @param pj Particle j.
159
160
 * @param a Current scale factor.
 * @param H Current Hubble parameter.
161
 */
162
__attribute__((always_inline)) INLINE static void runner_iact_gradient(
163
164
    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
    struct part *restrict pj, float a, float H) {
165

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

169
170
171
172
173
174
175
/**
 * @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.
 *
176
177
178
179
180
 * @param r2 Comoving squared distance between particle i and particle j.
 * @param dx Comoving distance vector between the particles (dx = pi->x -
 * pj->x).
 * @param hi Comoving smoothing-length of particle i.
 * @param hj Comoving smoothing-length of particle j.
181
182
 * @param pi Particle i.
 * @param pj Particle j.
183
184
 * @param a Current scale factor.
 * @param H Current Hubble parameter.
185
 */
186
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
187
188
    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
    struct part *restrict pj, float a, float H) {
189

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

193
194
195
196
197
/**
 * @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
198
199
200
 * very end (which is not done for particle j if mode is 0), both
 * runner_iact_force and runner_iact_nonsym_force call this method, with an
 * appropriate mode.
201
202
203
204
205
206
207
208
209
210
 *
 * This method calculates the surface area of the interface between particle i
 * and particle j, as well as the interface position and velocity. These are
 * then used to reconstruct and predict the primitive variables, which are then
 * fed to a Riemann solver that calculates a flux. This flux is used to update
 * the conserved variables of particle i or both particles.
 *
 * This method also calculates the maximal velocity used to calculate the time
 * step.
 *
211
212
213
214
215
 * @param r2 Comoving squared distance between particle i and particle j.
 * @param dx Comoving distance vector between the particles (dx = pi->x -
 * pj->x).
 * @param hi Comoving smoothing-length of particle i.
 * @param hj Comoving smoothing-length of particle j.
216
217
 * @param pi Particle i.
 * @param pj Particle j.
218
219
 * @param a Current scale factor.
 * @param H Current Hubble parameter.
220
 */
221
__attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
222
223
    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
    struct part *restrict pj, int mode, float a, float H) {
224
225
226

  float r = sqrtf(r2);
  float xi, xj;
227
228
  float hi_inv, hi_inv_dim;
  float hj_inv, hj_inv_dim;
229
230
231
232
233
234
235
236
  float wi, wj, wi_dx, wj_dx;
  int k, l;
  float A[3];
  float Anorm;
  float Bi[3][3];
  float Bj[3][3];
  float Vi, Vj;
  float xij_i[3], xfac, xijdotdx;
237
  float vmax, dvdotdx;
238
  float vi[3], vj[3], vij[3];
239
  float Wi[5], Wj[5];
240
241
242
243
244
245
246
247
  float n_unit[3];

  /* Initialize local variables */
  for (k = 0; k < 3; k++) {
    for (l = 0; l < 3; l++) {
      Bi[k][l] = pi->geometry.matrix_E[k][l];
      Bj[k][l] = pj->geometry.matrix_E[k][l];
    }
248
249
    vi[k] = pi->v[k]; /* particle velocities */
    vj[k] = pj->v[k];
250
251
252
253
254
255
256
257
258
259
260
261
262
  }
  Vi = pi->geometry.volume;
  Vj = pj->geometry.volume;
  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;
263

264
  /* calculate the maximal signal velocity */
265
  if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
266
267
268
269
270
    vmax =
        sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
  } else {
    vmax = 0.0f;
  }
271
272
  dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
            (Wi[3] - Wj[3]) * dx[2];
273
274
  dvdotdx = min(dvdotdx, (vi[0] - vj[0]) * dx[0] + (vi[1] - vj[1]) * dx[1] +
                             (vi[2] - vj[2]) * dx[2]);
275
  if (dvdotdx < 0.) {
276
    /* the magical factor 3 also appears in Gadget2 */
277
    vmax -= 3. * dvdotdx / r;
278
  }
279
  pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
280
  if (mode == 1) {
281
    pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
282
283
284
285
  }

  /* Compute kernel of pi. */
  hi_inv = 1.0 / hi;
286
  hi_inv_dim = pow_dimension(hi_inv);
287
288
289
290
291
  xi = r * hi_inv;
  kernel_deval(xi, &wi, &wi_dx);

  /* Compute kernel of pj. */
  hj_inv = 1.0 / hj;
292
  hj_inv_dim = pow_dimension(hj_inv);
293
294
295
  xj = r * hj_inv;
  kernel_deval(xj, &wj, &wj_dx);

296
  /* Compute h_dt. We are going to use an SPH-like estimate of div_v for that */
297
298
299
300
301
302
303
304
  float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
               (pi->v[2] - pj->v[2]) * dx[2];
  float ri = 1.0f / r;
  float hidp1 = pow_dimension_plus_one(hi_inv);
  float hjdp1 = pow_dimension_plus_one(hj_inv);
  float wi_dr = hidp1 * wi_dx;
  float wj_dr = hjdp1 * wj_dx;
  dvdr *= ri;
305
306
307
308
  if (pj->primitives.rho > 0.) {
    pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
  }
  if (mode == 1 && pi->primitives.rho > 0.) {
309
    pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
310
  }
311

312
313
314
  /* Compute area */
  /* eqn. (7) */
  Anorm = 0.0f;
315
316
317
318
319
320
321
322
323
324
  if (pi->density.wcorr > const_gizmo_min_wcorr &&
      pj->density.wcorr > const_gizmo_min_wcorr) {
    /* in principle, we use Vi and Vj as weights for the left and right
       contributions to the generalized surface vector.
       However, if Vi and Vj are very different (because they have very
       different
       smoothing lengths), then the expressions below are more stable. */
    float Xi = Vi;
    float Xj = Vj;
#ifdef GIZMO_VOLUME_CORRECTION
325
    if (fabsf(Vi - Vj) / min(Vi, Vj) > 1.5 * hydro_dimension) {
326
327
328
329
330
331
332
      Xi = (Vi * hj + Vj * hi) / (hi + hj);
      Xj = Xi;
    }
#endif
    for (k = 0; k < 3; k++) {
      /* we add a minus sign since dx is pi->x - pj->x */
      A[k] = -Xi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) *
333
                 wi * hi_inv_dim -
334
             Xj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) *
335
                 wj * hj_inv_dim;
336
337
338
339
340
341
342
343
344
      Anorm += A[k] * A[k];
    }
  } else {
    /* ill condition gradient matrix: revert to SPH face area */
    Anorm = -(hidp1 * Vi * Vi * wi_dx + hjdp1 * Vj * Vj * wj_dx) * ri;
    A[0] = -Anorm * dx[0];
    A[1] = -Anorm * dx[1];
    A[2] = -Anorm * dx[2];
    Anorm *= Anorm * r2;
345
346
  }

347
  if (Anorm == 0.) {
348
349
350
351
352
353
    /* if the interface has no area, nothing happens and we return */
    /* continuing results in dividing by zero and NaN's... */
    return;
  }

  Anorm = sqrtf(Anorm);
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371

#ifdef SWIFT_DEBUG_CHECKS
  /* For stability reasons, we do require A and dx to have opposite
     directions (basically meaning that the surface normal for the surface
     always points from particle i to particle j, as it would in a real
     moving-mesh code). If not, our scheme is no longer upwind and hence can
     become unstable. */
  float dA_dot_dx = A[0] * dx[0] + A[1] * dx[1] + A[2] * dx[2];
  /* In GIZMO, Phil Hopkins reverts to an SPH integration scheme if this
     happens. We curently just ignore this case and display a message. */
  const float rdim = pow_dimension(r);
  if (dA_dot_dx > 1.e-6 * rdim) {
    message("Ill conditioned gradient matrix (%g %g %g %g %g)!", dA_dot_dx,
            Anorm, Vi, Vj, r);
  }
#endif

  /* compute the normal vector of the interface */
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
  for (k = 0; k < 3; k++) n_unit[k] = A[k] / Anorm;

  /* Compute interface position (relative to pi, since we don't need the actual
   * position) */
  /* eqn. (8) */
  xfac = hi / (hi + hj);
  for (k = 0; k < 3; k++) xij_i[k] = -xfac * dx[k];

  /* Compute interface velocity */
  /* eqn. (9) */
  xijdotdx = xij_i[0] * dx[0] + xij_i[1] * dx[1] + xij_i[2] * dx[2];
  for (k = 0; k < 3; k++) vij[k] = vi[k] + (vi[k] - vj[k]) * xijdotdx / r2;

  /* complete calculation of position of interface */
  /* NOTE: dx is not necessarily just pi->x - pj->x but can also contain
           correction terms for periodicity. If we do the interpolation,
           we have to use xij w.r.t. the actual particle.
           => we need a separate xij for pi and pj... */
  /* tldr: we do not need the code below, but we do need the same code as above
     but then
     with i and j swapped */
  //    for ( k = 0 ; k < 3 ; k++ )
  //      xij[k] += pi->x[k];

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

405
  hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
406
407
408
409

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

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

413
414
415
416
417
418
  /* Multiply with the interface surface area */
  totflux[0] *= Anorm;
  totflux[1] *= Anorm;
  totflux[2] *= Anorm;
  totflux[3] *= Anorm;
  totflux[4] *= Anorm;
419

420
  /* Store mass flux */
421
  float mflux = totflux[0];
422
423
424
425
  pi->gravity.mflux[0] += mflux * dx[0];
  pi->gravity.mflux[1] += mflux * dx[1];
  pi->gravity.mflux[2] += mflux * dx[2];

426
427
  /* Update conserved variables */
  /* eqn. (16) */
428
429
430
431
432
  pi->conserved.flux.mass -= totflux[0];
  pi->conserved.flux.momentum[0] -= totflux[1];
  pi->conserved.flux.momentum[1] -= totflux[2];
  pi->conserved.flux.momentum[2] -= totflux[3];
  pi->conserved.flux.energy -= totflux[4];
433

434
#ifndef GIZMO_TOTAL_ENERGY
435
436
437
  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]);
438
439
440
441
  pi->conserved.flux.energy += totflux[1] * pi->primitives.v[0];
  pi->conserved.flux.energy += totflux[2] * pi->primitives.v[1];
  pi->conserved.flux.energy += totflux[3] * pi->primitives.v[2];
  pi->conserved.flux.energy -= totflux[0] * ekin;
442
#endif
443

444
445
446
447
448
449
  /* Note that this used to be much more complicated in early implementations of
   * the GIZMO scheme, as we wanted manifest conservation of conserved variables
   * and had to do symmetric flux exchanges. Now we don't care about manifest
   * conservation anymore and just assume the current fluxes are representative
   * for the flux over the entire time step. */
  if (mode == 1) {
450
    /* Store mass flux */
451
    mflux = totflux[0];
452
453
454
455
    pj->gravity.mflux[0] -= mflux * dx[0];
    pj->gravity.mflux[1] -= mflux * dx[1];
    pj->gravity.mflux[2] -= mflux * dx[2];

456
457
458
459
460
    pj->conserved.flux.mass += totflux[0];
    pj->conserved.flux.momentum[0] += totflux[1];
    pj->conserved.flux.momentum[1] += totflux[2];
    pj->conserved.flux.momentum[2] += totflux[3];
    pj->conserved.flux.energy += totflux[4];
461

462
#ifndef GIZMO_TOTAL_ENERGY
463
464
465
    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]);
466
467
468
469
    pj->conserved.flux.energy -= totflux[1] * pj->primitives.v[0];
    pj->conserved.flux.energy -= totflux[2] * pj->primitives.v[1];
    pj->conserved.flux.energy -= totflux[3] * pj->primitives.v[2];
    pj->conserved.flux.energy += totflux[0] * ekin;
470
#endif
471
472
473
  }
}

474
475
476
477
478
/**
 * @brief Flux calculation between particle i and particle j
 *
 * This method calls runner_iact_fluxes_common with mode 1.
 *
479
480
481
482
483
 * @param r2 Comoving squared distance between particle i and particle j.
 * @param dx Comoving distance vector between the particles (dx = pi->x -
 * pj->x).
 * @param hi Comoving smoothing-length of particle i.
 * @param hj Comoving smoothing-length of particle j.
484
485
 * @param pi Particle i.
 * @param pj Particle j.
486
487
 * @param a Current scale factor.
 * @param H Current Hubble parameter.
488
 */
489
__attribute__((always_inline)) INLINE static void runner_iact_force(
490
491
    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
    struct part *restrict pj, float a, float H) {
492

493
  runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1, a, H);
494
495
}

496
497
498
499
500
501
/**
 * @brief Flux calculation between particle i and particle j: non-symmetric
 * version
 *
 * This method calls runner_iact_fluxes_common with mode 0.
 *
502
503
504
505
506
 * @param r2 Comoving squared distance between particle i and particle j.
 * @param dx Comoving distance vector between the particles (dx = pi->x -
 * pj->x).
 * @param hi Comoving smoothing-length of particle i.
 * @param hj Comoving smoothing-length of particle j.
507
508
 * @param pi Particle i.
 * @param pj Particle j.
509
510
 * @param a Current scale factor.
 * @param H Current Hubble parameter.
511
 */
512
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
513
514
    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
    struct part *restrict pj, float a, float H) {
515

516
  runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0, a, H);
517
}
518
519

#endif /* SWIFT_GIZMO_HYDRO_IACT_H */