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
368
369
370
371
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
|
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#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))
// some prototypes
CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH);
static inline double SpatialDeterminantC(double gxx, double gxy,
double gxz, double gyy,
double gyz, double gzz);
static inline __attribute__((always_inline)) void prim2conC(double *w, double *dens, double *sx,
double *sy, double *sz, double *tau,
const double gxx, const double gxy,
const double gxz, const double gyy, const double gyz,
const double gzz, const double sdet,
const double rho, const double vx,
const double vy, const double vz,
const double eps, const double press);
static inline void GRHydro_SpeedOfSound(CCTK_ARGUMENTS);
void Primitive2ConservativeC(CCTK_ARGUMENTS);
CCTK_FCALL void CCTK_FNAME(Primitive2ConservativeCforF)(cGH ** p_cctkGH) {
Primitive2ConservativeC(*p_cctkGH);
}
CCTK_FCALL void CCTK_FNAME(GRHydro_SpeedOfSound)(cGH ** p_cctkGH) {
GRHydro_SpeedOfSound(*p_cctkGH);
}
void Primitive2ConservativeC(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
// save memory when multipatch is not used
CCTK_REAL * restrict g11, * restrict g12, * restrict g13;
CCTK_REAL * restrict g22, * restrict g23, * restrict g33;
if(GRHydro_UseGeneralCoordinates(cctkGH)) {
g11 = gaa;
g12 = gab;
g13 = gac;
g22 = gbb;
g23 = gbc;
g33 = gcc;
} else {
g11 = gxx;
g12 = gxy;
g13 = gxz;
g22 = gyy;
g23 = gyz;
g33 = gzz;
}
GRHydro_SpeedOfSound(CCTK_PASS_CTOC);
// if padding is used the the "extra" vector elements could contain junk
// which we do not know how to handle
assert(cctk_lsh[0] == cctk_ash[0]);
assert(cctk_lsh[1] == cctk_ash[1]);
assert(cctk_lsh[2] == cctk_ash[2]);
#pragma omp parallel for
for(int k = GRHydro_stencil-1; k < cctk_lsh[2]-GRHydro_stencil+1; k++)
for(int j = GRHydro_stencil-1; j < cctk_lsh[1]-GRHydro_stencil+1; j++)
#pragma ivdep // force compiler to vectorize the goddamn loop
for(int i = GRHydro_stencil-1; i < cctk_lsh[0]-GRHydro_stencil+1; i++) {
const int idx = CCTK_GFINDEX3D(cctkGH,i,j,k);
const int idxl = CCTK_GFINDEX3D(cctkGH,i-*xoffset,j-*yoffset,k-*zoffset);
const int idxr = CCTK_GFINDEX3D(cctkGH,i+*xoffset,j+*yoffset,k+*zoffset);
const double g11l = 0.5 * (g11[idx] + g11[idxl]);
const double g12l = 0.5 * (g12[idx] + g12[idxl]);
const double g13l = 0.5 * (g13[idx] + g13[idxl]);
const double g22l = 0.5 * (g22[idx] + g22[idxl]);
const double g23l = 0.5 * (g23[idx] + g23[idxl]);
const double g33l = 0.5 * (g33[idx] + g33[idxl]);
const double g11r = 0.5 * (g11[idx] + g11[idxr]);
const double g12r = 0.5 * (g12[idx] + g12[idxr]);
const double g13r = 0.5 * (g13[idx] + g13[idxr]);
const double g22r = 0.5 * (g22[idx] + g22[idxr]);
const double g23r = 0.5 * (g23[idx] + g23[idxr]);
const double g33r = 0.5 * (g33[idx] + g33[idxr]);
const double savg_detl =
sqrt(SpatialDeterminantC(g11l,g12l,g13l,g22l,g23l,g33l));
const double savg_detr =
sqrt(SpatialDeterminantC(g11r,g12r,g13r,g22r,g23r,g33r));
// minus call to p2c
prim2conC(&w_lorentzminus[idx], &densminus[idx], &sxminus[idx],
&syminus[idx], &szminus[idx], &tauminus[idx],
g11l,g12l,g13l,g22l,g23l,g33l,
savg_detl,rhominus[idx], velxminus[idx], velyminus[idx],
velzminus[idx], epsminus[idx], pressminus[idx]);
// plus call to p2c
prim2conC(&w_lorentzplus[idx], &densplus[idx], &sxplus[idx],
&syplus[idx], &szplus[idx], &tauplus[idx],
g11r,g12r,g13r,g22r,g23r,g33r,
savg_detr,rhoplus[idx], velxplus[idx], velyplus[idx],
velzplus[idx], epsplus[idx], pressplus[idx]);
}
} // end function Conservative2PrimitiveC
static inline void GRHydro_SpeedOfSound(CCTK_ARGUMENTS)
{
DECLARE_CCTK_PARAMETERS;
DECLARE_CCTK_ARGUMENTS;
// if padding is used the the "extra" vector elements could contain junk
// which we do not know how to handle
assert(cctk_lsh[0] == cctk_ash[0]);
assert(cctk_lsh[1] == cctk_ash[1]);
assert(cctk_lsh[2] == cctk_ash[2]);
// EOS calls (now GF-wide)
if(!*evolve_temper) {
// n needs to be computed using ash since ash is used when computing the
// linear index in CCTK_GFINDEX3D
int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2];
int *keyerr = malloc(sizeof(*keyerr)*n);
int anyerr = 0;
int keytemp = 0;
// don't need special error handling for analytic EOS
#pragma omp parallel for
for(int k=0;k<cctk_lsh[2];k++)
for(int j=0;j<cctk_lsh[1];j++) {
int i = CCTK_GFINDEX3D(cctkGH,0,j,k);
EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_lsh[0],
&(rhominus[i]),&(epsminus[i]),NULL,NULL,&(pressminus[i]),&(cs2minus[i]),
&(keyerr[i]),&anyerr);
EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_lsh[0],
&(rhoplus[i]),&(epsplus[i]),NULL,NULL,&(pressplus[i]),&(cs2plus[i]),
&(keyerr[i]),&anyerr);
}
free(keyerr);
} else {
if(reconstruct_temper) {
int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2];
int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1 + (transport_constraints && !*xoffset);
int *keyerr = malloc(sizeof(*keyerr)*n);
int keytemp = 1;
// ensure Y_e and temperature within bounds
#pragma omp parallel for
for(int i=0;i<n;i++) { // walks over slightly too many elements but cannot fail
Y_e_minus[i] = MAX(MIN(Y_e_minus[i],GRHydro_Y_e_max),GRHydro_Y_e_min);
Y_e_plus[i] = MAX(MIN(Y_e_plus[i],GRHydro_Y_e_max),GRHydro_Y_e_min);
tempminus[i] = MIN(MAX(tempminus[i],GRHydro_hot_atmo_temp),GRHydro_max_temp);
tempplus[i] = MIN(MAX(tempplus[i],GRHydro_hot_atmo_temp),GRHydro_max_temp);
keyerr[i] = 0;
}
// call the EOS with slices
#pragma omp parallel for
for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++)
for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) {
int anyerr = 0;
int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k);
EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx,
&rhominus[i],&epsminus[i],&tempminus[i],&Y_e_minus[i],
&pressminus[i],&cs2minus[i],
&keyerr[i],&anyerr);
if(anyerr) {
for(int ii=GRHydro_stencil-1;ii<cctk_lsh[0]-GRHydro_stencil+1;ii++) {
int idx = CCTK_GFINDEX3D(cctkGH,ii,j,k);
if(keyerr[idx]!=0) {
#pragma omp critical
{
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"rl: %d i,j,k,x,y,z: %d,%d,%d %15.6E %15.6E %15.6E, keyerr: %d",
*GRHydro_reflevel, ii, j, k, x[idx], y[idx], z[idx], keyerr[idx]);
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d",
*GRHydro_reflevel, rhominus[idx], tempminus[idx], Y_e_minus[idx], keyerr[idx]);
}
}
}
CCTK_ERROR("Aborting!");
} // if (anyerr)
} // loop
#pragma omp parallel for
for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++)
for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) {
int anyerr = 0;
int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k);
EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx,
&rhoplus[i],&epsplus[i],&tempplus[i],&Y_e_plus[i],
&pressplus[i],&cs2plus[i],
&keyerr[i],&anyerr);
if(anyerr) {
for(int ii=GRHydro_stencil-1;ii<cctk_lsh[0]-GRHydro_stencil+1;ii++) {
int idx = CCTK_GFINDEX3D(cctkGH,ii,j,k);
if(keyerr[idx]!=0) {
#pragma omp critical
{
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"rl: %d i,j,k,x,y,z: %d,%d,%d %15.6E %15.6E %15.6E, keyerr: %d",
*GRHydro_reflevel, ii, j, k, x[idx], y[idx], z[idx], keyerr[idx]);
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d",
*GRHydro_reflevel, rhoplus[idx], tempplus[idx], Y_e_plus[idx], keyerr[idx]);
}
}
}
CCTK_ERROR("Aborting!");
} // if (anyerr)
} // loop
free(keyerr);
} else {
// ******************** EPS RECONSTRUCTION BRANCH ******************
int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2];
int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1 + (transport_constraints && !*xoffset);
int *keyerr = malloc(sizeof(*keyerr)*n);
int keytemp = 0;
// ensure Y_e and temperature within bounds
#pragma omp parallel for
for(int i=0;i<n;i++) {
Y_e_minus[i] = MAX(MIN(Y_e_minus[i],GRHydro_Y_e_max),GRHydro_Y_e_min);
Y_e_plus[i] = MAX(MIN(Y_e_plus[i],GRHydro_Y_e_max),GRHydro_Y_e_min);
tempminus[i] = MIN(MAX(tempminus[i],GRHydro_hot_atmo_temp),GRHydro_max_temp);
tempplus[i] = MIN(MAX(tempplus[i],GRHydro_hot_atmo_temp),GRHydro_max_temp);
temperature[i] = MIN(MAX(temperature[i],GRHydro_hot_atmo_temp),GRHydro_max_temp);
keyerr[i] = 0;
}
// call the EOS with slices for minus states
#pragma omp parallel for
for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++)
for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) {
int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k);
int anyerr = 0;
EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx,
&rhominus[i],&epsminus[i],&tempminus[i],&Y_e_minus[i],
&pressminus[i],&cs2minus[i],
&keyerr[i],&anyerr);
if(anyerr) {
for(int ii=GRHydro_stencil-1;ii<cctk_lsh[0]-GRHydro_stencil+1;ii++) {
int idx = CCTK_GFINDEX3D(cctkGH,ii,j,k);
if(keyerr[idx]!=0) {
#pragma omp critical
{
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"rl: %d i,j,k,x,y,z: %d,%d,%d %15.6E %15.6E %15.6E, keyerr: %d",
*GRHydro_reflevel, ii, j, k, x[idx], y[idx], z[idx], keyerr[idx]);
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d",
*GRHydro_reflevel, rhominus[idx], tempminus[idx], Y_e_minus[idx], keyerr[idx]);
if(keyerr[idx] == 668) {
// This means the temperature came back negative.
// We'll try using piecewise constant for the temperature
tempminus[idx] = temperature[idx];
const int ln=1;
int lkeyerr[1];
int lanyerr = 0;
int lkeytemp = 1;
EOS_Omni_press_cs2(*GRHydro_eos_handle,lkeytemp,GRHydro_eos_rf_prec,ln,
&rhominus[i],&epsminus[i],&tempminus[i],
&Y_e_minus[i],&pressminus[i],&cs2minus[i],lkeyerr,&lanyerr);
if(lanyerr !=0) {
CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
"Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E",
lkeyerr[0],rhominus[i],tempminus[i],Y_e_minus[i]);
}
} else {
CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
"Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E",
keyerr[i],rhominus[i],tempminus[i],Y_e_minus[i]);
}
} // omp critical pragma
} // if keyerr
} // loop ii
} // if (anyerr)
} // big loop
// call the EOS with slices for plus states
#pragma omp parallel for
for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++)
for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) {
int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k);
int anyerr = 0;
EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx,
&rhoplus[i],&epsplus[i],&tempplus[i],&Y_e_plus[i],
&pressplus[i],&cs2plus[i],
&keyerr[i],&anyerr);
if(anyerr) {
for(int ii=GRHydro_stencil-1;ii<cctk_lsh[0]-GRHydro_stencil+1;ii++) {
int idx = CCTK_GFINDEX3D(cctkGH,ii,j,k);
if(keyerr[idx]!=0) {
#pragma omp critical
{
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"rl: %d i,j,k,x,y,z: %d,%d,%d %15.6E %15.6E %15.6E, keyerr: %d",
*GRHydro_reflevel, ii, j, k, x[idx], y[idx], z[idx], keyerr[idx]);
CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
"rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d",
*GRHydro_reflevel, rhoplus[idx], tempplus[idx], Y_e_plus[idx], keyerr[idx]);
if(keyerr[idx] == 668) {
// This means the temperature came back negative.
// We'll try using piecewise constant for the temperature
tempplus[idx] = temperature[idx];
const int ln=1;
int lkeyerr[1];
int lanyerr = 0;
int lkeytemp = 1;
EOS_Omni_press_cs2(*GRHydro_eos_handle,lkeytemp,GRHydro_eos_rf_prec,ln,
&rhoplus[i],&epsplus[i],&tempplus[i],
&Y_e_plus[i],&pressplus[i],&cs2plus[i],lkeyerr,&lanyerr);
if(lanyerr !=0) {
CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
"Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E",
lkeyerr[0],rhoplus[i],tempplus[i],Y_e_plus[i]);
}
} else {
CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
"Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E",
keyerr[i],rhoplus[i],tempplus[i],Y_e_plus[i]);
}
} // omp critical pragma
} // if keyerr
} // loop ii
} // if (anyerr)
} // big loop over plus states
free(keyerr);
} // end branch for no temp reconsturction
} // end of evolve temper branch
}
static inline double SpatialDeterminantC(double gxx, double gxy,
double gxz, double gyy,
double gyz, double gzz)
{
return -gxz*gxz*gyy + 2.0*gxy*gxz*gyz - gxx*gyz*gyz
- gxy*gxy*gzz + gxx*gyy*gzz;
}
static inline __attribute__((always_inline)) void prim2conC(double *w, double *dens, double *sx,
double *sy, double *sz, double *tau,
const double gxx, const double gxy,
const double gxz, const double gyy,
const double gyz,
const double gzz, const double sdet,
const double rho, const double vx,
const double vy, const double vz,
const double eps, const double press)
{
// local helpers
const double wtemp = 1.0 / sqrt(1.0 - (gxx*vx*vx + gyy*vy*vy +
gzz*vz*vz + 2.0*gxy*vx*vy
+ 2.0*gxz*vx*vz + 2.0*gyz*vy*vz));
const double vlowx = gxx*vx + gxy*vy + gxz*vz;
const double vlowy = gxy*vx + gyy*vy + gyz*vz;
const double vlowz = gxz*vx + gyz*vy + gzz*vz;
const double hrhow2 = (rho*(1.0+eps)+press)*(wtemp)*(wtemp);
const double denstemp = sdet*rho*(wtemp);
*w = wtemp;
*dens = denstemp;
*sx = sdet*hrhow2 * vlowx;
*sy = sdet*hrhow2 * vlowy;
*sz = sdet*hrhow2 * vlowz;
*tau = sdet*( hrhow2 - press) - denstemp;
}
|
