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

23
/* Includes. */
24
#include "const.h"
25
#include "kernel.h"
26
#include "part.h"
27
28
29
30
31
#include "vector.h"

/**
 * @brief SPH interaction functions following the Gadget-2 version of SPH.
 *
32
33
34
35
36
37
 * The interactions computed here are the ones presented in the Gadget-2 paper
 *and use the same
 * numerical coefficients as the Gadget-2 code. When used with the Spline-3
 *kernel, the results
 * should be equivalent to the ones obtained with Gadget-2 up to the rounding
 *errors and interactions
38
39
 * missed by the Gadget-2 tree-code neighbours search.
 *
40
41
 * The code uses internal energy instead of entropy as a thermodynamical
 *variable.
42
43
44
45
46
47
 */

/**
 * @brief Density loop
 */

48
49
50
__attribute__((always_inline)) INLINE static void runner_iact_density(
    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {

51
52
  float wi, wi_dx;
  float wj, wj_dx;
53
  float dv[3], curlvr[3];
54

55
  /* Get the masses. */
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
  const float mi = pj->mass;
  const float mj = pj->mass;

  /* Get r and r inverse. */
  const float r = sqrtf(r2);
  const float r_inv = 1.0f / r;

  /* Compute the kernel function for pi */
  const float hi_inv = 1.f / hi;
  const float ui = r * hi_inv;
  kernel_deval(ui, &wi, &wi_dx);

  /* Compute contribution to the density */
  pi->rho += mj * wi;
  pi->rho_dh -= mj * kernel_igamma * (3.f * wi + ui * wi_dx);
71

72
73
74
75
76
77
78
79
80
81
82
83
  /* Compute contribution to the number of neighbours */
  pi->density.wcount += wi;
  pi->density.wcount_dh -= ui * wi_dx;

  /* Compute the kernel function for pj */
  const float hj_inv = 1.f / hj;
  const float uj = r * hj_inv;
  kernel_deval(uj, &wj, &wj_dx);

  /* Compute contribution to the density */
  pj->rho += mi * wj;
  pj->rho_dh -= mi * kernel_igamma * (3.f * wj + uj * wj_dx);
84

85
86
87
  /* Compute contribution to the number of neighbours */
  pj->density.wcount += wj;
  pj->density.wcount_dh -= uj * wj_dx;
88

89
90
  const float faci = mj * wi_dx * r_inv;
  const float facj = mi * wj_dx * r_inv;
91

92
93
94
95
  /* Compute dv dot r */
  dv[0] = pi->v[0] - pj->v[0];
  dv[1] = pi->v[1] - pj->v[1];
  dv[2] = pi->v[2] - pj->v[2];
96
97
  const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];

98
99
  pi->div_v += faci * dvdr;
  pj->div_v += facj * dvdr;
100
101
102
103
104
105

  /* Compute dv cross r */
  curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
  curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
  curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];

106
107
108
  pi->rot_v[0] += faci * curlvr[0];
  pi->rot_v[1] += faci * curlvr[1];
  pi->rot_v[2] += faci * curlvr[2];
109

110
111
112
  pj->rot_v[0] += facj * curlvr[0];
  pj->rot_v[1] += facj * curlvr[1];
  pj->rot_v[2] += facj * curlvr[2];
113
114
}

115
116
117
118
/**
 * @brief Density loop (non-symmetric version)
 */

119
120
121
122
123
124
125
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {

  float wi, wi_dx;
  float dv[3], curlvr[3];

  /* Get the masses. */
126
  const float mj = pj->mass;
127
128

  /* Get r and r inverse. */
129
130
  const float r = sqrtf(r2);
  const float ri = 1.0f / r;
131

132
133
134
135
136
137
138
139
140
141
142
143
  /* Compute the kernel function */
  const float h_inv = 1.0f / hi;
  const float u = r * h_inv;
  kernel_deval(u, &wi, &wi_dx);

  /* Compute contribution to the density */
  pi->rho += mj * wi;
  pi->rho_dh -= mj * kernel_igamma * (3.f * wi + u * wi_dx);

  /* Compute contribution to the number of neighbours */
  pi->density.wcount += wi;
  pi->density.wcount_dh -= u * wi_dx;
144

Matthieu Schaller's avatar
Matthieu Schaller committed
145
146
  /* const float ih3 = h_inv * h_inv * h_inv; */
  /* const float ih4 = h_inv * h_inv * h_inv * h_inv; */
147

148
  const float fac = mj * wi_dx * ri;
149

150
151
152
153
154
155
156
  /* Compute dv dot r */
  dv[0] = pi->v[0] - pj->v[0];
  dv[1] = pi->v[1] - pj->v[1];
  dv[2] = pi->v[2] - pj->v[2];
  const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
  pi->div_v -= fac * dvdr;

Matthieu Schaller's avatar
Matthieu Schaller committed
157
  /* if(pi->id == 515050 && pj->id == 504849) */
158
159
  /*   message("Interacting with %lld. r=%e hi=%e u=%e W=%e dW/dx=%e dh_drho1=%e
   * dh_drho2=%e\n fac=%e dvdr=%e pj->v=[%.3e,%.3e,%.3e]", */
Matthieu Schaller's avatar
Matthieu Schaller committed
160
161
162
163
164
165
166
167
168
169
170
171
172
173
  /* 	    pj->id, */
  /* 	    r, */
  /* 	    hi, */
  /* 	    u, */
  /* 	    wi * ih3, */
  /* 	    wi_dx * ih4, */
  /* 	    -mj * (3.f * kernel_igamma * wi) * ih4, */
  /* 	    -mj * u * wi_dx * kernel_igamma * ih4, */
  /* 	    fac * ih4, */
  /* 	    dvdr, */
  /* 	    pj->v[0], */
  /* 	    pj->v[1], */
  /* 	    pj->v[2] */
  /* 	    ); */
174

175
176
177
178
179
  /* Compute dv cross r */
  curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
  curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
  curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];

180
181
  pi->rot_v[0] += fac * curlvr[0];
  pi->rot_v[1] += fac * curlvr[1];
182
  pi->rot_v[2] += fac * curlvr[2];
183
184
}

185
186
187
188
/**
 * @brief Force loop
 */

189
190
191
__attribute__((always_inline)) INLINE static void runner_iact_force(
    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {

192
193
194
  float wi, wj, wi_dx, wj_dx;

  const float fac_mu = 1.f; /* Will change with cosmological integration */
195

196
197
198
199
  const float r = sqrtf(r2);
  const float r_inv = 1.0f / r;

  /* Get some values in local variables. */
200
  // const float mi = pi->mass;
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
  const float mj = pj->mass;
  const float rhoi = pi->rho;
  const float rhoj = pj->rho;
  const float pressurei = pi->pressure;
  const float pressurej = pj->pressure;

  /* Get the kernel for hi. */
  const float hi_inv = 1.0f / hi;
  const float hi2_inv = hi_inv * hi_inv;
  const float ui = r * hi_inv;
  kernel_deval(ui, &wi, &wi_dx);
  const float wi_dr = hi2_inv * hi2_inv * wi_dx;

  /* Get the kernel for hj. */
  const float hj_inv = 1.0f / hj;
  const float hj2_inv = hj_inv * hj_inv;
  const float xj = r * hj_inv;
  kernel_deval(xj, &wj, &wj_dx);
  const float wj_dr = hj2_inv * hj2_inv * wj_dx;

  /* Compute gradient terms */
  const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
  const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;

  /* Compute sound speeds */
  const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
  const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
  float v_sig = ci + cj;
229

230
  /* Compute dv dot r. */
231
232
233
  const 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];
234
235
236
237
238
239
240

  /* Artificial viscosity term */
  float visc = 0.f;
  if (dvdr < 0.f) {
    const float mu_ij = fac_mu * dvdr * r_inv;
    v_sig -= 3.f * mu_ij;
    const float rho_ij = 0.5f * (rhoi + rhoj);
241
242
243
244
245
246
    const float balsara_i = fabsf(pi->div_v) / (fabsf(pi->div_v) + pi->curl_v +
                                                0.0001 * ci / fac_mu / hi);
    const float balsara_j = fabsf(pj->div_v) / (fabsf(pj->div_v) + pj->curl_v +
                                                0.0001 * cj / fac_mu / hj);
    visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij / rho_ij *
           (balsara_i + balsara_j);
247
248
249
250
  }

  /* Now, convolve with the kernel */
  const float visc_term = 0.5f * mj * visc * (wi_dr + wj_dr) * r_inv;
251
252
  const float sph_term =
      mj * (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
253
254
255
256

  /* Eventually got the acceleration */
  const float acc = visc_term + sph_term;

Matthieu Schaller's avatar
Matthieu Schaller committed
257
258
  /* //if(pi->id == 1000 && pj->id == 1100) */
  /* if(pi->id == 515050 && pj->id == 504849) */
259
260
  /*   message("Interacting with %lld. r=%e hi=%e hj=%e dWi/dx=%e dWj/dx=%3e
   * dvdr=%e visc=%e sph=%e", */
Matthieu Schaller's avatar
Matthieu Schaller committed
261
262
263
264
265
266
267
268
269
270
271
272
  /* 	    pj->id, */
  /* 	    r, */
  /* 	    2*hi, */
  /* 	    2*hj, */
  /* 	    wi_dr, */
  /* 	    wj_dr, */
  /* 	    dvdr, */
  /* 	    visc_term, */
  /* 	    sph_term */
  /* 	    ); */
  /* if(pi->id == 1100 && pj->id == 1000) */
  /*   message("oO"); */
273

274
275
276
277
278
279
280
281
282
283
  /* Use the force Luke ! */
  pi->a[0] -= acc * dx[0];
  pi->a[1] -= acc * dx[1];
  pi->a[2] -= acc * dx[2];

  pj->a[0] += acc * dx[0];
  pj->a[1] += acc * dx[1];
  pj->a[2] += acc * dx[2];

  /* Update the signal velocity. */
284
285
286
  pi->v_sig = fmaxf(pi->v_sig, v_sig);
  pj->v_sig = fmaxf(pj->v_sig, v_sig);

287
  /* Change in entropy */
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
288
289
  pi->entropy_dt += 0.5f * visc_term * dvdr;
  pj->entropy_dt -= 0.5f * visc_term * dvdr;
290
}
291
292
293
294
295

/**
 * @brief Force loop (non-symmetric version)
 */

296
297
298
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {

299
300
301
  float wi, wj, wi_dx, wj_dx;

  const float fac_mu = 1.f; /* Will change with cosmological integration */
302

303
304
305
306
  const float r = sqrtf(r2);
  const float r_inv = 1.0f / r;

  /* Get some values in local variables. */
307
  // const float mi = pi->mass;
308
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
  const float mj = pj->mass;
  const float rhoi = pi->rho;
  const float rhoj = pj->rho;
  const float pressurei = pi->pressure;
  const float pressurej = pj->pressure;

  /* Get the kernel for hi. */
  const float hi_inv = 1.0f / hi;
  const float hi2_inv = hi_inv * hi_inv;
  const float ui = r * hi_inv;
  kernel_deval(ui, &wi, &wi_dx);
  const float wi_dr = hi2_inv * hi2_inv * wi_dx;

  /* Get the kernel for hj. */
  const float hj_inv = 1.0f / hj;
  const float hj2_inv = hj_inv * hj_inv;
  const float xj = r * hj_inv;
  kernel_deval(xj, &wj, &wj_dx);
  const float wj_dr = hj2_inv * hj2_inv * wj_dx;

  /* Compute gradient terms */
  const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
  const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;

  /* Compute sound speeds */
  const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
  const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
  float v_sig = ci + cj;
336

337
  /* Compute dv dot r. */
338
339
340
  const 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];
341
342
343
344
345
346
347

  /* Artificial viscosity term */
  float visc = 0.f;
  if (dvdr < 0.f) {
    const float mu_ij = fac_mu * dvdr * r_inv;
    v_sig -= 3.f * mu_ij;
    const float rho_ij = 0.5f * (rhoi + rhoj);
348
349
350
351
352
353
    const float balsara_i = fabsf(pi->div_v) / (fabsf(pi->div_v) + pi->curl_v +
                                                0.0001 * ci / fac_mu / hi);
    const float balsara_j = fabsf(pj->div_v) / (fabsf(pj->div_v) + pj->curl_v +
                                                0.0001 * cj / fac_mu / hj);
    visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij / rho_ij *
           (balsara_i + balsara_j);
354
355
356
357
  }

  /* Now, convolve with the kernel */
  const float visc_term = 0.5f * mj * visc * (wi_dr + wj_dr) * r_inv;
358
359
  const float sph_term =
      mj * (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
360
361
362

  /* Eventually got the acceleration */
  const float acc = visc_term + sph_term;
363

364
365
366
367
368
369
  /* Use the force Luke ! */
  pi->a[0] -= acc * dx[0];
  pi->a[1] -= acc * dx[1];
  pi->a[2] -= acc * dx[2];

  /* Update the signal velocity. */
370
371
  pi->v_sig = fmaxf(pi->v_sig, v_sig);

372
  /* Change in entropy */
Matthieu Schaller's avatar
Done !    
Matthieu Schaller committed
373
  pi->entropy_dt += 0.5f * visc_term * dvdr;
374
}
375

376
#endif /* SWIFT_RUNNER_IACT_LEGACY_H */