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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
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
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
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
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
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
|
#include <cmath>
#include <algorithm>
#include <vector>
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#define MIN(a,b) (((a)<(b))?(a):(b))
#define MAX(a,b) (((a)>(b))?(a):(b))
#define MIN3(a,b,c) (MIN(a, MIN(b, c)))
#define MIN4(a,b,c,d) (MIN(a, MIN(b, MIN(c, d))))
#define MAX5(a,b,c,d,e) (MAX(a, MAX(b, MAX(c, MAX(d, e)))))
using namespace std;
/**
ePPM Reconstruction (see Reisswig et al. 2013) based on McCorquodale & Colella (2011).
Note: This routine differs from the Fortran routine (and the algorithm described in Reisswig et al. 2013)
in the following way: We do not treat the specific internal energy "epsilon" in any special way.
Rather, to avoid reconstruction to eps<0 for gamma-laws, we rely on
'eps<0'-error checks after reconstruction (a similar treatment is also necessary for other reconstruction methods such as WENO5 and MP5).
*/
/*
Cases that must be considered:
* basic hydro
* hydro + temperature + ye
* hydro + ye
* basic mhd
* mhd + temperature + ye
* mhd + ye
* mppm (not supported right now)
* not supporting trivial_rp
* with or without divergence cleaning
*/
static inline double approx_at_cell_interface(const double* const restrict a, const int i) {
return 7.0/12.0*(a[i]+a[i+1]) - 1.0/12.0*(a[i-1]+a[i+2]);
}
static inline double limit(const double* const restrict a, const double ah, const double C, const int i)
{
if ((MIN(a[i],a[i+1]) <= ah) && (ah <= MAX(a[i],a[i+1]))) {
return ah;
} else {
const double D2a = 3.0 * ((a[i] + a[i+1]) - 2.0*ah );
const double D2aL = ((a[i-1] + a[i+1]) - 2.0*a[i] );
const double D2aR = ((a[i] + a[i+2]) - 2.0*a[i+1]);
const double D2aLim = copysign(1.0, D2a)*MIN3(C*fabs(D2aL), C*fabs(D2aR), fabs(D2a)) * 1.0/3.0;
if (D2a*D2aR >= 0 && D2a*D2aL >= 0)
return 0.5*(a[i]+a[i+1]) - D2aLim;
else
return 0.5*(a[i]+a[i+1]);
}
return 0;
}
/// Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34.
/// This does not use the check for deviations from a cubic. Thus, it gets away with only 3 stencil points!
static inline void monotonize(double* const restrict aminus,
const double* const restrict a,
double* const restrict aplus,
const double C,
const int i)
{
double D2aLim = 0;
double rhi = 0;
const double daplus = aplus[i]-a[i];
const double daminus = a[i]-aminus[i];
if (daplus*daminus <= 0 || (a[i-2]-a[i])*(a[i]-a[i+2]) <= 0) {
const double D2a = - (12.0*a[i] - 6.0*(aminus[i]+aplus[i]));
const double D2aC = (a[i-1] + a[i+1]) - 2.0*a[i];
const double D2aL = (a[i-2] + a[i] ) - 2.0*a[i-1];
const double D2aR = (a[i] + a[i+2]) - 2.0*a[i+1];
if (copysign(1.0, D2a) == copysign(1.0, D2aC) && copysign(1.0, D2a) == copysign(1.0, D2aL) && copysign(1.0, D2a) == copysign(1.0, D2aR))
D2aLim = copysign(1.0, D2a) * MIN4(C*fabs(D2aL), C*fabs(D2aR), C*fabs(D2aC), fabs(D2a));
if (!(fabs(D2a) <= 1e-12*MAX5(fabs(a[i-2]), fabs(a[i-1]), fabs(a[i]), fabs(a[i+1]), fabs(a[i+2]))))
rhi = D2aLim / D2a;
if (! (rhi >= 1.0 - 1e-12)) {
if (daplus*daminus < 0) {
aplus[i] = a[i] + daplus * rhi;
aminus[i] = a[i] - daminus * rhi;
} else if (fabs(daminus) >= 2.0*fabs(daplus)) {
aminus[i] = a[i] - (2.0*(1.0-rhi)*daplus + rhi*daminus);
} else if (fabs(daplus) >= 2.0*fabs(daminus)) {
aplus[i] = a[i] + (2.0*(1.0-rhi)*daminus + rhi*daplus);
}
}
//trivial_rp(i-1) = .false.
//trivial_rp(i) = .false.
} else {
if (fabs(daplus) >= 2.0*fabs(daminus))
aplus[i] = a[i] + 2.0*daminus;
else if (fabs(daminus) >= 2.0*fabs(daplus))
aminus[i] = a[i] - 2.0*daplus;
//trivial_rp(i-1) = .false.
//trivial_rp(i) = .false.
}
return;
}
template<bool do_temp, bool do_ye, bool do_mhd,
bool dc_flag, bool do_ppm_detect>
void GRHydro_eppm1d_cxx(const int nx,
const double dx,
const double* const restrict rho,
const double* const restrict velx,
const double* const restrict vely,
const double* const restrict velz,
const double* const restrict eps,
const double* const restrict press,
const double* const restrict temp,
const double* const restrict ye,
const double* const restrict Bvcx,
const double* const restrict Bvcy,
const double* const restrict Bvcz,
const double* const restrict psidc,
double* const restrict rhominus,
double* const restrict velxminus,
double* const restrict velyminus,
double* const restrict velzminus,
double* const restrict epsminus,
double* const restrict tempminus,
double* const restrict yeminus,
double* const restrict Bvcxminus,
double* const restrict Bvcyminus,
double* const restrict Bvczminus,
double* const restrict psidcminus,
double* const restrict rhoplus,
double* const restrict velxplus,
double* const restrict velyplus,
double* const restrict velzplus,
double* const restrict epsplus,
double* const restrict tempplus,
double* const restrict yeplus,
double* const restrict Bvcxplus,
double* const restrict Bvcyplus,
double* const restrict Bvczplus,
double* const restrict psidcplus)
{
DECLARE_CCTK_PARAMETERS;
// We initialize "plus" \equiv a_j+1/2 with (16) via APPROX_AT_CELL_INTERFACE,
// then checking for (13) of Colella & Sekora 2008 and applying
// (18) and (19) if (13) is not satisfied. This is done with LIMIT.
for (int i=1; i < nx-2; ++i) {
rhoplus[i] = approx_at_cell_interface(rho, i);
rhoplus[i] = limit(rho, rhoplus[i], enhanced_ppm_C2, i);
rhominus[i+1] = rhoplus[i];
epsplus[i] = approx_at_cell_interface(eps, i);
epsplus[i] = limit(eps, epsplus[i], enhanced_ppm_C2, i);
epsminus[i+1] = epsplus[i];
velxplus[i] = approx_at_cell_interface(velx, i);
velxplus[i] = limit(velx, velxplus[i], enhanced_ppm_C2, i);
velxminus[i+1] = velxplus[i];
velyplus[i] = approx_at_cell_interface(vely, i);
velyplus[i] = limit(vely, velyplus[i], enhanced_ppm_C2, i);
velyminus[i+1] = velyplus[i];
velzplus[i] = approx_at_cell_interface(velz, i);
velzplus[i] = limit(velz, velzplus[i], enhanced_ppm_C2, i);
velzminus[i+1] = velzplus[i];
if(do_ye) {
if(do_temp) {
tempplus[i] = approx_at_cell_interface(temp, i);
tempplus[i] = limit(temp, tempplus[i], enhanced_ppm_C2, i);
tempminus[i+1] = tempplus[i];
}
yeplus[i] = approx_at_cell_interface(ye, i);
yeplus[i] = limit(ye, yeplus[i], enhanced_ppm_C2, i);
yeminus[i+1] = yeplus[i];
}
if(do_mhd) {
Bvcxplus[i] = approx_at_cell_interface(Bvcx, i);
Bvcxplus[i] = limit(Bvcx, Bvcxplus[i], enhanced_ppm_C2, i);
Bvcxminus[i+1] = Bvcxplus[i];
Bvcyplus[i] = approx_at_cell_interface(Bvcy, i);
Bvcyplus[i] = limit(Bvcy, Bvcyplus[i], enhanced_ppm_C2, i);
Bvcyminus[i+1] = Bvcyplus[i];
Bvczplus[i] = approx_at_cell_interface(Bvcz, i);
Bvczplus[i] = limit(Bvcz, Bvczplus[i], enhanced_ppm_C2, i);
Bvczminus[i+1] = Bvczplus[i];
if(dc_flag) {
psidcplus[i] = approx_at_cell_interface(psidc, i);
psidcplus[i] = limit(psidc, psidcplus[i], enhanced_ppm_C2, i);
psidcminus[i+1] = psidcplus[i];
}
}
}
// Monotonize
for(int i=GRHydro_stencil-1; i<nx-GRHydro_stencil+1; ++i) {
monotonize(rhominus,rho,rhoplus,enhanced_ppm_C2,i);
monotonize(epsminus,eps,epsplus,enhanced_ppm_C2,i);
monotonize(velxminus,velx,velxplus,enhanced_ppm_C2,i);
monotonize(velyminus,vely,velyplus,enhanced_ppm_C2,i);
monotonize(velzminus,velz,velzplus,enhanced_ppm_C2,i);
if(do_ye) {
if(do_temp) {
monotonize(tempminus,temp,tempplus,enhanced_ppm_C2,i);
}
monotonize(yeminus,ye,yeplus,enhanced_ppm_C2,i);
}
if(do_mhd) {
monotonize(Bvcxminus,Bvcx,Bvcxplus,enhanced_ppm_C2,i);
monotonize(Bvcyminus,Bvcy,Bvcyplus,enhanced_ppm_C2,i);
monotonize(Bvczminus,Bvcz,Bvczplus,enhanced_ppm_C2,i);
if(dc_flag) {
monotonize(psidcminus,psidc,psidcplus,enhanced_ppm_C2,i);
}
}
}
// Finally, apply flattening
vector<double> dpress(nx,0);
for(int i=1; i<nx-1; ++i)
dpress[i] = (press[i+1]-press[i-1]);
vector<double> tilde_flatten(nx,0);
for(int i=2; i<nx-2; ++i) {
const double dpress2 = press[i+2] - press[i-2];
const double dvel = velx[i+1] - velx[i-1];
double w=0.0;
if ( (fabs(dpress[i]) > ppm_epsilon * MIN(press[i-1],press[i+1]))
&& (dvel < 0.0) )
{
w = 1.0;
}
if (fabs(dpress2) < ppm_small)
{
tilde_flatten[i] = 1.0;
}
else
{
tilde_flatten[i] = MAX(0.0, 1.0 - w * MAX(0.0,
ppm_omega2 * (dpress[i] / dpress2 - ppm_omega1)));
}
}
for(int i=2; i<nx-2; ++i) {
const double flatten = tilde_flatten[i];
rhoplus[i] = flatten * rhoplus[i] + (1.0 - flatten) * rho[i];
rhominus[i] = flatten * rhominus[i] + (1.0 - flatten) * rho[i];
epsplus[i] = flatten * epsplus[i] + (1.0 - flatten) * eps[i];
epsminus[i] = flatten * epsminus[i] + (1.0 - flatten) * eps[i];
velxplus[i] = flatten * velxplus[i] + (1.0 - flatten) * velx[i];
velxminus[i] = flatten * velxminus[i] + (1.0 - flatten) * velx[i];
velyplus[i] = flatten * velyplus[i] + (1.0 - flatten) * vely[i];
velyminus[i] = flatten * velyminus[i] + (1.0 - flatten) * vely[i];
velzplus[i] = flatten * velzplus[i] + (1.0 - flatten) * velz[i];
velzminus[i] = flatten * velzminus[i] + (1.0 - flatten) * velz[i];
if(do_ye) {
if(do_temp) {
tempplus[i] = flatten * tempplus[i] + (1.0 - flatten) * temp[i];
tempminus[i] = flatten * tempminus[i] + (1.0 - flatten) * temp[i];
}
yeplus[i] = flatten * yeplus[i] + (1.0 - flatten) * ye[i];
yeminus[i] = flatten * yeminus[i] + (1.0 - flatten) * ye[i];
}
if(do_mhd) {
Bvcxplus[i] = flatten * Bvcxplus[i] + (1.0 - flatten) * Bvcx[i];
Bvcxminus[i] = flatten * Bvcxminus[i] + (1.0 - flatten) * Bvcx[i];
Bvcyplus[i] = flatten * Bvcyplus[i] + (1.0 - flatten) * Bvcy[i];
Bvcyminus[i] = flatten * Bvcyminus[i] + (1.0 - flatten) * Bvcy[i];
Bvczplus[i] = flatten * Bvczplus[i] + (1.0 - flatten) * Bvcz[i];
Bvczminus[i] = flatten * Bvczminus[i] + (1.0 - flatten) * Bvcz[i];
if(dc_flag) {
psidcplus[i] = flatten * psidcplus[i] + (1.0 - flatten) * psidc[i];
psidcminus[i] = flatten * psidcminus[i] + (1.0 - flatten) * psidc[i];
}
}
} // flattening
return;
}
//template template<bool do_temp, bool do_ye, bool do_mhd, bool dc_flag, bool do_ppm_detect>
#define INSTANTIATE(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect) \
template void GRHydro_eppm1d_cxx<do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect>(const int nx, \
const double dx, \
const double* restrict rho, \
const double* restrict velx, \
const double* restrict vely, \
const double* restrict velz, \
const double* restrict eps, \
const double* restrict press, \
const double* restrict temp, \
const double* restrict ye, \
const double* restrict Bvcx, \
const double* restrict Bvcy, \
const double* restrict Bvcz, \
const double* restrict psidc, \
double* restrict rhominus, \
double* restrict velxminus, \
double* restrict velyminus, \
double* restrict velzminus, \
double* restrict epsminus, \
double* restrict tempminus, \
double* restrict yeminus, \
double* restrict Bvcxminus, \
double* restrict Bvcyminus, \
double* restrict Bvczminus, \
double* restrict psidcminus, \
double* restrict rhoplus, \
double* restrict velxplus, \
double* restrict velyplus, \
double* restrict velzplus, \
double* restrict epsplus, \
double* restrict tempplus, \
double* restrict yeplus, \
double* restrict Bvcxplus, \
double* restrict Bvcyplus, \
double* restrict Bvczplus, \
double* restrict psidcplus);
//template template<bool do_temp, bool do_ye, bool do_mhd, bool dc_flag, bool do_ppm_detect>
// first stuff without MHD
INSTANTIATE(false,false,false,false,false)
INSTANTIATE(false,false,false,false,true)
// do_temp can only be true of do_ye is also true
INSTANTIATE(true,true,false,false,false)
INSTANTIATE(true,true,false,false,true)
// but do_ye can be true and do_temp can be false
INSTANTIATE(false,true,false,false,false)
INSTANTIATE(false,true,false,false,true)
// with MHD, but dc_flag false
INSTANTIATE(false,false,true,false,false)
INSTANTIATE(false,false,true,false,true)
// do_temp can only be true of do_ye is also true
INSTANTIATE(true,true,true,false,false)
INSTANTIATE(true,true,true,false,true)
// but do_ye can be true and do_temp can be false
INSTANTIATE(false,true,true,false,false)
INSTANTIATE(false,true,true,false,true)
// with MHD, dc_flag true
INSTANTIATE(false,false,true,true,false)
INSTANTIATE(false,false,true,true,true)
// do_temp can only be true of do_ye is also true
INSTANTIATE(true,true,true,true,false)
INSTANTIATE(true,true,true,true,true)
// but do_ye can be true and do_temp can be false
INSTANTIATE(false,true,true,true,false)
INSTANTIATE(false,true,true,true,true)
|
