-
Notifications
You must be signed in to change notification settings - Fork 38
/
GiRaFFE_NRPy_Main_Driver_staggered_new_way.py
704 lines (650 loc) · 42.5 KB
/
GiRaFFE_NRPy_Main_Driver_staggered_new_way.py
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
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
# Step 0: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os,sys
nrpy_dir_path = os.path.join("..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
from outputC import outCfunction, lhrh, add_to_Cfunction_dict, outC_function_outdir_dict, outC_function_dict, outC_function_prototype_dict, outC_function_master_list, outC_function_element # NRPy+: Core C code output module
import finite_difference as fin # NRPy+: Finite difference C code generation module
import NRPy_param_funcs as par # NRPy+: Parameter interface
import grid as gri # NRPy+: Functions having to do with numerical grids
import loop as lp # NRPy+: Generate C code loops
import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm # NRPy+: Reference metric support
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
thismodule = __name__
par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER",2)
out_dir = os.path.join("GiRaFFE_standalone_Ccodes")
cmd.mkdir(out_dir)
CoordSystem = "Cartesian"
par.set_parval_from_str("reference_metric::CoordSystem",CoordSystem)
rfm.reference_metric() # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc.
outCparams = "outCverbose=False,CSE_sorting=none"
xi_damping = par.Cparameters("REAL",thismodule,"xi_damping",0.1)
# Default Kreiss-Oliger dissipation strength
default_KO_strength = 0.1
diss_strength = par.Cparameters("REAL", thismodule, "diss_strength", default_KO_strength)
import GRHD.equations as GRHD # NRPy+: Generate general relativistic hydrodynamics equations
gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01",DIM=3)
betaU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","betaU",DIM=3)
alpha = gri.register_gridfunctions("AUXEVOL","alpha")
ixp.register_gridfunctions_for_single_rank1("EVOL","AD")
BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU")
ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BstaggerU")
ValenciavU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","ValenciavU")
gri.register_gridfunctions("EVOL","psi6Phi")
StildeD = ixp.register_gridfunctions_for_single_rank1("EVOL","StildeD")
gri.register_gridfunctions("AUXEVOL","psi6_temp")
gri.register_gridfunctions("AUXEVOL","psi6center")
gri.register_gridfunctions("AUXEVOL","cmax_x")
gri.register_gridfunctions("AUXEVOL","cmin_x")
gri.register_gridfunctions("AUXEVOL","cmax_y")
gri.register_gridfunctions("AUXEVOL","cmin_y")
gri.register_gridfunctions("AUXEVOL","cmax_z")
gri.register_gridfunctions("AUXEVOL","cmin_z")
# We will pass values of the gridfunction on the cell faces into the function. This requires us
# to declare them as C parameters in NRPy+. We will denote this with the _face infix/suffix.
alpha_face = gri.register_gridfunctions("AUXEVOL","alpha_face")
gamma_faceDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gamma_faceDD","sym01")
beta_faceU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","beta_faceU")
# We'll need some more gridfunctions, now, to represent the reconstructions of BU and ValenciavU
# on the right and left faces
Valenciav_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_rU",DIM=3)
B_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_rU",DIM=3)
Valenciav_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_lU",DIM=3)
B_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_lU",DIM=3)
ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Stilde_flux_HLLED")
ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_rrU",DIM=3)
ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_rlU",DIM=3)
ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_lrU",DIM=3)
ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_llU",DIM=3)
ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Bstagger_rU",DIM=3)
ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Bstagger_lU",DIM=3)
# Declare this symbol
sqrt4pi = par.Cparameters("REAL",thismodule,"sqrt4pi","sqrt(4.0*M_PI)")
import GiRaFFE_NRPy.GiRaFFE_NRPy_C2P_P2C as C2P_P2C
def add_to_Cfunction_dict__cons_to_prims(StildeD,BU,gammaDD,betaU,alpha, includes=None):
C2P_P2C.GiRaFFE_NRPy_C2P(StildeD,BU,gammaDD,betaU,alpha)
values_to_print = [
lhrh(lhs=gri.gfaccess("in_gfs","StildeD0"),rhs=C2P_P2C.outStildeD[0]),
lhrh(lhs=gri.gfaccess("in_gfs","StildeD1"),rhs=C2P_P2C.outStildeD[1]),
lhrh(lhs=gri.gfaccess("in_gfs","StildeD2"),rhs=C2P_P2C.outStildeD[2]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU0"),rhs=C2P_P2C.ValenciavU[0]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU1"),rhs=C2P_P2C.ValenciavU[1]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU2"),rhs=C2P_P2C.ValenciavU[2])
]
desc = "Apply fixes to \tilde{S}_i and recompute the velocity to match with current sheet prescription."
name = "GiRaFFE_NRPy_cons_to_prims"
params ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs,REAL *in_gfs"
body = fin.FD_outputC("returnstring",values_to_print,params=outCparams)
loopopts ="AllPoints,Read_xxs"
add_to_Cfunction_dict(
includes=includes,
desc=desc,
name=name, params=params,
body=body, loopopts=loopopts)
# TINYDOUBLE = par.Cparameters("REAL",thismodule,"TINYDOUBLE",1e-100)
def add_to_Cfunction_dict__prims_to_cons(gammaDD,betaU,alpha, ValenciavU,BU, sqrt4pi, includes=None):
C2P_P2C.GiRaFFE_NRPy_P2C(gammaDD,betaU,alpha, ValenciavU,BU, sqrt4pi)
values_to_print = [
lhrh(lhs=gri.gfaccess("in_gfs","StildeD0"),rhs=C2P_P2C.StildeD[0]),
lhrh(lhs=gri.gfaccess("in_gfs","StildeD1"),rhs=C2P_P2C.StildeD[1]),
lhrh(lhs=gri.gfaccess("in_gfs","StildeD2"),rhs=C2P_P2C.StildeD[2]),
]
desc = "Recompute StildeD after current sheet fix to Valencia 3-velocity to ensure consistency between conservative & primitive variables."
name = "GiRaFFE_NRPy_prims_to_cons"
params ="const paramstruct *params,REAL *auxevol_gfs,REAL *in_gfs"
body = fin.FD_outputC("returnstring",values_to_print,params=outCparams)
loopopts ="AllPoints"
add_to_Cfunction_dict(
includes=includes,
desc=desc,
name=name, params=params,
body=body, loopopts=loopopts)
main_evolution_prototype = "void GiRaFFE_NRPy_RHSs(const paramstruct *restrict params,REAL *restrict auxevol_gfs,REAL *restrict in_gfs,REAL *restrict rhs_gfs);"
main_evolution_func = """#include "NRPy_basic_defines.h"
#include "GiRaFFE_basic_defines.h"
#include "NRPy_function_prototypes.h"
#define WORKAROUND_ENABLED
#define MAXNUMINTERP 16
void workaround_Valencia_to_Drift_velocity(const paramstruct *params, REAL *vU0, const REAL *alpha, const REAL *betaU0, const REAL flux_dirn) {
#include "set_Cparameters.h"
// Converts Valencia 3-velocities to Drift 3-velocities for testing. The variable argument
// vu0 is any Valencia 3-velocity component or reconstruction thereof.
#pragma omp parallel for
for (int i2 = 2*(flux_dirn==3);i2 < Nxx_plus_2NGHOSTS2-1*(flux_dirn==3);i2++) for (int i1 = 2*(flux_dirn==2);i1 < Nxx_plus_2NGHOSTS1-1*(flux_dirn==2);i1++) for (int i0 = 2*(flux_dirn==1);i0 < Nxx_plus_2NGHOSTS0-1*(flux_dirn==1);i0++) {
int ii = IDX3S(i0,i1,i2);
// Here, we modify the velocity in place.
vU0[ii] = alpha[ii]*vU0[ii]-betaU0[ii];
}
}
void workaround_Drift_to_Valencia_velocity(const paramstruct *params, REAL *vU0, const REAL *alpha, const REAL *betaU0, const REAL flux_dirn) {
#include "set_Cparameters.h"
// Converts Drift 3-velocities to Valencia 3-velocities for testing. The variable argument
// vu0 is any drift (i.e. IllinoisGRMHD's definition) 3-velocity component or reconstruction thereof.
#pragma omp parallel for
for (int i2 = 2*(flux_dirn==3);i2 < Nxx_plus_2NGHOSTS2-1*(flux_dirn==3);i2++) for (int i1 = 2*(flux_dirn==2);i1 < Nxx_plus_2NGHOSTS1-1*(flux_dirn==2);i1++) for (int i0 = 2*(flux_dirn==1);i0 < Nxx_plus_2NGHOSTS0-1*(flux_dirn==1);i0++) {
int ii = IDX3S(i0,i1,i2);
// Here, we modify the velocity in place.
vU0[ii] = (vU0[ii]+betaU0[ii])/alpha[ii];
}
}
void GiRaFFE_NRPy_RHSs(const paramstruct *restrict params,REAL *restrict auxevol_gfs,REAL *restrict in_gfs,REAL *restrict rhs_gfs) {
#include "set_Cparameters.h"
// First thing's first: initialize the RHSs to zero!
#pragma omp parallel for
for(int ii=0;ii<Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2*NUM_EVOL_GFS;ii++) {
rhs_gfs[ii] = 0.0;
}
// Now, we set up a bunch of structs of pointers to properly guide the PPM algorithm.
// They also count the number of ghostzones available.
gf_and_gz_struct in_prims[NUM_RECONSTRUCT_GFS], out_prims_r[NUM_RECONSTRUCT_GFS], out_prims_l[NUM_RECONSTRUCT_GFS];
int which_prims_to_reconstruct[NUM_RECONSTRUCT_GFS],num_prims_to_reconstruct;
const int Nxxp2NG012 = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;
REAL *temporary = auxevol_gfs + Nxxp2NG012*PSI6_TEMPGF; // Using dedicated temporary variables for the staggered grid
REAL *psi6center = auxevol_gfs + Nxxp2NG012*PSI6CENTERGF; // Because the prescription requires more acrobatics.
// This sets pointers to the portion of auxevol_gfs containing the relevant gridfunction.
int ww=0;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAVU0GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU0GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU0GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAVU1GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU1GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU1GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAVU2GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU2GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU2GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*BU0GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU0GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU0GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*BU1GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU1GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU1GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*BU2GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU2GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU2GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*BSTAGGERU0GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*BSTAGGER_RU0GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*BSTAGGER_LU0GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*BSTAGGERU1GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*BSTAGGER_RU1GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*BSTAGGER_LU1GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*BSTAGGERU2GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*BSTAGGER_RU2GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*BSTAGGER_LU2GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU0GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RRU0GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RLU0GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU1GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RRU1GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RLU1GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU2GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RRU2GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RLU2GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU0GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LRU0GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LLU0GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU1GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LRU1GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LLU1GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU2GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LRU2GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LLU2GF;
ww++;
// Prims are defined AT ALL GRIDPOINTS, so we set the # of ghostzones to zero:
for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { in_prims[i].gz_lo[j]=0; in_prims[i].gz_hi[j]=0; }
// Left/right variables are not yet defined, yet we set the # of gz's to zero by default:
for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { out_prims_r[i].gz_lo[j]=0; out_prims_r[i].gz_hi[j]=0; }
for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { out_prims_l[i].gz_lo[j]=0; out_prims_l[i].gz_hi[j]=0; }
int flux_dirn;
flux_dirn=0;
interpolate_metric_gfs_to_cell_faces(params,auxevol_gfs,flux_dirn+1);
// ftilde = 0 in GRFFE, since P=rho=0.
/* There are two stories going on here:
* 1) Computation of \partial_x on RHS of \partial_t {mhd_st_{x,y,z}},
* via PPM reconstruction onto (i-1/2,j,k), so that
* \partial_y F = [ F(i+1/2,j,k) - F(i-1/2,j,k) ] / dx
* 2) Computation of \partial_t A_i, where A_i are *staggered* gridfunctions,
* where A_x is defined at (i,j+1/2,k+1/2), A_y at (i+1/2,j,k+1/2), etc.
* Ai_rhs = \partial_t A_i = \epsilon_{ijk} \psi^{6} (v^j B^k - v^j B^k),
* where \epsilon_{ijk} is the flat-space antisymmetric operator.
* 2A) Az_rhs is defined at (i+1/2,j+1/2,k), and it depends on {Bx,By,vx,vy},
* so the trick is to reconstruct {Bx,By,vx,vy} cleverly to get to these
* staggered points. For example:
* 2Aa) vx and vy are at (i,j,k), and we reconstruct them to (i-1/2,j,k) below. After
* this, we'll reconstruct again in the y-dir'n to get {vx,vy} at (i-1/2,j-1/2,k)
* 2Ab) By_stagger is at (i,j+1/2,k), and we reconstruct below to (i-1/2,j+1/2,k). */
ww=0;
which_prims_to_reconstruct[ww]=VX; ww++;
which_prims_to_reconstruct[ww]=VY; ww++;
which_prims_to_reconstruct[ww]=VZ; ww++;
//which_prims_to_reconstruct[ww]=BX_CENTER; ww++;
which_prims_to_reconstruct[ww]=BY_CENTER; ww++;
which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;
which_prims_to_reconstruct[ww]=BY_STAGGER;ww++;
num_prims_to_reconstruct=ww;
// This function is housed in the file: "reconstruct_set_of_prims_PPM_GRFFE.C"
reconstruct_set_of_prims_PPM_GRFFE_NRPy(params, auxevol_gfs, flux_dirn+1, num_prims_to_reconstruct,
which_prims_to_reconstruct, in_prims, out_prims_r, out_prims_l, temporary);
//Right and left face values of BI_CENTER are used in GRFFE__S_i__flux computation (first to compute b^a).
// Instead of reconstructing, we simply set B^x face values to be consistent with BX_STAGGER.
#pragma omp parallel for
for(int k=0;k<Nxx_plus_2NGHOSTS2;k++) for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) {
const int index=IDX3S(i,j,k), indexim1=IDX3S(i-1+(i==0),j,k); /* indexim1=0 when i=0 */
out_prims_r[BX_CENTER].gf[index]=out_prims_l[BX_CENTER].gf[index]=in_prims[BX_STAGGER].gf[indexim1]; }
// Then add fluxes to RHS for hydro variables {vx,vy,vz}:
// This function is housed in the file: "add_fluxes_and_source_terms_to_hydro_rhss.C"
calculate_StildeD0_source_term(params,auxevol_gfs,rhs_gfs);
calculate_Stilde_flux_D0(params,auxevol_gfs,rhs_gfs);
calculate_Stilde_rhsD(flux_dirn+1,params,auxevol_gfs,rhs_gfs);
// Note that we have already reconstructed vx and vy along the x-direction,
// at (i-1/2,j,k). That result is stored in v{x,y}{r,l}. Bx_stagger data
// are defined at (i+1/2,j,k).
// Next goal: reconstruct Bx, vx and vy at (i+1/2,j+1/2,k).
flux_dirn=1;
// ftilde = 0 in GRFFE, since P=rho=0.
// in_prims[{VXR,VXL,VYR,VYL}].gz_{lo,hi} ghostzones are set to all zeros, which
// is incorrect. We fix this below.
// [Note that this is a cheap operation, copying only 8 integers and a pointer.]
in_prims[VXR]=out_prims_r[VX];
in_prims[VXL]=out_prims_l[VX];
in_prims[VYR]=out_prims_r[VY];
in_prims[VYL]=out_prims_l[VY];
/* There are two stories going on here:
* 1) Computation of \partial_y on RHS of \partial_t {mhd_st_{x,y,z}},
* via PPM reconstruction onto (i,j-1/2,k), so that
* \partial_y F = [ F(i,j+1/2,k) - F(i,j-1/2,k) ] / dy
* 2) Computation of \partial_t A_i, where A_i are *staggered* gridfunctions,
* where A_x is defined at (i,j+1/2,k+1/2), A_y at (i+1/2,j,k+1/2), etc.
* Ai_rhs = \partial_t A_i = \epsilon_{ijk} \psi^{6} (v^j B^k - v^j B^k),
* where \epsilon_{ijk} is the flat-space antisymmetric operator.
* 2A) Az_rhs is defined at (i+1/2,j+1/2,k), and it depends on {Bx,By,vx,vy},
* so the trick is to reconstruct {Bx,By,vx,vy} cleverly to get to these
* staggered points. For example:
* 2Aa) VXR = [right-face of vx reconstructed along x-direction above] is at (i-1/2,j,k),
* and we reconstruct it to (i-1/2,j-1/2,k) below. Similarly for {VXL,VYR,VYL}
* 2Ab) Bx_stagger is at (i+1/2,j,k), and we reconstruct to (i+1/2,j-1/2,k) below
* 2Ac) By_stagger is at (i-1/2,j+1/2,k) already for Az_rhs, from the previous step.
* 2B) Ax_rhs is defined at (i,j+1/2,k+1/2), and it depends on {By,Bz,vy,vz}.
* Again the trick is to reconstruct these onto these staggered points.
* 2Ba) Bz_stagger is at (i,j,k+1/2), and we reconstruct to (i,j-1/2,k+1/2) below */
ww=0;
// NOTE! The order of variable reconstruction is important here,
// as we don't want to overwrite {vxr,vxl,vyr,vyl}!
which_prims_to_reconstruct[ww]=VXR; ww++;
which_prims_to_reconstruct[ww]=VYR; ww++;
which_prims_to_reconstruct[ww]=VXL; ww++;
which_prims_to_reconstruct[ww]=VYL; ww++;
num_prims_to_reconstruct=ww;
#ifdef WORKAROUND_ENABLED
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU0GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU0GF,flux_dirn+1);
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU1GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU1GF,flux_dirn+1);
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU0GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU0GF,flux_dirn+1);
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU1GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU1GF,flux_dirn+1);
#endif /*WORKAROUND_ENABLED*/
// This function is housed in the file: "reconstruct_set_of_prims_PPM_GRFFE.C"
reconstruct_set_of_prims_PPM_GRFFE_NRPy(params, auxevol_gfs, flux_dirn+1, num_prims_to_reconstruct,
which_prims_to_reconstruct, in_prims, out_prims_r, out_prims_l, temporary);
#ifdef WORKAROUND_ENABLED
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU0GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU0GF,flux_dirn+1);
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU1GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU1GF,flux_dirn+1);
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU0GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU0GF,flux_dirn+1);
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU1GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU1GF,flux_dirn+1);
#endif /*WORKAROUND_ENABLED*/
interpolate_metric_gfs_to_cell_faces(params,auxevol_gfs,flux_dirn+1);
ww=0;
// Reconstruct other primitives last!
which_prims_to_reconstruct[ww]=VX; ww++;
which_prims_to_reconstruct[ww]=VY; ww++;
which_prims_to_reconstruct[ww]=VZ; ww++;
which_prims_to_reconstruct[ww]=BX_CENTER; ww++;
//which_prims_to_reconstruct[ww]=BY_CENTER; ww++;
which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;
which_prims_to_reconstruct[ww]=BX_STAGGER;ww++;
which_prims_to_reconstruct[ww]=BZ_STAGGER;ww++;
num_prims_to_reconstruct=ww;
// This function is housed in the file: "reconstruct_set_of_prims_PPM_GRFFE.C"
reconstruct_set_of_prims_PPM_GRFFE_NRPy(params, auxevol_gfs, flux_dirn+1, num_prims_to_reconstruct,
which_prims_to_reconstruct, in_prims, out_prims_r, out_prims_l, temporary);
//Right and left face values of BI_CENTER are used in GRFFE__S_i__flux computation (first to compute b^a).
// Instead of reconstructing, we simply set B^y face values to be consistent with BY_STAGGER.
#pragma omp parallel for
for(int k=0;k<Nxx_plus_2NGHOSTS2;k++) for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) {
const int index=IDX3S(i,j,k), indexjm1=IDX3S(i,j-1+(j==0),k); /* indexjm1=0 when j=0 */
out_prims_r[BY_CENTER].gf[index]=out_prims_l[BY_CENTER].gf[index]=in_prims[BY_STAGGER].gf[indexjm1]; }
// Then add fluxes to RHS for hydro variables {vx,vy,vz}:
// This function is housed in the file: "add_fluxes_and_source_terms_to_hydro_rhss.C"
calculate_StildeD1_source_term(params,auxevol_gfs,rhs_gfs);
calculate_Stilde_flux_D1(params,auxevol_gfs,rhs_gfs);
calculate_Stilde_rhsD(flux_dirn+1,params,auxevol_gfs,rhs_gfs);
/*****************************************
* COMPUTING RHS OF A_z, BOOKKEEPING NOTE:
* We want to compute
* \partial_t A_z - [gauge terms] = \psi^{6} (v^x B^y - v^y B^x).
* A_z is defined at (i+1/2,j+1/2,k).
* ==========================
* Where defined | Variables
* (i-1/2,j-1/2,k)| {vxrr,vxrl,vxlr,vxll,vyrr,vyrl,vylr,vyll}
* (i+1/2,j-1/2,k)| {Bx_stagger_r,Bx_stagger_l} (see Table 1 in arXiv:1007.2848)
* (i-1/2,j+1/2,k)| {By_stagger_r,By_stagger_l} (see Table 1 in arXiv:1007.2848)
* (i,j,k) | {phi}
* ==========================
******************************************/
// Interpolates to i+1/2
#define IPH(METRICm1,METRICp0,METRICp1,METRICp2) (-0.0625*((METRICm1) + (METRICp2)) + 0.5625*((METRICp0) + (METRICp1)))
// Next compute sqrt(gamma)=psi^6 at (i+1/2,j+1/2,k):
// To do so, we first compute the sqrt of the metric determinant at all points:
#pragma omp parallel for
for(int k=0;k<Nxx_plus_2NGHOSTS2;k++) for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) {
const int index=IDX3S(i,j,k);
const REAL gxx = auxevol_gfs[IDX4ptS(GAMMADD00GF,index)];
const REAL gxy = auxevol_gfs[IDX4ptS(GAMMADD01GF,index)];
const REAL gxz = auxevol_gfs[IDX4ptS(GAMMADD02GF,index)];
const REAL gyy = auxevol_gfs[IDX4ptS(GAMMADD11GF,index)];
const REAL gyz = auxevol_gfs[IDX4ptS(GAMMADD12GF,index)];
const REAL gzz = auxevol_gfs[IDX4ptS(GAMMADD22GF,index)];
psi6center[index] = sqrt( gxx*gyy*gzz
- gxx*gyz*gyz
+2*gxy*gxz*gyz
- gyy*gxz*gxz
- gzz*gxy*gxy );
}
#pragma omp parallel for
for(int k=0;k<Nxx_plus_2NGHOSTS2;k++) for(int j=1;j<Nxx_plus_2NGHOSTS1-2;j++) for(int i=1;i<Nxx_plus_2NGHOSTS0-2;i++) {
temporary[IDX3S(i,j,k)]=
IPH(IPH(psi6center[IDX3S(i-1,j-1,k)],psi6center[IDX3S(i,j-1,k)],psi6center[IDX3S(i+1,j-1,k)],psi6center[IDX3S(i+2,j-1,k)]),
IPH(psi6center[IDX3S(i-1,j ,k)],psi6center[IDX3S(i,j ,k)],psi6center[IDX3S(i+1,j ,k)],psi6center[IDX3S(i+2,j ,k)]),
IPH(psi6center[IDX3S(i-1,j+1,k)],psi6center[IDX3S(i,j+1,k)],psi6center[IDX3S(i+1,j+1,k)],psi6center[IDX3S(i+2,j+1,k)]),
IPH(psi6center[IDX3S(i-1,j+2,k)],psi6center[IDX3S(i,j+2,k)],psi6center[IDX3S(i+1,j+2,k)],psi6center[IDX3S(i+2,j+2,k)]));
}
int A_directionz=3;
A_i_rhs_no_gauge_terms(A_directionz,params,out_prims_r,out_prims_l,temporary,
auxevol_gfs+Nxxp2NG012*CMAX_XGF,
auxevol_gfs+Nxxp2NG012*CMIN_XGF,
auxevol_gfs+Nxxp2NG012*CMAX_YGF,
auxevol_gfs+Nxxp2NG012*CMIN_YGF,
rhs_gfs+Nxxp2NG012*AD2GF);
// in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not correct, so we fix
// this below.
// [Note that this is a cheap operation, copying only 8 integers and a pointer.]
in_prims[VYR]=out_prims_r[VY];
in_prims[VYL]=out_prims_l[VY];
in_prims[VZR]=out_prims_r[VZ];
in_prims[VZL]=out_prims_l[VZ];
flux_dirn=2;
// ftilde = 0 in GRFFE, since P=rho=0.
/* There are two stories going on here:
* 1) Single reconstruction to (i,j,k-1/2) for {vx,vy,vz,Bx,By,Bz} to compute
* z-dir'n advection terms in \partial_t {mhd_st_{x,y,z}} at (i,j,k)
* 2) Multiple reconstructions for *staggered* gridfunctions A_i:
* \partial_t A_i = \epsilon_{ijk} \psi^{6} (v^j B^k - v^j B^k),
* where \epsilon_{ijk} is the flat-space antisymmetric operator.
* 2A) Ax_rhs is defined at (i,j+1/2,k+1/2), depends on v{y,z} and B{y,z}
* 2Aa) v{y,z}{r,l} are at (i,j-1/2,k), so we reconstruct here to (i,j-1/2,k-1/2)
* 2Ab) Bz_stagger{r,l} are at (i,j-1/2,k+1/2) already.
* 2Ac) By_stagger is at (i,j+1/2,k), and below we reconstruct its value at (i,j+1/2,k-1/2)
* 2B) Ay_rhs is defined at (i+1/2,j,k+1/2), depends on v{z,x} and B{z,x}.
* 2Ba) v{x,z} are reconstructed to (i,j,k-1/2). Later we'll reconstruct again to (i-1/2,j,k-1/2).
* 2Bb) Bz_stagger is at (i,j,k+1/2). Later we will reconstruct to (i-1/2,j,k+1/2).
* 2Bc) Bx_stagger is at (i+1/2,j,k), and below we reconstruct its value at (i+1/2,j,k-1/2)
*/
ww=0;
// NOTE! The order of variable reconstruction is important here,
// as we don't want to overwrite {vxr,vxl,vyr,vyl}!
which_prims_to_reconstruct[ww]=VYR; ww++;
which_prims_to_reconstruct[ww]=VZR; ww++;
which_prims_to_reconstruct[ww]=VYL; ww++;
which_prims_to_reconstruct[ww]=VZL; ww++;
num_prims_to_reconstruct=ww;
#ifdef WORKAROUND_ENABLED
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU1GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU1GF,flux_dirn+1);
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU2GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU2GF,flux_dirn+1);
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU1GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU1GF,flux_dirn+1);
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU2GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU2GF,flux_dirn+1);
#endif /*WORKAROUND_ENABLED*/
// This function is housed in the file: "reconstruct_set_of_prims_PPM_GRFFE.C"
reconstruct_set_of_prims_PPM_GRFFE_NRPy(params, auxevol_gfs, flux_dirn+1, num_prims_to_reconstruct,
which_prims_to_reconstruct, in_prims, out_prims_r, out_prims_l, temporary);
#ifdef WORKAROUND_ENABLED
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU1GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU1GF,flux_dirn+1);
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU2GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU2GF,flux_dirn+1);
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU1GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU1GF,flux_dirn+1);
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU2GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU2GF,flux_dirn+1);
#endif /*WORKAROUND_ENABLED*/
interpolate_metric_gfs_to_cell_faces(params,auxevol_gfs,flux_dirn+1);
// Reconstruct other primitives last!
ww=0;
which_prims_to_reconstruct[ww]=VX; ww++;
which_prims_to_reconstruct[ww]=VY; ww++;
which_prims_to_reconstruct[ww]=VZ; ww++;
which_prims_to_reconstruct[ww]=BX_CENTER; ww++;
which_prims_to_reconstruct[ww]=BY_CENTER; ww++;
//which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;
which_prims_to_reconstruct[ww]=BX_STAGGER; ww++;
which_prims_to_reconstruct[ww]=BY_STAGGER; ww++;
num_prims_to_reconstruct=ww;
// This function is housed in the file: "reconstruct_set_of_prims_PPM_GRFFE.C"
reconstruct_set_of_prims_PPM_GRFFE_NRPy(params, auxevol_gfs, flux_dirn+1, num_prims_to_reconstruct,
which_prims_to_reconstruct, in_prims, out_prims_r, out_prims_l, temporary);
//Right and left face values of BI_CENTER are used in GRFFE__S_i__flux computation (first to compute b^a).
// Instead of reconstructing, we simply set B^z face values to be consistent with BZ_STAGGER.
#pragma omp parallel for
for(int k=0;k<Nxx_plus_2NGHOSTS2;k++) for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) {
const int index=IDX3S(i,j,k), indexkm1=IDX3S(i,j,k-1+(k==0)); /* indexkm1=0 when k=0 */
out_prims_r[BZ_CENTER].gf[index]=out_prims_l[BZ_CENTER].gf[index]=in_prims[BZ_STAGGER].gf[indexkm1]; }
// Then add fluxes to RHS for hydro variables {vx,vy,vz}:
// This function is housed in the file: "add_fluxes_and_source_terms_to_hydro_rhss.C"
calculate_StildeD2_source_term(params,auxevol_gfs,rhs_gfs);
calculate_Stilde_flux_D2(params,auxevol_gfs,rhs_gfs);
calculate_Stilde_rhsD(flux_dirn+1,params,auxevol_gfs,rhs_gfs);
// in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not set correcty.
// We fix this below.
// [Note that this is a cheap operation, copying only 8 integers and a pointer.]
in_prims[VXR]=out_prims_r[VX];
in_prims[VZR]=out_prims_r[VZ];
in_prims[VXL]=out_prims_l[VX];
in_prims[VZL]=out_prims_l[VZ];
// FIXME: lines above seem to be inconsistent with lines below.... Possible bug, not major enough to affect evolutions though.
in_prims[VZR].gz_lo[1]=in_prims[VZR].gz_hi[1]=0;
in_prims[VXR].gz_lo[1]=in_prims[VXR].gz_hi[1]=0;
in_prims[VZL].gz_lo[1]=in_prims[VZL].gz_hi[1]=0;
in_prims[VXL].gz_lo[1]=in_prims[VXL].gz_hi[1]=0;
/*****************************************
* COMPUTING RHS OF A_x, BOOKKEEPING NOTE:
* We want to compute
* \partial_t A_x - [gauge terms] = \psi^{6} (v^y B^z - v^z B^y).
* A_x is defined at (i,j+1/2,k+1/2).
* ==========================
* Where defined | Variables
* (i,j-1/2,k-1/2)| {vyrr,vyrl,vylr,vyll,vzrr,vzrl,vzlr,vzll}
* (i,j+1/2,k-1/2)| {By_stagger_r,By_stagger_l} (see Table 1 in arXiv:1007.2848)
* (i,j-1/2,k+1/2)| {Bz_stagger_r,Bz_stagger_l} (see Table 1 in arXiv:1007.2848)
* (i,j,k) | {phi}
* ==========================
******************************************/
// Next compute phi at (i,j+1/2,k+1/2):
#pragma omp parallel for
for(int k=1;k<Nxx_plus_2NGHOSTS2-2;k++) for(int j=1;j<Nxx_plus_2NGHOSTS1-2;j++) for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) {
temporary[IDX3S(i,j,k)]=
IPH(IPH(psi6center[IDX3S(i,j-1,k-1)],psi6center[IDX3S(i,j,k-1)],psi6center[IDX3S(i,j+1,k-1)],psi6center[IDX3S(i,j+2,k-1)]),
IPH(psi6center[IDX3S(i,j-1,k )],psi6center[IDX3S(i,j,k )],psi6center[IDX3S(i,j+1,k )],psi6center[IDX3S(i,j+2,k )]),
IPH(psi6center[IDX3S(i,j-1,k+1)],psi6center[IDX3S(i,j,k+1)],psi6center[IDX3S(i,j+1,k+1)],psi6center[IDX3S(i,j+2,k+1)]),
IPH(psi6center[IDX3S(i,j-1,k+2)],psi6center[IDX3S(i,j,k+2)],psi6center[IDX3S(i,j+1,k+2)],psi6center[IDX3S(i,j+2,k+2)]));
}
int A_directionx=1;
A_i_rhs_no_gauge_terms(A_directionx,params,out_prims_r,out_prims_l,temporary,
auxevol_gfs+Nxxp2NG012*CMAX_YGF,
auxevol_gfs+Nxxp2NG012*CMIN_YGF,
auxevol_gfs+Nxxp2NG012*CMAX_ZGF,
auxevol_gfs+Nxxp2NG012*CMIN_ZGF,
rhs_gfs+Nxxp2NG012*AD0GF);
// We reprise flux_dirn=0 to finish up computations of Ai_rhs's!
flux_dirn=0;
// ftilde = 0 in GRFFE, since P=rho=0.
ww=0;
// NOTE! The order of variable reconstruction is important here,
// as we don't want to overwrite {vxr,vxl,vyr,vyl}!
which_prims_to_reconstruct[ww]=VXR; ww++;
which_prims_to_reconstruct[ww]=VZR; ww++;
which_prims_to_reconstruct[ww]=VXL; ww++;
which_prims_to_reconstruct[ww]=VZL; ww++;
which_prims_to_reconstruct[ww]=BZ_STAGGER;ww++;
num_prims_to_reconstruct=ww;
#ifdef WORKAROUND_ENABLED
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU0GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU0GF,flux_dirn+1);
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU2GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU2GF,flux_dirn+1);
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU0GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU0GF,flux_dirn+1);
workaround_Valencia_to_Drift_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU2GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU2GF,flux_dirn+1);
#endif /*WORKAROUND_ENABLED*/
// This function is housed in the file: "reconstruct_set_of_prims_PPM_GRFFE.C"
reconstruct_set_of_prims_PPM_GRFFE_NRPy(params, auxevol_gfs, flux_dirn+1, num_prims_to_reconstruct,
which_prims_to_reconstruct, in_prims, out_prims_r, out_prims_l, temporary);
#ifdef WORKAROUND_ENABLED
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU0GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU0GF,flux_dirn+1);
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_RU2GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU2GF,flux_dirn+1);
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU0GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU0GF,flux_dirn+1);
workaround_Drift_to_Valencia_velocity(params,auxevol_gfs+Nxxp2NG012*VALENCIAV_LU2GF,auxevol_gfs+Nxxp2NG012*ALPHA_FACEGF,auxevol_gfs+Nxxp2NG012*BETA_FACEU2GF,flux_dirn+1);
#endif /*WORKAROUND_ENABLED*/
/*****************************************
* COMPUTING RHS OF A_y, BOOKKEEPING NOTE:
* We want to compute
* \partial_t A_y - [gauge terms] = \psi^{6} (v^z B^x - v^x B^z).
* A_y is defined at (i+1/2,j,k+1/2).
* ==========================
* Where defined | Variables
* (i-1/2,j,k-1/2)| {vyrr,vyrl,vylr,vyll,vzrr,vzrl,vzlr,vzll}
* (i+1/2,j,k-1/2)| {By_stagger_r,By_stagger_l} (see Table 1 in arXiv:1007.2848)
* (i-1/2,j,k+1/2)| {Bz_stagger_r,Bz_stagger_l} (see Table 1 in arXiv:1007.2848)
* (i,j,k) | {phi}
* ==========================
******************************************/
// Next compute phi at (i+1/2,j,k+1/2):
#pragma omp parallel for
for(int k=1;k<Nxx_plus_2NGHOSTS2-2;k++) for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) for(int i=1;i<Nxx_plus_2NGHOSTS0-2;i++) {
temporary[IDX3S(i,j,k)]=
IPH(IPH(psi6center[IDX3S(i-1,j,k-1)],psi6center[IDX3S(i,j,k-1)],psi6center[IDX3S(i+1,j,k-1)],psi6center[IDX3S(i+2,j,k-1)]),
IPH(psi6center[IDX3S(i-1,j,k )],psi6center[IDX3S(i,j,k )],psi6center[IDX3S(i+1,j,k )],psi6center[IDX3S(i+2,j,k )]),
IPH(psi6center[IDX3S(i-1,j,k+1)],psi6center[IDX3S(i,j,k+1)],psi6center[IDX3S(i+1,j,k+1)],psi6center[IDX3S(i+2,j,k+1)]),
IPH(psi6center[IDX3S(i-1,j,k+2)],psi6center[IDX3S(i,j,k+2)],psi6center[IDX3S(i+1,j,k+2)],psi6center[IDX3S(i+2,j,k+2)]));
}
int A_directiony=2;
A_i_rhs_no_gauge_terms(A_directiony,params,out_prims_r,out_prims_l,temporary,
auxevol_gfs+Nxxp2NG012*CMAX_ZGF,
auxevol_gfs+Nxxp2NG012*CMIN_ZGF,
auxevol_gfs+Nxxp2NG012*CMAX_XGF,
auxevol_gfs+Nxxp2NG012*CMIN_XGF,
rhs_gfs+Nxxp2NG012*AD1GF);
// Next compute psi6phi_rhs, and add gauge terms to A_i_rhs terms!
// Note that in the following function, we don't bother with reconstruction, instead interpolating.
// We need A^i, but only have A_i. So we add gtupij to the list of input variables.
REAL *interp_vars[MAXNUMINTERP];
ww=0;
interp_vars[ww]=auxevol_gfs+Nxxp2NG012*BETAU0GF; ww++;
interp_vars[ww]=auxevol_gfs+Nxxp2NG012*BETAU1GF; ww++;
interp_vars[ww]=auxevol_gfs+Nxxp2NG012*BETAU2GF; ww++;
interp_vars[ww]=auxevol_gfs+Nxxp2NG012*GAMMADD00GF; ww++;
interp_vars[ww]=auxevol_gfs+Nxxp2NG012*GAMMADD01GF; ww++;
interp_vars[ww]=auxevol_gfs+Nxxp2NG012*GAMMADD02GF; ww++;
interp_vars[ww]=auxevol_gfs+Nxxp2NG012*GAMMADD11GF; ww++;
interp_vars[ww]=auxevol_gfs+Nxxp2NG012*GAMMADD12GF; ww++;
interp_vars[ww]=auxevol_gfs+Nxxp2NG012*GAMMADD22GF; ww++;
interp_vars[ww]=temporary;ww++;
interp_vars[ww]=auxevol_gfs+Nxxp2NG012*ALPHAGF; ww++;
interp_vars[ww]=in_gfs+Nxxp2NG012*AD0GF; ww++;
interp_vars[ww]=in_gfs+Nxxp2NG012*AD1GF; ww++;
interp_vars[ww]=in_gfs+Nxxp2NG012*AD2GF; ww++;
const int max_num_interp_variables=ww;
// if(max_num_interp_variables>MAXNUMINTERP) {CCTK_VError(VERR_DEF_PARAMS,"Error: Didn't allocate enough space for interp_vars[]."); }
// We are FINISHED with v{x,y,z}{r,l} and P{r,l} so we use these 8 gridfunctions' worth of space as temp storage.
Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs(params,interp_vars,
in_gfs+Nxxp2NG012*PSI6PHIGF,
auxevol_gfs+Nxxp2NG012*VALENCIAV_RU0GF, // WARNING:
auxevol_gfs+Nxxp2NG012*VALENCIAV_RU1GF, // ALL VARIABLES
auxevol_gfs+Nxxp2NG012*VALENCIAV_RU2GF, // ON THESE LINES
auxevol_gfs+Nxxp2NG012*VALENCIAV_LU0GF, // ARE OVERWRITTEN
auxevol_gfs+Nxxp2NG012*VALENCIAV_LU1GF, // FOR TEMP STORAGE
auxevol_gfs+Nxxp2NG012*VALENCIAV_LU2GF, // .
auxevol_gfs+Nxxp2NG012*VALENCIAV_RRU0GF, // .
auxevol_gfs+Nxxp2NG012*VALENCIAV_RLU0GF, // .
rhs_gfs+Nxxp2NG012*PSI6PHIGF,
rhs_gfs+Nxxp2NG012*AD0GF,
rhs_gfs+Nxxp2NG012*AD1GF,
rhs_gfs+Nxxp2NG012*AD2GF);
/*#pragma omp parallel for
for(int k=0;k<Nxx_plus_2NGHOSTS2;k++) for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) {
REAL x = xx[0][i];
REAL y = xx[1][j];
REAL z = xx[2][k];
if(sqrt(x*x+y*y+z*z)<min_radius_inside_of_which_conserv_to_prims_FFE_and_FFE_evolution_is_DISABLED) {
rhs_gfs[IDX4S(STILDED0GF,i,j,k)] = 0.0;
rhs_gfs[IDX4S(STILDED1GF,i,j,k)] = 0.0;
rhs_gfs[IDX4S(STILDED2GF,i,j,k)] = 0.0;
rhs_gfs[IDX4S(PSI6PHIGF,i,j,k)] = 0.0;
rhs_gfs[IDX4S(AD0GF,i,j,k)] = 0.0;
rhs_gfs[IDX4S(AD1GF,i,j,k)] = 0.0;
rhs_gfs[IDX4S(AD2GF,i,j,k)] = 0.0;
}
}*/
}"""
post_step_prototype = "void GiRaFFE_NRPy_post_step(const paramstruct *restrict params,REAL *xx[3],REAL *restrict auxevol_gfs,REAL *restrict evol_gfs,const int n);"
post_step_func = """#include "NRPy_basic_defines.h"
#include "GiRaFFE_basic_defines.h"
#include "NRPy_function_prototypes.h"
void GiRaFFE_NRPy_post_step(const paramstruct *restrict params,REAL *xx[3],REAL *restrict auxevol_gfs,REAL *restrict evol_gfs,const int n) {
#include "set_Cparameters.h"
// First, apply BCs to AD and psi6Phi. Then calculate BU from AD
apply_bcs_potential(params,evol_gfs);
const int Nxxp2NG012 = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;
GiRaFFE_compute_B_and_Bstagger_from_A(params,
auxevol_gfs+Nxxp2NG012*GAMMADD00GF,
auxevol_gfs+Nxxp2NG012*GAMMADD01GF,
auxevol_gfs+Nxxp2NG012*GAMMADD02GF,
auxevol_gfs+Nxxp2NG012*GAMMADD11GF,
auxevol_gfs+Nxxp2NG012*GAMMADD12GF,
auxevol_gfs+Nxxp2NG012*GAMMADD22GF,
auxevol_gfs + Nxxp2NG012*PSI6_TEMPGF, /* Temporary storage */
evol_gfs+Nxxp2NG012*AD0GF,
evol_gfs+Nxxp2NG012*AD1GF,
evol_gfs+Nxxp2NG012*AD2GF,
auxevol_gfs+Nxxp2NG012*BU0GF,
auxevol_gfs+Nxxp2NG012*BU1GF,
auxevol_gfs+Nxxp2NG012*BU2GF,
auxevol_gfs+Nxxp2NG012*BSTAGGERU0GF,
auxevol_gfs+Nxxp2NG012*BSTAGGERU1GF,
auxevol_gfs+Nxxp2NG012*BSTAGGERU2GF);
//override_BU_with_old_GiRaFFE(params,auxevol_gfs,n);
// Apply fixes to StildeD, then recompute the velocity at the new timestep.
// Apply the current sheet prescription to the velocities
GiRaFFE_NRPy_cons_to_prims(params,xx,auxevol_gfs,evol_gfs);
// Then, recompute StildeD to be consistent with the new velocities
//GiRaFFE_NRPy_prims_to_cons(params,auxevol_gfs,evol_gfs);
// Finally, apply outflow boundary conditions to the velocities.
apply_bcs_velocity(params,auxevol_gfs);
}"""
def add_to_Cfunction_dict__driver_function():
outC_function_master_list.append(outC_function_element("empty", "empty", "empty", "empty", "GiRaFFE_NRPy_RHSs", "empty",
"empty", "empty", "empty", "empty",
"empty", "empty"))
outC_function_outdir_dict["GiRaFFE_NRPy_RHSs"] = "default"
outC_function_dict["GiRaFFE_NRPy_RHSs"] = main_evolution_func
outC_function_prototype_dict["GiRaFFE_NRPy_RHSs"] = main_evolution_prototype
outC_function_master_list.append(outC_function_element("empty", "empty", "empty", "empty", "GiRaFFE_NRPy_post_step", "empty",
"empty", "empty", "empty", "empty",
"empty", "empty"))
outC_function_outdir_dict["GiRaFFE_NRPy_post_step"] = "default"
outC_function_dict["GiRaFFE_NRPy_post_step"] = post_step_func
outC_function_prototype_dict["GiRaFFE_NRPy_post_step"] = post_step_prototype