-
Notifications
You must be signed in to change notification settings - Fork 7
/
advance.c
2706 lines (1975 loc) · 89.7 KB
/
advance.c
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
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#include "decs.h"
// static declarations
static int compute_dt_fromsource(struct of_geom *ptrgeom, struct of_state *state, FTYPE *U, FTYPE *pr, FTYPE *dUevolve, FTYPE *dUgeomevolveUU, FTYPE *dtij, FTYPE *gravitydt);
static int dUtodt(struct of_geom *ptrgeom, struct of_state *state, FTYPE *pr, FTYPE *dUgeom, FTYPE *dUriemann, FTYPE *dUgeomgravity, FTYPE *accdt, FTYPE *gravitydt);
static int check_point_vs_average(int timeorder, int numtimeorders, PFTYPE *lpflag, FTYPE *pb, FTYPE *pf, FTYPE *upoint, FTYPE *uavg, struct of_geom *ptrgeom, struct of_newtonstats *newtonstats);
static int prepare_globaldt(
int truestep,
FTYPE ndt1,FTYPE ndt2,FTYPE ndt3,
FTYPE accdt,int accdti,int accdtj,int accdtk,
FTYPE gravitydt,int gravitydti,int gravitydtj,int gravitydtk,
FTYPE *ndt);
static void flux2dUavg(int i, int j, int k, FTYPE (*F1)[NSTORE2][NSTORE3][NPR],FTYPE (*F2)[NSTORE2][NSTORE3][NPR],FTYPE (*F3)[NSTORE2][NSTORE3][NPR],FTYPE *dUavg1,FTYPE *dUavg2,FTYPE *dUavg3);
static void dUtoU(int i, int j, int k, int loc, FTYPE *dUgeom, FTYPE *dUriemann, FTYPE *CUf, FTYPE *CUnew, FTYPE *Ui, FTYPE *uf, FTYPE *ucum);
static void dUtoU_check(int i, int j, int k, int loc, int pl, FTYPE *dUgeom, FTYPE *dUriemann, FTYPE *CUf, FTYPE *CUnew, FTYPE *Ui, FTYPE *Uf, FTYPE *ucum);
static int asym_compute_1(FTYPE (*prim)[NSTORE2][NSTORE3][NPR]);
static int asym_compute_2(FTYPE (*prim)[NSTORE2][NSTORE3][NPR]);
static FTYPE fractional_diff( FTYPE a, FTYPE b );
// AVG_2_POINT functions:
static void debugeno_compute(FTYPE (*p)[NSTORE2][NSTORE3][NPR],FTYPE (*U)[NSTORE2][NSTORE3][NPR],FTYPE (*debugvars)[NSTORE2][NSTORE3][NUMENODEBUGS]);
static void show_fluxes(int i, int j, int k, int loc, int pl,FTYPE (*F1)[NSTORE2][NSTORE3][NPR],FTYPE (*F2)[NSTORE2][NSTORE3][NPR],FTYPE (*F3)[NSTORE2][NSTORE3][NPR]);
static int advance_standard(int truestep,int stage, FTYPE (*pi)[NSTORE2][NSTORE3][NPR],FTYPE (*pb)[NSTORE2][NSTORE3][NPR], FTYPE (*pf)[NSTORE2][NSTORE3][NPR],
FTYPE (*pstag)[NSTORE2][NSTORE3][NPR],
FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],
FTYPE (*F1)[NSTORE2][NSTORE3][NPR],FTYPE (*F2)[NSTORE2][NSTORE3][NPR],FTYPE (*F3)[NSTORE2][NSTORE3][NPR],
FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],
FTYPE (*ui)[NSTORE2][NSTORE3][NPR], FTYPE (*uf)[NSTORE2][NSTORE3][NPR], FTYPE (*ucum)[NSTORE2][NSTORE3][NPR],
FTYPE *CUf,FTYPE *CUnew,SFTYPE fluxdt, SFTYPE boundtime, SFTYPE fluxtime, int stagenow, int numstages, FTYPE *ndt);
static int advance_finitevolume(int truestep,int stage, FTYPE (*pi)[NSTORE2][NSTORE3][NPR],FTYPE (*pb)[NSTORE2][NSTORE3][NPR], FTYPE (*pf)[NSTORE2][NSTORE3][NPR],
FTYPE (*pstag)[NSTORE2][NSTORE3][NPR],
FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],
FTYPE (*F1)[NSTORE2][NSTORE3][NPR],FTYPE (*F2)[NSTORE2][NSTORE3][NPR],FTYPE (*F3)[NSTORE2][NSTORE3][NPR],
FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],
FTYPE (*ui)[NSTORE2][NSTORE3][NPR],FTYPE (*uf)[NSTORE2][NSTORE3][NPR], FTYPE (*ucum)[NSTORE2][NSTORE3][NPR],
FTYPE *CUf,FTYPE *CUnew, SFTYPE fluxdt, SFTYPE boundtime, SFTYPE fluxtime, int stagenow, int numstages, FTYPE *ndt);
// things to do before any interpolation or advance step
// includes pre-computed things for interpolation and advance that (e.g.) aren't required to perform for each interpolation or advance call or portion of a call
void pre_interpolate_and_advance(FTYPE (*pb)[NSTORE2][NSTORE3][NPR])
{
/////////////////////////////////////
//
// Compute and Store (globally) the get_state() data for the CENT position to avoid computing later and for merged-higher-order method
//
/////////////////////////////////////
#if(STOREFLUXSTATE||STORESHOCKINDICATOR)
// NOTE: This is done before advance since always needed, and only needed once for all dimensions, and don't here instead of inside advance() since needed during fluxcalc() that is called first before any use of get_geometry() that we would use to put this call with
compute_and_store_fluxstatecent(pb);
// now flux_compute() and other flux-position-related things will obtain get_state() data for p_l and p_r from global arrays
#endif
}
// pi: initial values at t=t0 to compute Ui
// pb: values used to compute flux/source
// pf: solution using flux(pb) from pi's Ui -> Uf
// pi, pb, and pf can all be the same since
// 1) pb used first on a stencil, not modified, to compute fluxes
// 2) pf=pi is assigned by value at each zone
// 3) pf is modified using Utoprim at each zone using pb for sources (to correspond to fluxes which used pb)
//
// So in the end only pf is modified at each zone, so the loop changing p at previous (i,j) location doesn't affect the any new location in (i,j)
int advance(int truestep, int stage, FTYPE (*pi)[NSTORE2][NSTORE3][NPR],FTYPE (*pb)[NSTORE2][NSTORE3][NPR], FTYPE (*pf)[NSTORE2][NSTORE3][NPR],
FTYPE (*pstag)[NSTORE2][NSTORE3][NPR],
FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],
FTYPE (*F1)[NSTORE2][NSTORE3][NPR],FTYPE (*F2)[NSTORE2][NSTORE3][NPR],FTYPE (*F3)[NSTORE2][NSTORE3][NPR],
FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],
FTYPE (*ui)[NSTORE2][NSTORE3][NPR],FTYPE (*uf)[NSTORE2][NSTORE3][NPR], FTYPE (*ucum)[NSTORE2][NSTORE3][NPR],
FTYPE *CUf, FTYPE *CUnew, SFTYPE fluxdt, SFTYPE boundtime, SFTYPE fluxtime, int timeorder, int numtimeorders, FTYPE *ndt)
{
////////////////
//
// setup state and interpolation stuff for interpolation and advance calls
//
////////////////
pre_interpolate_and_advance(pb);
////////////////
//
// advance
//
////////////////
if(DOENOFLUX==ENOFINITEVOLUME){
MYFUN(advance_finitevolume(truestep,stage,pi,pb,pf,pstag,pl_ct, pr_ct, F1, F2, F3, vpot,ui,uf,ucum,CUf,CUnew,fluxdt,boundtime,fluxtime,timeorder,numtimeorders,ndt),"advance.c:advance()", "advance_finitevolume()", 1);
}
else if((DOENOFLUX==NOENOFLUX)||(DOENOFLUX==ENOFLUXRECON)||(DOENOFLUX==ENOFLUXSPLIT)){
// new standard preserves conserved quantities even when metric changes
MYFUN(advance_standard(truestep,stage,pi,pb,pf,pstag,pl_ct, pr_ct, F1, F2, F3, vpot,ui,uf,ucum,CUf,CUnew,fluxdt,boundtime,fluxtime,timeorder,numtimeorders,ndt),"advance.c:advance()", "advance_standard()", 1);
}
else{
dualfprintf(fail_file,"No such DOENOFLUX=%d\n",DOENOFLUX);
myexit(1);
}
return(0);
}
// Notes:
// loop range defined with +SHIFT? so consistent with requirement by IF3DSPCTHENMPITRANSFERATPOLE or consistent with setting upper face flux and using it, which is only used for FLUXBSTAG for evolving field on faces. This forces field on upper face to be evolved as required in some cases.
// This is bit excessive for non-face quantities, but fake partial evolution of centered value at "N" is ok as long as don't hit NaN's that slow down code.
// this method guarantees conservation of non-sourced conserved quantities when metric is time-dependent
// this method has updated field staggered method
// Note that when dt==0.0, assume no fluxing, just take ucum -> ui -> {uf,ucum} and invert. Used with metric update.
static int advance_standard(
int truestep,
int stage,
FTYPE (*pi)[NSTORE2][NSTORE3][NPR],
FTYPE (*pb)[NSTORE2][NSTORE3][NPR],
FTYPE (*pf)[NSTORE2][NSTORE3][NPR],
FTYPE (*pstag)[NSTORE2][NSTORE3][NPR],
FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],
FTYPE (*F1)[NSTORE2][NSTORE3][NPR],FTYPE (*F2)[NSTORE2][NSTORE3][NPR],FTYPE (*F3)[NSTORE2][NSTORE3][NPR],
FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],
FTYPE (*ui)[NSTORE2][NSTORE3][NPR],
FTYPE (*uf)[NSTORE2][NSTORE3][NPR],
FTYPE (*ucum)[NSTORE2][NSTORE3][NPR],
FTYPE *CUf,
FTYPE *CUnew,
SFTYPE fluxdt,
SFTYPE boundtime,
SFTYPE fluxtime,
int timeorder, int numtimeorders,
FTYPE *ndt)
{
FTYPE ndt1, ndt2, ndt3;
FTYPE dUtot;
FTYPE idx1,idx2;
SFTYPE dt4diag;
static SFTYPE dt4diag_willbe=0;
int finalstep,initialstep;
FTYPE accdt, accdt_ij;
int accdti,accdtj,accdtk;
FTYPE gravitydt, gravitydt_ij;
int gravitydti,gravitydtj,gravitydtk;
// FTYPE (*dUriemannarray)[NSTORE2][NSTORE3][NPR];
FTYPE (*ucumformetric)[NSTORE2][NSTORE3][NPR];
int enerregion;
int *localenerpos;
int jj;
FTYPE (*utoinvert)[NSTORE2][NSTORE3][NPR];
FTYPE (*tempucum)[NSTORE2][NSTORE3][NPR];
FTYPE (*useducum)[NSTORE2][NSTORE3][NPR];
FTYPE (*myupoint)[NSTORE2][NSTORE3][NPR];
int whichpltoavg[NPR];
int ifnotavgthencopy[NPR];
int is,ie,js,je,ks,ke;
int doingextrashiftforstag;
if(FLUXB==FLUXCTSTAG){
// fill in tempucum with changes that are not controlled well in space, but later fill in ucum from this in controlled way
tempucum=GLOBALPOINT(utemparray);
useducum=tempucum; // unless changed
}
else{
// no special shifting of indices occurs that requires loop shifting
tempucum=ucum;
useducum=ucum;
}
ucumformetric=GLOBALPOINT(ucumformetric);// temporary space for ucum for metric that is same time as "pb", so not updated yet or is ui
/////////////////////////////////////////////
//
// Setup function tasks
//
////////////////////////////////////////////
// for ZLOOP:
// avoid looping over region outside active+ghost grid
// good because somewhat general and avoid bad inversions, etc.
enerregion=TRUEGLOBALWITHBNDENERREGION;
localenerpos=enerposreg[enerregion];
accdt=BIG; // initially no limit to dt due to acceleration
accdti=accdtj=accdtk=-100;
gravitydt=BIG; // initially no limit to dt due to time derivatives in gravity
gravitydti=gravitydtj=gravitydtk=-100;
#if(ASYMDIAGCHECK)
dualfprintf(fail_file,"BEGINNING steppart=%d nstep=%ld\n",steppart,nstep);
dualfprintf(fail_file,"pi in advance\n");
asym_compute_1(pi);
dualfprintf(fail_file,"pb in advance\n");
asym_compute_1(pb);
#endif
/////////
//
// set initialstep and finalstep to tell some procedures and diagnostic functions if should be accounting or not
//
/////////
if(timeorder==0){
initialstep=1;
}
else{
initialstep=0;
}
if(timeorder==numtimeorders-1){
finalstep=1;
}
else{
finalstep=0;
}
/////////
//
// set dt4diag for source diagnostics
//
/////////
if(timeorder==numtimeorders-1 && (nstep%DIAGSOURCECOMPSTEP==0) ){
// every 4 full steps as well as on final step of full step (otherwise diag_source_comp() too expensive)
dt4diag=dt4diag_willbe;
dt4diag_willbe=0;
}
else{
dt4diag_willbe+=dt;
dt4diag=-1.0;
}
/////////////////////////////////////////////
//
// Setup loops [+1 extra compared to normal comp region if doing FLUXCTSTAG]
//
////////////////////////////////////////////
get_flux_startendindices(Uconsevolveloop,&is,&ie,&js,&je,&ks,&ke);
/////////////////////////////////////////////
//
// Set initial ui, temporary space, so ok that COMPZLOOP() goes over +1 in FLUXB==FLUXCTSTAG case
//
////////////////////////////////////////////
if(timeorder==0){
// last timestep's final ucum is stored into ucumformetric and ui as initial term in ucum
// copy ucum -> {ucumformetric,ui}
if(doingmetricsubstep()) copy_3dnpr_2ptrs(is,ie,js,je,ks,ke,ucum,ucumformetric,ui);
else copy_3dnpr(is,ie,js,je,ks,ke,ucum,ui); // only need to setup ui then
}
else{
// preserve this time's value of ucum for the metric (for timeorder==0 ucumformetric is assigned from ui)
// copy ucum -> ucumformetric
if(doingmetricsubstep()) copy_3dnpr(is,ie,js,je,ks,ke,ucum,ucumformetric);
}
/////////////////////////////////////////////
//
// Compute flux
//
////////////////////////////////////////////
#if(PRODUCTION==0)
trifprintf( "#0f");
#endif
if(truestep){ // only do if not just passing through
if(1){
// NORMAL:
ndt1=ndt2=ndt3=BIG;
// pb used here on a stencil, so if pb=pf or pb=pi in pointers, shouldn't change pi or pf yet -- don't currently
MYFUN(fluxcalc(stage,initialstep,finalstep,pb,pstag,pl_ct, pr_ct, vpot,F1,F2,F3,CUf,CUnew,fluxdt,fluxtime,&ndt1,&ndt2,&ndt3),"advance.c:advance_standard()", "fluxcalcall", 1);
}
#if(0)// DEBUG:
if(1){
ndt1donor=ndt2donor=ndt3donor=BIG;
MYFUN(fluxcalc_donor(stage,pb,pstag,pl_ct, pr_ct, vpot,GLOBALPOINT(F1EM),GLOBALPOINT(F2EM),GLOBALPOINT(F3EM),CUf,CUnew,fluxdt,fluxtime,&ndt1donor,&ndt2donor,&ndt3donor),"advance.c:advance_standard()", "fluxcalcall", 1);
}
// DEBUG:
if(1){
int i,j,k,pl,pliter;
FULLLOOP{
PLOOP(pliter,pl){
if(N1>1) MACP0A1(F1,i,j,k,pl)=GLOBALMACP0A1(F1EM,i,j,k,pl);
if(N2>1) MACP0A1(F2,i,j,k,pl)=GLOBALMACP0A1(F2EM,i,j,k,pl);
if(N3>1) MACP0A1(F3,i,j,k,pl)=GLOBALMACP0A1(F3EM,i,j,k,pl);
}
}
ndt1=ndt1donor;
ndt2=ndt2donor;
ndt3=ndt3donor;
}
#endif
#if(0)// DEBUG:
if(1){
int i,j,k,pliter,pl;
dualfprintf(fail_file, "BADLOOPCOMPARE: nstep=%ld steppart=%d\n",nstep,steppart);
dualfprintf(fail_file, "ndt1orig=%21.15g ndt1new=%21.15g ndt2orig=%21.15g ndt2new=%21.15g\n",ndt1,ndt1donor,ndt2,ndt2donor);
COMPFULLLOOP{
if(i==390 && j==1 && k==0){
dualfprintf(fail_file, "i=%d j=%d k=%d\n",i,j,k);
PLOOP(pliter,pl) dualfprintf(fail_file," pl=%d F1orig=%21.15g F1new=%21.15g :: F2orig=%21.15g F2new=%21.15g \n",pl,MACP0A1(F1,i,j,k,pl),GLOBALMACP0A1(F1EM,i,j,k,pl),MACP0A1(F2,i,j,k,pl),GLOBALMACP0A1(F2EM,i,j,k,pl));
}
}
}
#endif
#if(0)// DEBUG:
if(1){
int i,j,k,pliter,pl;
FTYPE diff1,diff2;
FTYPE sum1,sum2;
dualfprintf(fail_file, "BADLOOPCOMPARE: nstep=%ld steppart=%d\n",nstep,steppart);
dualfprintf(fail_file, "ndt1orig=%21.15g ndt1new=%21.15g ndt2orig=%21.15g ndt2new=%21.15g\n",ndt1,ndt1donor,ndt2,ndt2donor);
COMPFULLLOOP{
PLOOP(pliter,pl){
diff1=fabs(MACP0A1(F1,i,j,k,pl)-GLOBALMACP0A1(F1EM,i,j,k,pl));
sum1=fabs(MACP0A1(F1,i,j,k,pl))+fabs(GLOBALMACP0A1(F1EM,i,j,k,pl))+SMALL;
diff2=fabs(MACP0A1(F2,i,j,k,pl)-GLOBALMACP0A1(F2EM,i,j,k,pl));
sum2=fabs(MACP0A1(F2,i,j,k,pl))+fabs(GLOBALMACP0A1(F2EM,i,j,k,pl))+SMALL;
if(diff1/sum1>100.0*NUMEPSILON || diff2/sum2>100.0*NUMEPSILON){
dualfprintf(fail_file, "i=%d j=%d k=%d\n",i,j,k);
dualfprintf(fail_file," pl=%d diff1/sum1=%21.15g F1orig=%21.15g F1new=%21.15g :: diff2/sum2=%21.15g F2orig=%21.15g F2new=%21.15g \n",pl,diff1/sum1,MACP0A1(F1,i,j,k,pl),GLOBALMACP0A1(F1EM,i,j,k,pl),diff2/sum2,MACP0A1(F2,i,j,k,pl),GLOBALMACP0A1(F2EM,i,j,k,pl));
}
}
}
}
#endif
}// end if not just passing through
#if(PRODUCTION==0)
trifprintf( "1f");
#endif
// from here on, pi/pb/pf are only used a zone at a time rather than on a stencil
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
//
// now update get flux update [loop should go over normal computational region +1 on outer edges so automatically set field if staggered. Since we only set tempucum so far, ucum in that +1 region is unaffected]
//
/////////////////////////////////////////////////////
if(truestep){
// initialize uf and ucum if very first time here since ucum is cumulative (+=) [now tempucum is cumulative]
// copy 0 -> {uf,tempucum}
if(timeorder==0) init_3dnpr_2ptrs(is, ie, js, je, ks, ke,0.0, uf,tempucum);
#if(WHICHEOM==WITHNOGDET && (NOGDETB1==1 || NOGDETB2==1 || NOGDETB3==1) )
// for FLUXB==FLUXCTSTAG, assume no source term for field
if(FLUXB==FLUXCTSTAG){
dualfprintf(fail_file,"Not setup for field source term if staggered field\n");
myexit(176023);
}
#endif
#pragma omp parallel OPENMPGLOBALPRIVATEFORSTATEANDGEOM
{
int pl,pliter,i,j,k;
struct of_geom geomdontuse;
struct of_geom *ptrgeom=&geomdontuse;
FTYPE Uitemp[NPR];
FTYPE dUgeom[NPR],dUriemann[NPR],dUriemann1[NPR],dUriemann2[NPR],dUriemann3[NPR],dUcomp[NUMSOURCES][NPR];
struct of_state qdontuse;
struct of_state *qptr=&qdontuse;
struct of_state qdontuse2;
struct of_state *qptr2=&qdontuse2; // different qptr since call normal and special get_state()
OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUP(is,ie,js,je,ks,ke);
#pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
OPENMP3DLOOPBLOCK{
OPENMP3DLOOPBLOCK2IJK(i,j,k);
// set geometry for centered zone to be updated
get_geometry(i, j, k, CENT, ptrgeom);
// find Ui(pi)
// force use of primitive to set Ui since otherwise wherever corrected/changed primitive (in fixup, etc.) then would have to change conserved quantity, while same since both are point values
// only field for staggered method is special point value at faces that needs to come from conserved quantity
MYFUN(get_state(MAC(pi,i,j,k), ptrgeom, qptr),"step_ch.c:advance()", "get_state()", 1);
MYFUN(primtoU(UEVOLVE,MAC(pi,i,j,k), qptr, ptrgeom, Uitemp),"step_ch.c:advance()", "primtoU()", 1);
if(FLUXB==FLUXCTSTAG || DOENOFLUX != NOENOFLUX ){
// then field version of ui[] is stored as "conserved" value at FACE, not CENT
PLOOPNOB1(pl) MACP0A1(ui,i,j,k,pl)=Uitemp[pl]; // CENT
//PLOOPBONLY(pl) MACP0A1(ui,i,j,k,pl) is itself // FACE (see step_ch.c's setup_rktimestep and know that ui=unew for before first substep)
PLOOPNOB2(pl) MACP0A1(ui,i,j,k,pl)=Uitemp[pl]; // CENT
}
else{
PLOOP(pliter,pl) MACP0A1(ui,i,j,k,pl)=Uitemp[pl]; // all at CENT
}
// dUriemann is actually average quantity, but we treat is as a point quantity at the zone center
flux2dUavg(i,j,k,F1,F2,F3,dUriemann1,dUriemann2,dUriemann3);
PLOOP(pliter,pl) dUriemann[pl]=dUriemann1[pl]+dUriemann2[pl]+dUriemann3[pl]; // this addition is one type of avg->point mistake
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
//
// now update get source update (only affects stress-energy tensor in general)
// [Loop goes over up to +1 normal computational region for getting new staggered U if doing FLUXCTSTAG]
//
/////////////////////////////////////////////////////
// get state since both source() and dUtodt() need same state
// From pb, so different than state for Ui(pi)
MYFUN(get_stateforsource(MAC(pb,i,j,k), ptrgeom, &qptr2) ,"advance.c:()", "get_state() dir=0", 1);
// note that uf and ucum are initialized inside setup_rktimestep() before first substep
// find dU(pb)
// source() doesn't actually use CUf[2]=dt right now
MYFUN(source(MAC(pb,i,j,k), ptrgeom, qptr2, MAC(ui,i,j,k), dUriemann, dUcomp, dUgeom),"step_ch.c:advance()", "source", 1);
// assumes final dUcomp is nonzero and representative of source term over this timestep
#if(SPLITNPR)
// don't update metric if only doing B1-B3
if(advancepassnumber==-1 || advancepassnumber==1)
#endif
{
#if(ACCURATESOURCEDIAG)
if(DODIAGS){
diag_source_comp(ptrgeom,dUcomp,fluxdt);
diag_source_all(ptrgeom,dUgeom,fluxdt);
}
#else
if(DODIAGS){
diag_source_comp(ptrgeom,dUcomp,dt4diag);
diag_source_all(ptrgeom,dUgeom,dt4diag);
}
#endif
}
// Get update
dUtoU(i,j,k,ptrgeom->p,dUgeom, dUriemann, CUf, CUnew, MAC(ui,i,j,k), MAC(uf,i,j,k), MAC(tempucum,i,j,k));
// get timestep limit from acceleration
#if(LIMITDTWITHSOURCETERM)
#if(SPLITNPR)
// don't do dUtodt if only doing B1-B3
if(advancepassnumber==-1 || advancepassnumber==1)
#endif
{
// geometry is post-metric update, but should still give good estimate of future dt
dUtodt(ptrgeom, qptr2, MAC(pb,i,j,k), dUgeom, dUriemann, dUcomp[GEOMSOURCE], &accdt_ij, &gravitydt_ij);
#pragma omp critical
{
if(accdt_ij<accdt){
accdt=accdt_ij;
accdti=i;
accdtj=j;
accdtk=k;
}
if(gravitydt_ij<gravitydt){
gravitydt=gravitydt_ij;
gravitydti=i;
gravitydtj=j;
gravitydtk=k;
}
}// end critical region
}// end block that may mean: if not only doing B1-B3
#endif // end if doing LIMITDTWITHSOURCETERM
#if(FLUXDUMP)
fluxdumpdt=dt; // store current dt so when dump fluxdump use that dt instead of updated dt
fluxdumprealnstep=realnstep;
// DEBUG - DIAG:
PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,0*NPR + pl)=dUgeom[pl];
if(N1>1) PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,1*NPR + pl)=dUriemann1[pl];
else PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,1*NPR + pl)=0.0;
if(N2>1) PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,2*NPR + pl)=dUriemann2[pl];
else PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,2*NPR + pl)=0.0;
if(N3>1) PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,3*NPR + pl)=dUriemann3[pl];
else PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,3*NPR + pl)=0.0;
#endif
} // end COMPZLOOP :: end looping to obtain dUriemann and full unew update
}// end parallel block
} // end if truestep
else{
// then nothing to do since nothing changed
// previously updated dU and got new ucum as fed into metric, but now metric has its own ucummetric so all that is not needed
// SUPERGODMARK: no, I guess I don't recall what was being done for metric and why when passing through with dt==0.0
// just need to copy ui -> {uf,tempucum} to duplicate behavior of dUtoU()
copy_3dnpr_2ptrs(is,ie,js,je,ks,ke,ui,uf,tempucum);
}
#if(PRODUCTION==0)
trifprintf( "#0m");
#endif
/////////////////////////////////////////////
//
// EVOLVE (update/compute) METRIC HERE
// In general computes stress-energy tensor (T) from pb and T is then used to compute metric
// Done here instead of after flux since horizon_flux() updates flux through boundary that would change metric
// want metric to be at same time as everything else done with pb so RK method is consistent
//
// uses unew that's NOT updated yet
/////////////////////////////////////////////
if(truestep){
if(doingmetricsubstep()){
#if(SPLITNPR)
// don't update metric if only doing B1-B3
if(advancepassnumber==-1 || advancepassnumber==1)
#endif
{
compute_new_metric_substep(CUf,CUnew,pb,ucumformetric); // CHANGINGMARK : Not sure if CUnew here is correct
}
}// end if doing metric substepping
}
////////////////////////////////
//
// compute flux diagnostics (accurately using all substeps)
//
///////////////////////////////
if(truestep){
// must come after metric changes that can change where flux surface is since otherwise when flux surface changes, we won't track this substep's flux through the new surface but the old surface (which could even be at r=0 and have no flux)
// if using unew, then since metric update above uses old unew, need present flux at new horizon surface
#if(SPLITNPR)
// don't update metric if only doing B1-B3
if(advancepassnumber==-1 || advancepassnumber==1)
#endif
{
#if(ACCURATEDIAG)
diag_flux_pureflux(pb,F1, F2, F3, fluxdt); // fluxdt is true dt for this flux as added in dUtoU() as part of ucum
#endif
}
}
#if(PRODUCTION==0)
trifprintf( "#0s");
#endif
// Do comp full loop in case subset of full grid so that boundary cells of subset domain are defined without having to copy over those cells in bounds.c or initbase.gridsectioning.c
/////////////
//
// Utoprim as initial conditions : can't assume want these to be same in the end, so assign
//
// Since final step often sets pointers of pi=pf, in order to use arbitrary guess, must set here once done with pi,pb,pf.
//
////////////
// copy pb->pf
copy_3dnpr_fullloop(pb,pf);
///////////////////////////////////////
//
// Copy over tempucum -> ucum per pl to account for staggered field
//
// And choose which RK-quantity to invert
//
///////////////////////////////////////
if(finalstep){
if(FLUXB==FLUXCTSTAG){
// copy over new ucum in only desired locations irrespective of where tempucum was updated
// copy tempucum -> ucum
copy_tempucum_finalucum(Uconsevolveloop,tempucum,ucum);
}
utoinvert = ucum;
useducum=ucum;
}
else{
// tempucum cumulates for now
utoinvert = uf;
// tempucum cumulates for now
useducum=tempucum;
}
////////////////////////////
//
// split ZLOOP above and below to allow staggered field method
// [loop only over "centered" cells]
//
////////////////////////////
if(FLUXB==FLUXCTSTAG){
// if using staggered grid for magnetic field, then need to convert ucum to pstag to new pb/pf
// GODMARK: Use of globals
myupoint=GLOBALPOINT(upointglobal);
// first copy over all quantities as point, which is true except for fields if FLUXRECON active
// copy utoinvert -> myupoint
// copy all pl's since myupoint eventually used to invert rest of non-field quantities
copy_tempucum_finalucum(Uconsevolveloop,utoinvert,myupoint);
if(extrazones4emf && dofluxreconevolvepointfield==0){
//bound_uavg(STAGEM1,utoinvert); // DEBUG
// uses utoinvert and gets reaveraged field into myupoint
field_integrate_fluxrecon(stage, pb, utoinvert, myupoint);
}
// first pb entry is used for shock indicator, second is filled with new field
// myupoint goes in as staggered point value of magnetic flux and returned as centered point value of magnetic flux
// first pb entry is used for shock indicator, second is filled with new field
// myupoint goes in as staggered point value of magnetic flux and returned as centered point value of magnetic flux
if(1){ // must manually override to use below donor version
interpolate_ustag2fieldcent(stage, boundtime, timeorder, numtimeorders, pb, pstag, myupoint, pf);
}
else{
// interpolate_ustag2fieldcent_donor(stage, boundtime, timeorder, numtimeorders, pb, pstag, myupoint, pf);
}
////////////////////
// now myupoint contains centered point conserved quantities ready for inversion
////////////////////
}
else{
// utoinvert never reassigned from global a_utoinvert assignment since if here not doing FLUXCTSTAG
myupoint=utoinvert;
}
////////////////////////////
//
// INVERT [loop only over "centered" cells]
//
////////////////////////////
// get loop range
get_inversion_startendindices(Uconsevolveloop,&is,&ie,&js,&je,&ks,&ke);
#pragma omp parallel OPENMPGLOBALPRIVATEFORINVERSION
{
int pl,pliter,i,j,k;
struct of_geom geomdontuse;
struct of_geom *ptrgeom=&geomdontuse;
FTYPE prbefore[NPR];
struct of_newtonstats newtonstats; // not pointer
OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUP(is,ie,js,je,ks,ke);
// initialize counters
newtonstats.nstroke=newtonstats.lntries=0;
#pragma omp for schedule(OPENMPVARYENDTIMESCHEDULE(),OPENMPCHUNKSIZE(blocksize)) reduction(+: nstroke)
OPENMP3DLOOPBLOCK{
OPENMP3DLOOPBLOCK2IJK(i,j,k);
// set geometry for centered zone to be updated
get_geometry(i, j, k, CENT, ptrgeom);
// invert U->p
if(finalstep){ // last call, so ucum is cooked and ready to eat!
// store guess for diss_compute before changed by normal inversion
PALLLOOP(pl) prbefore[pl]=MACP0A1(pf,i,j,k,pl);
MYFUN(Utoprimgen(finalstep,EVOLVEUTOPRIM,UEVOLVE,MAC(myupoint,i,j,k), ptrgeom, MAC(pf,i,j,k),&newtonstats),"step_ch.c:advance()", "Utoprimgen", 1);
nstroke+=newtonstats.nstroke; newtonstats.nstroke=newtonstats.lntries=0;
#if(DODISS||DODISSVSR)
// then see what entropy inversion would give
diss_compute(EVOLVEUTOPRIM,UEVOLVE,MAC(myupoint,i,j,k),ptrgeom,prbefore,MAC(pf,i,j,k),&newtonstats);
#endif
}
else{ // otherwise still iterating on primitives
MYFUN(Utoprimgen(finalstep,EVOLVEUTOPRIM,UEVOLVE,MAC(myupoint,i,j,k), ptrgeom, MAC(pf,i,j,k),&newtonstats),"step_ch.c:advance()", "Utoprimgen", 1);
nstroke+=newtonstats.nstroke; newtonstats.nstroke=newtonstats.lntries=0;
}
#if(SPLITNPR)
// don't update metric if only doing B1-B3
if(advancepassnumber==-1 || advancepassnumber==1)
#endif
{
// immediate local (i.e. 1-zone) fix
#if(FIXUPZONES==FIXUP1ZONE)
// SUPERGODMARK: Below should differentiate beteween whether want negative densities fixed or not, but right now fixup1zone() does all
if((STEPOVERNEGU==0)||(STEPOVERNEGRHO==0)||(STEPOVERNEGRHOU==0)||(finalstep)){
MYFUN(fixup1zone(MAC(pf,i,j,k),MAC(useducum,i,j,k), ptrgeom,finalstep),"fixup.c:fixup()", "fixup1zone()", 1);
}
#endif
}
}// end COMPZLOOP
}// end parallel section
#if(ASYMDIAGCHECK)
dualfprintf(fail_file,"useducum in advance\n");
asym_compute_2(useducum);
dualfprintf(fail_file,"ENDING steppart=%d nstep=%ld\n",steppart,nstep);
#endif
/////////////////////////////////
//
// If not fixing up primitives after inversion immediately, then fix up all zones at once afterwards
//
/////////////////////////////////
#if(SPLITNPR)
// don't update metric if only doing B1-B3
if(advancepassnumber==-1 || advancepassnumber==1)
#endif
{
#if(FIXUPZONES==FIXUPALLZONES)
fixup(stage,pf,useducum,finalstep);
#endif
}
/////////////////////////////////
//
// Determine next timestep from waves, fluxes, and source updates
//
/////////////////////////////////
prepare_globaldt(truestep,ndt1,ndt2,ndt3,accdt,accdti,accdtj,accdtk,gravitydt,gravitydti,gravitydtj,gravitydtk,ndt);
#if(PRODUCTION==0)
trifprintf( "2f");
#endif
return (0);
}
// finite volume method NOT SETUP FOR CONSISTENT METRIC EVOLUTION YET -- EASY, JUST NOT DOING IT YET -- FOLLOW ABOVE AS EXAMPLE OF WHAT TO DO
// also not setup for staggered field method
static int advance_finitevolume(
int truestep,
int stage,
FTYPE (*pi)[NSTORE2][NSTORE3][NPR],
FTYPE (*pb)[NSTORE2][NSTORE3][NPR],
FTYPE (*pf)[NSTORE2][NSTORE3][NPR],
FTYPE (*pstag)[NSTORE2][NSTORE3][NPR],
FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],
FTYPE (*F1)[NSTORE2][NSTORE3][NPR],FTYPE (*F2)[NSTORE2][NSTORE3][NPR],FTYPE (*F3)[NSTORE2][NSTORE3][NPR],
FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],
FTYPE (*ui)[NSTORE2][NSTORE3][NPR],
FTYPE (*uf)[NSTORE2][NSTORE3][NPR],
FTYPE (*ucum)[NSTORE2][NSTORE3][NPR],
FTYPE *CUf,
FTYPE *CUnew,
SFTYPE fluxdt,
SFTYPE boundtime,
SFTYPE fluxtime,
int timeorder, int numtimeorders,
FTYPE *ndt)
{
int sc;
FTYPE ndt1, ndt2, ndt3;
FTYPE dUtot;
FTYPE idx1,idx2;
SFTYPE dt4diag;
static SFTYPE dt4diag_willbe=0;
int finalstep,initialstep;
FTYPE accdt, accdt_ij;
int accdti,accdtj,accdtk;
FTYPE gravitydt, gravitydt_ij;
int gravitydti,gravitydtj,gravitydtk;
int enerregion;
int *localenerpos;
int jj;
FTYPE (*utoinvert)[NSTORE2][NSTORE3][NPR];
FTYPE (*tempucum)[NSTORE2][NSTORE3][NPR];
FTYPE (*useducum)[NSTORE2][NSTORE3][NPR];
FTYPE (*myupoint)[NSTORE2][NSTORE3][NPR];
FTYPE (*ucumformetric)[NSTORE2][NSTORE3][NPR];
FTYPE (*dUgeomarray)[NSTORE2][NSTORE3][NPR];
// staggered field function:
int whichpltoavg[NPR];
int ifnotavgthencopy[NPR];
int docons,dosource;
int locpl[NPR];
int is,ie,js,je,ks,ke;
int doingextrashiftforstag;
if(FLUXB==FLUXCTSTAG){
tempucum=GLOBALPOINT(utemparray);
useducum=tempucum; // unless changed
}
else{
// no special shifting of indices occurs that requires loop shifting
tempucum=ucum;
useducum=ucum;
}
ucumformetric=GLOBALPOINT(ucumformetric);// temporary space for ucum for metric that is same time as "pb", so not updated yet or is ui
dUgeomarray=GLOBALPOINT(dUgeomarray);// temporary space for dUgeomarray
{// begin block
int pl,pliter;
// any cons
docons=0;
PLOOP(pliter,pl) docons+=do_conserved_integration[pl];
docons*=(DOENOFLUX == ENOFINITEVOLUME);
// any source
dosource=0;
PLOOP(pliter,pl) dosource+=do_source_integration[pl];
dosource*=(DOENOFLUX == ENOFINITEVOLUME);
}// end block
/////////////////////////////////////////////
//
// Setup loops and dt's
//
////////////////////////////////////////////
// for CZLOOP:
// avoid looping over region outside active+ghost grid
// good because somewhat general and avoid bad inversions, etc.
enerregion=TRUEGLOBALWITHBNDENERREGION;
localenerpos=enerposreg[enerregion];
accdt=BIG; // initially no limit to dt due to acceleration
accdti=accdtj=accdtk=-100;
gravitydt=BIG; // initially no limit to dt due to time-changes in gravity
gravitydti=gravitydtj=gravitydtk=-100;
/////////
//
// set initialstep and finalstep to tell some procedures and diagnostic functions if should be accounting or not
//
/////////
if(timeorder==0){
initialstep=1;
}
else{
initialstep=0;
}
if(timeorder==numtimeorders-1){
finalstep=1;