29 #if(ANALYTICMEMORY==0)
30 dualfprintf(
fail_file,
"ANALYTICMEMORY==0 but called bound_x1dn_analytic\n");
46 pl=
B1;
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
49 PBOUNDLOOP(pliter,pl)
if(pl!=
B1)
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
79 #if(ANALYTICMEMORY==0)
80 dualfprintf(
fail_file,
"ANALYTICMEMORY==0 but called bound_x1up_analytic\n");
94 pl=
B1;
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
97 PBOUNDLOOP(pliter,pl)
if(pl!=
B1)
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
104 PBOUNDLOOP(pliter,pl)
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(panalytic,i,j,k,pl);
126 #if(ANALYTICMEMORY==0)
127 dualfprintf(
fail_file,
"ANALYTICMEMORY==0 but called bound_x2dn_analytic\n");
144 pl=
B2;
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
147 PBOUNDLOOP(pliter,pl)
if(pl!=
B2)
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
154 PBOUNDLOOP(pliter,pl)
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(panalytic,i,j,k,pl);
177 #if(ANALYTICMEMORY==0)
178 dualfprintf(
fail_file,
"ANALYTICMEMORY==0 but called bound_x2up_analytic\n");
192 pl=
B2;
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
195 PBOUNDLOOP(pliter,pl)
if(pl!=
B2)
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
202 PBOUNDLOOP(pliter,pl)
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(panalytic,i,j,k,pl);
225 #if(ANALYTICMEMORY==0)
226 dualfprintf(
fail_file,
"ANALYTICMEMORY==0 but called bound_x3dn_analytic\n");
243 pl=
B3;
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
246 PBOUNDLOOP(pliter,pl)
if(pl!=
B3)
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
253 PBOUNDLOOP(pliter,pl)
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(panalytic,i,j,k,pl);
276 #if(ANALYTICMEMORY==0)
277 dualfprintf(
fail_file,
"ANALYTICMEMORY==0 but called bound_x3up_analytic\n");
291 pl=
B3;
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
294 PBOUNDLOOP(pliter,pl)
if(pl!=
B3)
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
301 PBOUNDLOOP(pliter,pl)
MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(panalytic,i,j,k,pl);
314 int bound_x1dn_outflow_simple(
315 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
330 #pragma omp parallel // assume don't require EOS
333 struct of_geom geom[NPR],rgeom[NPR];
341 struct of_geom geomdontuse[NPR];
343 struct of_geom rgeomdontuse[NPR];
349 ptrgeom[pl]=&(geomdontuse[pl]);
350 ptrrgeom[pl]=&(rgeomdontuse[pl]);
366 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
394 inflow_check_4vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[
U1], 0) ;
395 #elif(WHICHVEL==VEL3)
397 inflow_check_3vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[
U1], 0) ;
402 failreturn=check_pr(
MAC(prim,i,j,k),
MAC(prim,i,j,k),ptrgeom[
U1],-3);
409 #elif(WHICHVEL==VELREL4)
411 inflow_check_rel4vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[
U1],0) ;
428 dualfprintf(
fail_file,
"Shouldn't be here in bounds: dir=%d\n",
X1DN);
440 int bound_x1up_outflow_simple(
441 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
447 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
457 #pragma omp parallel // assume don't require EOS
468 struct of_geom geomdontuse[NPR];
470 struct of_geom rgeomdontuse[NPR];
475 ptrgeom[pl]=&(geomdontuse[pl]);
476 ptrrgeom[pl]=&(rgeomdontuse[pl]);
493 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
519 inflow_check_4vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[
U1],0) ;
520 #elif(WHICHVEL==VEL3)
522 inflow_check_3vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[
U1],0) ;
527 failreturn=check_pr(
MAC(prim,i,j,k),
MAC(prim,i,j,k),ptrgeom[
U1],-3);
534 #elif(WHICHVEL==VELREL4)
536 inflow_check_rel4vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[
U1],0) ;
552 dualfprintf(
fail_file,
"Shouldn't be here in bounds: dir=%d\n",
X1UP);
579 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
585 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
594 #pragma omp parallel // assume don't require EOS
597 struct of_geom geom[NPR],rgeom[NPR];
606 struct of_geom geomdontuse[NPR];
608 struct of_geom rgeomdontuse[NPR];
614 ptrgeom[pl]=&(geomdontuse[pl]);
615 ptrrgeom[pl]=&(rgeomdontuse[pl]);
627 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
638 #if(HORIZONEXTRAP==0)
644 #elif(HORIZONEXTRAP==1)
652 if(pl>=
RHO && pl<=
UU){
653 MACP0A1(prim,i,j,k,pl) =
MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]) ;
659 for(pl=
U2;pl<=
U3;pl++){
666 MACP0A1(prim,i,j,k,pl) =
MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]);
667 for(pl=
B2;pl<=
B3;pl++){
673 if(pl>=
B3+1 && pl<NPRBOUND){
674 MACP0A1(prim,i,j,k,pl) =
MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]);
678 #elif(HORIZONEXTRAP==2)
680 rescale(
DORESCALE,1,
MAC(prim,ri,rj,rk),ptrrgeom[0],prescale);
685 rescale(
UNRESCALE,1,
MAC(prim,i,j,k),ptrgeom[0],prescale);
687 #elif(HORIZONEXTRAP==3)
688 extrapfunc(
X1DN,j,k,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,
inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
701 inflow_check_4vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[U1], 0) ;
702 #elif(WHICHVEL==VEL3)
704 inflow_check_3vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[U1], 0) ;
709 failreturn=check_pr(
MAC(prim,i,j,k),
MAC(prim,i,j,k),ptrgeom[U1],-3);
716 #elif(WHICHVEL==VELREL4)
718 inflow_check_rel4vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
735 dualfprintf(
fail_file,
"Shouldn't be here in bounds: horizonextrap version: dir=%d\n",
X1DN);
759 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
765 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
775 #pragma omp parallel // assume don't require EOS
786 struct of_geom geomdontuse[NPR];
788 struct of_geom rgeomdontuse[NPR];
793 ptrgeom[pl]=&(geomdontuse[pl]);
794 ptrrgeom[pl]=&(rgeomdontuse[pl]);
807 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
821 #elif(OUTEREXTRAP==1)
828 if(pl>=
RHO && pl<=
UU){
829 MACP0A1(prim,i,j,k,pl) =
MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]) ;
834 MACP0A1(prim,i,j,k,pl) =
MACP0A1(prim,ri,rj,rk,pl) * (1. - 2*(i-ri)*
dx[1]) ;
835 for(pl=
U2;pl<=
U3;pl++){
842 MACP0A1(prim,i,j,k,pl) =
MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]) ;
843 for(pl=
B2;pl<=
B3;pl++){
849 if(pl>=
B3+1 && pl<NPRBOUND){
850 MACP0A1(prim,i,j,k,pl) =
MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]) ;
854 #elif(OUTEREXTRAP==2)
856 rescale(
DORESCALE,1,
MAC(prim,ri,rj,rk),ptrrgeom[0],prescale);
861 rescale(
UNRESCALE,1,
MAC(prim,i,j,k),ptrgeom[0],prescale);
863 #elif(OUTEREXTRAP==3)
864 extrapfunc(
X1UP,j,k,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,
inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
878 inflow_check_4vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
879 #elif(WHICHVEL==VEL3)
881 inflow_check_3vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
886 failreturn=check_pr(
MAC(prim,i,j,k),
MAC(prim,i,j,k),ptrgeom[U1],-3);
893 #elif(WHICHVEL==VELREL4)
895 inflow_check_rel4vel(1,
MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
911 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
928 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
934 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
942 #pragma omp parallel // assume don't require EOS
966 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
992 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1000 if(pl==
U1 || pl==URAD1 || pl==
B1){
1009 if(pl==
U1 || pl==URAD1 || pl==
B1){
1010 MACP0A1(prim,i,j,k,pl) *= -1.;
1019 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
1026 bound_x1dn_r0singfixinterior(boundstage,finalstep, boundtime, whichdir, boundvartype, dirprim,ispstag, prim,inboundloop,outboundloop,innormalloop,outnormalloop,
inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1039 int bound_x2dn_outflow_simple(
1040 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
1046 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
1055 #pragma omp parallel // assume don't require EOS
1057 int i,
j,k,pl,pliter;
1058 int localj1,globalj1;
1059 int localj2,globalj2;
1061 struct of_geom geom[NPR],rgeom[NPR];
1067 FTYPE prescale[NPR];
1069 struct of_geom geomdontuse[NPR];
1071 struct of_geom rgeomdontuse[NPR];
1072 struct of_geom *ptrrgeom[NPR];
1077 ptrgeom[pl]=&(geomdontuse[pl]);
1078 ptrrgeom[pl]=&(rgeomdontuse[pl]);
1090 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1120 inflow_check_4vel(2,
MAC(prim,i,j,k),NULL,ptrgeom[
U2], 0) ;
1121 #elif(WHICHVEL==VEL3)
1123 inflow_check_3vel(2,
MAC(prim,i,j,k),NULL,ptrgeom[
U2], 0) ;
1128 failreturn=check_pr(
MAC(prim,i,j,k),
MAC(prim,i,j,k),ptrgeom[
U2],-3);
1135 #elif(WHICHVEL==VELREL4)
1137 inflow_check_rel4vel(2,
MAC(prim,i,j,k),NULL,ptrgeom[
U2],0) ;
1154 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
1166 int bound_x2up_outflow_simple(
1167 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
1173 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
1183 #pragma omp parallel // assume don't require EOS
1186 int i,
j,k,pl,pliter;
1187 int localj1,globalj1;
1188 int localj2,globalj2;
1195 FTYPE prescale[NPR];
1197 struct of_geom geomdontuse[NPR];
1199 struct of_geom rgeomdontuse[NPR];
1200 struct of_geom *ptrrgeom[NPR];
1204 ptrgeom[pl]=&(geomdontuse[pl]);
1205 ptrrgeom[pl]=&(rgeomdontuse[pl]);
1218 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1242 inflow_check_4vel(2,
MAC(prim,i,j,k),NULL,ptrgeom[
U2],0) ;
1243 #elif(WHICHVEL==VEL3)
1245 inflow_check_3vel(2,
MAC(prim,i,j,k),NULL,ptrgeom[
U2],0) ;
1250 failreturn=check_pr(
MAC(prim,i,j,k),
MAC(prim,i,j,k),ptrgeom[
U2],-3);
1257 #elif(WHICHVEL==VELREL4)
1259 inflow_check_rel4vel(2,
MAC(prim,i,j,k),NULL,ptrgeom[
U2],0) ;
1275 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
1307 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
1313 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
1322 #pragma omp parallel // assume don't require EOS
1324 int i,
j,k,pl,pliter;
1330 FTYPE prescale[NPR];
1340 if(whichcall==1 &&
ncpux3==1){
1345 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1350 rk=(k+(int)
N3/2)%
N3;
1369 #if(DEBUGINOUTLOOPS)
1370 dualfprintf(
fail_file,
"INNER pole1: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,i,j,k);
1372 dualfprintf(
fail_file,
"INNER pole1: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1400 #if(DEBUGINOUTLOOPS)
1401 dualfprintf(
fail_file,
"INNER pole2: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,i,j,k);
1403 dualfprintf(
fail_file,
"INNER pole2: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1417 if(
POLEDEATH ||
POLEGAMMADEATH)
poledeath(
X2DN,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,
inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1418 if(
POLESMOOTH)
polesmooth(
X2DN,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,
inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1430 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
1451 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
1457 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
1466 #pragma omp parallel // assume don't require EOS
1468 int i,
j,k,pl,pliter;
1474 FTYPE prescale[NPR];
1487 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1511 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1519 if(pl==
U2 || pl==URAD2 || pl==
B2){
1545 if(pl==
U2 || pl==URAD2 || pl==
B2)
MACP0A1(prim,i,j,k,pl) *= -1.0;
1555 if(
POLEDEATH ||
POLEGAMMADEATH)
poledeath(
X2DN,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,
inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1556 if(
POLESMOOTH)
polesmooth(
X2DN,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,
inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1565 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
1583 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
1589 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
1600 #pragma omp parallel // assume don't require EOSs
1602 int i,
j,k,pl,pliter;
1608 FTYPE prescale[NPR];
1619 if(whichcall==1 &&
ncpux3==1){
1624 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1630 rk=(k+(int)
N3/2)%
N3;
1635 else rj=2*(
N2-1)-j+1;
1649 #if(DEBUGINOUTLOOPS)
1650 dualfprintf(
fail_file,
"OUTER pole1: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,i,j,k);
1652 dualfprintf(
fail_file,
"OUTER pole1: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1679 #if(DEBUGINOUTLOOPS)
1680 dualfprintf(
fail_file,
"OUTER pole2: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,i,j,k);
1682 dualfprintf(
fail_file,
"OUTER pole2: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1695 if(
POLEDEATH ||
POLEGAMMADEATH)
poledeath(
X2UP,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,
inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1696 if(
POLESMOOTH)
polesmooth(
X2UP,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,
inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1703 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
1719 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
1725 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
1736 #pragma omp parallel // assume don't require EOS
1738 int i,
j,k,pl,pliter;
1744 FTYPE prescale[NPR];
1759 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1770 else rj=2*(
N2-1)-j+1;
1782 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1789 if(pl==
U2 || pl==URAD2 || pl==
B2){
1815 if(pl==
U2 || pl==URAD2 || pl==
B2)
MACP0A1(prim,i,j,k,pl) *= -1.0;
1825 if(
POLEDEATH ||
POLEGAMMADEATH)
poledeath(
X2UP,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,
inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1826 if(
POLESMOOTH)
polesmooth(
X2UP,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,
inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1835 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
1850 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
1856 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
1865 #pragma omp parallel // assume don't require EOS
1867 int i,
j,k,pl,pliter;
1873 FTYPE prescale[NPR];
1887 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1900 #if(DEBUGINOUTLOOPS)
1901 dualfprintf(
fail_file,
"INNER X1: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherri=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,ri+
SHIFT1+i,i,j,k);
1903 dualfprintf(
fail_file,
"INNER X1: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1916 #if(DEBUGINOUTLOOPS)
1917 dualfprintf(
fail_file,
"OUTER X1: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherri=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,ri+(i-
N1),rk,i,j,k);
1919 dualfprintf(
fail_file,
"INNER X1: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1930 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
1942 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
1948 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
1957 #pragma omp parallel // assume don't require EOS
1959 int i,
j,k,pl,pliter;
1965 FTYPE prescale[NPR];
1979 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1992 #if(DEBUGINOUTLOOPS)
1993 dualfprintf(
fail_file,
"INNER X2: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherrj=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,rj+
SHIFT2+j,i,j,k);
1996 dualfprintf(
fail_file,
"INNER X2: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1997 myexit(29533246346);
2010 #if(DEBUGINOUTLOOPS)
2011 dualfprintf(
fail_file,
"OUTER X2: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherrj=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,rj+(j-
N2),i,j,k);
2014 dualfprintf(
fail_file,
"INNER X2: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
2015 myexit(3294863498634);
2026 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
2042 int bound_x3dn_outflow_simple(
2043 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
2049 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
2058 #pragma omp parallel // assume don't require EOS
2060 int i,
j,k,pl,pliter;
2061 int localj1,globalj1;
2062 int localk2,globalk2;
2064 struct of_geom geom[NPR],rgeom[NPR];
2070 FTYPE prescale[NPR];
2072 struct of_geom geomdontuse[NPR];
2074 struct of_geom rgeomdontuse[NPR];
2075 struct of_geom *ptrrgeom[NPR];
2080 ptrgeom[pl]=&(geomdontuse[pl]);
2081 ptrrgeom[pl]=&(rgeomdontuse[pl]);
2093 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2120 inflow_check_4vel(3,
MAC(prim,i,j,k),NULL,ptrgeom[
U3], 0) ;
2121 #elif(WHICHVEL==VEL3)
2123 inflow_check_3vel(3,
MAC(prim,i,j,k),NULL,ptrgeom[
U3], 0) ;
2128 failreturn=check_pr(
MAC(prim,i,j,k),
MAC(prim,i,j,k),ptrgeom[
U3],-3);
2135 #elif(WHICHVEL==VELREL4)
2137 inflow_check_rel4vel(3,
MAC(prim,i,j,k),NULL,ptrgeom[
U3],0) ;
2154 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
2166 int bound_x3up_outflow_simple(
2167 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
2173 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
2183 #pragma omp parallel // assume don't require EOS
2186 int i,
j,k,pl,pliter;
2187 int localk1,globalk1;
2188 int localk2,globalk2;
2195 FTYPE prescale[NPR];
2197 struct of_geom geomdontuse[NPR];
2199 struct of_geom rgeomdontuse[NPR];
2200 struct of_geom *ptrrgeom[NPR];
2204 ptrgeom[pl]=&(geomdontuse[pl]);
2205 ptrrgeom[pl]=&(rgeomdontuse[pl]);
2218 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2242 inflow_check_4vel(3,
MAC(prim,i,j,k),NULL,ptrgeom[
U3],0) ;
2243 #elif(WHICHVEL==VEL3)
2245 inflow_check_3vel(3,
MAC(prim,i,j,k),NULL,ptrgeom[
U3],0) ;
2250 failreturn=check_pr(
MAC(prim,i,j,k),
MAC(prim,i,j,k),ptrgeom[
U3],-3);
2257 #elif(WHICHVEL==VELREL4)
2259 inflow_check_rel4vel(3,
MAC(prim,i,j,k),NULL,ptrgeom[
U3],0) ;
2275 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
2306 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
2312 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
2321 #pragma omp parallel // assume don't require EOS
2323 int i,
j,k,pl,pliter;
2329 FTYPE prescale[NPR];
2343 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2356 #if(DEBUGINOUTLOOPS)
2357 dualfprintf(
fail_file,
"INNER X3: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherrk=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,rk+
SHIFT3+k,i,j,k);
2359 dualfprintf(
fail_file,
"INNER X3: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
2372 #if(DEBUGINOUTLOOPS)
2373 dualfprintf(
fail_file,
"OUTER X3: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherrk=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,rk+(k-
N3),i,j,k);
2375 dualfprintf(
fail_file,
"INNER X3: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
2386 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
2401 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
2407 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
2413 int i,
j,k,pl,pliter;
2414 int locali1,globali1;
2415 int locali2,globali2;
2422 FTYPE prescale[NPR];
2451 else if(globali1>
endpos[1]+N1BND-1) locali1=
N1+N1BND-1;
2452 else locali1=horizoni-
N1BND;
2460 if(globali2<
startpos[1]) locali2=0;
2461 else if(globali2>
endpos[1]-1) locali2=
N1-1;
2486 #if(UTOPRIMADJUST==UTOPRIMAVG)
2495 if(
MACP0A1(prim,i,j,k,
U1)>0)
MACP0A1(prim,i,j,k,
U1)=-0.5*(fabs(
MACP0A1(prim,ri1,rj,rk,
U1))+fabs(
MACP0A1(prim,ri1+1,rj,rk,
U1)));
2504 if(
MACP0A1(prim,i,j,k,
U1)>0)
MACP0A1(prim,i,j,k,
U1)=-0.5*(fabs(
MACP0A1(prim,ri1,rj,rk,
U1))+fabs(
MACP0A1(prim,ri1+1,rj,rk,
U1)));
2515 dualfprintf(
fail_file,
"Shouldn't be here in bounds\n");
2532 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
2538 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
2544 int i,
j,k,pl,pliter;
2545 int locali1,globali1;
2546 int locali2,globali2;
2553 FTYPE prescale[NPR];
2569 if(!finite(
MACP0A1(prim,i,j,k,pl))){
2571 dualfprintf(
fail_file,
"whichdir=%d ispstag=%d trigger=%d :: BC didn't set properly: #1: i=%d j=%d k=%d ti=%d tj=%d tk=%d pl=%d\n",whichdir,ispstag,trigger,i,j,k,
startpos[1]+i,
startpos[2]+j,
startpos[3]+k,pl);
2581 if(!finite(
MACP0A1(prim,i,j,k,pl))){
2583 dualfprintf(
fail_file,
"whichdir=%d ispstag=%d trigger=%d :: BC didn't set properly: #1: i=%d j=%d k=%d ti=%d tj=%d tk=%d pl=%d\n",whichdir,ispstag,trigger,i,j,k,
startpos[1]+i,
startpos[2]+j,
startpos[3]+k,pl);
2606 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
2612 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
2620 int iloopstart,iloopend,iloopstep;
2623 FTYPE Bd3,Bd3ri2,Bd3ri3;
2625 FTYPE Bu1,Bu2,gcon03,gcon13,gcon23,gcon33;
2626 FTYPE gcov01,gcov02,gcov11,gcov12,gcov21,gcov22,gcov03,gcov13,gcov23;
2631 FTYPE primtemp[NPR];
2640 FTYPE Dqp,Dqm,Dqc,dqMdot;
2645 FTYPE dqgdetpl[NPR];
2646 FTYPE dqorthopl[NPR];
2647 FTYPE dqlogdensity[NPR];
2652 FTYPE xdqfrac,ydqfrac;
2653 FTYPE ftemp2,linearvalue,expvalue;
2657 struct of_geom geomdontuse[NPR];
2659 struct of_geom rgeomdontuse[NPR];
2660 struct of_geom *ptrrgeom[NPR];
2661 struct of_geom ri2geomdontuse[NPR];
2662 struct of_geom *ptrri2geom[NPR];
2663 struct of_geom ri3geomdontuse[NPR];
2664 struct of_geom *ptrri3geom[NPR];
2669 ptrgeom[pl]=&(geomdontuse[pl]);
2670 ptrrgeom[pl]=&(rgeomdontuse[pl]);
2671 ptrri2geom[pl]=&(ri2geomdontuse[pl]);
2672 ptrri3geom[pl]=&(ri3geomdontuse[pl]);
2691 else if(boundary==
X1UP){
2702 dualfprintf(
fail_file,
"extrapfunc not setup for this boundary = %d\n",boundary);
2715 get_geometry(ri2, rj, rk, dirprim[pl], ptrri2geom[pl]);
2716 get_geometry(ri3, rj, rk, dirprim[pl], ptrri3geom[pl]);
2718 bl_coord_ijk_2(ri,rj,rk,dirprim[pl],Xr[pl],Vr[pl]);
2729 Bd3ri2=0.0;
SLOOPA(jj) Bd3ri2 +=
MACP0A1(prim,ri2,rj,rk,
B1+jj-1)*(ptrri2geom[
B1+jj-1]->gcov[
GIND(3,jj)]);
2730 Bd3ri3=0.0;
SLOOPA(jj) Bd3ri3 +=
MACP0A1(prim,ri3,rj,rk,
B1+jj-1)*(ptrri3geom[
B1+jj-1]->gcov[
GIND(3,jj)]);
2749 Mdot = (ptrrgeom[
U1]->gdet)*
MACP0A1(prim,ri,rj,rk,
RHO)*uconref[1];
2752 ucon_calc(
MAC(prim,ri2,rj,rk),ptrri2geom[U1],uconrefri2,othersrefri2);
2753 Mdotri2 = (ptrri2geom[
U1]->gdet)*
MACP0A1(prim,ri2,rj,rk,
RHO)*uconrefri2[1];
2756 ucon_calc(
MAC(prim,ri3,rj,rk),ptrri3geom[U1],uconrefri3,othersrefri3);
2757 Mdotri3 = (ptrri3geom[
U1]->gdet)*
MACP0A1(prim,ri3,rj,rk,
RHO)*uconrefri3[1];
2760 Dqp=Mdotri3-Mdotri2;
2768 ucon_calc(&
MACP0A1(prim,ri,rj,rk,URAD1-U1),ptrrgeom[URAD1],uradconref,othersradref);
2769 ucon_calc(&
MACP0A1(prim,ri2,rj,rk,URAD1-U1),ptrri2geom[URAD1],uradconrefri2,othersradrefri2);
2770 ucon_calc(&
MACP0A1(prim,ri3,rj,rk,URAD1-U1),ptrri3geom[URAD1],uradconrefri3,othersradrefri3);
2813 Dqp=(ptrri3geom[pl]->gdet)*
MACP0A1(prim,ri3,rj,rk,pl)-(ptrri2geom[pl]->gdet)*
MACP0A1(prim,ri2,rj,rk,pl);
2814 Dqm=(ptrri2geom[pl]->gdet)*
MACP0A1(prim,ri2,rj,rk,pl)-(ptrrgeom[pl]->gdet)*
MACP0A1(prim,ri,rj,rk,pl);
2815 Dqc=(ptrri3geom[pl]->gdet)*
MACP0A1(prim,ri3,rj,rk,pl)-(ptrrgeom[pl]->gdet)*
MACP0A1(prim,ri,rj,rk,pl);
2829 if(pl>=
U1 || pl<=
U3) ddir=pl-
U1+1;
2830 if(pl>=URAD1 || pl<=URAD3) ddir=pl-URAD1+1;
2831 if(pl>=
B1 || pl<=
B3) ddir=pl-
B1+1;
2832 Dqp=sqrt(fabs(ptrri3geom[pl]->
gcov[
GIND(ddir,ddir)]))*
MACP0A1(prim,ri3,rj,rk,pl)-sqrt(fabs(ptrri2geom[pl]->
gcov[
GIND(ddir,ddir)]))*
MACP0A1(prim,ri2,rj,rk,pl);
2833 Dqm=sqrt(fabs(ptrri2geom[pl]->
gcov[
GIND(ddir,ddir)]))*
MACP0A1(prim,ri2,rj,rk,pl)-sqrt(fabs(ptrrgeom[pl]->
gcov[
GIND(ddir,ddir)]))*
MACP0A1(prim,ri,rj,rk,pl);
2834 Dqc=sqrt(fabs(ptrri3geom[pl]->
gcov[
GIND(ddir,ddir)]))*
MACP0A1(prim,ri3,rj,rk,pl)-sqrt(fabs(ptrrgeom[pl]->
gcov[
GIND(ddir,ddir)]))*
MACP0A1(prim,ri,rj,rk,pl);
2848 Dqp=uconrefri3[jj]-uconrefri2[jj];
2849 Dqm=uconrefri2[jj]-uconref[jj];
2850 Dqc=uconrefri3[jj]-uconref[jj];
2857 xfrac = fabs(uconref[
TT]-uconrefri2[
TT])/uconref[
TT];
2858 #define XFRAC1 (0.1) // 0.1 was chosen from experience with BH+torus and what the u^t value was doing at the pole near the horizon -- JCM found once d(u^t)>0.1 then started to grow unbound. Certainly related to resolution so assume resolving of horizon is generally similar to the 64^2 case with Rout=40 and exp grid
2859 #define YFRAC1 (0.0)
2860 #define XFRAC2 (0.2)
2861 #define YFRAC2 (1.0)
2862 #define YFRAC12(frac) (YFRAC1 + (YFRAC2-YFRAC1)/(XFRAC2-XFRAC1)*(frac-XFRAC1))
2869 #define XDQFRAC1 (0.1)
2870 #define YDQFRAC1 (0.0)
2871 #define XDQFRAC2 (0.2)
2872 #define YDQFRAC2 (1.0)
2873 #define YDQFRAC12(frac) (YDQFRAC1 + (YDQFRAC2-YDQFRAC1)/(XDQFRAC2-XDQFRAC1)*(frac-XDQFRAC1))
2876 xdqfrac = fabs(uconref[TT]-uconrefri2[TT])/uconref[
TT];
2885 dqucon[jj] = dqucon[jj]*(1.0-ydqfrac);
2886 dq[
UU+jj] = dq[
UU+jj]*(1.0-ydqfrac);
2887 dqlogdensity[
UU+jj] = dqlogdensity[
UU+jj]*(1.0-ydqfrac);
2888 dqMdot = dqMdot*(1.0-ydqfrac);
2889 dqgdetpl[
UU+jj] = dqgdetpl[
UU+jj]*(1.0-ydqfrac);
2929 bl_coord_ijk_2(i,j,k,dirprim[pl],X[pl],V[pl]);
2957 #define POWEREXTRAFRAC (0.1)
2958 #define POWERRATIO (1.0+POWEREXTRAFRAC) // ratio of densities not allow to go bigger than this when using dqlogdensity
2965 if(pl!=
B1 && pl!=
B2 && pl!=
B3 && pl!=URAD1 && pl!=URAD2 && pl!=URAD3 &&
SCALARPL(pl)==0){
2968 ftemp = exp(-signdq*dqlogdensity[pl]) -
POWERRATIO;
2971 if(ftemp>=0.0) mydqlog=log(
POWERRATIO)/(-signdq);
2972 else mydqlog = dqlogdensity[pl];
2976 MACP0A1(prim,i,j,k,pl) = exp(log(
SMALL+fabs(
MACP0A1(prim,ri,rj,rk,pl))) + mydqlog*(i-ri));
3023 myMdot = Mdot + dqMdot*(i-ri);
3024 ucon[1] = myMdot/((ptrgeom[
U1]->gdet)*
MACP0A1(prim,i,j,k,
RHO));
3027 ucon[2] = uconref[2] + dqucon[2]*(i-ri);
3029 ucon[3] = uconref[3] + dqucon[3]*(i-ri);
3035 primtemp[
U1] = ucon[1];
3036 primtemp[
U2] = ucon[2];
3037 primtemp[
U3] = ucon[3];
3054 primtemp[
U3] = uconref[3]/uconref[
TT]*ucontouse[
TT];
3065 ucon2pr(
WHICHVEL,ucontouse,ptrgeom[U1],
MAC(prim,i,j,k));
3078 #define DO_CONSERVE_D_INBOUNDS 0 // CHANGINGMARK
3082 #if(DO_CONSERVE_D_INBOUNDS)
3090 if(pl==U1 || pl==
U2 || pl==
U3 || pl==URAD1 || pl==URAD2 || pl==URAD3){
3093 ftemp =
MACP0A1(prim,ri,rj,rk,pl)*(1.0-yfrac) +
MACP0A1(prim,ri2,rj,rk,pl)*yfrac;
3095 MACP0A1(prim,i,j,k,pl) = (ftemp + dq[pl]*(i-ri));
3108 MACP0A1(prim,ri,rj,rk,pl) =
MACP0A1(prim,ri,rj,rk,pl)*(1.0-yfrac) + (
MACP0A1(prim,ri2,rj,rk,pl) + dq[pl]*(ri-(ri2)))*yfrac;
3113 #if(DO_CONSERVE_D_INBOUNDS)
3115 ucon_calc(
MAC(prim,ri,rj,rk),ptrrgeom[U1],uconrefnew, othersrefnew);
3126 if(pl==U1 || pl==
U2 || pl==
U3 || pl==URAD1 || pl==URAD2 || pl==URAD3){
3129 ftemp =
MACP0A1(prim,ri,rj,rk,pl)*(1.0-yfrac) +
MACP0A1(prim,ri3,rj,rk,pl)*yfrac;
3131 MACP0A1(prim,i,j,k,pl) = ftemp + dq[pl]*(i-ri);
3137 // linear extrap for velocities
3139 if(pl==U1 || pl==
U2 || pl==
U3 || pl==URAD1 || pl==URAD2 || pl==URAD3){
3141 ftemp =
MACP0A1(prim,ri,rj,rk,pl);
3143 MACP0A1(prim,i,j,k,pl) = ftemp + dq[pl]*(i-ri);
3148 if(pl==URAD1 || pl==URAD2 || pl==URAD3){
3149 ftemp =
MACP0A1(prim,ri,rj,rk,pl);
3151 MACP0A1(prim,i,j,k,pl) = ftemp;
3160 if(boundary==
X1DN && (uconref[1]>0.0 || ucon[1]>0.0))
MACP0A1(prim,i,j,k,U1)=0.0;
3161 if(boundary==
X1UP && (uconref[1]<0.0 || ucon[1]<0.0))
MACP0A1(prim,i,j,k,U1)=0.0;
3162 if(boundary==
X2DN && (uconref[2]>0.0 || ucon[2]>0.0))
MACP0A1(prim,i,j,k,
U2)=0.0;
3163 if(boundary==
X2UP && (uconref[2]<0.0 || ucon[2]<0.0))
MACP0A1(prim,i,j,k,
U2)=0.0;
3164 if(boundary==
X3DN && (uconref[3]>0.0 || ucon[3]>0.0))
MACP0A1(prim,i,j,k,
U3)=0.0;
3165 if(boundary==
X3UP && (uconref[3]<0.0 || ucon[3]<0.0))
MACP0A1(prim,i,j,k,
U3)=0.0;
3168 ucon_calc(&
MACP0A1(prim,ri,rj,rk,URAD1-U1),ptrgeom[URAD1],uradcon,othersrad);
3169 if(boundary==
X1DN && (uradconref[1]>0.0 || uradcon[1]>0.0))
MACP0A1(prim,i,j,k,URAD1)=0.0;
3170 if(boundary==
X1UP && (uradconref[1]<0.0 || uradcon[1]<0.0))
MACP0A1(prim,i,j,k,URAD1)=0.0;
3171 if(boundary==
X2DN && (uradconref[2]>0.0 || uradcon[2]>0.0))
MACP0A1(prim,i,j,k,URAD2)=0.0;
3172 if(boundary==
X2UP && (uradconref[2]<0.0 || uradcon[2]<0.0))
MACP0A1(prim,i,j,k,URAD2)=0.0;
3173 if(boundary==
X3DN && (uradconref[3]>0.0 || uradcon[3]>0.0))
MACP0A1(prim,i,j,k,URAD3)=0.0;
3174 if(boundary==
X3UP && (uradconref[3]<0.0 || uradcon[3]<0.0))
MACP0A1(prim,i,j,k,URAD3)=0.0;
3202 for(pl=
B1;pl<=
B2;pl++){
3205 MACP0A1(prim,i,j,k,pl) = (
MACP0A1(prim,ri,rj,rk,pl)*sqrt(fabs(ptrrgeom[pl]->
gcov[
GIND(ddirl,ddirl)])) + 0.0*dqorthopl[pl]*(i-ri))/sqrt(fabs(ptrgeom[pl]->
gcov[
GIND(ddirl,ddirl)]));
3213 for(pl=
B1;pl<=
B2;pl++){
3215 MACP0A1(prim,i,j,k,pl) = (
MACP0A1(prim,ri,rj,rk,pl)*(ptrrgeom[pl]->gdet) + dqgdetpl[pl]*(i-ri))*igdetnosing;
3243 gcon03=ptrgeom[
B3]->gcon[
GIND(0,3)];
3244 gcon13=ptrgeom[
B3]->gcon[
GIND(1,3)];
3245 gcon23=ptrgeom[
B3]->gcon[
GIND(2,3)];
3246 gcon33=ptrgeom[
B3]->gcon[
GIND(3,3)];
3248 gcov01=ptrgeom[
B3]->gcov[
GIND(0,1)];
3249 gcov02=ptrgeom[
B3]->gcov[
GIND(0,2)];
3250 gcov11=ptrgeom[
B3]->gcov[
GIND(1,1)];
3251 gcov12=gcov21=ptrgeom[
B3]->gcov[
GIND(1,2)];
3252 gcov22=ptrgeom[
B3]->gcov[
GIND(2,2)];
3253 gcov03=ptrgeom[
B3]->gcov[
GIND(0,3)];
3254 gcov13=ptrgeom[
B3]->gcov[
GIND(1,3)];
3255 gcov23=ptrgeom[
B3]->gcov[
GIND(2,3)];
3258 myBd3=Bd3 + dqBd3*(i-ri);
3259 ftemp=(1.0 - gcon03*gcov03 - gcon13*gcov13 - gcon23*gcov23);
3260 igdetnosing=
sign(ftemp)/(fabs(ftemp)+
SMALL);
3261 pl=
B3;
MACP0A1(prim,i,j,k,pl) = (myBd3*gcon33 + Bu1*gcon03*gcov01 + Bu2*gcon03*gcov02 + Bu1*gcon13*gcov11 + Bu2*gcon13*gcov12 + Bu1*gcon23*gcov21 + Bu2*gcon23*gcov22)*igdetnosing;
3271 MACP0A1(prim,i,j,k,pl) = (
MACP0A1(prim,ri,rj,rk,pl)*(ptrrgeom[pl]->gdet) + dqgdetpl[pl]*(i-ri))*igdetnosing;
3291 #define UTHETAPOLEDEATH ((PRIMTOINTERP_3VELREL_GAMMAREL_DXDXP==VARTOINTERP)?(1):(0)) //only interpolate u^\theta instead of u^2 here when doing the same in interp
3293 #if(0 == UTHETAPOLEDEATH)
3294 # define MACP0A1mod(prim,ri,rj,rk,pl) MACP0A1(prim,ri,rj,rk,pl)
3305 if ( dir < U2 || dir >
B3 ) {
3306 dualfprintf(
fail_file,
"dir cannot be < U2 or > B3 in MACP0A1mod()\n" );
3311 dir = 1 + (pl-
U1)%3;
3315 return(
MACP0A1(prim,ri,rj,rk,pl) * dxdxp[dir][dir] );
3331 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
3337 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
3357 int rjtest,rjstag0,deathstagjs0,deathstagje0,rjstagtest,deathstagjstest,deathstagjetest;
3360 int deathjs0,deathje0;
3361 int deathjs,deathje;
3362 int deathstagjs,deathstagje;
3363 int gammadeathjs,gammadeathje;
3368 int lowrho,lowuu,highu1;
3369 int deathjstest,deathjetest;
3370 FTYPE gamma,gammarj0,gammarjtest;
3371 FTYPE qsq,qsqrj0,qsqrjtest;
3377 FTYPE gammavaluelimit,gammaradvaluelimit;
3378 int doavginradius[NPR];
3380 struct of_geom geomdontuse[NPR];
3382 struct of_geom rgeomdontuse[NPR];
3383 struct of_geom *ptrrgeom[NPR];
3387 struct of_geom geompoledontuse;
3388 struct of_geom *ptrgeompole=&geompoledontuse;
3392 int testplexist[MAXNPR]={0};
3396 if(testplexist[pl]!=1){
3418 ptrgeom[pl]=&(geomdontuse[pl]);
3419 ptrrgeom[pl]=&(rgeomdontuse[pl]);
3444 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3456 int poledeathreal,polegammadeathreal;
3462 #define RADIUSOUTERDEATHMORE (BIG) //(OUTERDEATHRADIUS) // i.e. inactive
3466 poledeathreal=
N2BND;
3467 polegammadeathreal=
N2BND;
3477 rj0 = poledeathreal;
3482 deathjs0 = 0-poledeathreal;
3483 deathje0 = 0+poledeathreal-1;
3485 deathjstest = deathjs0-DEATHEXPANDAMOUNT;
3486 deathjetest = deathje0+DEATHEXPANDAMOUNT;
3493 deathstagjs0 = 0-poledeathreal+1;
3494 deathstagje0 = 0+poledeathreal-1;
3496 rjstagtest = rjtest;
3497 deathstagjstest = deathstagjs0-DEATHEXPANDAMOUNT;
3498 deathstagjetest = deathstagje0+DEATHEXPANDAMOUNT;
3503 gammadeathjs=0-polegammadeathreal;
3504 gammadeathje=0+polegammadeathreal-1;
3519 else if(whichx2==
X2UP){
3523 poledeathreal=
N2BND;
3524 polegammadeathreal=
N2BND;
3531 rj0=
N2-1-poledeathreal;
3536 deathjs0 =
N2-1+1-poledeathreal;
3537 deathje0 =
N2-1+poledeathreal;
3539 deathjstest = deathjs0-DEATHEXPANDAMOUNT;
3540 deathjetest = deathje0+DEATHEXPANDAMOUNT;
3545 if(dirprim[
B2]==
FACE2) rjstag0=
N2-poledeathreal;
3546 else if(dirprim[
B2]==
CENT) rjstag0=rj0;
3547 deathstagjs0 =
N2+1-poledeathreal;
3548 deathstagje0 =
N2-1+poledeathreal;
3550 rjstagtest = rjtest;
3551 deathstagjstest = deathstagjs0-DEATHEXPANDAMOUNT;
3552 deathstagjetest = deathstagje0+DEATHEXPANDAMOUNT;
3558 gammadeathjs =
N2-1+1-polegammadeathreal;
3559 gammadeathje =
N2-1+polegammadeathreal;
3576 FTYPE Rhorref=rhor_calc(0);
3589 int dontskippoledeath;
3590 dontskippoledeath=0;
3595 dontskippoledeath=1;
3683 gamma_calc(
MAC(prim,ri,rj0,rk), ptrrgeom[U1],&gammarj0,&qsqrj0);
3685 get_geometry(ri, rjtest, rk, dirprim[U1], ptrrgeom[U1]);
3686 gamma_calc(
MAC(prim,ri,rjtest,rk), ptrrgeom[U1],&gammarjtest,&qsqrjtest);
3690 for (j = deathjs0; j <= deathje0; j++) {
3697 if(fabs(
MACP0A1(prim,i,j,k,pl))< fabs(
MACP0A1(prim,ri,rj0,rk,pl))/
POLEDENSITYDROPFACTOR || fabs(
MACP0A1(prim,i,j,k,pl))< fabs(
MACP0A1(prim,ri,rjtest,rk,pl))/
POLEDENSITYDROPFACTOR){
3701 if(fabs(
MACP0A1(prim,i,j,k,pl))< fabs(
MACP0A1(prim,ri,rj0,rk,pl)/
POLEDENSITYDROPFACTOR) || fabs(
MACP0A1(prim,i,j,k,pl))< fabs(
MACP0A1(prim,ri,rjtest,rk,pl))/
POLEDENSITYDROPFACTOR){
3708 gamma_calc(
MAC(prim,i,j,k), ptrrgeom[U1],&gamma,&qsq);
3724 if(lowuu || lowrho || highu1){
3726 deathjs=deathjstest;
3727 deathje=deathjetest;
3730 deathstagjs=deathstagjstest;
3731 deathstagje=deathstagjetest;
3740 deathstagjs=deathstagjs0;
3741 deathstagje=deathstagje0;
3746 deathjs=deathjstest;
3747 deathje=deathjetest;
3750 deathstagjs=deathstagjstest;
3751 deathstagje=deathstagjetest;
3766 get_geometry(ri, rjstag, rk, dirprim[pl], ptrrgeom[pl]);
3782 #define DEATHLOOP2(whichx2) for (j = (whichx2==X2UP ? MIN(deathstagjs,deathjs) : MAX(deathstagje,deathje)) ; (whichx2==X2UP ? (j <= MAX(deathstagje,deathje)) : (j >= MIN(deathstagjs,deathjs)) ) ; (whichx2==X2UP ? j++ : j--) )
3813 bl_coord_ijk_2(i,j,k,dirprim[pl],X[pl],V[pl]);
3817 doavginradius[pl]=1;
3819 else doavginradius[pl]=0;
3833 if(j>=deathjs && j<=deathje){
3840 if(!(pl==
U1 || pl==
U3 || pl==URAD1 || pl==URAD3))
continue;
3843 myvalue=
MACP0A1(prim,i,j,k,pl);
3845 if(doavginradius[pl]){
3846 avgresult=
THIRD*(
MACP0A1(prim,rim1,rj,rk,pl)+
MACP0A1(prim,ri,rj,rk,pl)+
MACP0A1(prim,rip1,rj,rk,pl));
3849 avgresult=
MACP0A1(prim,ri,rj,rk,pl);
3852 if(1||fabs(avgresult)<myvalue){
3853 MACP0A1(prim,i,j,k,pl) = avgresult;
3865 if(!(pl==
RHO || pl==
UU || pl==ENTROPY || pl==URAD0))
continue;
3868 myvalue=
MACP0A1(prim,i,j,k,pl);
3870 if(doavginradius[pl]){
3871 avgresult=
THIRD*(
MACP0A1(prim,rim1,rj,rk,pl)+
MACP0A1(prim,ri,rj,rk,pl)+
MACP0A1(prim,rip1,rj,rk,pl));
3874 avgresult=
MACP0A1(prim,ri,rj,rk,pl);
3877 if(1||fabs(avgresult)>myvalue){
3878 MACP0A1(prim,i,j,k,pl) = avgresult;
3893 if (
bsq_calc(
MAC(prim,ri,rj,rk), ptrrgeom[pl], &bsqr) >= 1)
FAILSTATEMENT(
"bounds.tools.c:poledeath()",
"bsq_calc()", 1);
3894 sigmar=bsqr/(2.0*fabs(
MACP0A1(prim,ri,rj,rk,pl)));
3898 if (
bsq_calc(
MAC(prim,i,j,k), ptrgeom[pl], &bsq) >= 1)
FAILSTATEMENT(
"bounds.tools.c:poledeath()",
"bsq_calc()", 2);
3899 sigma=bsq/(2.0*fabs(
MACP0A1(prim,i,j,k,pl)));
3902 MACP0A1(prim,i,j,k,pl) = fabs(
MACP0A1(prim,i,j,k,pl)*(sigma/sigmar));
3921 if(j>=deathjs && j<=deathje){
3933 dirprim[
B2]==
FACE2 && j>=deathstagjs && j<=deathstagje ||
3934 dirprim[
B2]==
CENT && j>=deathjs && j<=deathje
3940 #if(POLEINTERPTYPE==0)
3946 #elif(POLEINTERPTYPE==1 || POLEINTERPTYPE==2)
3950 if(doavginradius[pl]) ftemp=
THIRD*(
MACP0A1(prim,rim1,rjstag,rk,pl) +
MACP0A1(prim,ri,rjstag,rk,pl) +
MACP0A1(prim,rip1,rjstag,rk,pl));
3951 else ftemp=
MACP0A1(prim,ri,rjstag,rk,pl);
3952 MACP0A1(prim,i,j,k,pl) = ftemp + (X[pl][2]-Xr[pl][2])*(ftemp-0.0)/(Xr[pl][2]-X0[2]);
3955 #elif(POLEINTERPTYPE==3)
3974 if(j>=deathjs && j<=deathje){
3982 #if(POLEINTERPTYPE==0 || POLEINTERPTYPE==1)
3984 if(doavginradius[pl]) ftemp=
THIRD*(
MACP0A1(prim,rim1,rj,rk,pl) +
MACP0A1(prim,ri,rj,rk,pl) +
MACP0A1(prim,rip1,rj,rk,pl));
3985 else ftemp=
MACP0A1(prim,ri,rj,rk,pl);
3989 #elif(POLEINTERPTYPE==2)
3992 if(doavginradius[pl]) ftemp=
THIRD*(
MACP0A1(prim,rim1,rj,rk,pl) +
MACP0A1(prim,ri,rj,rk,pl) +
MACP0A1(prim,rip1,rj,rk,pl));
3993 else ftemp=
MACP0A1(prim,ri,rj,rk,pl);
3997 #elif(POLEINTERPTYPE==3)
3999 if(doavginradius[pl]) ftemp=
THIRD*(
MACP0A1(prim,rim1,rj,rk,pl) +
MACP0A1(prim,ri,rj,rk,pl) +
MACP0A1(prim,rip1,rj,rk,pl));
4000 else ftemp=
MACP0A1(prim,ri,rj,rk,pl);
4019 if((pl==
RHO || pl==
UU || pl==
U1 || pl==
U2 || pl==
U3 || pl==ENTROPY || pl==
B1 || pl==
B2 || pl==
B3 || pl==URAD0 || pl==URAD1 || pl==URAD2 || pl==URAD3))
continue;
4021 if(doavginradius[pl]) ftemp=
THIRD*(
MACP0A1(prim,rim1,rj,rk,pl) +
MACP0A1(prim,ri,rj,rk,pl) +
MACP0A1(prim,rip1,rj,rk,pl));
4022 else ftemp=
MACP0A1(prim,ri,rj,rk,pl);
4032 if(madechange && candodiag){
4045 fixupptrgeom=ptrgeom[pl];
4049 int doingmhdfixup=1;
4050 diag_fixup(modcons,prdiag,pr,ucons,fixupptrgeom,finalstep,doingmhdfixup,
COUNTBOUND1);
4092 bl_coord_ijk_2(i,j,k,dirprim[pl],X[pl],V[pl]);
4096 doavginradius[pl]=1;
4098 else doavginradius[pl]=0;
4106 if(j>=deathjs && j<=deathje){
4110 for(iteru2=0;iteru2<=1;iteru2++){
4120 if(URAD0<0)
continue;
4138 if(j>=jother) signD=+1.0;
4143 FTYPE rhovjhere,rhovjother,rhovDiff;
4147 rhovDiff=signD*(rhovjhere-rhovjother);
4155 FTYPE dU2 =-signD*rhovDiff*0.5/
MACP0A1(prim,i,j,k,plrho);
4159 if( (fabs(U2jhere)>fabs(dU2))&&(fabs(U2jother)>fabs(dU2)) ){
4161 MACP0A1(prim,i,j,k,pl) += dU2;
4162 MACP0A1(prim,i,jother,k,pl) -= dU2;
4167 if(fabs(U2jhere)<fabs(U2jother)){
4169 MACP0A1(prim,i,j,k,pl) *= -1.0;
4175 MACP0A1(prim,i,jother,k,pl) *= -1.0;
4197 #if(POLEINTERPTYPE==0)
4203 #elif(POLEINTERPTYPE==1 || POLEINTERPTYPE==2)
4206 if(doavginradius[pl]) ftemp=
THIRD*(
MACP0A1mod(prim,rim1,rj,rk,pl) +
MACP0A1mod(prim,ri,rj,rk,pl) +
MACP0A1mod(prim,rip1,rj,rk,pl));
4208 MACP0A1(prim,i,j,k,pl) = ftemp + (X[pl][2]-Xr[pl][2])*(ftemp-0.0)/(Xr[pl][2]-X0[2]);
4211 #elif(POLEINTERPTYPE==3)
4221 if(whichx2==
X2DN && ftemp>0.0){
4225 if(doavginradius[pl]){
4226 for(jj=0;jj<=rj+DEATHEXPANDAMOUNT;jj++) ftemp=
MIN(ftemp,
MACP0A1mod(prim,rip1,jj,rk,pl));
4227 for(jj=0;jj<=rj+DEATHEXPANDAMOUNT;jj++) ftemp=
MIN(ftemp,
MACP0A1mod(prim,rim1,jj,rk,pl));
4233 MACP0A1(prim,i,j,k,pl) = ftemp + (X[pl][2]-Xr[pl][2])*(ftemp-0.0)/(Xr[pl][2]-X0[2]);
4236 else if(whichx2==
X2UP && ftemp<0.0){
4239 if(doavginradius[pl]){
4240 for(jj=
N2-1;jj>=rj-DEATHEXPANDAMOUNT;jj--) ftemp=
MAX(ftemp,
MACP0A1mod(prim,rip1,jj,rk,pl));
4241 for(jj=
N2-1;jj>=rj-DEATHEXPANDAMOUNT;jj--) ftemp=
MAX(ftemp,
MACP0A1mod(prim,rim1,jj,rk,pl));
4247 MACP0A1(prim,i,j,k,pl) = ftemp + (X[pl][2]-Xr[pl][2])*(ftemp-0.0)/(Xr[pl][2]-X0[2]);
4252 if(doavginradius[pl]) ftemp=
THIRD*(
MACP0A1mod(prim,rim1,rj,rk,pl) +
MACP0A1mod(prim,ri,rj,rk,pl) +
MACP0A1mod(prim,rip1,rj,rk,pl));
4255 MACP0A1(prim,i,j,k,pl) = ftemp + (X[pl][2]-Xr[pl][2])*(ftemp-0.0)/(Xr[pl][2]-X0[2]);
4259 #endif // endif POLEINTERPTYPE==3
4263 #if( UTHETAPOLEDEATH )
4267 MACP0A1(prim,i,j,k,pl) /= dxdxp[1 + (pl-
U1)%3][1 + (pl-
U1)%3];
4273 MACP0A1(prim,i,j,k,pl) /= dxdxp[1 + (pl-URAD1)%3][1 + (pl-URAD1)%3];
4284 if(madechange&&candodiag){
4297 fixupptrgeom=ptrgeom[pl];
4302 int doingmhdfixup=1;
4303 diag_fixup(modcons,prdiag,pr,ucons,fixupptrgeom,finalstep,doingmhdfixup,
COUNTBOUND1);
4331 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
4345 int poledeathreal,polegammadeathreal;
4352 poledeathreal=
N2BND;
4353 polegammadeathreal=
N2BND;
4363 rj0 = poledeathreal;
4368 deathjs0 = 0-poledeathreal;
4369 deathje0 = 0+poledeathreal-1;
4371 deathjstest = deathjs0-DEATHEXPANDAMOUNT;
4372 deathjetest = deathje0+DEATHEXPANDAMOUNT;
4379 deathstagjs0 = 0-poledeathreal+1;
4380 deathstagje0 = 0+poledeathreal-1;
4382 rjstagtest = rjtest;
4383 deathstagjstest = deathstagjs0-DEATHEXPANDAMOUNT;
4384 deathstagjetest = deathstagje0+DEATHEXPANDAMOUNT;
4389 gammadeathjs=0-polegammadeathreal;
4390 gammadeathje=0+polegammadeathreal-1;
4406 else if(whichx2==
X2UP){
4411 poledeathreal=
N2BND;
4412 polegammadeathreal=
N2BND;
4419 rj0=
N2-1-poledeathreal;
4424 deathjs0 =
N2-1+1-poledeathreal;
4425 deathje0 =
N2-1+poledeathreal;
4427 deathjstest = deathjs0-DEATHEXPANDAMOUNT;
4428 deathjetest = deathje0+DEATHEXPANDAMOUNT;
4433 if(dirprim[
B2]==
FACE2) rjstag0=
N2-poledeathreal;
4434 else if(dirprim[
B2]==
CENT) rjstag0=rj0;
4435 deathstagjs0 =
N2+1-poledeathreal;
4436 deathstagje0 =
N2-1+poledeathreal;
4438 rjstagtest = rjtest;
4439 deathstagjstest = deathstagjs0-DEATHEXPANDAMOUNT;
4440 deathstagjetest = deathstagje0+DEATHEXPANDAMOUNT;
4446 gammadeathjs =
N2-1+1-polegammadeathreal;
4447 gammadeathje =
N2-1+polegammadeathreal;
4466 FTYPE Rhorref=rhor_calc(0);
4475 for (j = gammadeathjs; j <= gammadeathje; j++) {
4494 dualfprintf(
fail_file,
"BNDNAN: ispstag=%d t=%21.15g nstep=%ld steppart=%d :: i=%d j=%d k=%d pl2=%d prim=%21.15g\n",ispstag,
t,
nstep,
steppart,i,j,k,pl2,
MACP0A1(prim,i,j,k,pl2));
4515 bl_coord_ijk_2(i,j,k,dirprim[pl],X[pl],V[pl]);
4536 else if(ucon[
RR]<0.0){
4544 if(1||ucon[
RR]>0.0){
4573 if(uconrad[
RR]>0.0){
4577 else if(uconrad[
RR]<0.0){
4584 if(1||uconrad[
RR]>0.0){
4609 limit_gamma(0,gammavaluelimit,gammaradvaluelimit,
MAC(prim,i,j,k),NULL,ptrgeom[pl],-1);
4628 if(madechange&&candodiag){
4642 fixupptrgeom=ptrgeom[pl2];
4647 int doingmhdfixup=1;
4648 diag_fixup(modcons,prdiag,pr,ucons,fixupptrgeom,finalstep,doingmhdfixup,
COUNTBOUND2);
4672 #define DEBUGPOLESMOOTH 0
4678 int boundstage,
int finalstep,
SFTYPE boundtime,
int whichdir,
int boundvartype,
int *dirprim,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
4684 int riin,
int riout,
int rjin,
int rjout,
int rkin,
int rkout,
4691 int j0, dj, stopj, rj;
4697 static FTYPE *fullpr;
4698 FTYPE cartavgpr[NPR], spcavgpr[NPR],spcpr[NPR];
4699 static int firsttime=1;
4714 if(ispstag)
return(0);
4717 dualfprintf(
fail_file,
"SUPERWARNING: Must enable special3dspc==1 for it to be useful when N3>1\n");
4721 dualfprintf(
fail_file,
"SUPERWARNING: Must enable dofull2pi==1 and have full 2\\pi doain when using polesmooth() for it to be useful when N3>1\n");
4725 #if(DEBUGPOLESMOOTH)
4743 #define MAPFULLPR(i,fakej,k,pl) ( (fakej)*(N3BND+N3+N3BND)*(N1BND+N1+N1BND)*NPR + ((k)+N3BND)*(N1BND+N1+N1BND)*NPR + ((i)+N1BND)*NPR + (pl) )
4745 #define N1TOTFULLPR (N1BND+N1+N1BND)
4746 #define N2TOTFULLPR (ncpux3>1&&USEMPI&&N3>1 ? 1 : 2) // need 2 j positions if single cpu
4747 #define N3TOTFULLPR (N3BND+totalsize[3]+N3BND)
4752 dualfprintf(
fail_file,
"Cannot allocate fullpr\n");
4766 if (whichx2==
X2DN) {
4776 rj =
N2-1-POLESMOOTH;
4792 static MPI_Request *srequest;
4793 static MPI_Request *rrequest;
4798 srequest=(MPI_Request *)malloc(
sizeof(MPI_Request)*
ncpux3);
4800 dualfprintf(
fail_file,
"Cannot allocate srequest\n");
4803 rrequest=(MPI_Request *)malloc(
sizeof(MPI_Request)*
ncpux3);
4805 dualfprintf(
fail_file,
"Cannot allocate rrequest\n");
4821 int count=
N3*
N1M*NPR;
4823 for(posk=0;posk<
ncpux3;posk++){
4827 int recvtag=TAGSTARTBOUNDMPIPOLESMOOTH + originmyid +
numprocs*mycpupos[3];
4830 MPI_Irecv(&fullpr[
MAPFULLPR(-
N1BND,fakej,posk*
N3,0)] , count , MPI_FTYPE , MPIid[originmyid] , recvtag , MPI_COMM_GRMHD , &rrequest[posk]);
4831 #if(DEBUGPOLESMOOTH)
4832 dualfprintf(
fail_file,
"MPI_Irecv: posk=%d : %d %d\n",posk,originmyid,recvtag);
4870 pr =
MAC(prim,i,rj,k);
4874 PBOUNDLOOP(pliter,pl)
if(pl<U1 || pl>
B3) spcpr[pl] = pr[pl];
4890 bl_coord_ijk_2(i,rj,k,
CENT,X,V);
4898 fullpr[
MAPFULLPR(i,fakej,
mycpupos[3]*
N3+k,
B1)] = -(r*sin(ph)*sin(th)*spcpr[
U3]) + cos(ph)*sin(th)*spcpr[
U1] + r*cos(ph)*cos(th)*spcpr[
U2];
4899 fullpr[
MAPFULLPR(i,fakej,
mycpupos[3]*
N3+k,
B2)] = r*cos(ph)*sin(th)*spcpr[
U3] + sin(ph)*sin(th)*spcpr[
U1] + r*cos(th)*sin(ph)*spcpr[
U2];
4902 #if(DEBUGPOLESMOOTH)
4903 PBOUNDLOOP(pliter,pl) dualfprintf(
fail_file,
"Got hereINITIAL: t=%g i=%d j=%d k=%d : pl=%d pr=%g spc=%g cart=%g i/dxdxp11=%g %g\n",
t,i,j,k,pl,pr[pl],spcpr[pl],fullpr[
MAPFULLPR(i,fakej,
mycpupos[3]*
N3+k,pl)],idxdxp[
RR][
RR],dxdxp[RR][RR]);
4920 MPI_Status mpistatus;
4933 int count=
N3*
N1M*NPR;
4936 for(posk=0;posk<
ncpux3;posk++){
4940 int recvtag=TAGSTARTBOUNDMPIPOLESMOOTH + originmyid +
numprocs*mycpupos[3];
4943 int sendtag=TAGSTARTBOUNDMPIPOLESMOOTH + myid +
numprocs*posk;
4945 if(MPIFLOWCONTROL==1){
4952 ¬hingsend,0,MPI_INT,
4956 ¬hingrecv,0,MPI_INT,
4960 MPI_COMM_GRMHD,MPI_STATUS_IGNORE);
4966 MPI_Isend(&fullpr[
MAPFULLPR(-
N1BND,fakej,mycpupos[3]*
N3,0)] , count , MPI_FTYPE , MPIid[destmyid] , sendtag , MPI_COMM_GRMHD , &srequest[posk]);
4967 #if(DEBUGPOLESMOOTH)
4968 dualfprintf(
fail_file,
"MPI_Isend: posk=%d : %d %d\n",posk,destmyid,sendtag);
4979 for(posk=0;posk<
ncpux3;posk++){
4981 MPI_Wait(&rrequest[posk],&mpistatus);
4982 #if(DEBUGPOLESMOOTH)
4983 dualfprintf(
fail_file,
"MPI_Wait: posk=%d\n",posk);
4990 for(posk=0;posk<
ncpux3;posk++){
4992 MPI_Wait(&srequest[posk],&mpistatus);
4993 #if(DEBUGPOLESMOOTH)
4994 dualfprintf(
fail_file,
"MPI_Wait: posk=%d\n",posk);
5004 #if(DEBUGPOLESMOOTH)
5021 #if(DEBUGPOLESMOOTH)
5022 dualfprintf(
fail_file,
"Got here t=%g i=%d\n",
t,i);
5042 for(k=0;k<totalsize[3];k++){
5044 #if(DEBUGPOLESMOOTH)
5045 dualfprintf(
fail_file,
"Got here t=%g i=%d k=%d\n",
t,i,k);
5052 PBOUNDLOOP(pliter,pl)
if(pl<U1 || pl>B3) cartavgpr[pl] += fullpr[
MAPFULLPR(i,fakej,k,pl)];
5070 cartavgpr[pl] /= (
FTYPE)totalsize[3];
5071 spcavgpr[pl] /= (
FTYPE)totalsize[3];
5079 cartavgpr[
U1] = 0.0;
5080 cartavgpr[
U2] = 0.0;
5083 #if(DEBUGPOLESMOOTH)
5084 PBOUNDLOOP(pliter,pl) dualfprintf(
fail_file,
"Got hereCARTSPC: t=%g i=%d pl=%d %g %g\n",
t,i,pl,cartavgpr[pl],spcavgpr[pl]);
5095 for (j=j0; j != stopj; j+=dj) {
5101 bl_coord_ijk_2(i,j,k,
CENT,X,V);
5107 PBOUNDLOOP(pliter,pl)
if(pl<U1 || pl>
B3) spcpr[pl] = spcavgpr[pl];
5111 spcpr[
U1] = cos(ph)*sin(th)*cartavgpr[
U1] + sin(ph)*sin(th)*cartavgpr[
U2] + cos(th)*cartavgpr[
U3];
5112 spcpr[
U2] = pow(r,-1)*(cos(ph)*cos(th)*cartavgpr[
U1] + cos(th)*sin(ph)*cartavgpr[
U2] - sin(th)*cartavgpr[
U3]);
5113 spcpr[
U3] = (1./sin(th))*pow(r,-1)*(-sin(ph)*cartavgpr[
U1] + cos(ph)*cartavgpr[
U2]);
5116 spcpr[
U3] += spcavgpr[
U3];
5129 pr =
MAC(prim,i,j,k);
5132 PBOUNDLOOP(pliter,pl)
if(pl<U1 || pl>B3) pr[pl] = spcpr[pl];
5162 #if(DEBUGPOLESMOOTH)
5163 PBOUNDLOOP(pliter,pl) dualfprintf(
fail_file,
"Got hereFINAL: t=%g i=%d j=%d k=%d : pl=%d pr=%g\n",
t,i,j,k,pl,pr[pl]);
5194 if(ADJUSTFLUXCT==0){
5218 if( i < -
MAXBND )
continue;
5222 if(dir==
X1UP &&
BCtype[dir]==FIXEDUSEPANALYTIC){
5235 #if(0) // below not currently used
5248 if( j < -
MAXBND )
continue;
5268 if( k < -
MAXBND )
continue;
5315 if( i < -
MAXBND )
continue;
5338 if( i < -
MAXBND )
continue;
5390 if( i < -
MAXBND )
continue;
5431 void check_spc_singularities_user(
void)
5438 int numlocs=
NPG,indloc,loc;
5443 if(
nstep==0) doprintout=1;
5461 for(indloc=0;indloc<numlocs;indloc++){
5463 if(indloc==0) loc=
FACE1;
5464 else if(indloc==1) loc=
CORN3;
5465 else if(indloc==2) loc=
CORN2;
5466 else if(indloc==3) loc=
CORNT;
5472 singfound+=(localgdet[0]!=0.0);
5473 #if(WHICHEOM!=WITHGDET)
5478 if(doprintout) dualfprintf(
fail_file,
"Detected singularity at i=%d j=%d k=%d loc=%d with gdet=%21.15g so resetting it to 0.0\n",i,j,k,loc,localgdet[0]);
5488 #if(WHICHEOM!=WITHGDET)
5524 for(poleloop=0;poleloop<=1;poleloop++){
5538 for(indloc=0;indloc<numlocs;indloc++){
5540 if(indloc==0) loc=
FACE2;
5541 else if(indloc==1) loc=
CORN3;
5542 else if(indloc==2) loc=
CORN1;
5543 else if(indloc==3) loc=
CORNT;
5549 singfound+=(fabs(localgdet[0])>0.0);
5550 #if(WHICHEOM!=WITHGDET)
5555 nonsingfound+=(fabs(localgdet[0])<100*
NUMEPSILON);
5556 #if(WHICHEOM!=WITHGDET)
5560 if(needzero && singfound){
5561 if(doprintout) dualfprintf(
fail_file,
"Detected singularity at i=%d j=%d k=%d loc=%d with gdet=%21.15g so resetting it to 0.0\n",i,j,k,loc,localgdet[0]);
5571 #if(WHICHEOM!=WITHGDET)
5576 else if(needzero==0 && nonsingfound){
5577 dualfprintf(
fail_file,
"Detected |\\detg|<=NUMEPSILON at i=%d j=%d k=%d loc=%d with gdet=%21.15g: Must be non-zero\n",i,j,k,loc,localgdet[0]);
5591 int debugspecial3dspc(
int which,
int whichdir,
int ispstag,
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR])
5593 int i,
j,k,pliter,pl;
5596 dualfprintf(
fail_file,
"nstep=%ld steppart=%d which=%d whichdir=%d ispstag=%d\n",
nstep,
steppart,which,whichdir,ispstag);
5597 for(pl=0;pl<=7;pl++){
5598 if(ispstag==1 && (pl<5 || pl>7))
continue;
5602 dualfprintf(
fail_file,
"%02s | i=\n ",
"j");
5604 if(ispstag==1 && (pl==5 && i==-
N1BND)) dualfprintf(
fail_file,
"%021s ",
"NA");
5605 else dualfprintf(
fail_file,
"%19s%02d ",
" ",i);
5611 if(ispstag==1 && (pl==5 && i==-
N1BND || pl==6 && j==-
N2BND)) dualfprintf(
fail_file,
"%21s ",
"NA");
5633 struct of_geom *ptrgeom=&geomdontuse;
5641 ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
5645 bl_coord_ijk_2(i,j,k,
CENT,X, V) ;
5660 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
5661 ufix[ENTROPY] = ufix[
UU];
5662 ufix[ENTROPY] =
MAX(0.0001,
MIN(ufix[ENTROPY],1.0));
5667 prfix[URAD0] =
MIN(
MIN(prfix[
RHO],prfix[URAD0]),prfix[
UU]);
5669 ufix[URAD0]=
MAX(-prfix[URAD0],ufix[URAD0]);
5674 limit_gamma(0,OUTERDEATHGAMMAMAX,OUTERDEATHGAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
5688 if(prfix[URAD1]<0.0){
5707 struct of_geom *ptrgeom=&geomdontuse;
5722 ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
5742 else if(AVOIDCS==2){
5747 bl_coord_ijk_2(i,j,k,
CENT,X, V) ;
5750 FTYPE Rhorlocal=1.0+sqrt(1.0-
a*
a);
5756 ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
5763 else if(AVOIDCS==3){
5771 FTYPE Rhorlocal=1.0+sqrt(1.0-
a*
a);
5773 FTYPE maxbsq=0.0,mybsq=0.0;
5782 prfix=&
MACP0A1(prim,fii,fjj,fkk,0);
5783 ufix=&GLOBALMACP0A1(unewglobal,fii,fjj,fkk,0);
5794 maxbsq=
MAX(maxbsq,bsq);
5796 if(ii==0 && jj==0 && kk==0) mybsq=bsq;
5802 if(mybsq*10.0<maxbsq){
5804 dualfprintf(
fail_file,
"i=%d j=%d k=%d mybsq=%g maxbsq=%g\n",i,j,k,mybsq,maxbsq);
5817 prfix=&
MACP0A1(prim,fii,fjj,fkk,0);
5818 ufix=&GLOBALMACP0A1(unewglobal,fii,fjj,fkk,0);
5846 #define OUTERDEATHDAMP 0
5848 if(DAMPGAMMA || OUTERDEATHDAMP || DAMPECOV){
5852 ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
5856 bl_coord_ijk_2(i,j,k,
CENT,X, V) ;
5865 gamma_calc(prfix,ptrgeom,&gamma,&qsq);
5866 FTYPE gammanew=((gamma-1.0)*0.1)+1.0;
5867 limit_gamma(0,gammanew,OUTERDEATHGAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
5890 FTYPE prfloor[NPR],prceiling[NPR];
5891 set_density_floors(ptrgeom,prfix,prfloor,prceiling);
5894 prfix[
UU]=prfloor[
UU];
5896 limit_gamma(0,OUTERDEATHGAMMAMAX,OUTERDEATHGAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
5907 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
5908 ufix[ENTROPY] = ufix[
UU];
5909 ufix[ENTROPY] =
MAX(0.0001,
MIN(ufix[ENTROPY],1.0));
5914 prfix[URAD0] =
MIN(
MIN(prfix[
RHO],prfix[URAD0]),prfix[
UU]);
5916 ufix[URAD0]=
MAX(-prfix[URAD0],ufix[URAD0]);
5921 limit_gamma(0,OUTERDEATHGAMMAMAX,OUTERDEATHGAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
5934 FTYPE localrbrr0=
MAX(100.0,localrbr/2.0);
5936 FTYPE switch0 = 0.5+1.0/M_PI*atan((V[1]-localrbr)/localrbrr0);
5937 FTYPE switch2 = 1.0-switch0;
5942 if(V[1]<100.0) switch2=1.0;
5958 if(prfix[URAD1]<0.0){
5960 prfix[URAD1]*=switch2;
5962 ufix[URAD1]*=switch2;
5973 FTYPE prfixtry[NPR],ufixtry[NPR];
5974 FTYPE prfixtryb[NPR],ufixtryb[NPR];
5980 ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
5982 PLOOP(pliter,pl) prfixtry[pl]=prfix[pl];
5984 PLOOP(pliter,pl) ufixtry[pl]=ufix[pl];
5987 PLOOP(pliter,pl) prfixtryb[pl]=prfix[pl];
5989 PLOOP(pliter,pl) ufixtryb[pl]=ufix[pl];
5994 bl_coord_ijk_2(i,j,k,CENT,X, V) ;
5999 #define BEGINDEATHTRANSITION (OUTERDEATHRADIUS*0.8)
6001 if(V[1]>OUTERDEATHRADIUS|| V[1]>BEGINDEATHTRANSITION){
6006 FTYPE localrbrdr=localrbr/3.0;
6008 #define cr(x) (exp(-1.0/(x)))
6009 #define tr(x) (cr(x)/(cr(x) + cr(1.0-(x))))
6010 #define trans(x,L,R) ((x)<=(L) ? 0.0 : ((x)>=(R) ? 1.0 : tr(((x)-(L))/((R)-(L)))) )
6012 FTYPE switch0,switch2;
6013 FTYPE switch0b,switch2b;
6017 switch0 = 0.5+1.0/M_PI*atan((V[1]-localrbr)/localrbrdr);
6018 switch2 = 1.0-switch0;
6020 switch0b = 0.5+1.0/M_PI*atan((V[1]-localrbr)/localrbrdr);
6021 switch2b = 1.0-switch0;
6024 rinner=
MAX(rbr-localrbrdr,BEGINDEATHTRANSITION);
6025 switch0 =
trans(V[1],rinner,rbr+localrbrdr);
6026 switch2 = 1.0-switch0;
6028 #define FACTORNEXTSWTICH (2.0)
6031 switch0b =
trans(V[1],rbr+localrbrdr,FACTORNEXTSWTICH*rbr+localrbrdr);
6032 switch2b = 1.0-switch0b;
6035 if(V[1]<BEGINDEATHTRANSITION){
6059 prfixtry[
UU]=
MIN(1.0*prfixtry[
RHO],prfixtry[
UU]);
6062 ufixtry[
UU]=
MAX(-prfixtry[
UU],ufixtry[UU]);
6063 ufixtry[ENTROPY] = ufixtry[
UU];
6064 ufixtry[ENTROPY] =
MAX(0.0001,
MIN(ufixtry[ENTROPY],1.0));
6070 prfixtry[URAD0] =
MIN(10.0*prfixtry[
RHO],prfixtry[URAD0]);
6072 ufixtry[URAD0]=
MAX(-prfixtry[URAD0],ufixtry[URAD0]);
6077 limit_gamma(0,OUTERDEATHGAMMAMAX,OUTERDEATHGAMMAMAXRAD,prfixtry,NULL,ptrgeom,-1);
6082 if(prfixtry[
U1]<0.0){
6094 if(prfixtry[URAD1]<0.0){
6095 prfixtry[URAD1]=0.0;
6106 prfixtryb[pl] = prfixtry[pl];
6108 ufixtryb[pl] = ufixtry[pl];
6117 prfixtryb[URAD0] =
MIN(prfixtryb[
RHO],prfixtryb[URAD0]);
6123 ufixtryb[URAD0]=
MAX(-prfixtryb[URAD0],ufixtryb[URAD0]);
6130 prfix[pl] = prfixtryb[pl]*switch0b + prfixtry[pl]*switch0*switch2b + prfix[pl]*switch2*switch2b;
6132 ufix[pl] = ufixtryb[pl]*switch0b + ufixtry[pl]*switch0*switch2b + ufix[pl]*switch2*switch2b;
6148 void debugfixup(
void)
6157 prfix=&GLOBALMACP0A1(pglobal,i,j,k,0);
6158 ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
6162 bl_coord_ijk_2(i,j,k,CENT,X, V) ;
6164 if(V[1]>300.0 && V[2]<0.075){
6169 prfix[URAD0] =
MIN(
MIN(prfix[
RHO],prfix[URAD0]),prfix[
UU]);
6172 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
6173 ufix[ENTROPY] = ufix[
UU];
6174 ufix[ENTROPY] =
MAX(0.0001,
MIN(ufix[ENTROPY],1.0));
6175 ufix[URAD0]=
MAX(-prfix[URAD0],ufix[URAD0]);
6179 SLOOPA(jjj) prfix[URAD1+jjj-1] = 0.0;
6183 SLOOPA(jjj) ufix[URAD1+jjj-1]=0.0;
6192 void debugfixupnonmad(
void)
6202 prfix=&GLOBALMACP0A1(pglobal,i,j,k,0);
6203 ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
6207 bl_coord_ijk_2(i,j,k,CENT,X, V) ;
6214 if(V[1]>1E3 && V[2]>M_PI*0.5*1.1){
6220 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
6221 ufix[ENTROPY] = ufix[
UU];
6222 ufix[ENTROPY] =
MAX(0.0001,
MIN(ufix[ENTROPY],1.0));
6225 prfix[URAD0] =
MIN(
MIN(prfix[
RHO],prfix[URAD0]),prfix[
UU]);
6227 ufix[URAD0]=
MAX(-prfix[URAD0],ufix[URAD0]);
6231 SLOOPA(jjj) prfix[URAD1+jjj-1] = 0.0;
6235 SLOOPA(jjj) ufix[URAD1+jjj-1]=0.0;
6244 void debugfixupmad(
void)
6249 struct of_geom *ptrgeom=&geomdontuse;
6254 prfix=&GLOBALMACP0A1(pglobal,i,j,k,0);
6255 ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
6259 bl_coord_ijk_2(i,j,k,CENT,X, V) ;
6266 limit_gamma(0,1.5,1.5,prfix,NULL,ptrgeom,-1);
6269 if(V[1]>60.0 && bsq/prfix[
RHO]>1.0){
6273 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
6274 ufix[ENTROPY] = ufix[
UU];
6275 ufix[ENTROPY] =
MAX(0.0001,
MIN(ufix[ENTROPY],1.0));
6278 prfix[URAD0] =
MIN(
MIN(prfix[
RHO],prfix[URAD0]),prfix[
UU]);
6280 ufix[URAD0]=
MAX(-prfix[URAD0],ufix[URAD0]);
6283 limit_gamma(0,1.5,
GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6287 if(V[1]>900.0 && (V[2]>M_PI*0.5*1.1 || V[2]<M_PI*0.5*0.9) ){
6289 prfix[
RHO] = 1.5E-9*pow(V[1]/500.0,-1.5);
6293 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
6294 ufix[ENTROPY] = ufix[
UU];
6295 ufix[ENTROPY] =
MAX(0.0001,
MIN(ufix[ENTROPY],1.0));
6299 prfix[URAD0] = 1
E-9*pow(V[1]/900.0,-1.0);
6300 prfix[URAD0] =
MIN(
MIN(prfix[
RHO],prfix[URAD0]),prfix[
UU]);
6302 ufix[URAD0]=
MAX(-prfix[URAD0],ufix[URAD0]);
6305 limit_gamma(0,1.5,
GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6314 void debugfixupalt1(
void)
6319 struct of_geom *ptrgeom=&geomdontuse;
6324 prfix=&GLOBALMACP0A1(pglobal,i,j,k,0);
6325 ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
6329 bl_coord_ijk_2(i,j,k,CENT,X, V) ;
6336 limit_gamma(0,1.5,
GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6340 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
6346 if(V[1]>60.0 && bsq/prfix[
RHO]>1.0){
6350 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
6351 ufix[ENTROPY] = ufix[
UU];
6352 ufix[ENTROPY] =
MAX(0.0001,
MIN(ufix[ENTROPY],1.0));
6355 prfix[URAD0] =
MIN(
MIN(prfix[
RHO],prfix[URAD0]),prfix[
UU]);
6357 ufix[URAD0]=
MAX(-prfix[URAD0],ufix[URAD0]);
6360 limit_gamma(0,1.5,
GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6372 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
6373 ufix[ENTROPY] = ufix[
UU];
6374 ufix[ENTROPY] =
MAX(0.0001,
MIN(ufix[ENTROPY],1.0));
6378 prfix[URAD0] =
MIN(
MIN(prfix[
RHO],prfix[URAD0]),prfix[
UU]);
6380 ufix[URAD0]=
MAX(-prfix[URAD0],ufix[URAD0]);
6383 limit_gamma(0,1.5,
GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6386 if(ufix[
U1]<0.0) ufix[
U1]=0.0;
6387 ufix[
U2]=ufix[
U3]=0.0;
6388 if(ufix[URAD1]<0.0) ufix[URAD1]=0.0;
6389 ufix[URAD2]=ufix[URAD3]=0.0;
6401 void debugfixupaltdeath(
void)
6406 struct of_geom *ptrgeom=&geomdontuse;
6411 prfix=&GLOBALMACP0A1(pglobal,i,j,k,0);
6412 ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
6416 bl_coord_ijk_2(i,j,k,CENT,X, V) ;
6423 limit_gamma(0,1.5,
GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6427 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
6433 if(V[1]>60.0 && bsq/prfix[
RHO]>1.0){
6437 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
6438 ufix[ENTROPY] = ufix[
UU];
6439 ufix[ENTROPY] =
MAX(0.0001,
MIN(ufix[ENTROPY],1.0));
6442 prfix[URAD0] =
MIN(
MIN(prfix[
RHO],prfix[URAD0]),prfix[
UU]);
6444 ufix[URAD0]=
MAX(-prfix[URAD0],ufix[URAD0]);
6447 limit_gamma(0,1.5,
GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6460 ufix[
UU]=
MAX(-prfix[
UU],ufix[UU]);
6461 ufix[ENTROPY] = ufix[
UU];
6462 ufix[ENTROPY] =
MAX(0.0001,
MIN(ufix[ENTROPY],1.0));
6466 prfix[URAD0] =
MIN(
MIN(prfix[
RHO],prfix[URAD0]),prfix[
UU]);
6468 ufix[URAD0]=
MAX(-prfix[URAD0],ufix[URAD0]);
6472 limit_gamma(0,6.0,
GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6475 if(ufix[
U1]<0.0) ufix[
U1]=0.0;
6478 if(ufix[URAD1]<0.0) ufix[URAD1]=0.0;