Skip to content

Commit f950115

Browse files
pseudotensorpseudotensor
committed
Updated fixes for kaz+gravity.
1 parent d0aba75 commit f950115

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

52 files changed

+2719
-1071
lines changed

‎advance.c

+116-77
Original file line numberDiff line numberDiff line change
@@ -136,6 +136,7 @@ int advance(int stage, FTYPE (*pi)[NSTORE2][NSTORE3][NPR],FTYPE (*pb)[NSTORE2][N
136136

137137
// this method guarantees conservation of non-sourced conserved quantities when metric is time-dependent
138138
// this method has updated field staggered method
139+
// Note that when dt==0.0, assume no fluxing, just take ucum -> ui -> {uf,ucum} and invert. Used with metric update.
139140
static int advance_standard(
140141
int stage,
141142
FTYPE (*pi)[NSTORE2][NSTORE3][NPR],
@@ -158,8 +159,8 @@ static int advance_standard(
158159
FTYPE ndt1, ndt2, ndt3;
159160
FTYPE dUtot;
160161
FTYPE idx1,idx2;
161-
SFTYPE dt4diag,dt4diag_forcomp;
162-
static SFTYPE dt4diag_forcomp_willbe=0;
162+
SFTYPE dt4diag;
163+
static SFTYPE dt4diag_willbe=0;
163164
int finalstep;
164165
FTYPE accdt, accdt_ij;
165166
int accdti,accdtj,accdtk;
@@ -234,26 +235,30 @@ static int advance_standard(
234235

235236
/////////
236237
//
237-
// tells diagnostics functions if should be accounting or not
238+
// set finalstep to tell some procedures and diagnostic functions if should be accounting or not
238239
//
239240
/////////
240241
if(timeorder==numtimeorders-1){
241-
dt4diag=dt;
242242
finalstep=1;
243243
}
244244
else{
245-
dt4diag=-1.0;
246245
finalstep=0;
247246
}
248247

248+
249+
/////////
250+
//
251+
// set dt4diag for source diagnostics
252+
//
253+
/////////
249254
if(timeorder==numtimeorders-1 && (nstep%DIAGSOURCECOMPSTEP==0) ){
250255
// every 4 full steps as well as on final step of full step (otherwise diag_source_comp() too expensive)
251-
dt4diag_forcomp=dt4diag_forcomp_willbe;
252-
dt4diag_forcomp_willbe=0;
256+
dt4diag=dt4diag_willbe;
257+
dt4diag_willbe=0;
253258
}
254259
else{
255-
dt4diag_forcomp_willbe+=dt;
256-
dt4diag_forcomp=-1.0;
260+
dt4diag_willbe+=dt;
261+
dt4diag=-1.0;
257262
}
258263

259264

@@ -298,71 +303,74 @@ static int advance_standard(
298303

299304

300305

301-
if(1){
302-
// NORMAL:
303-
ndt1=ndt2=ndt3=BIG;
304-
// 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
305-
MYFUN(fluxcalc(stage,pb,pstag,pl_ct, pr_ct, vpot,F1,F2,F3,CUf[2],fluxdt,&ndt1,&ndt2,&ndt3),"advance.c:advance_standard()", "fluxcalcall", 1);
306-
}
306+
if(dt!=0.0){ // only do if not just passing through
307+
308+
if(1){
309+
// NORMAL:
310+
ndt1=ndt2=ndt3=BIG;
311+
// 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
312+
MYFUN(fluxcalc(stage,pb,pstag,pl_ct, pr_ct, vpot,F1,F2,F3,CUf[2],fluxdt,&ndt1,&ndt2,&ndt3),"advance.c:advance_standard()", "fluxcalcall", 1);
313+
}
307314

308315
#if(0)// DEBUG:
309-
if(1){
310-
ndt1donor=ndt2donor=ndt3donor=BIG;
311-
MYFUN(fluxcalc_donor(stage,pb,pstag,pl_ct, pr_ct, vpot,GLOBALPOINT(F1EM),GLOBALPOINT(F2EM),GLOBALPOINT(F3EM),CUf[2],fluxdt,&ndt1donor,&ndt2donor,&ndt3donor),"advance.c:advance_standard()", "fluxcalcall", 1);
312-
}
313-
// DEBUG:
314-
if(1){
315-
int i,j,k,pl,pliter;
316-
FULLLOOP{
317-
PLOOP(pliter,pl){
318-
if(N1>1) MACP0A1(F1,i,j,k,pl)=GLOBALMACP0A1(F1EM,i,j,k,pl);
319-
if(N2>1) MACP0A1(F2,i,j,k,pl)=GLOBALMACP0A1(F2EM,i,j,k,pl);
320-
if(N3>1) MACP0A1(F3,i,j,k,pl)=GLOBALMACP0A1(F3EM,i,j,k,pl);
316+
if(1){
317+
ndt1donor=ndt2donor=ndt3donor=BIG;
318+
MYFUN(fluxcalc_donor(stage,pb,pstag,pl_ct, pr_ct, vpot,GLOBALPOINT(F1EM),GLOBALPOINT(F2EM),GLOBALPOINT(F3EM),CUf[2],fluxdt,&ndt1donor,&ndt2donor,&ndt3donor),"advance.c:advance_standard()", "fluxcalcall", 1);
319+
}
320+
// DEBUG:
321+
if(1){
322+
int i,j,k,pl,pliter;
323+
FULLLOOP{
324+
PLOOP(pliter,pl){
325+
if(N1>1) MACP0A1(F1,i,j,k,pl)=GLOBALMACP0A1(F1EM,i,j,k,pl);
326+
if(N2>1) MACP0A1(F2,i,j,k,pl)=GLOBALMACP0A1(F2EM,i,j,k,pl);
327+
if(N3>1) MACP0A1(F3,i,j,k,pl)=GLOBALMACP0A1(F3EM,i,j,k,pl);
328+
}
321329
}
330+
ndt1=ndt1donor;
331+
ndt2=ndt2donor;
332+
ndt3=ndt3donor;
322333
}
323-
ndt1=ndt1donor;
324-
ndt2=ndt2donor;
325-
ndt3=ndt3donor;
326-
}
327334
#endif
328335

329336

330337
#if(0)// DEBUG:
331-
if(1){
332-
int i,j,k,pliter,pl;
333-
dualfprintf(fail_file, "BADLOOPCOMPARE: nstep=%ld steppart=%d\n",nstep,steppart);
334-
dualfprintf(fail_file, "ndt1orig=%21.15g ndt1new=%21.15g ndt2orig=%21.15g ndt2new=%21.15g\n",ndt1,ndt1donor,ndt2,ndt2donor);
335-
COMPFULLLOOP{
336-
if(i==390 && j==1 && k==0){
337-
dualfprintf(fail_file, "i=%d j=%d k=%d\n",i,j,k);
338-
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));
338+
if(1){
339+
int i,j,k,pliter,pl;
340+
dualfprintf(fail_file, "BADLOOPCOMPARE: nstep=%ld steppart=%d\n",nstep,steppart);
341+
dualfprintf(fail_file, "ndt1orig=%21.15g ndt1new=%21.15g ndt2orig=%21.15g ndt2new=%21.15g\n",ndt1,ndt1donor,ndt2,ndt2donor);
342+
COMPFULLLOOP{
343+
if(i==390 && j==1 && k==0){
344+
dualfprintf(fail_file, "i=%d j=%d k=%d\n",i,j,k);
345+
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));
346+
}
339347
}
340348
}
341-
}
342349
#endif
343350

344351
#if(0)// DEBUG:
345-
if(1){
346-
int i,j,k,pliter,pl;
347-
FTYPE diff1,diff2;
348-
FTYPE sum1,sum2;
349-
dualfprintf(fail_file, "BADLOOPCOMPARE: nstep=%ld steppart=%d\n",nstep,steppart);
350-
dualfprintf(fail_file, "ndt1orig=%21.15g ndt1new=%21.15g ndt2orig=%21.15g ndt2new=%21.15g\n",ndt1,ndt1donor,ndt2,ndt2donor);
351-
COMPFULLLOOP{
352-
PLOOP(pliter,pl){
353-
diff1=fabs(MACP0A1(F1,i,j,k,pl)-GLOBALMACP0A1(F1EM,i,j,k,pl));
354-
sum1=fabs(MACP0A1(F1,i,j,k,pl))+fabs(GLOBALMACP0A1(F1EM,i,j,k,pl))+SMALL;
355-
diff2=fabs(MACP0A1(F2,i,j,k,pl)-GLOBALMACP0A1(F2EM,i,j,k,pl));
356-
sum2=fabs(MACP0A1(F2,i,j,k,pl))+fabs(GLOBALMACP0A1(F2EM,i,j,k,pl))+SMALL;
357-
if(diff1/sum1>100.0*NUMEPSILON || diff2/sum2>100.0*NUMEPSILON){
358-
dualfprintf(fail_file, "i=%d j=%d k=%d\n",i,j,k);
359-
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));
352+
if(1){
353+
int i,j,k,pliter,pl;
354+
FTYPE diff1,diff2;
355+
FTYPE sum1,sum2;
356+
dualfprintf(fail_file, "BADLOOPCOMPARE: nstep=%ld steppart=%d\n",nstep,steppart);
357+
dualfprintf(fail_file, "ndt1orig=%21.15g ndt1new=%21.15g ndt2orig=%21.15g ndt2new=%21.15g\n",ndt1,ndt1donor,ndt2,ndt2donor);
358+
COMPFULLLOOP{
359+
PLOOP(pliter,pl){
360+
diff1=fabs(MACP0A1(F1,i,j,k,pl)-GLOBALMACP0A1(F1EM,i,j,k,pl));
361+
sum1=fabs(MACP0A1(F1,i,j,k,pl))+fabs(GLOBALMACP0A1(F1EM,i,j,k,pl))+SMALL;
362+
diff2=fabs(MACP0A1(F2,i,j,k,pl)-GLOBALMACP0A1(F2EM,i,j,k,pl));
363+
sum2=fabs(MACP0A1(F2,i,j,k,pl))+fabs(GLOBALMACP0A1(F2EM,i,j,k,pl))+SMALL;
364+
if(diff1/sum1>100.0*NUMEPSILON || diff2/sum2>100.0*NUMEPSILON){
365+
dualfprintf(fail_file, "i=%d j=%d k=%d\n",i,j,k);
366+
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));
367+
}
360368
}
361369
}
362370
}
363-
}
364371
#endif
365372

373+
}// end if not just passing through
366374

367375

368376

@@ -390,13 +398,16 @@ static int advance_standard(
390398
/////////////////////////////////////////////////////
391399

392400

393-
// initialize uf and ucum if very first time here since ucum is cumulative (+=) [now tempucum is cumulative]
394-
// copy 0 -> {uf,tempucum}
395-
if(timeorder==0) init_3dnpr_2ptrs(is, ie, js, je, ks, ke,0.0, uf,tempucum);
396-
397401

398402
if(dt!=0.0){
399403

404+
405+
// initialize uf and ucum if very first time here since ucum is cumulative (+=) [now tempucum is cumulative]
406+
// copy 0 -> {uf,tempucum}
407+
if(timeorder==0) init_3dnpr_2ptrs(is, ie, js, je, ks, ke,0.0, uf,tempucum);
408+
409+
410+
400411
#if(WHICHEOM==WITHNOGDET && (NOGDETB1==1 || NOGDETB2==1 || NOGDETB3==1) )
401412
// for FLUXB==FLUXCTSTAG, assume no source term for field
402413
if(FLUXB==FLUXCTSTAG){
@@ -487,7 +498,7 @@ static int advance_standard(
487498
}
488499
#else
489500
if(DODIAGS){
490-
diag_source_comp(ptrgeom,dUcomp,dt4diag_forcomp);
501+
diag_source_comp(ptrgeom,dUcomp,dt4diag);
491502
diag_source_all(ptrgeom,dUgeom,dt4diag);
492503
}
493504
#endif
@@ -557,6 +568,9 @@ static int advance_standard(
557568
// then nothing to do since nothing changed
558569
// previously updated dU and got new ucum as fed into metric, but now metric has its own ucummetric so all that is not needed
559570
// SUPERGODMARK: no, I guess I don't recall what was being done for metric and why when passing through with dt==0.0
571+
572+
// just need to copy ui -> {uf,tempucum} to duplicate behavior of dUtoU()
573+
copy_3dnpr_2ptrs(is,ie,js,je,ks,ke,ui,uf,tempucum);
560574
}
561575

562576

@@ -876,8 +890,8 @@ static int advance_finitevolume(int stage,
876890
FTYPE ndt1, ndt2, ndt3;
877891
FTYPE dUtot;
878892
FTYPE idx1,idx2;
879-
SFTYPE dt4diag,dt4diag_forcomp;
880-
static SFTYPE dt4diag_forcomp_willbe=0;
893+
SFTYPE dt4diag;
894+
static SFTYPE dt4diag_willbe=0;
881895
int finalstep;
882896
FTYPE accdt, accdt_ij;
883897
int accdti,accdtj,accdtk;
@@ -953,26 +967,29 @@ static int advance_finitevolume(int stage,
953967

954968
/////////
955969
//
956-
// tells diagnostics functions if should be accounting or not
970+
// set finalstep to tell some procedures and diagnostic functions if should be accounting or not
957971
//
958972
/////////
959973
if(timeorder==numtimeorders-1){
960-
dt4diag=dt;
961974
finalstep=1;
962975
}
963976
else{
964-
dt4diag=-1.0;
965977
finalstep=0;
966978
}
967979

980+
/////////
981+
//
982+
// set dt4diag for source diagnostics
983+
//
984+
/////////
968985
if(timeorder==numtimeorders-1 && (nstep%DIAGSOURCECOMPSTEP==0) ){
969986
// every 4 full steps as well as on final step of full step (otherwise diag_source_comp() too expensive)
970-
dt4diag_forcomp=dt4diag_forcomp_willbe;
971-
dt4diag_forcomp_willbe=0;
987+
dt4diag=dt4diag_willbe;
988+
dt4diag_willbe=0;
972989
}
973990
else{
974-
dt4diag_forcomp_willbe+=dt;
975-
dt4diag_forcomp=-1.0;
991+
dt4diag_willbe+=dt;
992+
dt4diag=-1.0;
976993
}
977994

978995

@@ -1105,6 +1122,27 @@ static int advance_finitevolume(int stage,
11051122
}// end COMPZLOOP
11061123

11071124

1125+
// Store diagnostics related to component form of sources. Done here since don't integrate-up compnents of source. So only point accurate
1126+
#if(SPLITNPR)
1127+
// don't update metric if only doing B1-B3
1128+
if(advancepassnumber==-1 || advancepassnumber==1)
1129+
#endif
1130+
{
1131+
#pragma omp critical // since diagnostics store in same global cumulative variables
1132+
{
1133+
#if(ACCURATESOURCEDIAG)
1134+
if(DODIAGS){
1135+
diag_source_comp(ptrgeom,dUcomp,fluxdt);
1136+
}
1137+
#else
1138+
if(DODIAGS){
1139+
diag_source_comp(ptrgeom,dUcomp,dt4diag);
1140+
}
1141+
#endif
1142+
}
1143+
}
1144+
1145+
11081146
// volume integrate dUgeom
11091147
// c2a_1 c2a_2 c2a_3
11101148
if(dosource){
@@ -1175,6 +1213,8 @@ static int advance_finitevolume(int stage,
11751213
// get source term (volume integrated)
11761214
PLOOP(pliter,pl) dUgeom[pl]=MACP0A1(dUgeomarray,i,j,k,pl);
11771215

1216+
1217+
// store diagnostics related to source. Totals, so use volume-integrated source.
11781218
#if(SPLITNPR)
11791219
// don't update metric if only doing B1-B3
11801220
if(advancepassnumber==-1 || advancepassnumber==1)
@@ -1183,14 +1223,13 @@ static int advance_finitevolume(int stage,
11831223
#pragma omp critical // since diagnostics store in same global cumulative variables
11841224
{
11851225
#if(ACCURATESOURCEDIAG)
1186-
if(DODIAGS){
1187-
// do diagnostics for volume integrated source term
1188-
diag_source_all(ptrgeom,dUgeom,fluxdt);
1189-
}
1226+
if(DODIAGS){
1227+
diag_source_all(ptrgeom,dUgeom,fluxdt);
1228+
}
11901229
#else
1191-
if(DODIAGS){
1192-
diag_source_all(ptrgeom,dUgeom,dt4diag);
1193-
}
1230+
if(DODIAGS){
1231+
diag_source_all(ptrgeom,dUgeom,dt4diag);
1232+
}
11941233
#endif
11951234
}
11961235
}

‎bounds.tools.c

+6-6
Original file line numberDiff line numberDiff line change
@@ -1593,9 +1593,9 @@ int extrapfunc(int boundary, int j,int k,
15931593
//////////////////////////
15941594
PBOUNDLOOP(pliter,pl){
15951595
// determine MINM slope for extrapolation
1596-
Dqp=log(fabs(MACP0A1(prim,ri3,rj,rk,pl)))-log(fabs(MACP0A1(prim,ri2,rj,rk,pl)));
1597-
Dqm=log(fabs(MACP0A1(prim,ri2,rj,rk,pl)))-log(fabs(MACP0A1(prim,ri,rj,rk,pl)));
1598-
Dqc=log(fabs(MACP0A1(prim,ri3,rj,rk,pl)))-log(fabs(MACP0A1(prim,ri,rj,rk,pl)));
1596+
Dqp=log(SMALL+fabs(MACP0A1(prim,ri3,rj,rk,pl)))-log(SMALL+fabs(MACP0A1(prim,ri2,rj,rk,pl)));
1597+
Dqm=log(SMALL+fabs(MACP0A1(prim,ri2,rj,rk,pl)))-log(SMALL+fabs(MACP0A1(prim,ri,rj,rk,pl)));
1598+
Dqc=log(SMALL+fabs(MACP0A1(prim,ri3,rj,rk,pl)))-log(SMALL+fabs(MACP0A1(prim,ri,rj,rk,pl)));
15991599
dqlogdensity[pl] = signdq*MINMOD(MINMOD(Dqp,Dqm),Dqc);
16001600
}
16011601

@@ -1746,7 +1746,7 @@ int extrapfunc(int boundary, int j,int k,
17461746

17471747
// log extrap (very speculative and can cause problems if used alone when (say) density is super low on ri+1 and relatively high on ri, then i will be super huge
17481748
// should use this when values that go into slope are much different, or equally when dqlogdensity is large
1749-
MACP0A1(prim,i,j,k,pl) = exp(log(fabs(MACP0A1(prim,ri,rj,rk,pl))) + mydqlog*(i-ri));
1749+
MACP0A1(prim,i,j,k,pl) = exp(log(SMALL+fabs(MACP0A1(prim,ri,rj,rk,pl))) + mydqlog*(i-ri));
17501750

17511751
// DEBUG:
17521752
// dualfprintf(fail_file,"i=%d j=%d pl=%d ftemp=%21.15g linearvalue=%21.15g expvalue=%21.15g final=%21.15g\n",i,j,pl,ftemp,linearvalue,expvalue,MACP0A1(prim,i,j,k,pl));
@@ -1763,8 +1763,8 @@ int extrapfunc(int boundary, int j,int k,
17631763
pl=UU; MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) + dq[pl]*(i-ri);
17641764
#else
17651765
// log extrap
1766-
pl=RHO; MACP0A1(prim,i,j,k,pl) = exp(log(fabs(MACP0A1(prim,ri,rj,rk,pl))) + dqlogdensity[pl]*(i-ri));
1767-
pl=UU; MACP0A1(prim,i,j,k,pl) = exp(log(fabs(MACP0A1(prim,ri,rj,rk,pl))) + dqlogdensity[pl]*(i-ri));
1766+
pl=RHO; MACP0A1(prim,i,j,k,pl) = exp(log(SMALL+fabs(MACP0A1(prim,ri,rj,rk,pl))) + dqlogdensity[pl]*(i-ri));
1767+
pl=UU; MACP0A1(prim,i,j,k,pl) = exp(log(SMALL+fabs(MACP0A1(prim,ri,rj,rk,pl))) + dqlogdensity[pl]*(i-ri));
17681768
#endif
17691769

17701770
#endif

0 commit comments

Comments
 (0)