59 #define MAXPASSPARMS 10
65 #define USER_THETAROTMETRIC (0.0) // WALD
66 #define USER_THETAROTPRIMITIVES (0.0) // probably want to choose 0, so initial conditions are as if no tilt // WALD -> make same as USER_THETAROTMETRIC
76 #define BLANDFORDQUAD 5
77 #define TOROIDALFIELD 6
81 #define FIELDJONMAD 10
84 #define SPLITMONOPOLE 13
341 static int firsttime=1;
348 if((fstar=fopen(
"star.txt",
"rt"))==NULL){
349 dualfprintf(
fail_file,
"Couldn't open star.txt, assume values not used.\n");
352 logfprintf(
"opened star.txt and got contents\n");
358 MPI_Barrier(MPI_COMM_WORLD);
368 int set_fieldfrompotential(
int *fieldfrompotential)
393 int fieldfrompotential[
NDIM];
395 set_fieldfrompotential(fieldfrompotential);
397 funreturn=user1_init_conservatives(fieldfrompotential, prim,pstag, Utemp, U);
398 if(funreturn!=0)
return(funreturn);
411 if(funreturn!=0)
return(funreturn);
417 trifprintf(
"Constants\n");
427 #define DIMVARLIST GGG,CCCTRUE,MSUNCM,MPERSUN,LBAR,TBAR,VBAR,RHOBAR,MBAR,ENBAR,UBAR,TEMPBAR,ARAD_CODE_DEF,XFACT,YFACT,ZFACT,MUMEAN,MUMEAN,OPACITYBAR,MASSCM,KORAL2HARMRHO(1.0),TEMPMIN
428 #if(REALTYPE==FLOATYPE || REALTYPE==DOUBLETYPE)
429 #define DIMTYPELIST "%21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g\n"
430 #elif(REALTYPE==LONGDOUBLETYPE)
431 #define DIMTYPELIST "%26.20Lg %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g %26.20g\n"
433 #error "no such type."
437 dimfile=fopen(
"dimensions.txt",
"wt");
443 dualfprintf(
fail_file,
"Could not open dimensions.txt\n");
488 if(funreturn!=0)
return(funreturn);
700 #define NTUBE 31 // harder near t=0 at discontinuity
887 dualfprintf(
fail_file,
"Must choose MCOORD (currently %d) to be spherical polar grid type for RADBEAM2D\n",
MCOORD);
913 tf = 10.0*(M_PI*0.5)*3.0;
917 tf = 10.0*(M_PI*0.5)*6.0;
921 tf = 10.0*(M_PI*0.5)*15.0;
925 tf = 10.0*(M_PI*0.5)*40.0;
957 dualfprintf(
fail_file,
"Must choose MCOORD (currently %d) to be spherical polar grid type for RADBEAM2DKSVERT\n",
MCOORD);
1016 dualfprintf(
fail_file,
"Must choose MCOORD (currently %d) to be spherical polar grid type for ATMSTATIC\n",
MCOORD);
1067 dualfprintf(
fail_file,
"Must choose MCOORD (currently %d) to be spherical polar grid type for RADATM\n",
MCOORD);
1740 dualfprintf(
fail_file,
"RADWAVE_RHOZERO=%g, RADWAVE_KK=%g, RADWAVE_UINT=%g, RADWAVE_TEMP=%g, ARAD_CODE=%g, RADWAVE_ERAD=%21.15g\n",
1744 if((out=fopen(
"radtestparams.dat",
"wt"))==NULL){
1745 dualfprintf(
fail_file,
"Couldn't write radtestparams.dat file\n");
1749 fprintf(out,
"#%20s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s %21s\n",
1779 "RADWAVE_WAVETYPE");
1780 fprintf(out,
"%21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g\n",
2020 dualfprintf(
fail_file,
"Must choose MCOORD (currently %d) to be spherical polar grid type for RADBONDI\n",
MCOORD);
2140 int set_fieldtype(
void);
2351 dualfprintf(
fail_file,
"Must choose MCOORD (currently %d) to be spherical polar grid type for RADNT,\n",
MCOORD);
2355 dualfprintf(
fail_file,
"Must choose MCOORD (currently %d) to be CYLMINKMETRIC for RADCYLBEAM.\n",
MCOORD);
2359 dualfprintf(
fail_file,
"Must choose MCOORD (currently %d) to be CARTMINKMETRIC2 for RADCYLBEAMCART.\n",
MCOORD);
2925 #if(WHICHPROBLEM==KOMIPROBLEM)
3135 dualfprintf(
fail_file,
"WARNING: Have boundary oustide horizon in KS coords\n");
3298 dualfprintf(
fail_file,
"RADDONUT setup for 128x64 with that grid\n");
3470 funreturn=user1_init_atmosphere(whichvel, whichcoord,i, j, k, pr);
3471 if(funreturn!=0)
return(funreturn);
3495 int init_grid_post_set_grid(
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*pstag)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*ucons)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*vpot)[
NSTORE1+
SHIFTSTORE1][
NSTORE2+
SHIFTSTORE2][
NSTORE3+
SHIFTSTORE3],
FTYPE (*Bhat)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*panalytic)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*pstaganalytic)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*vpotanalytic)[
NSTORE1+
SHIFTSTORE1][
NSTORE2+
SHIFTSTORE2][
NSTORE3+
SHIFTSTORE3],
FTYPE (*Bhatanalytic)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*F1)[
NSTORE2][
NSTORE3][NPR+
NSPECIAL],
FTYPE (*F2)[
NSTORE2][
NSTORE3][NPR+
NSPECIAL],
FTYPE (*F3)[
NSTORE2][
NSTORE3][NPR+
NSPECIAL],
FTYPE (*Atemp)[
NSTORE1+
SHIFTSTORE1][
NSTORE2+
SHIFTSTORE2][
NSTORE3+
SHIFTSTORE3])
3499 extern void check_spc_singularities_user(
void);
3510 int set_fieldtype(
void);
3594 trifprintf(
"RADDOT: myX: %g %g %g\n",myX[1],myX[2],myX[3]);
3597 extern void icoord_round(
FTYPE *X,
int loc,
int *i,
int *j,
int *k);
3610 trifprintf(
"BEGIN check_rmin\n");
3613 trifprintf(
"END check_rmin\n");
3618 trifprintf(
"BEGIN check_spc_singularities_user\n");
3620 dualfprintf(
fail_file,
"WARNING: check_spc_singularities_user() oddly stalls sometimes...\n");
3621 check_spc_singularities_user();
3622 trifprintf(
"END check_spc_singularities_user\n");
3623 dualfprintf(
fail_file,
"WARNING: done with check_spc_singularities_user(), but it sometimes stalls or goes very very slow for no apparently good reason. E.g., on NAUTILUS with -O0, very slow checks. But just putting dualfprintf before and after the above call leads to quick finish.\n");
3633 int init_primitives(
FTYPE (*prim)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*pstag)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*ucons)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*vpot)[
NSTORE1+
SHIFTSTORE1][
NSTORE2+
SHIFTSTORE2][
NSTORE3+
SHIFTSTORE3],
FTYPE (*Bhat)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*panalytic)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*pstaganalytic)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*vpotanalytic)[
NSTORE1+
SHIFTSTORE1][
NSTORE2+
SHIFTSTORE2][
NSTORE3+
SHIFTSTORE3],
FTYPE (*Bhatanalytic)[
NSTORE2][
NSTORE3][NPR],
FTYPE (*F1)[
NSTORE2][
NSTORE3][NPR+
NSPECIAL],
FTYPE (*F2)[
NSTORE2][
NSTORE3][NPR+
NSPECIAL],
FTYPE (*F3)[
NSTORE2][
NSTORE3][NPR+
NSPECIAL],
FTYPE (*Atemp)[
NSTORE1+
SHIFTSTORE1][
NSTORE2+
SHIFTSTORE2][
NSTORE3+
SHIFTSTORE3])
3646 funreturn=
user1_init_primitives(inittype, prim, pstag, ucons, vpot, Bhat, panalytic, pstaganalytic, vpotanalytic, Bhatanalytic, F1, F2, F3,Atemp);
3647 if(funreturn!=0)
return(funreturn);
3678 int set_fieldtype(
void);
3684 int whichinversion=0;
3685 fieldprim(whichmethod, whichinversion, whichvel, whichcoord, i, j, k, pr);
3694 if(YFL1>=0) pr[YFL1] =
SMALL + rhofloor;
3695 if(YFL2>=0) pr[YFL2] =
SMALL + rhofloor*vfloor*vfloor;
3696 if(YFL3>=0) pr[YFL3] =
SMALL + rhofloor*vfloor;
3697 if(YFL4>=0) pr[YFL4] =
SMALL + enfloor;
3698 if(YFL5>=0) pr[YFL5] =
SMALL + enfloor*vfloor;
3783 FTYPE pradffortho[NPR];
3785 pradffortho[PRAD1] = 0;
3786 pradffortho[PRAD2] = 0;
3787 pradffortho[PRAD3] = 0;
3793 prad_fforlab(whichvel, whichcoord,
FF2LAB, i,j,k,
CENT,NULL,pradffortho,pr, pr);
3820 FTYPE Trad,Tgas,ERAD,uint;
3830 rsq=(xx)*(xx)+(yy)*(yy)+(zz)*(zz);
3846 Trad=T_AMB*(1.+BLOBP*exp(-rsq/(BLOBW*BLOBW)));
3897 FTYPE rho,mx,my,mz,m,ERAD,uint,E0,Fx,Fy,Fz,pLTE;
3898 FTYPE rho0,Tgas0,
ur,Tgas,Trad,r,rcm,prad,pgas,vx,ut,ux;
3910 if(
NTUBE==1) {uint = 3.e-5 / (
gamideal - 1.); ERAD=1.e-8; Fx=1.e-2*ERAD;ux=0.015;}
3911 if(
NTUBE==2) {uint = 4.e-3 / (
gamideal - 1.);ERAD=2.e-5; Fx=1.e-2*ERAD;ux=0.25;}
3912 if(
NTUBE==3 ||
NTUBE==31) {uint = 60. / (
gamideal - 1.);ERAD=2.; Fx=1.e-2*ERAD;ux=10.;}
3913 if(
NTUBE==4 ||
NTUBE==41) {uint = 6.e-3 / (
gamideal - 1.);ERAD=0.18; Fx=1.e-2*ERAD;ux=0.69;}
3914 if(
NTUBE==5) {uint = 60. / (
gamideal - 1.);ERAD=2.; Fx=1.e-2*ERAD;ux=1.25;}
3917 if(
NTUBE==1) {rho=2.4;uint = 1.61e-4/ (
gamideal - 1.); ERAD=2.51e-7; Fx=1.e-2*ERAD;ux=6.25e-3;}
3918 if(
NTUBE==2) {rho=3.11;uint = 0.04512 / (
gamideal - 1.);ERAD=3.46e-3; Fx=1.e-2*ERAD;ux=0.0804;}
3919 if(
NTUBE==3 ||
NTUBE==31) {rho=8.0;uint = 2.34e3 / (
gamideal - 1.);ERAD=1.14e3; Fx=1.e-2*ERAD;ux=1.25;}
3920 if(
NTUBE==4 ||
NTUBE==41) {rho=3.65;uint =3.59e-2 / (
gamideal - 1.);ERAD=1.30; Fx=1.e-2*ERAD;ux=0.189;}
3921 if(
NTUBE==5) {rho=1.0;uint = 60. / (
gamideal - 1.);ERAD=2.; Fx=1.e-2*ERAD;ux=1.10;}
3954 FTYPE pradffortho[NPR];
3955 pradffortho[PRAD0] = ERAD;
3956 pradffortho[PRAD1] = Fx;
3957 pradffortho[PRAD2] = Fy;
3958 pradffortho[PRAD3] = Fz;
3962 prad_fforlab(whichvel, whichcoord,
FF2LAB, i,j,k,
CENT,NULL,pradffortho,pr, pr);
3986 FTYPE Trad,Tgas,ERAD;
3994 FTYPE rho,uint,Fx,Fy,Fz,pLTE;
3999 rsq=xx*xx+yy*yy+zz*zz;
4000 rho=(RHOBLOB-RHOAMB)*exp(-sqrt(rsq)/(BLOBW*BLOBW))+RHOAMB;
4001 Tgas=TAMB*RHOAMB/rho;
4045 FTYPE pradffortho[NPR];
4046 pradffortho[PRAD0] = ERAD;
4047 pradffortho[PRAD1] = Fx;
4048 pradffortho[PRAD2] = Fy;
4049 pradffortho[PRAD3] = Fz;
4052 prad_fforlab(whichvel, whichcoord,
FF2LAB, i,j,k,
CENT,NULL,pradffortho,pr, pr);
4140 Vr=sqrt(2./r)*(1.-2./r);
4144 struct of_geom geomrealdontuse;
4145 struct of_geom *ptrgeomreal=&geomrealdontuse;
4146 gset(getprim,*whichcoord,i,j,k,ptrgeomreal);
4148 FTYPE W=1./sqrt(1.-Vr*Vr*ptrgeomreal->gcov[
GIND(1,1)]);
4179 pr[PRAD0] = ERADAMB;
4185 FTYPE pradffortho[NPR];
4186 pradffortho[PRAD0] = ERADAMB;
4187 pradffortho[PRAD1] = Fx;
4188 pradffortho[PRAD2] = Fy;
4189 pradffortho[PRAD3] = Fz;
4192 prad_fforlab(whichvel, whichcoord,
FF2LAB, i,j,k,
CENT,NULL,pradffortho,pr, pr);
4224 uint=(4*r*u0 - 4*r*u0*
gamideal - 2*r*rho0 + 2*r0*rho0 + r*r0*rho0*
Log(-2 + r) - r*r0*rho0*
Log(r) - r*r0*rho0*
Log(-2 + r0) + r*r0*rho0*
Log(r0))/(4*r - 4*r*
gamideal);
4226 FTYPE uradx,urady,uradz;
4227 uradx=urady=uradz=0.;
4248 pr[PRAD0] = uint*1
E-20;
4280 #define WHICHRADATM 0 // 0,1,2,3
4301 FTYPE kappaesperrho=calc_kappaes_user(1,0, 0,0,0);
4319 f = (
FTYPE)kappaesperrho*FLUXLEFT*MINX*MINX;
4336 FTYPE Fx=FLUXLEFT*(MINX/xx)*(MINX/xx);
4347 dualfprintf(
fail_file,
"i=%d f=%g p0=%g KKK=%g C3=%g rho=%g uint=%g Fx=%g ERAD=%g : kappaesperrho=%g \n",i,f,p0,KKK,C3,rho,uint,Fx,ERAD , kappaesperrho);
4349 dualfprintf(
fail_file,
"i=%d f=%g p0=%g KKK=%g C3=%g rho=%g uint=%g Fx=%g ERAD=%g : kappaesperrho=%g \n",i,f,p0*
UBAR,KKK*UBAR/pow(
RHOBAR,
gamideal),C3,rho*
RHOBAR,uint*UBAR,Fx*
ENBAR/
TBAR/
LBAR/
LBAR,ERAD*UBAR , kappaesperrho*
OPACITYBAR);
4381 FTYPE pradffortho[NPR];
4382 pradffortho[PRAD0] = ERAD;
4383 pradffortho[PRAD1] = Fx;
4384 pradffortho[PRAD2] = Fy;
4385 pradffortho[PRAD3] = Fz;
4388 prad_fforlab(whichvel, whichcoord,
FF2LAB, i,j,k,
CENT,NULL,pradffortho, pr, pr);
4394 struct of_geom geomrealdontuse;
4395 struct of_geom *ptrgeomreal=&geomrealdontuse;
4396 gset(getprim,*whichcoord,i,j,k,ptrgeomreal);
4398 dualfprintf(
fail_file,
"AFTER: i=%d rho=%g uint=%g vx=%g ERAD=%g uradx=%g\n",i,pr[
RHO]*
RHOBAR,pr[
UU]*
UBAR,pr[
U1]*sqrt(ptrgeomreal->gcov[
GIND(1,1)])*
VBAR,pr[URAD0]*UBAR,pr[URAD1]*sqrt(ptrgeomreal->gcov[
GIND(1,1)])*
VBAR);
4464 FTYPE rho,ERAD,uint,Fx,Fy,Fz;
4569 FTYPE pradffortho[NPR];
4570 pradffortho[PRAD0] = ERAD;
4571 pradffortho[PRAD1] = Fx;
4572 pradffortho[PRAD2] = Fy;
4573 pradffortho[PRAD3] = Fz;
4577 prad_fforlab(whichvel, whichcoord,
FF2LAB, i,j,k,
CENT,NULL, pradffortho, pr, pr);
4604 FTYPE pleft[NPR], pright[NPR], P;
4621 pleft[
UU] = P/(
gam-1);
4625 pright[
U2] = 0.3923;
4631 pright[
RHO] = 25.48;
4632 pright[
UU] = P/(
gam-1);
4645 pleft[
UU] = P/(
gam-1);
4649 pright[
U2] = -0.6822;
4655 pright[
RHO] = 3.323;
4656 pright[
UU] = P/(
gam-1);
4669 pleft[
UU] = P/(
gam-1);
4672 pright[
U1] = -0.212;
4673 pright[
U2] = -0.590;
4679 pright[
RHO] = 0.562;
4680 pright[
UU] = P/(
gam-1);
4692 pleft[
RHO] = 1.78e-3;
4693 pleft[
UU] = P/(
gam-1);
4704 pright[
UU] = P/(
gam-1);
4717 pleft[
UU] = P/(
gam-1);
4724 pright[
B2] = -6.857;
4728 pright[
UU] = P/(
gam-1);
4741 pleft[
UU] = P/(
gam-1);
4752 pright[
UU] = P/(
gam-1);
4765 pleft[
UU] = P/(
gam-1);
4776 pright[
UU] = P/(
gam-1);
4789 pleft[
UU] = P/(
gam-1);
4800 pright[
UU] = P/(
gam-1);
4807 if(
WHICHKOMI==6){ xcl=-0.025; xcr=0.0; }
4808 if(x<=xcl)
PALLLOOP(pl) pr[pl] = pleft[pl];
4809 else if(x>xcr)
PALLLOOP(pl) pr[pl] = pright[pl];
4817 if(
WHICHKOMI==6){ xcl=-0.025; xcr=0.0; }
4819 if(x<=xcl) phi0=0.0;
4820 else if((x>xcl)&&(x<xcr)) phi0=0.5*M_PI*(x-xcl)/(xcr-xcl);
4821 else if(x>=xcr) phi0=0.5*M_PI;
4824 if(x<=xcl) phi1=0.0;
4825 else if((x>xcl)&&(x<xcr)) phi1=M_PI*(x-xcl)/(xcr-xcl);
4826 else if(x>=xcr) phi1=M_PI;
4830 PALLLOOP(pl) pr[pl] = pleft[pl] + (pright[pl]-pleft[pl])*sin(phi0);
4835 FTYPE gamma = sqrt(1.0 + fabs(usq));
4838 FTYPE gammaleft = sqrt(1.0 + fabs(usqleft));
4842 (-1.*pleft[
B1]*pleft[
U1] - 1.*pleft[
B2]*pleft[
U2])*
4843 (pleft[B1]*pleft[
U1] + pleft[
B2]*pleft[
U2]) +
4844 Power(pleft[B1]/gammaleft + (pleft[
U1]*
4845 (pleft[B1]*pleft[
U1] + pleft[
B2]*pleft[U2]))/gammaleft,2) +
4846 Power(pleft[
B2]/gammaleft + (pleft[U2]*
4847 (pleft[B1]*pleft[
U1] + pleft[
B2]*pleft[U2]))/gammaleft,2);
4849 dualfprintf(
fail_file,
"x=%g mybsqconst=%g\n",x,mybsqconst);
4873 if(disc1>=0.0){
PALLLOOP(pl)
if(pl==B3) pr[pl] = +sqrt(disc1); }
4874 else if(disc2>=0){
PALLLOOP(pl)
if(pl==B3) pr[pl] = -sqrt(disc2); }
4875 else PALLLOOP(pl){
if(pl==B3) pr[pl] = 0.0; }
4910 pr[
U1]=pr[
U2]=pr[
U3]=0.0;
4930 struct of_geom *ptrgeom=&geomdontuse;
4931 gset(getprim,*whichcoord,i,j,k,ptrgeom);
4956 if(x-x0<=-dx0) B[2]=1.0;
4957 else if((x-x0>-dx0)&&(x-x0<dx0)) B[2]=1.0-(0.3/(2.0*dx0))*((x-x0)+dx0);
4958 else if(x-x0>=dx0) B[2]=0.7;
4969 EBtopr(E,B,ptrgeom,pr);
4976 computeKK(pr,ptrgeom,&KK);
4978 dualfprintf(
fail_file,
"i=%d KK=%21.15g\n",i,KK);
4994 if(x-x0<=-dx0) bcon[2]=1.0;
4995 else if((x-x0>-dx0)&&(x-x0<dx0)) bcon[2]=1.0-(0.3/(2.0*dx0))*((x-x0)+dx0);
4996 else if(x-x0>=dx0) bcon[2]=0.7;
5007 vbtopr(vcon,bcon,ptrgeom,pr);
5009 computeKK(pr,ptrgeom,&KK);
5011 dualfprintf(
fail_file,
"i=%d KK=%21.15g\n",i,KK);
5019 bcon[1]=bcon[2]=1.0;
5022 if(x-x0<=-0.1) bcon[3]=1.0;
5023 else if((x-x0>-0.1)&&(x-x0<0.1)) bcon[3]=1.0+3.0/2.0*((x-x0)+0.1);
5024 else if(x-x0>=0.1) bcon[3]=1.3;
5026 econ[2]=econ[3]=0.0;
5032 FTYPE CONSTECON1=1.0;
5033 econ[1]=sqrt(-CONSTECON1 + bcon[3]*bcon[3]);
5044 EBvetatopr(econ, bcon, vcon, ptrgeom, pr);
5047 computeKK(pr,ptrgeom,&KK);
5049 dualfprintf(
fail_file,
"i=%d KK=%21.15g\n",i,KK);
5061 if(x-x0<=-0.1) phi0=0.0;
5062 else if((x-x0>-0.1)&&(x-x0<0.1)) phi0=5.0/2.0*M_PI*((x-x0)+0.1);
5063 else if(x-x0>=0.1) phi0=M_PI*0.5;
5065 bcon[2]=2.0*cos(phi0);
5066 bcon[3]=2.0*sin(phi0);
5075 vbtopr(vcon,bcon,ptrgeom,pr);
5077 computeKK(pr,ptrgeom,&KK);
5079 dualfprintf(
fail_file,
"i=%d KK=%21.15g\n",i,KK);
5111 EBtopr(E,B,ptrgeom,pr);
5114 computeKK(pr,ptrgeom,&KK);
5116 dualfprintf(
fail_file,
"i=%d KK=%21.15g\n",i,KK);
5153 EBtopr(E,B,ptrgeom,pr);
5156 computeKK(pr,ptrgeom,&KK);
5158 dualfprintf(
fail_file,
"i=%d KK=%21.15g\n",i,KK);
5184 else if((x-x0>=-0.1)&&(x-x0<=0.1)){
5186 B[2]=1.0+(x-x0+0.1)*(-2.0/0.2);
5187 B[3]=1.0+(x-x0+0.1)*(-2.0/0.2);
5199 EBtopr(E,B,ptrgeom,pr);
5202 computeKK(pr,ptrgeom,&KK);
5204 dualfprintf(
fail_file,
"i=%d KK=%21.15g\n",i,KK);
5226 else if(x>=0.2+0.5){
5230 B[3]=1.0+0.15*(1.0+sin(5.0*M_PI*(x-0.1-0.5)));
5239 EBtopr(E,B,ptrgeom,pr);
5242 computeKK(pr,ptrgeom,&KK);
5244 dualfprintf(
fail_file,
"i=%d KK=%21.15g\n",i,KK);
5271 EBtopr(E,B,ptrgeom,pr);
5274 computeKK(pr,ptrgeom,&KK);
5276 dualfprintf(
fail_file,
"i=%d KK=%21.15g\n",i,KK);
5321 struct of_geom geomrealdontuse;
5322 struct of_geom *ptrgeomreal=&geomrealdontuse;
5323 gset(getprim,*whichcoord,i,j,k,ptrgeomreal);
5326 FTYPE rho,ERAD,uint;
5327 FTYPE rho0,Tgas0,
ur,Tgas,Trad,r,rcm,prad,pgas,vx,ut;
5341 ut=sqrt((-1.-ur*ur*ptrgeomreal->gcov[
GIND(1,1)])/ptrgeomreal->gcov[
GIND(0,0)]);
5344 Tgas=Tgas0*pow(rho/rho0,
gam-1.);
5387 FTYPE pradffortho[NPR];
5388 pradffortho[PRAD0] = ERAD;
5389 pradffortho[PRAD1] = Fx;
5390 pradffortho[PRAD2] = Fy;
5391 pradffortho[PRAD3] = Fz;
5395 dualfprintf(
fail_file,
"i=%d rho=%g uint=%g ERAD=%g\n",i,rho,uint,ERAD);
5402 prad_fforlab(whichvel, whichcoord,
FF2LAB, i,j,k,
CENT,ptrgeomreal, pradffortho, pr, pr);
5405 dualfprintf(
fail_file,
"AFTER: i=%d rho=%g uint=%g vx=%g ERAD=%g uradx=%g\n",i,pr[
RHO]*
RHOBAR,pr[
UU]*
UBAR,pr[
U1]*sqrt(ptrgeomreal->gcov[
GIND(1,1)])*
VBAR,pr[URAD0]*UBAR,pr[URAD1]*sqrt(ptrgeomreal->gcov[
GIND(1,1)])*
VBAR);
5456 FTYPE pradffortho[NPR];
5458 pradffortho[PRAD1] = 0;
5459 pradffortho[PRAD2] = 0;
5460 pradffortho[PRAD3] = 0;
5473 prad_fforlab(whichvel, whichcoord,
FF2LAB, i,j,k,
CENT,NULL, pradffortho, pr, pr);
5494 struct of_geom geomraddontuse;
5495 struct of_geom *ptrgeomrad=&geomraddontuse;
5496 gset(getprim,*whichcoord,i,j,k,ptrgeomrad);
5504 int set_fieldtype(
void);
5517 set_zamo_velocity(*whichvel,ptrgeomrad,pr);
5538 FTYPE pradffortho[NPR],pradfforthoatm[NPR];
5544 pradfforthoatm[PRAD1]=0;
5545 pradfforthoatm[PRAD2]=0;
5546 pradfforthoatm[PRAD3]=0;
5549 pradffortho[PRAD0]=pradfforthoatm[PRAD0];
5550 pradffortho[PRAD1]=pradfforthoatm[PRAD1];
5551 pradffortho[PRAD2]=pradfforthoatm[PRAD2];
5552 pradffortho[PRAD3]=pradfforthoatm[PRAD3];
5561 pr[PRAD0]=pradffortho[PRAD0];
5562 pr[PRAD1]=pradffortho[PRAD1];
5563 pr[PRAD2]=pradffortho[PRAD2];
5564 pr[PRAD3]=pradffortho[PRAD3];
5572 pradffortho[PRAD0]=pr[PRAD0];
5573 pradffortho[PRAD1]=pr[PRAD1];
5574 pradffortho[PRAD2]=pr[PRAD2];
5575 pradffortho[PRAD3]=pr[PRAD3];
5593 prad_fforlab(whichvel, whichcoord,
FF2LAB, i,j,k,
CENT,ptrgeomrad, pradffortho, pr, pr);
5599 if(debugfail>=2 && returndonut==0 && pradfforthoatm[PRAD0]>pradffortho[PRAD0]){
5600 dualfprintf(
fail_file,
"WARNING: Torus radiation pressure below atmosphere.\n");
5601 dualfprintf(
fail_file,
"CHECKPOST: ijk=%d %d %d : %g %g %g %g\n",i,j,k,pr[PRAD0],pr[PRAD1],pr[PRAD2],pr[PRAD3]);
5608 FTYPE prreport[NPR];
5609 set_zamo_velocity(*whichvel,ptrgeomrad,prreport);
5610 dualfprintf(
fail_file,
"ZAMO: ijk=%d %d %d : %g %g %g : fluid: %g %g %g\n",i,j,k,prreport[
U1],prreport[
U2],prreport[
U3],pr[U1],pr[U2],pr[U3]);
5612 DLOOP(jj,kk) dualfprintf(
fail_file,
"gn%d%d=%26.20g\n",jj+1,kk+1,ptrgeomrad->gcon[
GIND(jj,kk)]);
5613 DLOOP(jj,kk) dualfprintf(
fail_file,
"gv%d%d=%26.20g\n",jj+1,kk+1,ptrgeomrad->gcov[
GIND(jj,kk)]);
5625 set_zamo_velocity(*whichvel,ptrgeomrad,&pr[URAD1-
U1]);
5692 #define fsx(x) ((x)>0 ? exp(-1.0/(x)) : 0.0)
5693 #define gsx(x) (fsx(x)/(fsx(x) + fsx(1.0-(x))))
5694 #define stepfunction(x) (1.0 - gsx((x)-0.5))
5695 #define stepfunctionab(x,a,b) (((a)-(b))*stepfunction(x) + (b))
5696 #define stepfunctionab2(x,a,b) (((a)-(b))*stepfunction(x)*stepfunction(x) + (b))
5716 E0=E0perRho0*Rho0*rhojet;
5726 E0=E0perRho0*Rho0*rhojet;
5739 pr[
U2]=vz0*stepfunctionab(r,1,0.001);
5749 pr[URAD0]=Ehat0*stepfunctionab(r,Ehatjet/Ehat0,Ehatstar/Ehat0);
5759 FTYPE ptot=1.1*(4.0/3.0-1.0)*Ehatstar;
5760 pr[
UU]=( ptot - (4.0/3.0-1.0)*pr[URAD0])/(
gam-1.0);
5763 pr[
UU]=0.1*pr[URAD0];
5789 FTYPE Ehat0=0.1*rhojet;
5791 FTYPE Ehatjet=Ehat0*1
E-2;
5813 FTYPE ptot=1.5*(4.0/3.0-1.0)*Ehatjet;
5814 pr[
UU]=( ptot - (4.0/3.0-1.0)*pr[URAD0])/(
gam-1.0);
5834 FTYPE rhojet=0.1*100.0;
5837 FTYPE Ehat0=0.1*rhojet;
5839 FTYPE Ehatjet=Ehat0*1
E-2;
5889 FTYPE Rho0=1
E-5*100.0*100.0/79.4*40.0;
5912 pr[
UU]=u_rho0_T_simple(i,j,k,
CENT,pr[
RHO],Tstar);
5916 static int firsttime=1;
5920 FTYPE ptot=((
gam-1.0)*pr[
UU] + (4.0/3.0-1.0)*pr[URAD0]);
5924 FTYPE gamptot=(
gam*(
gam-1.0)*pr[
UU] + (4.0/3.0)*(4.0/3.0-1.0)*pr[URAD0]);
5926 dualfprintf(
fail_file,
"STAR: ptrue=%21.15g ptot=%21.15g ptrue/rho=%21.15g cs2=%21.15g cs=%21.15g\n",ptrue,ptot,ptrue/pr[RHO],cs2,sqrt(cs2));
5943 static int firsttime=1;
5944 if(myid==0&&firsttime==1&&i==0&&j==0&&k==0){
5947 if((fstar=fopen(
"star.txt",
"wt"))==NULL){
5948 dualfprintf(
fail_file,
"Couldn't open star.txt\n");
5979 int i=(*ptrptrgeom)->i;
5980 int j=(*ptrptrgeom)->j;
5981 int k=(*ptrptrgeom)->k;
5982 int loc=(*ptrptrgeom)->p;
5986 PLOOP(pliter,pl) ppback[pl]=pp[pl];
5988 int whichvelback=*whichvel;
5992 int whichcoordreal=
MCOORD;
5995 struct of_geom *ptrgeombl=&geombldontuse;
5996 gset(0,whichcoordreal,i,j,k,ptrgeombl);
6026 FTYPE kappa,kappaes,chi;
6033 calc_kappa_user(pp[RHO],B,Tg,Tg,1.0,V[1],V[2],V[3])
6039 FTYPE pptemp[NPR],E1=0,E2=0;
6043 struct of_geom *ptrgeomt=&geomtdontuse;
6049 Xtemp1[1]=1.01*X[1];
6053 gset_X(getprim,whichcoordback,i,j,k,
NOWHERE,Xtemp1,ptrgeomt);
6054 PLOOP(pliter,pl) pptemp[pl]=ppback[pl];
6056 anretlocal=
make_nonrt2rt_solution(&whichvelback, &whichcoordback, opticallythick, pptemp,Xtemp1,Vtemp1,&ptrgeomt);
6057 if(anretlocal<0) anretmin=-1;
6067 gset_X(getprim,whichcoordback,i,j,k,
NOWHERE,Xtemp2,ptrgeomt);
6068 PLOOP(pliter,pl) pptemp[pl]=ppback[pl];
6070 anretlocal=
make_nonrt2rt_solution(&whichvelback, &whichcoordback, opticallythick, pptemp,Xtemp2,Vtemp2,&ptrgeomt);
6071 if(anretlocal<0) anretmin=-1;
6078 Fx=-
THIRD*(E2-E1)/((Vtemp2[1]-Vtemp1[1])*sqrt(fabs(ptrgeombl->gcov[
GIND(1,1)])))/chi;
6086 Xtemp1[2]=1.01*X[2];
6089 gset_X(getprim,whichcoordback,i,j,k,
NOWHERE,Xtemp1,ptrgeomt);
6090 PLOOP(pliter,pl) pptemp[pl]=ppback[pl];
6092 anretlocal=
make_nonrt2rt_solution(&whichvelback, &whichcoordback, opticallythick, pptemp,Xtemp1,Vtemp1,&ptrgeomt);
6093 if(anretlocal<0) anretmin=-1;
6100 Xtemp2[2]=0.99*X[2];
6103 gset_X(getprim,whichcoordback,i,j,k,
NOWHERE,Xtemp2,ptrgeomt);
6104 PLOOP(pliter,pl) pptemp[pl]=ppback[pl];
6106 anretlocal=
make_nonrt2rt_solution(&whichvelback, &whichcoordback, opticallythick, pptemp,Xtemp2,Vtemp2,&ptrgeomt);
6107 if(anretlocal<0) anretmin=-1;
6114 Fy=-
THIRD*(E2-E1)/((Vtemp2[2]-Vtemp1[2])*sqrt(fabs(ptrgeombl->gcov[
GIND(2,2)])))/chi;
6124 Xtemp1[3]=1.01*X[3];
6126 gset_X(getprim,whichcoordback,i,j,k,
NOWHERE,Xtemp1,ptrgeomt);
6127 PLOOP(pliter,pl) pptemp[pl]=ppback[pl];
6129 anretlocal=
make_nonrt2rt_solution(&whichvelback, &whichcoordback, opticallythick, pptemp,Xtemp1,Vtemp1,&ptrgeomt);
6130 if(anretlocal<0) anretmin=-1;
6138 Xtemp2[3]=0.99*X[3];
6140 gset_X(getprim,whichcoordback,i,j,k,
NOWHERE,Xtemp2,ptrgeomt);
6141 PLOOP(pliter,pl) pptemp[pl]=ppback[pl];
6143 anretlocal=
make_nonrt2rt_solution(&whichvelback, &whichcoordback, opticallythick, pptemp,Xtemp2,Vtemp2,&ptrgeomt);
6144 if(anretlocal<0) anretmin=-1;
6150 Fz=-
THIRD*(E2-E1)/((Vtemp2[3]-Vtemp1[3])*sqrt(fabs(ptrgeombl->gcov[
GIND(3,3)])))/chi;
6170 FTYPE Fl=sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
6175 FTYPE Flnew=sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
6176 if(
PRODUCTION==0) dualfprintf(
fail_file,
"limited: Fl=%g E=%g Fx=%g Fy=%g Fz=%g : Flnew=%g\n",Fl,E,Fx,Fy,Fz,Flnew);
6193 if(anret==0)
return(anret);
6205 int i=(*ptrptrgeom)->i;
6206 int j=(*ptrptrgeom)->j;
6207 int k=(*ptrptrgeom)->k;
6208 int p=(*ptrptrgeom)->p;
6246 #define NUMTgasITERS 4 // for now, was 4.
6251 FTYPE kappaabs,kappaes,kappatot;
6256 kappaabs=calc_kappa_user(rho,B,Tgas,Tgas,1.0,V[1],V[2],V[3]);
6257 kappaes=calc_kappaes_user(rho,Tgas,V[1],V[2],V[3]);
6258 kappatot=kappaabs+kappaes;
6269 dl = r*
MAX(dl,
dx[2]*dxdxp[2][2]);
6271 FTYPE tautot,tauabs;
6272 tautot = kappatot*dl+
SMALL;
6273 tauabs = kappaabs*dl+
SMALL;
6274 FTYPE gtau = (tautot/2.0 + 1.0/sqrt(3.0))/(tautot/2.0 + 1.0/sqrt(3.0) + 1.0/(3.0*tauabs));
6291 Tgas=-
Sqrt((-4.*
Power(0.6666666666666666,0.3333333333333333)*P)/naw1 + naw1/(
Power(2.,0.3333333333333333)*
Power(3.,0.6666666666666666)*aaa))/2. +
Sqrt((4.*
Power(0.6666666666666666,0.3333333333333333)*P)/naw1 - naw1/(
Power(2.,0.3333333333333333)*
Power(3.,0.6666666666666666)*aaa) + (2.*bbb)/(aaa*
Sqrt((-4.*
Power(0.6666666666666666,0.3333333333333333)*P)/naw1 + naw1/(
Power(2.,0.3333333333333333)*
Power(3.,0.6666666666666666)*aaa))))/2.;
6295 if(fabs(rho*Tgas + aaa*pow(Tgas,4.0)*gtau - P)>1E-5){
6296 dualfprintf(
fail_file,
"2: NOT right equation: rho=%26.21g Tgas=%26.21g aaa=%26.21g Ptried=%26.21g P=%26.21g\n",rho,Tgas,aaa,rho*Tgas + aaa*pow(Tgas,4.),P);
6359 FTYPE Vphi=0.0,Vr=0.0,Vh=0.0;
6360 FTYPE D,W,uT=1.0,uphi,uPhi,rho,ucon[
NDIM],uint,
E,Fx,Fy,Fz;
6366 PLOOP(pliter,pl) ppback[pl]=pp[pl];
6371 FTYPE Rcyl=fabs(r*sin(th));
6372 int i=(*ptrptrgeom)->i;
6373 int j=(*ptrptrgeom)->j;
6374 int k=(*ptrptrgeom)->k;
6375 int p=(*ptrptrgeom)->p;
6380 if(opticallythick==1) gamtorus=4.0/3.0;
6392 if(*whichcoord!=mycoords){
6395 gset_X(getprim,mycoords,i,j,k,
NOWHERE,X,*ptrptrgeom);
6397 *whichcoord=mycoords;
6421 FTYPE xr=pow(r,xrpow);
6422 FTYPE xr0=pow(Risco,xrpow);
6423 FTYPE xr1=pow(1.0,xrpow);
6434 FTYPE z = r*cos(th) ;
6435 FTYPE zisco = Risco*cos(th) ;
6450 FTYPE NN=1.0/(gamtorus-1.0);
6460 FTYPE rhopowisco=7.0;
6461 rho = rho*pow(r/Rtrans,rhopowisco);
6466 if(rho<ppback[
RHO]){
6477 FTYPE omegakep=1./(pow(r,1.5) +
a);
6478 FTYPE omegakepisco=1./(pow(Risco,1.5) +
a);
6479 FTYPE omega,omegaisco=omegakepisco;
6481 FTYPE sigma=r*r+
a*
a*cos(th)*cos(th);
6482 FTYPE gppvsr=sigma +
a*
a*(1+2*r/sigma)*sin(th)*sin(th);
6483 FTYPE sigmaisco=Risco*Risco+
a*
a*cos(th)*cos(th);
6484 FTYPE gppvsrisco=sigmaisco +
a*
a*(1+2*Risco/sigmaisco)*sin(th)*sin(th);
6489 omega=omegaisco*gppvsrisco/gppvsr;
6494 FTYPE gtt=-(1.0-2.0*r/sigma);
6495 FTYPE gpp=sin(th)*sin(th)*(sigma+
a*
a*(1.0+2.0*r/sigma)*sin(th)*sin(th));
6496 uT = 1.0/pow(-gtt-uPhi*uPhi*gpp,0.5);
6498 FTYPE AA=pow(r*r+a*a,2.0)-a*a*Delta*sin(th)*sin(th);
6499 FTYPE gttbl=-(1.0-2.0*r/sigma);
6500 FTYPE gppbl=AA*sin(th)*sin(th)/sigma;
6501 FTYPE uTbl = 1.0/pow(-gttbl-uPhi*uPhi*gppbl,0.5);
6551 else if(r<1.1*
Rhor){
6583 uint = rho*pow(H/R,2.0)*pow(r*omega,2.0)/(gamtorus*(gamtorus-1.0))*uintfact;
6584 FTYPE uintisco = rhoisco*pow(Hisco/Risco,2.0)*pow(Risco*omegaisco,2.0)/(gamtorus*(gamtorus-1.0))*uintfact;
6588 FTYPE pint=(gamtorus-1.0)*uint;
6589 FTYPE KK= pint*pow(rho,-gamtorus);
6590 FTYPE pintisco=(gamtorus-1.0)*uintisco;
6591 FTYPE KKisco= pintisco*pow(rhoisco,-gamtorus);
6597 pint = pintisco*pow(rhoisco,-gamtorus)*pow(rho,gamtorus);
6600 pt = uint * (gamtorus-1.0);
6618 if(*whichcoord!=mycoords){
6621 gset_X(getprim,mycoords,i,j,k,
NOWHERE,X,*ptrptrgeom);
6623 *whichcoord=mycoords;
6642 FTYPE SMALLESTH=1E-15;
6648 FTYPE xr=pow(r,xrpow);
6650 FTYPE xr1=pow(1.0,xrpow);
6671 if(r<Risco) nz=nzisco;
6675 FTYPE zisco = Risco*cos(th) ;
6677 FTYPE Sisco=1./(HSisco*HSisco*nzisco) ;
6679 FTYPE csisco = Hisco*nzisco ;
6684 rho = rho0 * (
MAX(0.0,xr-xr1)/xr) ;
6687 FTYPE rhopowisco=10.0;
6688 rho = (rho0isco * (
MAX(0.0,xr-xr1)/xr)) * pow(r/Risco,rhopowisco);
6690 uint = rho*cs*cs/(gamtorus - 1.) ;
6692 uPhi = 1./(pow(r,1.5) +
a) ;
6693 FTYPE uPhiisco = 1./(pow(Risco,1.5) +
a) ;
6698 uint *= (1.0 +
randfact * (ranc(0,0) - 0.5));
6701 FTYPE Vphiisco=uPhiisco;
6705 pt = uint * (gamtorus-1.0);
6708 if(rho<ppback[
RHO]){
6719 FTYPE sigma=r*r+
a*
a*cos(th)*cos(th);
6720 FTYPE gpp=sigma +
a*
a*(1+2*r/sigma)*sin(th)*sin(th);
6721 FTYPE sigmaisco=Risco*Risco+
a*
a*cos(th)*cos(th);
6722 FTYPE gppisco=sigmaisco +
a*
a*(1+2*Risco/sigmaisco)*sin(th)*sin(th);
6723 Vphi=Vphiisco*gppisco/gpp;
6726 if(1||(r<1.3*
Rhor || r<Risco) && usingback==0){
6751 if(*whichcoord!=mycoords){
6754 gset_X(getprim,mycoords,i,j,k,
NOWHERE,X,*ptrptrgeom);
6756 *whichcoord=mycoords;
6760 FTYPE ut=-1./sqrt(podpierd);
6773 FTYPE eps=(h-1.)/gamtorus;
6776 rho=pow(eps*(gamtorus-1.)/
RADNT_KKK,1./(gamtorus-1.));
6778 pt = uint * (gamtorus-1.0);
6780 uT=((*ptrptrgeom)->gcon[
GIND(0,0)])*ut+((*ptrptrgeom)->gcon[
GIND(0,3)])*uphi;
6781 uPhi=((*ptrptrgeom)->gcon[
GIND(3,3)])*uphi+((*ptrptrgeom)->gcon[
GIND(0,3)])*ut;
6789 if(rho<ppback[
RHO]){
6805 if(*whichcoord!=mycoords){
6808 gset_X(getprim,mycoords,i,j,k,
NOWHERE,X,*ptrptrgeom);
6810 *whichcoord=mycoords;
6818 FTYPE ll=l0*pow(Rcyl/r0,aa);
6820 FTYPE nn=1.0/(gamtorus-1.0);
6832 FTYPE phieff=phi + 0.5*pow(ll/Rcyl,2.0)/(1.0-aa);
6833 FTYPE phieffr0=phir0 + 0.5*pow(l0/r0,2.0)/(1.0-aa);
6834 FTYPE phieffeq=phieq + 0.5*pow(ll/Rcyl,2.0)/(1.0-aa);
6835 FTYPE rhot,rhotarg,rhoteq,rhotargeq;
6838 FTYPE epsilon0=1.45E-3;
6839 rhotarg=1.0 - (1.0/(gamtorus*epsilon0*epsilon0*phir0*phir0))*(phieff-phieffr0)/(nn+1.0);
6840 rhotargeq=1.0 - (1.0/(gamtorus*epsilon0*epsilon0*phir0*phir0))*(phieffeq-phieffr0)/(nn+1.0);
6841 rhot = rhoc*pow(rhotarg,nn);
6842 rhoteq = rhoc*pow(rhotargeq,nn);
6843 pt = rhoc*gamtorus*epsilon0*epsilon0*phir0*phir0*pow(rhot/rhoc,1.0+1.0/nn);
6851 FTYPE vs0=HoR0*vphi0;
6852 rhotarg=1.0 - (gamtorus/(vs0*vs0))*(phieff-phieffr0)/(nn+1.0);
6853 rhotargeq=1.0 - (gamtorus/(vs0*vs0))*(phieffeq-phieffr0)/(nn+1.0);
6854 rhot = rhoc*pow(rhotarg,nn);
6855 rhoteq = rhoc*pow(rhotargeq,nn);
6856 pt = rhoc*vs0*vs0/(gamtorus)*pow(rhot/rhoc,1.0+1.0/nn);
6857 FTYPE Eth0=vs0*vs0/(gamtorus*fabs(phir0));
6860 if(r<3.0 || rhotarg<0.0 || phieff<phieffr0 || rhot<0.0 || rhot>rhoc || phieffeq>0.0 || rhoteq>rhoc || rhoteq<0.0 ){
6865 uint=pt/(gamtorus-1.0);
6868 if(rho<ppback[
RHO]){
6879 dualfprintf(
fail_file,
"donutohsuga: r(%d)=%g th(%d)=%g Rcyl=%g : l0=%g ll=%g phi=%g phieff=%g phir0=%g phieffr0=%g rhot=%g rhotarg=%g pt=%g Vphi=%g\n",r,(*ptrptrgeom)->i,th,(*ptrptrgeom)->j,Rcyl,l0,ll,phi,phieff,phir0,phieffr0,rhot,rhotarg,pt,Vphi);
6895 pt = pressure_rho0_u_simple(i,j,k,
CENT,rho,uint);
6904 gset(getprim,*whichcoord,i,j,k,*ptrptrgeom);
6905 set_zamo_velocity(*whichvel,*ptrptrgeom,ppback);
6962 FTYPE gamtorus=4.0/3.0;
6964 *ptptr=(gamtorus-1.0)*pp[
UU];
6984 bl_coord_ijk_2(i,j,k,loc,X,V);
6985 struct of_geom geomrealdontuse;
6986 struct of_geom *ptrgeomrad=&geomrealdontuse;
6989 FTYPE *pr=&GLOBALMACP0A1(pglobal,i,j,k,0);
6991 FTYPE pradffortho[NPR];
6994 pradffortho[PRAD0]=pr[PRAD0];
6995 pradffortho[PRAD1]=pr[PRAD1];
6996 pradffortho[PRAD2]=pr[PRAD2];
6997 pradffortho[PRAD3]=pr[PRAD3];
6999 prad_fforlab(&whichvel, &whichcoord,
FF2LAB, i,j,k,loc,ptrgeomrad, pradffortho, pr, pr);
7014 (3.*
a*
a - 4.*
a*sqrt(R) + R*R)/
7015 pow(R*(
a + pow(R,1.5)),2)
7029 return(1. - POTENTIALorPRESSURE * exp((0.95*rin - R)*softer)) ;
7038 return(pow( 1.-pow(0.95*rin/R,0.5) , 1.)) ;
7042 int set_fieldtype(
void)
7095 fneg=1.0-pow(cos(th),thpower);
7096 fpos=1.0+pow(cos(th),thpower);
7097 gpara=0.5*(myr*fneg + 2.0*fpos*(1.0-log(fpos)));
7099 gpara=gpara-2.0*(1.0-log(2.0));
7109 FTYPE rshift,myr,rpower,myz,myR,myvert;
7110 FTYPE thother,thpower,gparalow,gparahigh,mygpara;
7119 myr=pow(r+rshift,rpower);
7122 myvert = (th>M_PI*0.5) ? (myr*sin(th)) : (myr*sin(-th));
7126 gparahigh=
setgpara(myr,thother,thpower);
7127 mygpara=(th<0.5*M_PI) ? gparalow : gparahigh;
7134 aphi=mygpara*cos(th);
7154 SFTYPE rho_av, p_av,u_av,q;
7159 int set_fieldtype(
void);
7160 int FIELDTYPE=set_fieldtype();
7187 int ibound=0,jbound=0,kbound=0;
7192 int iout=0,jout=0,kout=0;
7200 if(iout||jout||kout){
7203 else if(ibound||jbound||kbound){
7207 p_av = pressure_rho0_u_simple(i,j,k,
CENT,pr[
RHO],pr[
UU]);
7226 p_av = pressure_rho0_u_simple(i,j,k,loc,rho_av,u_av);
7229 else if(loc==
CORN2){
7232 p_av = pressure_rho0_u_simple(i,j,k,loc,rho_av,u_av);
7235 else if(loc==
CORN1){
7238 p_av = pressure_rho0_u_simple(i,j,k,loc,rho_av,u_av);
7242 dualfprintf(
fail_file,
"No such setup for loc=%d\n",loc);
7262 q = q*(p_av /
ptotmax - FRACAPHICUT);
7290 vpot += -(pow(r,rpow)*pow(sin(th),hpow)*sin(FIELDROT)*sin(ph));
7311 vpot += pow(r,rpow)*pow(sin(th),hpow)*(cos(FIELDROT) - cos(ph)*cot(th)*sin(FIELDROT));
7316 vpot += (1.0-cos(th));
7320 vpot +=
MAX(pow(r,4.0)*pow(rho_av,2.0)*1E40-0.02,0.0)*pow(sin(th),4.0);
7323 #define JONMADHPOW (4.0)
7324 #define JONMADR0 (0.0)
7325 #define JONMADROUT (300.0)
7335 if(V[2]<1
E-5 || V[2]>M_PI-1
E-5){
7349 q = rho_av /
rhomax - 0.2;
7351 if (q > 0.) vpot += q;
7355 q = rho_av /
rhomax - 0.2;
7357 if (q > 0.) vpot += q;
7377 q = 1
E-30*(p_av /
ptotmax*1E30 - FRACAPHICUT);
7379 if(rho_av/
rhomax-FRACAPHICUT<0) q=0;
7384 #define QPOWER (1.0)
7386 #define POWERNU (2.0) // 2.5 for toroidal field SUPERMAD paper
7390 FTYPE fact1,fact2,SSS,TTT;
7393 if(r<
rin) fact1=0.0;
7396 fact2=sin(log(r/SSS)/TTT);
7399 if (q > 0.) vpot += fact1*fact2;
7401 if(
PRODUCTION==0) dualfprintf(
fail_file,
"ijk=%d %d %d : ptotmax=%g p_av=%g q=%g %g %g rin=%g vpot=%g : rho_av=%g rhomax=%g\n",i,j,k,
ptotmax,p_av,q,fact1,fact2,
rin,vpot,rho_av,
rhomax);
7433 struct of_geom geomrealdontuse;
7434 struct of_geom *ptrgeomreal=&geomrealdontuse;
7436 gset_genloc(getprim,*whichcoord,i,j,k,loc,ptrgeomreal);
7438 lower_vec(mcon,ptrgeomreal,mcov);
7439 lower_vec(kcon,ptrgeomreal,kcov);
7443 FTYPE Rhorlocal=rhor_calc(0);
7464 FTYPE rV=Vmetric[1];
7465 FTYPE thV=Vmetric[2];
7466 FTYPE phV=Vmetric[3];
7475 delta= rV*rV - 2.0*
MBH*rV +
a*
a;
7476 sigma = rV*rV + a*a*cos(thV)*cos(thV);
7484 Acovblnonrot[0] = -1.*a*B0z + (a*B0z*
MBH*rV*(1. +
Power(
Cos(thV),2)))/sigma +
7487 Acovblnonrot[1] = -1.*B0x*(-1.*
MBH + rV)*
Cos(thV)*
Sin(psi)*
Sin(thV);
7494 Acovblnonrot[3] = -1.*B0x*
Cos(thV)*(delta*
Cos(psi) +
7504 DLOOPA(jj) Acovksnonrot[jj]=Acovblnonrot[jj];
7510 DLOOPA(jj) Acovks[jj]=Acovksnonrot[jj];
7542 static int fieldprim(
int whichmethod,
int whichinversion,
int *whichvel,
int*whichcoord,
int ii,
int jj,
int kk,
FTYPE *pr)
7552 int set_fieldtype(
void);
7553 int FIELDTYPE=set_fieldtype();
7575 struct of_geom *ptrgeom=&geomdontuse;
7586 gset(getprim,*whichcoord,ii,jj,kk,ptrgeom);
7599 FTYPE Rhorlocal=rhor_calc(0);
7603 if(th<M_PI*0.5) pr[
B1]=
B0WALD/ptrgeom->gdet;
7614 int doit=ii==
N1-1 && jj==
N2/4-1 && kk==0;
7629 extern void lower_A(
FTYPE (*Acon)[NDIM],
struct of_geom *geom,
FTYPE (*Acov)[NDIM]);
7637 PLOOP(pliter,pl) prold[pl]=pr[pl];
7653 MtoF(3,Fcov,ptrgeom,Mcon);
7657 DLOOP(j,k) Mcon[
j][k]=0.0;
7658 SLOOPA(j) Mcon[
j][0] = Bcon[
j] = prold[B1+j-1];
7659 SLOOPA(j) Mcon[0][
j] = -Mcon[
j][0];
7665 FTYPE Ed1=Mcon[2][3];
7666 FTYPE Ed2=Mcon[3][1];
7667 FTYPE Ed3=Mcon[1][2];
7679 Bcon[3] = Mcon[3][0] = 0.0;
7682 Bcon[3] = Mcon[3][0] =(Bu1*gv01*gv03 + Bu2*gv02*gv03 - Bu1*gv00*gv13 - Ed3*gv02*gv13 +
7683 Ed2*gv03*gv13 - Bu2*gv00*gv23 + Ed3*gv01*gv23 - Ed1*gv03*gv23 -
7684 Ed2*gv01*gv33 + Ed1*gv02*gv33 + myBd3)/denom;
7689 Mcon[0][3] = -Mcon[3][0];
7694 MtoF(0,Mcon,ptrgeom,Fcov);
7698 indices_2221(Fcon,Fud,ptrgeom);
7699 indices_2212(Fcon,Fdu,ptrgeom);
7701 indices_2221(Mcon,Mud,ptrgeom);
7702 indices_2212(Mcon,Mdu,ptrgeom);
7703 lower_A(Mcon,ptrgeom,Mcov);
7715 DLOOP(j,k) Fsq += Fcon[
j][k]*Fcov[
j][k];
7716 DLOOP(j,k) Tud[
j][k] = 0.0;
7718 DLOOPA(ll) Tud[
j][k] += Fud[
j][ll]*Fdu[k][ll];
7719 Tud[
j][k] += - 0.25*
delta(j,k)*Fsq;
7721 if(doit)
DLOOP(j,k) dualfprintf(
fail_file,
"Tud[%d][%d]=%g\n",j,k,Tud[j][k]);
7727 alpha = 1./sqrt(-ptrgeom->gcon[
GIND(0,0)]);
7734 raise_vec(etacov,ptrgeom,etacon);
7742 DLOOP(j,k) Betacon[
j]+=etacov[k]*Mcon[k][
j];
7743 lower_vec(Betacon,ptrgeom,Betacov);
7752 DLOOP(j,k) Eetacov[
j]+=etacon[k]*Fcov[
j][k];
7753 raise_vec(Eetacov,ptrgeom,Eetacon);
7754 if(doit) dualfprintf(
fail_file,
"etacon=%g %g %g %g Eetacov[3]=%21.15g Fcov03=%g\n",etacon[0],etacon[1],etacon[2],etacon[3],Eetacov[3],Fcov[0][3]);
7760 EBtoT(Eetacon,Betacon,ptrgeom,TfromEB);
7761 if(doit)
DLOOP(j,k) dualfprintf(
fail_file,
"TudfromEB[%d][%d]=%g\n",j,k,TfromEB[j][k]);
7765 if(whichinversion==0){
7769 SLOOPA(j) U[B1+j-1] = Bcon[
j] = Mcon[
j][0];
7773 if(doit) dualfprintf(
fail_file,
"CONS: %g %g %g %g : %g %g %g: F=%g %g\n",U[
UU],U[
U1],U[
U2],U[
U3],Bcon[1],Bcon[2],Bcon[3],Fcov[0][3],Fcov[3][0]);
7775 struct of_newtonstats newtonstats; setnewtonstatsdefault(&newtonstats);
7784 if(doit) dualfprintf(
fail_file,
"BEFORE Utoprimgen()\n");
7785 if(doit)
PLOOP(pliter,pl) dualfprintf(
fail_file,
"oldpr[%d]=%g\n",pl,pr[pl]);
7786 Utoprimgen(0,0,0,0,1,&eomtypelocal,
CAPTYPEBASIC,0,1,
EVOLVEUTOPRIM,
UNOTHING,U,qptr,ptrgeom,0,pr,pr,&newtonstats);
7787 PLOOP(pliter,pl)
if(pl>=URAD1 && pl<=URAD3) pr[pl]=prold[pl];
7788 if(doit) dualfprintf(
fail_file,
"AFTER Utoprim()\n");
7789 if(doit)
PLOOP(pliter,pl) dualfprintf(
fail_file,
"newpr[%d]=%g\n",pl,pr[pl]);
7793 if(pl==
RHO || pl==UU || pl>
B3 || pl>=URAD1 && pl<=URAD3) pr[pl]=prold[pl];
7799 FTYPE Ueomtype[NPR],Ueomtypeabs[NPR];
7802 Utoprimgen(0,0,0,0,1,&eomtypelocal,
CAPTYPEBASIC,0,1,
EVOLVEUTOPRIM,
UNOTHING,Ueomtype,qptr4,ptrgeom,0,pr,pr,&newtonstats);
7804 if(doit)
PLOOP(pliter,pl) dualfprintf(
fail_file,
"newprMHD[%d]=%g\n",pl,pr[pl]);
7811 else if(whichinversion==1){
7819 if(doit) dualfprintf(
fail_file,
"BEFORE EBtopr()\n");
7820 if(doit) dualfprintf(
fail_file,
"%g %g %g %g : %g %g %g %g\n",Eetacon[0],Eetacon[1],Eetacon[2],Eetacon[3],Betacon[0],Betacon[1],Betacon[2],Betacon[3]);
7821 EBtopr(Eetacon,Betacon,ptrgeom,pr);
7822 if(doit) dualfprintf(
fail_file,
"AFTER EBtopr()\n");
7823 if(doit)
PLOOP(pliter,pl) dualfprintf(
fail_file,
"newpr[%d]=%g\n",pl,pr[pl]);
7828 PLOOP(pliter,pl)
if(pl>=URAD1 && pl<=URAD3) pr[pl]=prold[pl];
7842 fdd03 = ptrgeom->gdet * (qptr2->
ucon[1]*qptr2->
bcon[2] - qptr2->
ucon[2]*qptr2->
bcon[1]) ;
7843 if(doit) dualfprintf(
fail_file,
"fdd03=%g\n",fdd03);
7847 if(doit) dualfprintf(
fail_file,
"mhdflux=%g %g %g\n",mhdflux[1][0],mhdflux[2][0],mhdflux[3][0]);
7850 if(doit) dualfprintf(
fail_file,
"mhdfluxma=%g %g %g\n",mhdfluxma[1][0],mhdfluxma[2][0],mhdfluxma[3][0]);
7853 if(doit) dualfprintf(
fail_file,
"mhdfluxem=%g %g %g\n",mhdfluxem[1][0],mhdfluxem[2][0],mhdfluxem[3][0]);
7860 int return1=pr2ucon(
WHICHVEL,pr,ptrgeom,ucontemp);
7863 for(pl=
RHO;pl<=
U3;pl++) pr[pl]=prold[pl];
7868 ucon2pr(*whichvel,ucontemp,ptrgeom,pr);
7873 for(pl=
RHO;pl<=
U3;pl++) pr[pl]=prold[pl] + (pr[pl]-prold[pl])/(
Rhor-0.0)*(r-0.0);
7894 FTYPE xx=r*sin(th)*cos(ph);
7895 FTYPE yy=r*sin(th)*sin(ph);
7900 lambdatrans[
TT][
TT]=1.0;
7901 SLOOPA(j) lambdatrans[
TT][
j] = lambdatrans[
j][
TT] = 0.0;
7903 ilambdatrans[
TT][
TT]=1.0;
7904 SLOOPA(j) ilambdatrans[
TT][
j] = ilambdatrans[
j][
TT] = 0.0;
7908 lambdatrans[1][
RR] = sin(th)*cos(ph);
7909 lambdatrans[1][
TH] = cos(th)*cos(ph);
7910 lambdatrans[1][
PH] = -sin(ph);
7912 lambdatrans[2][
RR] = sin(th)*sin(ph);
7913 lambdatrans[2][
TH] = cos(th)*sin(ph);
7914 lambdatrans[2][
PH] = cos(ph);
7916 lambdatrans[3][
RR] = cos(th);
7917 lambdatrans[3][
TH] = -sin(th);
7918 lambdatrans[3][
PH] = 0.0;
7921 ilambdatrans[1][
RR] = sin(th)*cos(ph);
7922 ilambdatrans[1][
TH] = cos(th)*cos(ph);
7923 ilambdatrans[1][
PH] = -sin(ph);
7925 ilambdatrans[2][
RR] = sin(th)*sin(ph);
7926 ilambdatrans[2][
TH] = cos(th)*sin(ph);
7927 ilambdatrans[2][
PH] = cos(ph);
7929 ilambdatrans[3][
RR] = cos(th);
7930 ilambdatrans[3][
TH] = -sin(th);
7931 ilambdatrans[3][
PH] = 0.0;
7945 tempcomp[k] += ilambdatrans[
j][k]*finalvec[
j];
7952 finalvec[
RR]=pr[
U1] + finalvec[1]/sqrt(ptrgeom->gcov[
GIND(1,1)]);
7953 finalvec[
TH]=pr[
U2] + finalvec[2]/sqrt(ptrgeom->gcov[
GIND(2,2)]);
7954 finalvec[
PH]=pr[
U3] + finalvec[3]/sqrt(ptrgeom->gcov[
GIND(3,3)]);
7956 pr[
U1] = finalvec[
RR];
7957 pr[
U2] = finalvec[
TH];
7958 pr[
U3] = finalvec[
PH];
7964 for(k=
U1;k<=
U1+3;k++){
7983 pr[
B1]=Econ[1]*geom.g/etacov[
TT];
7997 if(fabs(faradaytest[j][k]-Fcov[j][k])>1E-10){
7998 dualfprintf(
fail_file,
"1 %d %d : %21.15g %21.15g\n",ii,jj,faradaytest[j][k],Fcov[j][k]);
8001 if(fabs(faradaytest[0][3])>1E-10) dualfprintf(
fail_file,
"1 Fcov=%21.15g faraday=%21.15g\n",Fcov[0][3],faradaytest[0][3]);
8013 int set_fieldtype(
void);
8014 int FIELDTYPE=set_fieldtype();
8021 int whichinversion=0;
8022 fieldprim(whichmethod, whichinversion, &whichvel, &whichcoord, i, j, k, pr);
8052 int i,
j,k,pliter,pl;
8054 struct of_geom *ptrgeom=&geomdontuse;
8060 FTYPE Bu1,Bu2,gcon03,gcon13,gcon23,gcon33;
8061 FTYPE gcov01,gcov02,gcov11,gcov12,gcov21,gcov22,gcov03,gcov13,gcov23;
8066 gcon03=ptrgeom->gcon[
GIND(0,3)];
8067 gcon13=ptrgeom->gcon[
GIND(1,3)];
8068 gcon23=ptrgeom->gcon[
GIND(2,3)];
8069 gcon33=ptrgeom->gcon[
GIND(3,3)];
8071 gcov01=ptrgeom->gcov[
GIND(0,1)];
8072 gcov02=ptrgeom->gcov[
GIND(0,2)];
8073 gcov11=ptrgeom->gcov[
GIND(1,1)];
8074 gcov12=gcov21=ptrgeom->gcov[
GIND(1,2)];
8075 gcov22=ptrgeom->gcov[
GIND(2,2)];
8076 gcov03=ptrgeom->gcov[
GIND(0,3)];
8077 gcov13=ptrgeom->gcov[
GIND(1,3)];
8078 gcov23=ptrgeom->gcov[
GIND(2,3)];
8082 FTYPE ftemp=(1.0 - gcon03*gcov03 - gcon13*gcov13 - gcon23*gcov23);
8084 pl=
B3; pr[pl] = (myBd3*gcon33 + Bu1*gcon03*gcov01 + Bu2*gcon03*gcov02 + Bu1*gcon13*gcov11 + Bu2*gcon13*gcov12 + Bu1*gcon23*gcov21 + Bu2*gcon23*gcov22)*igdetnosing;
8097 if(r>3.0) ucon[1]=ucon[2]=ucon[3]=0.0;
8102 extern void filterffde(
int i,
int j,
int k,
FTYPE *pr);
8105 filterffde(i,j,k,
GLOBALMAC(pglobal,i,j,k));
8112 #define FCOVDERTYPE DIFFGAMMIE
8115 #if((REALTYPE==DOUBLETYPE)||(REALTYPE==FLOATTYPE))
8116 #define FCOVDXDELTA 1E-5
8117 #elif(REALTYPE==LONGDOUBLETYPE)
8118 #define FCOVDXDELTA 1E-6
8128 FTYPE mcovhj,mcovlj,kcovhj,kcovlj;
8129 FTYPE mcovhk,mcovlk,kcovhk,kcovlk;
8164 for(k=0;k<
NDIM;k++){
8165 for(j=0;j<
NDIM;j++){
8167 for(l=0;l<
NDIM;l++) Xlk[l]=Xhk[l]=Xlj[l]=Xhj[l]=X[l];
8222 +(mcovhj - mcovlj) / (Vhk[k] - Vlk[k])
8224 +(kcovhj - kcovlj) / (Vhk[k] - Vlk[k])
8230 -(mcovhk - mcovlk) / (Vhj[j] - Vlj[j])
8232 -(kcovhk - kcovlk) / (Vhj[j] - Vlj[j])
8243 for(k=0;k<
NDIM;k++)
for(j=0;j<
NDIM;j++){
8247 FTYPE ans1; dfridr(mcov_func_mcoord,ptrgeom,X,0,j,k,&ans1);
8248 FTYPE ans2; dfridr(mcov_func_mcoord,ptrgeom,X,0,k,j,&ans2);
8249 FTYPE ans3; dfridr(kcov_func_mcoord,ptrgeom,X,0,j,k,&ans3);
8250 FTYPE ans4; dfridr(kcov_func_mcoord,ptrgeom,X,0,k,j,&ans4);
8259 dualfprintf(
fail_file,
"NOT SETUP FOR NUMREC\n");
8278 gset_X(getprim, whichcoord, ptrgeom->i, ptrgeom->j, ptrgeom->k, ptrgeom->p, X, ptrgeom);
8288 lower_vec(mcon,ptrgeom,mcov);
8302 gset_X(getprim, whichcoord, ptrgeom->i, ptrgeom->j, ptrgeom->k, ptrgeom->p, X, ptrgeom);
8310 lower_vec(kcon,ptrgeom,kcov);
8339 for(k=0;k<
NDIM;k++){
8342 for(j=0;j<
NDIM;j++){
8344 for(l=0;l<
NDIM;l++) Xl[l]=Xh[l]=X[l];
8366 Jcon[k] += (Fconh - Fconl) / (Vh[j] - Vl[j]) ;
8373 for(k=0;k<
NDIM;k++){
8375 for(j=0;j<
NDIM;j++){
8376 FTYPE ans1; dfridr(Fcon_func_mcoord,ptrgeom,X,0,j,k,&ans1);
8393 DLOOP(j,k) idxdxp[
j][k]=0.0;
8402 Jcon[k] += Jconorig[
j] * idxdxp[k][
j];
8424 gset_X(getprim, whichcoord, ptrgeom->i, ptrgeom->j, ptrgeom->k, ptrgeom->p, X, ptrgeom);
8433 raise_A(Fcov, ptrgeom, Fcon) ;
8436 return(ptrgeom->gdet*Fcon[ii][jj]);
8451 int fieldfrompotential[
NDIM];
8454 if(funreturn!=0)
return(funreturn);
8493 if(funreturn!=0)
return(funreturn);
8506 int set_fieldtype(
void);
8507 int FIELDTYPE=set_fieldtype();
8527 funreturn=
user1_get_maxes(eqslice, parms,prim, bsq_max, ptot_max, beta_min);
8528 if(funreturn!=0)
return(funreturn);
8539 int set_fieldtype(
void);
8540 int FIELDTYPE=set_fieldtype();
8543 dualfprintf(
fail_file,
"DID NORM FIELD\n");
8546 if(funreturn!=0)
return(funreturn);
8565 #define TAUADJUSTATM (10.0) // timescale for boundary to adjust to using preset inflow
8566 int set_atmosphere(
int whichcond,
int whichvel,
struct of_geom *ptrgeom,
FTYPE *pr)
8573 funreturn=user1_set_atmosphere(atmospheretype, whichcond, whichvel, ptrgeom, pr);
8574 if(funreturn!=0)
return(funreturn);
8596 prceiling[PRAD0]=
BIG;
8603 funreturn=set_density_floors_default(ptrgeom, pr, prfloor, prceiling);
8608 bl_coord_ijk(ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p,V);
8613 if(pr[
RHO]<rholimit) pr[
RHO]=rholimit;
8615 if(funreturn!=0)
return(funreturn);
8635 prceiling[PRAD0]=
BIG;
8642 funreturn=set_density_floors_default_alt(ptrgeom, q, pr, U, bsq, prfloor, prceiling);
8645 bl_coord_ijk(ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p,V);
8650 if(pr[
RHO]<rholimit) pr[
RHO]=rholimit;
8654 if(funreturn!=0)
return(funreturn);
8666 int theproblem_set_enerregiondef(
int forceupdate,
int timeorder,
int numtimeorders,
long int thenstep,
FTYPE thetime,
int (*enerregiondef)[
NDIM] )
8709 int theproblem_set_enerregionupdate(
int forceupdate,
int timeorder,
int numtimeorders,
long int thenstep,
FTYPE thetime,
int *updateeverynumsteps,
int *everynumsteps)
8725 *updateeverynumsteps=1;
8729 *updateeverynumsteps=1;
8737 *everynumsteps = *updateeverynumsteps*100;
8745 int theproblem_set_myid(
void)
8783 set_boundloop(boundvartype, inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi, &riin, &riout, &rjin, &rjout, &rkin, &rkout, dosetbc);
8799 struct
of_geom *ptrgeom=&geomdontuse;
8812 FTYPE interptime=(fluxtime<1.0 ? (fluxtime/t0) : 1.0);
8813 if(interptime>1.0) interptime=1.0;
8814 if(interptime<0.0) interptime=0.0;
8839 prflux[URAD0]=Ehatstar*(1.0-interptime) + interptime*Ehatstar*10.0;
8842 FTYPE kappa=calc_kappaes_user(rhostar,1.0,0,0,0);
8843 vr0=sqrt(1.0/(3.0*kappa*t0));
8844 static
int firsttime=1;
8845 if(firsttime) dualfprintf(
fail_file,"vr0=%21.15g\n",vr0);
8847 prflux[URAD1]=vr0/sqrt(fabs(ptrgeom->
gcov[
GIND(1,1)]));
8873 MACP0A1(F1,i,j,k,URAD1)=fluxrad[URAD1];
8874 MACP0A1(F1,i,j,k,URAD0)=fluxrad[URAD0];
8893 int inboundloop[
NDIM];
8894 int outboundloop[
NDIM];
8895 int innormalloop[
NDIM];
8896 int outnormalloop[
NDIM];
8908 set_boundloop(boundvartype, inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi, &riin, &riout, &rjin, &rjout, &rkin, &rkout, dosetbc);
8924 struct
of_geom *ptrgeom=&geomdontuse;
8938 FTYPE interptime=(fluxtime<1.0 ? (fluxtime/t0) : 1.0);
8939 if(interptime>1.0) interptime=1.0;
8940 if(interptime<0.0) interptime=0.0;
8956 primtoflux(UEVOLVE, prflux, &q, 2, ptrgeom, flux, NULL);
8976 struct of_geom *ptrgeom=&geomdontuse;
9009 else ujet=ustar*1
E-2;
9013 FTYPE rhojet; rhojet=rhostar*1E-2;
9021 FTYPE gammacore; gammacore = 7.0;
9022 FTYPE uzcore; uzcore=sqrt(gammacore*gammacore-1.0);
9023 FTYPE rhocore; rhocore=rhojet;
9025 FTYPE efluxcore; efluxcore=rhocore*uzcore*uzcore;
9029 FTYPE uz=uzcore*stepfunctionab2(V[1],1,0);
9033 FTYPE rho; rho = eflux/(uz*uz);
9034 if(rho>rhostar) rho=rhostar;
9040 FTYPE Tjet = 1.62*Tstar;
9042 Temp = stepfunctionab2(V[1],Tjet,Tstar);
9044 ug = Temp*rho/(
gam-1.0);
9055 FTYPE Tgas=compute_temp_simple(i,j,k,loc,prflux[RHO],prflux[UU]);
9058 prflux[
U2]=uz/sqrt(fabs(ptrgeom->gcov[GIND(2,2)]));
9063 prflux[URAD2]=prflux[
U2];
9070 static int firsttime=1;
9074 FTYPE ptot=((
gam-1.0)*prflux[
UU] + (4.0/3.0-1.0)*prflux[URAD0]);
9078 FTYPE gamptot=(
gam*(
gam-1.0)*prflux[
UU] + (4.0/3.0)*(4.0/3.0-1.0)*prflux[URAD0]);
9080 dualfprintf(fail_file,
"JET: ptrue=%21.15g ptot=%21.15g cs2=%21.15g cs=%21.15g\n",ptrue,ptot,cs2,sqrt(cs2));
9082 dualfprintf(fail_file,
"Tgasinj=%21.15g u/rho=%21.15g prad0/rho=%21.15g\n",Tgas,prflux[UU]/prflux[RHO],prflux[URAD0]/prflux[RHO]);
9083 dualfprintf(fail_file,
"TEST: %21.15g\n",0.99*
stepfunctionab(1E-2,1.0,0.0));
9084 dualfprintf(fail_file,
"TESTUCON0: %21.15g %21.15g %21.15g %21.15g\n",ucon[
TT],ucon[1]*sqrt(ptrgeom->gcov[GIND(1,1)]),ucon[2]*sqrt(ptrgeom->gcov[GIND(2,2)]),ucon[3]*sqrt(ptrgeom->gcov[GIND(3,3)]));
9085 dualfprintf(fail_file,
"Need: pradstar/rhojet=%21.15g > %21.15g and taujet=%21.15g >>1\n",pradstar/prflux[RHO],ucon[TT]*ucon[TT]*calc_kappaes_user(prflux[RHO],Temp,0,0,0)*width*width/(2.0*
Rout_array[2]),calc_kappaes_user(prflux[RHO],Temp,0,0,0)*width);
9087 FTYPE ljet,ljet1,ljet2,ljet3;
9088 ljet1 = prflux[
RHO]*ucon[
TT]*ucon[
TT]*M_PI*width*width;
9089 ljet2 = (prflux[
UU] + (
gam-1.0)*prflux[
UU])*ucon[TT]*ucon[TT]*M_PI*width*width;
9090 ljet3 = (prflux[URAD0] + (4.0/3.0-1.0)*prflux[URAD0])*ucon[TT]*ucon[TT]*M_PI*width*width;
9091 ljet=ljet1+ljet2+ljet3;
9094 FTYPE lrad,lradalt1,lradalt2;
9095 lrad = (pradstar/calc_kappaes_user(prflux[RHO],Temp,0,0,0))*(2.0*M_PI*
Rout_array[2]);
9096 lradalt1 = (pradstar/calc_kappaes_user(rhostar,Tstar,0,0,0))*(2.0*M_PI*Rout_array[2]);
9097 lradalt2 = (
ARAD_CODE*pow(Tstar,4.0)/4.0)*(2.0*M_PI*width*Rout_array[2]);
9099 dualfprintf(fail_file,
"Need: ljet=%21.15g < lrad=%21.15g lradalt1=%21.15g lradalt2=%21.15g\n",ljet,lrad,lradalt1,lradalt2);
9100 dualfprintf(fail_file,
"ljet1=%21.15g ljet2=%21.15g ljet3=%21.15g\n",ljet1,ljet2,ljet3);
9133 if(mycpupos[1]==ncpux1-1){
9134 dualfprintf(fail_file,
"did modify\n");
9137 MACP1A1(fluxvec,2,
N1,j,k,B1)=
MACP1A1(fluxvec,1,
N1,j,k,
B2)=0.0;
9206 #define KAPPA_ES_FERMICORR(rhocode,Tcode) (1.0/(1.0+2.7E11*((rhocode)*RHOBAR)/prpow((Tcode)*TEMPBAR,2.0))) // Buchler and Yueh 1976 (just Fermi part). Fewer electrons when near Fermi fluid limit.
9207 #define KAPPA_ES_KNCORR(rhocode,Tcode) (1.0/(1.0+prpow((Tcode)*TEMPBAR/4.5E8,0.86))) // Buchler and Yueh 1976 . Klein-Nishina for thermal electrons.
9209 #define KAPPA_ES_CODE(rhocode,Tcode) (0.2*(1.0+XFACT)*KAPPA_ES_FERMICORR(rhocode,Tcode)*KAPPA_ES_KNCORR(rhocode,Tcode)/OPACITYBAR)
9210 #define KAPPA_ES_BASIC_CODE(rhocode,Tcode) (0.2*(1.0+XFACT)/OPACITYBAR)
9220 #define KAPPA_FORCOMPT_RELCORREP(ep) ((1.0 + 3.683*(ep)+4.0*(ep)*(ep))/(1.0 + (ep))) // Sadowski et al. (2014) Eq 26 and 27.
9221 #define KAPPA_FORCOMPT_RELCORR(rhocode,Tcode) (KAPPA_FORCOMPT_RELCORREP(K_BOLTZ*(Tcode)*TEMPBAR/(MELE*CCCTRUE*CCCTRUE)))
9222 #define KAPPA_FORCOMPT_CODE(rhocode,Tcode) (0.2*(1.0+XFACT)*KAPPA_ES_FERMICORR(rhocode,Tcode)*KAPPA_FORCOMPT_RELCORR(rhocode,Tcode)/OPACITYBAR)
9234 #define KAPPA_ZETA(Tgcode,Trcode) ((TEMPMIN+Trcode)/(TEMPMIN+Tgcode))
9236 #define KAPPA_FF_CODE(rhocode,Tgcode,Trcode) (4.0E22*(1.0+XFACT)*(1.0-ZFACT)*((rhocode)*RHOBAR)*prpow((Tgcode)*TEMPBAR,-0.5)*prpow((Trcode)*TEMPBAR,-3.0)*prlog(1.0+1.6*KAPPA_ZETA(Tgcode,Trcode))*(1.0+4.4E-10*(Tgcode*TEMPBAR))/OPACITYBAR) // ASSUMPTION: Thermal ele and no pairs. See Rybicki & Lightman Eq~5.25 and McKinney & Uzdensky (2012) . For Tr,Tg split, see Ramesh notes.
9254 #define KAPPAN_FF_CODE(rhocode,Tgcode,Trcode) KAPPA_FF_CODE(rhocode,Tgcode,Trcode)
9258 #define KAPPA_BF_CODE(rhocode,Tgcode,Trcode) (3.0E25*(ZFACT)*(1.0+XFACT+0.75*YFACT)*((rhocode)*RHOBAR)*prpow((Tgcode)*TEMPBAR,-0.5)*prpow((Trcode)*TEMPBAR,-3.0)*prlog(1.0+1.6*KAPPA_ZETA(Tgcode,Trcode))/OPACITYBAR) // ASSUMPTION: Number of electrons similar to for solar abundances for 1+X+(3/4)Y term. For Tr,Tg split, see Ramesh notes.
9259 #define KAPPA_CHIANTIBF_CODE(rhocode,Tgcode,Trcode) (4.0E34*((rhocode*RHOBAR))*(ZFACT/ZSOLAR)*YELE*prpow((Tgcode)*TEMPBAR,-1.7)*prpow((Trcode)*TEMPBAR,-3.0)/OPACITYBAR) // *XFACT literally from Fig 34.1 in Draine book, but for solar n_H\sim n_b\sim 1/cm^3 only
9260 #define KAPPA_HN_CODE(rhocode,Tgcode,Trcode) (1.1E-25*prpow(ZFACT,0.5)*prpow((rhocode)*RHOBAR,0.5)*prpow((Tgcode)*TEMPBAR,7.7)/OPACITYBAR) // other sources cite 2.5E-31 (Z/0.02)(rho)^(1/2)(T)^9
9261 #define KAPPA_MOL_CODE(rhocode,Tgcode,Trcode) (0.1*ZFACT/OPACITYBAR)
9263 #define KAPPA_GENFF_CODE(rhocode,Tgcode,Trcode) (1.0/(1.0/(KAPPA_MOL_CODE(rhocode,Tgcode,Trcode)+KAPPA_HN_CODE(rhocode,Tgcode,Trcode)) + 1.0/(KAPPA_CHIANTIBF_CODE(rhocode,Tgcode,Trcode)+KAPPA_BF_CODE(rhocode,Tgcode,Trcode)+KAPPA_FF_CODE(rhocode,Tgcode,Trcode)))) // for 1.3E3K \le T \le 1E9K or higher. Numerically better to have kappa bottom out at low T so no diverent opacity as T->0
9300 #if(WHICHPROBLEM==FLATNESS)
9306 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA*KAPPA_FF_CODE(rho,Tg,Tr))
9308 #define KAPPAESUSER(rho,T) (rho*KAPPAES*KAPPA_ES_BASIC_CODE(rho,T))
9317 #if(WHICHPROBLEM==RADPULSE || WHICHPROBLEM==RADPULSEPLANAR || WHICHPROBLEM==RADPULSE3D)
9322 #if(WHICHPROBLEM==RADPULSEPLANAR)
9341 #define KAPPAES (1E3) // takes VERY long time with sub-cycling, but works.
9348 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA)
9349 #define KAPPAESUSER(rho,T) (rho*KAPPAES)
9352 #else // PULSE and PULSE3D
9356 #define KAPPAES (SMALL)
9359 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA*KAPPA_FF_CODE(rho,Tg,Tr))
9361 #define KAPPAESUSER(rho,T) (rho*KAPPAES*KAPPA_ES_BASIC_CODE(rho,T))
9374 #if(WHICHPROBLEM==RADBEAMFLAT)
9381 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA*KAPPA_FF_CODE(rho,Tg,Tr))
9383 #define KAPPAESUSER(rho,T) (rho*KAPPAES*KAPPA_ES_BASIC_CODE(rho,T))
9392 #if(WHICHPROBLEM==RADTUBE)
9394 #define KAPPAESUSER(rho,T) (0.0)
9397 #define KAPPAUSER(rho,B,Tg,Tr) (0.4*rho)
9399 #define KAPPAUSER(rho,B,Tg,Tr) (0.2*rho)
9401 #define KAPPAUSER(rho,B,Tg,Tr) (0.3*rho)
9403 #define KAPPAUSER(rho,B,Tg,Tr) (25*rho)
9405 #define KAPPAUSER(rho,B,Tg,Tr) (0.08*rho)
9407 #define KAPPAUSER(rho,B,Tg,Tr) (0.7*rho)
9409 #define KAPPAUSER(rho,B,Tg,Tr) (1000*rho)
9419 #if(WHICHPROBLEM==RADSHADOW || WHICHPROBLEM==RADDBLSHADOW)
9422 #define KAPPAUSER(rho,B,Tg,Tr) (rho*1.0) // paper
9423 #define KAPPAESUSER(rho,T) (rho*0.0)
9437 #if(WHICHPROBLEM==RADBEAM2D || WHICHPROBLEM==RADBEAM2DKS || WHICHPROBLEM==RADBEAM2DKSVERT)
9444 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA*KAPPA_FF_CODE(rho,Tg,Tr))
9446 #define KAPPAESUSER(rho,T) (rho*KAPPAES*KAPPA_ES_BASIC_CODE(rho,T))
9452 #if(WHICHPROBLEM==ATMSTATIC)
9459 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA*KAPPA_FF_CODE(rho,Tg,Tr))
9461 #define KAPPAESUSER(rho,T) (rho*KAPPAES*KAPPA_ES_BASIC_CODE(rho,T))
9467 #if(WHICHPROBLEM==RADATM)
9471 #define KAPPAES 1. // only scattering
9474 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA*KAPPA_FF_CODE(rho,Tg,Tr))
9476 #define KAPPAESUSER(rho,T) (rho*KAPPAES*KAPPA_ES_BASIC_CODE(rho,T))
9482 #if(WHICHPROBLEM==RADWALL)
9489 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA*KAPPA_FF_CODE(rho,Tg,Tr))
9491 #define KAPPAESUSER(rho,T) (rho*KAPPAES*KAPPA_ES_BASIC_CODE(rho,T))
9499 #if(WHICHPROBLEM==RADWAVE)
9501 #define KAPPAUSER(rho,B,Tg,Tr) (rho*RADWAVE_KAPPA)
9502 #define KAPPAESUSER(rho,T) (rho*RADWAVE_KAPPAES)
9508 #if(WHICHPROBLEM==KOMIPROBLEM)
9510 #define KAPPAUSER(rho,B,Tg,Tr) (0.)
9511 #define KAPPAESUSER(rho,T) (0.)
9516 #if(WHICHPROBLEM==RADBONDI)
9523 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA*KAPPA_FF_CODE(rho,Tg,Tr))
9525 #define KAPPAESUSER(rho,T) (rho*KAPPAES*KAPPA_ES_BASIC_CODE(rho,T))
9531 #if(WHICHPROBLEM==RADDOT)
9538 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA*KAPPA_FF_CODE(rho,Tg,Tr))
9540 #define KAPPAESUSER(rho,T) (rho*KAPPAES*KAPPA_ES_BASIC_CODE(rho,T))
9545 #if(WHICHPROBLEM==RADNT || WHICHPROBLEM==RADFLATDISK)
9547 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA_ES_CODE(rho,Tg)/1E14*0.1) // wierd use of kappa_{es} in koral
9549 #define KAPPAESUSER(rho,T) (0.0)
9553 #if(WHICHPROBLEM==RADDONUT)
9567 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA*(KAPPA_GENFF_CODE(SMALL+rho,Tg+TEMPMIN,Tr+TEMPMIN)))
9570 #define KAPPANUSER(rho,B,Tg,Tr) (rho*KAPPA*(KAPPA_GENFF_CODE(SMALL+rho,Tg+TEMPMIN,Tr+TEMPMIN)))
9573 #define KAPPAESUSER(rho,T) (rho*KAPPAES*KAPPA_ES_CODE(rho,T))
9578 #if(WHICHPROBLEM==RADCYLBEAM || WHICHPROBLEM==RADCYLBEAMCART)
9580 #define KAPPAUSER(rho,B,Tg,Tr) (rho*KAPPA_ES_BASIC_CODE(rho,Tg)/1E14*0.0) // note 0.0
9581 #define KAPPAESUSER(rho,T) (0.0)
9585 #if(WHICHPROBLEM==RADCYLJET)
9587 #define KAPPAUSER(rho,B,Tg,Tr) (rho*(KAPPA_FF_CODE(SMALL+rho,Tg+TEMPMIN,Tr+TEMPMIN))) // SMALL
9588 #define KAPPAESUSER(rho,Tg) (rho*KAPPA_ES_BASIC_CODE(rho,Tg)/100.0)
9595 #define KAPPANUSER(rho,B,Tg,Tr) KAPPAUSER(rho,B,Tg,Tr)
9605 #define ISKAPPAEABS 0
9606 #define ISKAPPANABS 1
9607 #define ISKAPPAEEMIT 2
9608 #define ISKAPPANEMIT 3
9611 #define E 2.718281828459045
9632 FTYPE kappa,kappaemit,kappan,kappanemit,kappaes;
9638 else if(which==
ISKAPPAES)
return(kappaes);
9652 #if(WHICHFIT==ISFITORIG)
9653 dualfprintf(fail_file,
"Shouldn't be here\n");
9662 FTYPE nereal=3.0110683499999995e23*rhoreal*(1.0 +
XFACT);
9667 FTYPE xi = Trreal/Tereal;
9670 FTYPE thetaesq=thetae*thetae;
9671 FTYPE thetaecubed=thetaesq*thetae;
9681 FTYPE kappaffreal=1.2E24*pow(Tereal,-7.0/2.0)*pow(rhoreal,2.0)*(1.0+
XFACT)*(1.0-
ZFACT);
9684 FTYPE Reilow = 1.0 + 1.7644624334086192*
Power(thetae,1.34);
9685 FTYPE Reihigh = 1.4019447514099515*
Power(thetae,0.5)*(1.5 +
Log(0.48 + 1.123*thetae));
9686 FTYPE Rei = (thetae>=1.0 ? Reihigh : Reilow);
9691 FTYPE kappanffreal=kappaffreal;
9692 FTYPE kappaemitffreal=kappaffreal;
9693 FTYPE kappanemitffreal=kappaffreal;
9698 FTYPE aa = 0.3564568198863157 - 0.19982886985727735*
Power(1. - 1.*varexpf,0.5646962005387126) + 0.18848660385247928*
Power(varexpf,13.940235230123522);
9699 FTYPE bb = 3.0641759806549147 + 0.25543546368991477*
Power(1. - 1.*varexpf,0.3128511829901932) + 0.07219025685305303*
Power(varexpf,1.3638629108321003);
9700 FTYPE cc = 5.99474041542733 - 1.4436038489430345*
Power(1. - 1.*varexpf,0.12828451539717775) - 1.4051865724895833*
Power(varexpf,3.0775239339815585);
9705 kappaffreal *= aa*pow(xi,-bb)*log(1.0+cc*xi);
9711 FTYPE aan = 4. - 2.06*
Power(1. - 1.*varexpf,1.) + 21.*
Power(varexpf,5.);
9712 FTYPE bbn = 3.1503757694129613 + 0.0008942404015805927*
Power(1. - 1.*varexpf,10.174621780103456) - 0.41243605382204196*
Power(varexpf,59.06672963853068);
9713 FTYPE ccn = 4.552732466653472e-7 + 2.3916452081662603*
Power(1. - 1.*varexpf,0.5522198226571239) + 5.27043093753271*
Power(varexpf,69.17066973904404);
9718 kappanffreal *= aan*pow(xi,-bbn)*log(1.0+ccn*xi);
9724 FTYPE aae = 0.5316463286647214;
9725 FTYPE bbe = 3.140128803410833;
9726 FTYPE cce = 4.522392868247536;
9731 kappaemitffreal *= aae*pow(xi,-bbe)*log(1.0+cce*xi);
9736 FTYPE bbne = 2.6669560547974855;
9742 kappanemitffreal *= aane*pow(xi,-bbne)*log(1.0+ccne*xi);
9756 FTYPE Reelow =1.706666666666667*thetae*(1 + 1.1*thetae +
Power(thetae,2) - 1.063803438337589*
Power(thetae,2.5));
9757 FTYPE Reehigh = 2.489326395711546*
Power(thetae,0.5)*(1.28 +
Log(1.123*thetae));
9758 FTYPE Ree = (thetae>=1.0 ? Reehigh : Reelow);
9762 FTYPE kappaffeereal = kappaffreal*factorffee;
9763 FTYPE kappanffeereal = kappanffreal*factorffee;
9764 FTYPE kappaemitffeereal = kappaemitffreal*factorffee;
9765 FTYPE kappanemitffeereal = kappanemitffreal*factorffee;
9776 FTYPE kappabfreal = kappaffreal*factorbf;
9777 FTYPE kappanbfreal = kappanffreal*factorbf;
9778 FTYPE kappaemitbfreal = kappaemitffreal*factorbf;
9779 FTYPE kappanemitbfreal = kappanemitffreal*factorbf;
9794 FTYPE kappachiantireal = kappaffreal*factorchianti;
9795 FTYPE kappanchiantireal = kappanffreal*factorchianti;
9796 FTYPE kappaemitchiantireal = kappaemitffreal*factorchianti;
9797 FTYPE kappanemitchiantireal = kappanemitffreal*factorchianti;
9808 FTYPE factorhm=2.75739E-48*pow(Tereal,11.2)*pow(
ZFACT,0.5)/(pow(rhoreal,0.5)*(1.0+
XFACT)*(1.0-
ZFACT)*Rei);
9809 FTYPE kappahmreal = kappaffreal*factorhm;
9810 FTYPE kappanhmreal = kappanffreal*factorhm;
9811 FTYPE kappaemithmreal = kappaemitffreal*factorhm;
9812 FTYPE kappanemithmreal = kappanemitffreal*factorhm;
9822 FTYPE factorchiantiopal=3
E-13*pow(Tereal,1.6)*pow(rhoreal,-0.4);
9823 FTYPE kappachiantiopalreal = kappaffreal*factorchianti*factorchiantiopal;
9824 FTYPE kappanchiantiopalreal = kappanffreal*factorchianti*factorchiantiopal;
9825 FTYPE kappaemitchiantiopalreal = kappaemitffreal*factorchianti*factorchiantiopal;
9826 FTYPE kappanemitchiantiopalreal = kappanemitffreal*factorchianti*factorchiantiopal;
9836 FTYPE factorhmopal=1E4*pow(Tereal,-1.2);
9837 FTYPE kappahmopalreal = kappaffreal*factorhm*factorhmopal;
9838 FTYPE kappanhmopalreal = kappanffreal*factorhm*factorhmopal;
9839 FTYPE kappaemithmopalreal = kappaemitffreal*factorhm*factorhmopal;
9840 FTYPE kappanemithmopalreal = kappanemitffreal*factorhm*factorhmopal;
9850 FTYPE kappanmolreal = kappamolreal;
9851 FTYPE kappaemitmolreal = kappamolreal;
9852 FTYPE kappanemitmolreal = kappamolreal;
9860 FTYPE kappalowdensityreal = 1.0/ ( 1.0/(kappamolreal + kappahmreal) + 1.0/(kappachiantireal + kappaffreal + kappaffeereal + kappabfreal) );
9861 FTYPE kappanlowdensityreal = 1.0/ ( 1.0/(kappanmolreal + kappanhmreal) + 1.0/(kappanchiantireal + kappanffreal + kappanffeereal + kappanbfreal) );
9862 FTYPE kappaemitlowdensityreal = 1.0/ ( 1.0/(kappaemitmolreal + kappaemithmreal) + 1.0/(kappaemitchiantireal + kappaemitffreal + kappaemitffeereal + kappaemitbfreal) );
9863 FTYPE kappanemitlowdensityreal = 1.0/ ( 1.0/(kappanemitmolreal + kappanemithmreal) + 1.0/(kappanemitchiantireal + kappanemitffreal + kappanemitffeereal + kappanemitbfreal) );
9871 FTYPE kappadensityreal = 1.0/ ( 1.0/(kappamolreal + kappahmopalreal) + 1.0/(kappachiantiopalreal) + 1.0/(kappachiantireal + kappaffreal + kappaffeereal + kappabfreal) );
9872 FTYPE kappandensityreal = 1.0/ ( 1.0/(kappanmolreal + kappanhmopalreal) + 1.0/(kappanchiantiopalreal) + 1.0/(kappanchiantireal + kappanffreal + kappanffeereal + kappanbfreal) );
9873 FTYPE kappaemitdensityreal = 1.0/ ( 1.0/(kappaemitmolreal + kappaemithmopalreal) + 1.0/(kappaemitchiantiopalreal) + 1.0/(kappaemitchiantireal + kappaemitffreal + kappaemitffeereal + kappaemitbfreal) );
9874 FTYPE kappanemitdensityreal = 1.0/ ( 1.0/(kappanemitmolreal + kappanemithmopalreal) + 1.0/(kappanemitchiantiopalreal) + 1.0/(kappanemitchiantireal + kappanemitffreal + kappanemitffeereal + kappanemitbfreal) );
9888 FTYPE kappasyreal=2.13E-11*nereal/Breal*pow(Tereal/1E10,-5.0);
9891 FTYPE q0 = 1.0/ (1.0+pow(3.3*thetae,-5.0));
9892 FTYPE q1 = 1.0 + pow(3.3*thetae,-5.0);
9893 FTYPE qsyn = exp(log(q0) + ( 0.5 + (1.0/M_PI)*atan(3.0+log(phi)) )*(log(q1) - log(q0)) );
9894 kappasyreal *= qsyn;
9896 FTYPE kappansyreal=kappasyreal;
9897 FTYPE kappaemitsyreal=kappasyreal;
9898 FTYPE kappanemitsyreal=kappasyreal;
9903 FTYPE aas=0.001 - 0.00030000000000000003*
Power(1. - 1.*varexpf,0.1) + 9.999*
Power(varexpf,100);
9904 FTYPE bbs=0.34915123451867575 + 0.6639854085243424*
Power(1. - 1.*varexpf,0.21551972527431834) + 0.6969780620639547*
Power(varexpf,109.95164261078826);
9905 FTYPE ccs=3.5 - 0.5*
Power(1. - 1.*varexpf,0.1) + 1.7999999999999998*
Power(varexpf,10);
9906 FTYPE dds=80.*
Power(2.718281828459045,1.611809565095832*
Power(varexpf,2) + 1.701470224678968*
Power(varexpf,19.123258595041275));
9907 FTYPE ees=2.9814806065542103 + 1.0252733145171038*
Power(1. - 1.*varexpf,0.03215194092829952);
9909 kappasyreal *= 1.0/( 1.0/(aas*pow(phi,-bbs)*log(1.0+ccs*phi)) + 1.0/(dds*pow(phi,-ees)) );
9913 FTYPE aans=0.0012 + 0.00005000000000000013*
Power(1. - 1.*varexpf,0.05) - 0.0009099999999999999*
Power(varexpf,7);
9914 FTYPE bbns=2.65 + 0.040000000000000036*
Power(1. - 1.*varexpf,0.02) - 0.5499999999999998*
Power(varexpf,30);
9915 FTYPE ccns=1.31 - 0.06000000000000005*
Power(1. - 1.*varexpf,0.02) + 1.69*
Power(varexpf,8);
9916 FTYPE ddns=
Power(2.718281828459045,1.5686159179138452 + 1.151292546497023*
Power(varexpf,2) + 1.151292546497023*
Power(varexpf,4) + 2.9834230661458117*
Power(varexpf,195.52338999877892));
9917 FTYPE eens=2.974042326402226 + 0.23373841071210855*
Power(1. - 1.*varexpf,11.45996831663565) - 0.1009858740118661*
Power(varexpf,201.41910724608374);
9919 kappansyreal *= 1.0/( 1.0/(aans*pow(phi,-bbns)*log(1.0+ccns*phi)) + 1.0/(ddns*pow(phi,-eens)) );
9925 FTYPE aaes=0.3188964065440192;
9926 FTYPE bbes=1.0768900596611533;
9927 FTYPE cces=0.0026066604097452917;
9928 FTYPE ddes=1.2627623470514082;
9929 FTYPE eees=2.987395175467846;
9931 kappaemitsyreal *= 1.0/( 1.0/(aaes*pow(phi,-bbes)*log(1.0+cces*phi)) + 1.0/(ddes*pow(phi,-eees)) );
9935 FTYPE aanes=0.00029169186404094216;
9936 FTYPE bbnes=1.3218673308302118;
9937 FTYPE ccnes=1.1634612365295334;
9938 FTYPE ddnes=120.01337286010565;
9939 FTYPE eenes=2.3252442378326856;
9941 kappanemitsyreal *= 1.0/( 1.0/(aanes*pow(phi,-bbnes)*log(1.0+ccnes*phi)) + 1.0/(ddnes*pow(phi,-eenes)) );
9955 FTYPE kappadcreal=7.36E-46*nereal*Trreal*Trreal*varexpf;
9958 FTYPE pdc = pow(1.0+thetae,-3.0);
9961 FTYPE kappandcreal=kappadcreal;
9962 FTYPE kappaemitdcreal=kappadcreal;
9963 FTYPE kappanemitdcreal=kappadcreal;
9968 FTYPE aadc=3.10482346702441e-8 + 4.162489998316388*
Power(1. - 1.*varexpf,1.6896115787181434) + 6.702210386421933*
Power(varexpf,0.941782674568028);
9969 FTYPE bbdc=0.04202835820213653 - 0.03337036710130055*
Power(1. - 1.*varexpf,0.4693413976733784) - 0.002097439704449637*
Power(varexpf,0.021742118711173992);
9970 FTYPE ccdc=3.8020287401719655 + 0.2009643395406986*
Power(1. - 1.*varexpf,0.25767159028663056) - 0.18040020940821044*
Power(varexpf,33.009903857413704);
9971 FTYPE dddc=0.11793169411613985 - 0.06262831233321572*
Power(1. - 1.*varexpf,0.34961114748511685) + 0.016889330731907487*
Power(varexpf,35.3723390749425);
9973 kappadcreal *= 1.0/( 1.0/aadc + 1.0/(bbdc*pow(thetar,-ccdc)) + 1.0/(dddc*pow(thetar,-ccdc/3.0)) );
9977 FTYPE aandc=87.4586627372423 - 76.35909344479292*
Power(1. - 1.*varexpf,0.13586634790456995) + 29.385274489218205*
Power(varexpf,285.27285622653653);
9978 FTYPE bbndc=1.1604225625247175 - 1.1192436863350592*
Power(1. - 1.*varexpf,0.1339258484550137) + 0.1955857773358507*
Power(varexpf,18.1030753003189);
9979 FTYPE ccndc=3.9257505990037376 + 0.04266044504755229*
Power(1. - 1.*varexpf,181.71820529516145) - 0.8002705331174753*
Power(varexpf,21.623589714867247);
9980 FTYPE ddndc=2.8625119194776856 - 2.716652341873481*
Power(1. - 1.*varexpf,0.1055549863784506) + 1.8741117433895793*
Power(varexpf,309.10043138151474);
9982 kappandcreal *= 1.0/( 1.0/aandc + 1.0/(bbndc*pow(thetar,-ccndc)) + 1.0/(ddndc*pow(thetar,-ccndc/3.0)) );
9988 FTYPE aaedc=6.335417556116487 - 0.058933985431348646*
Power(1. - 1.*varexpf,10.71634135441363) + 0.48765467682596686*
Power(varexpf,1.7536859250821957);
9989 FTYPE bbedc=0.008747549563345837 + 0.014220493179111593*
Power(1. - 1.*varexpf,0.36109285953630127) + 0.028227296616735897*
Power(varexpf,1.556532326060757);
9990 FTYPE ccedc=3.7824716575816524 + 0.18426457897750437*
Power(1. - 1.*varexpf,0.3661610552743282) - 0.1602799611377681*
Power(varexpf,15.388944611235807);
9991 FTYPE ddedc=0.11936934944026895 - 0.025640260828114658*
Power(1. - 1.*varexpf,0.39821885979760463) + 0.014988119208161788*
Power(varexpf,26.289513268487767);
9993 kappaemitdcreal *= 1.0/( 1.0/aaedc + 1.0/(bbedc*pow(thetar,-ccedc)) + 1.0/(ddedc*pow(thetar,-ccedc/3.0)) );
9997 FTYPE aanedc=197.97520843920014 - 94.77219201589888*
Power(1. - 1.*varexpf,0.9245008022924389) - 81.66905945089535*
Power(varexpf,1.010185534936844);
9998 FTYPE bbnedc=4.798508367755025e-11 + 1.0490039263100823*
Power(1. - 1.*varexpf,0.24861894212727034) + 1.3060425347854707*
Power(varexpf,1.1248777879043648);
9999 FTYPE ccnedc=3.4389272678970677 + 0.441710530722887*
Power(1. - 1.*varexpf,0.3614676129393117) - 0.4179301735940477*
Power(varexpf,14.840910987577235);
10000 FTYPE ddnedc=3.375250751641187 - 1.3820659504372945*
Power(1. - 1.*varexpf,0.31579132247596853) + 1.3742704814428137*
Power(varexpf,31.28257475560334);
10002 kappanemitdcreal *= 1.0/( 1.0/aanedc + 1.0/(bbnedc*pow(thetar,-ccnedc)) + 1.0/(ddnedc*pow(thetar,-ccnedc/3.0)) );
10006 FTYPE aadcpl=6.8308477566235427614037721740062696378080011960029;
10007 FTYPE bbdcpl=0.03737660642072567319753249796331010429215268940195;
10008 FTYPE ccdcpl=3.6263522982307391011876292625424350179419071859114;
10009 FTYPE dddcpl=0.13396823846337876662165099344371722644585939843797;
10012 FTYPE aandcpl=116.24082637910697352499285024420803245395311274652;
10013 FTYPE bbndcpl=1.3436648321999717296470441601823998776061527255987;
10014 FTYPE ccndcpl=3.0309933617509186060061141161349078396542823722783;
10015 FTYPE ddndcpl=4.7183006844597535352187951503958715513547062753028;
10025 FTYPE kappareal = kappadensityreal + kappasyreal + kappadcreal;
10026 FTYPE kappanreal = kappandensityreal + kappansyreal + kappandcreal;
10027 FTYPE kappaemitreal = kappaemitdensityreal + kappaemitsyreal + kappaemitdcreal;
10028 FTYPE kappanemitreal = kappanemitdensityreal + kappanemitsyreal + kappanemitdcreal;
10039 FTYPE kappa_es_fermicorr = (1.0/(1.0+2.7E11*
prpow(rhoreal,2.0)/
prpow(Tereal,2.0)));
10042 FTYPE kappa_es_kncorr = (1.0/(1.0+
prpow(Tereal/4.5E8,0.86)));
10046 FTYPE kappaesreal = 0.2*(1.0+
XFACT)*(rhoreal)*kappa_es_fermicorr*kappa_es_kncorr;
10057 *kappa=kappareal*overopacitybaralt;
10058 *kappan=kappanreal*overopacitybaralt;
10059 *kappaemit=kappaemitreal*overopacitybaralt;
10060 *kappanemit=kappanemitreal*overopacitybaralt;
10061 *kappaes=kappaesreal*overopacitybaralt;
10091 FTYPE kappa_forcompt_relcorr = (1.0 + 3.683*thetae+4.0*thetae*thetae)/(1.0 + thetae);
10094 FTYPE kappa_es_fermicorr = 1.0/(1.0+2.7E11*rhoreal/
prpow(Tereal,2.0));
10096 FTYPE kappa_forcompt_code = 0.2*(1.0+
XFACT)*kappa_es_fermicorr*kappa_forcompt_relcorr/
OPACITYBAR;
10100 FTYPE preterm3 = -4.0*rho0*kappa_forcompt_code*(Tgas - Tradff)*(TEMPBAR/
TEMPELE)*Ruu;
10120 #if(WHICHPROBLEM==RADDONUT)
10137 #if(WHICHPROBLEM==RADDONUT)
10152 #if(WHICHPROBLEM==RADDONUT)
10167 #if(WHICHPROBLEM==RADDONUT)
10179 #if(WHICHPROBLEM==RADDONUT)
10199 #define JET6LIKEUSERCOORD 0
10200 #define UNIHALFUSERCOORD 1
10203 #define WHICHUSERCOORD JET6LIKEUSERCOORD
10207 static FTYPE npow,
r1jet,
njet1,
njet,
r0jet,
rsjet,
Qjet,
ntheta,
htheta,
rsjet2,
r0jet2,
rsjet3,
r0jet3,
rs,
r0,
npow2,
cpow2,
rbr,
x1br,
h0,
cpow3;
10306 x1br = log( rbr -
R0 ) /
npow;
10314 fprintf(out,
"%21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g\n",npow,r1jet,njet,r0jet,rsjet,Qjet,ntheta,htheta,rsjet2,r0jet2,rsjet3,r0jet3,rs,r0,npow2,cpow2,rbr,x1br,cpow3);
10322 fscanf(in,
HEADER19IN,&npow,&r1jet,&njet,&r0jet,&rsjet,&Qjet,&ntheta,&htheta,&rsjet2,&r0jet2,&rsjet3,&r0jet3,&rs,&r0,&npow2,&cpow2,&rbr,&x1br,&cpow3);
10327 void read_coord_parms_mpi_user(
int defcoordlocal)
10331 MPI_Bcast(&npow, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10332 MPI_Bcast(&r1jet, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10333 MPI_Bcast(&njet, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10334 MPI_Bcast(&r0jet, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10335 MPI_Bcast(&rsjet, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10336 MPI_Bcast(&Qjet, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10337 MPI_Bcast(&ntheta, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10338 MPI_Bcast(&htheta, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10339 MPI_Bcast(&rsjet2, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10340 MPI_Bcast(&r0jet2, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10341 MPI_Bcast(&rsjet3, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10342 MPI_Bcast(&r0jet3, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10343 MPI_Bcast(&rs, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10344 MPI_Bcast(&r0, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10345 MPI_Bcast(&npow2, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10346 MPI_Bcast(&cpow2, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10347 MPI_Bcast(&rbr, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10348 MPI_Bcast(&x1br, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10349 MPI_Bcast(&cpow3, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
10362 #if(0) // no change in exponentiation
10364 V[1] =
R0+exp(pow(X[1],npow)) ;
10365 #elif(WHICHUSERCOORD==UNIHALFUSERCOORD)
10368 theexp = npow*X[1];
10371 FTYPE gconst2=gconst1*.000001;
10372 V[1] =
R0 + gconst1*X[1] + gconst2*exp(theexp);
10375 #elif(WHICHUSERCOORD==JET6LIKEUSERCOORD)
10378 #define cr(x) (exp(-1.0/(x)))
10379 #define tr(x) (cr(x)/(cr(x) + cr(1.0-(x))))
10380 #define trans(x,L,R) ((x)<=(L) ? 0.0 : ((x)>=(R) ? 1.0 : tr(((x)-(L))/((R)-(L)))) )
10381 #define transR(x,center,width) ( 0.5*(1.0-tanh(+((x)-(center))/(width))))
10382 #define transL(x,center,width) ( 0.5*(1.0-tanh(-((x)-(center))/(width))))
10383 #define transM(x,center,width) ( exp(-pow(((x)-(center))/((width)*0.5),2.0) ) )
10385 #define plateau(x,L,R,W) (trans(x,(L)-0.5*(W),(L)+0.5*(W))*(1.0-trans(x,(R)-0.5*(W),(R)+0.5*(W))))
10389 FTYPE theexp = npow*X[1];
10390 if( X[1] > x1br ) {
10391 theexp += cpow2 * pow(X[1]-x1br,npow2);
10393 V[1] =
R0+exp(theexp);
10397 #define line1r(x,w) (Rout)
10398 #define line2r(x,w) (Routeq)
10399 #define line3r(x,w) (Rout))
10403 #define thetajon(x,w,xp1,xp2) (line1r(x,w)*(1.0-trans(x,xp1,xp2)) + line2r(x,w)*trans(x,xp1,xp2))
10412 #define lineeqr(x,w) (R0 + exp(npow*(X[1]-mysx1)*0.70 + npow*mysx1 ) )
10413 #define linepoler(x,w) (R0 + exp(npow*X[1]))
10414 #define thetaLr(x,wp,weq,xp1,xp2) ( linepoler(x,wp)*(1.0-trans(x,xp1,xp2)) + lineeqr(x,weq)*trans(x,xp1,xp2) )
10415 #define thetajon2(x,wp,weq,xp1,xp2) ( x<0.5 ? thetaLr(x,wp,weq,xp1,xp2) : thetaLr(1.0-x,wp,weq,xp1,xp2) )
10417 FTYPE Routvsx2,R0vsx2;
10418 FTYPE Routeq = 100.0;
10436 Routvsx2 =
thetajon2(X[2], 0.0, 0.0, xpole, xeq);
10440 R0vsx2 = R0eq +
R0*((Routvsx2-Routeq)/(
Rout-Routeq));
10452 Routvsx2 =
thetajon(X[2],0.0,0.25,0.75);
10468 V[1] =
R0+exp(npow*X[1]);
10473 V[1] =
Rin + (V[1]-
Rin)*(1.0 - AA*sin(X[2]*M_PI));
10486 FTYPE r0rbr=rbr/2.0;
10487 FTYPE switch0 = 0.5+1.0/M_PI*atan((V[1]-rbr)/r0rbr);
10489 FTYPE V1 =
R0+exp(npow*X[1]);
10490 FTYPE V2 =
R0+exp(npow*X[1] + cpow2 * pow(cpow3*(X[1]-x1br*1.0),npow2));
10492 V[1] = V1*(1.0-switch0) + V2*switch0;
10503 FTYPE theta1,theta2,arctan2;
10508 FTYPE myhslope=2.0-Qjet*pow(V[1]/r1jet,-njet*(0.5+1.0/M_PI*atan(V[1]/r0jet-rsjet/r0jet)));
10509 theta1 = M_PI * X[2] + ((1. - myhslope) * 0.5) * mysin(2. * M_PI * X[2]);
10522 FTYPE localrbrr0=100.0;
10524 FTYPE switch0 = 0.5+1.0/M_PI*atan((V[1]-localrbr)/localrbrr0);
10525 FTYPE switch2 = 1.0-switch0;
10529 FTYPE myhslope1=h0 + pow( (V[1]-rsjet3)/r0jet3 , njet);
10530 FTYPE myhslope2=h0 + pow( (localrbr-rsjet3)/r0jet3 , njet);
10531 FTYPE myhslope = myhslope1*switch2 + myhslope2*switch0;
10535 myhslope=myhslope1;
10539 if(X[2]>1.0) myx2=2.0-X[2];
10540 else if(X[2]<0.0) myx2=-X[2];
10543 FTYPE th2 = 0.5*M_PI*(1.0 + atan(myhslope*(myx2-0.5))/atan(myhslope*0.5));
10545 if(X[2]>1.0) th2=2.0*M_PI-th2;
10546 else if(X[2]<0.0) th2=-th2;
10550 myhslope1=2.0-Qjet*pow(V[1]/r1jet,-njet*(0.5+1.0/M_PI*atan(V[1]/r0jet-rsjet/r0jet)));
10551 myhslope2=2.0-Qjet*pow(localrbr/r1jet,-njet*(0.5+1.0/M_PI*atan(localrbr/r0jet-rsjet/r0jet)));
10552 myhslope = myhslope1*switch2 + myhslope2*switch0;
10554 FTYPE th0 = M_PI * X[2] + ((1. - myhslope) * 0.5) * mysin(2. * M_PI * X[2]);
10558 switch0 = 0.5+1.0/M_PI*atan((V[1]-rs)/r0);
10559 switch2 = 1.0-switch0;
10562 theta1 = th0*switch2 + th2*switch0;
10570 theta2 = M_PI*0.5*(htheta*(2.0*X[2]-1.0)+(1.0-
htheta)*pow(2.0*X[2]-1.0,ntheta)+1.0);
10573 arctan2 = 0.5 + 1.0/M_PI*(atan( (V[1]-rsjet2)/r0jet2) );
10576 V[2] = theta2 + arctan2*(theta1-theta2);
10584 FTYPE fracpole=(1.0-fraceq)/2.0;
10585 FTYPE x2p1=0.0+fracpole;
10586 FTYPE x2p2=1.0-fracpole;
10594 FTYPE switchh0 = 0.5+0.5*tanh((X[2]-x2p1)/swide);
10595 FTYPE switchh2 = 0.5+0.5*tanh((X[2]-x2p2)/swide);
10598 FTYPE theq = M_PI * X[2] + ((1. - eqh) * 0.5) * mysin(2. * M_PI * X[2]);
10600 FTYPE thup=switchh0*theq + (1.0-switchh0)*theta1;
10602 V[2]=switchh2*theta1 + (1.0-switchh2)*thup;
10610 FTYPE transwidth=0.06;
10612 FTYPE transR=0.5*(1.0-tanh(+(X[2] - xcent)/transwidth));
10613 FTYPE transL=0.5*(1.0-tanh(-(X[2] - xcent)/transwidth));
10617 FTYPE wpar=pow( (V[1])/1 , -njet);
10622 V[2] = line1*transR + line2*transL;
10632 #define line1(x,w) ((x)*(w))
10633 #define line2(x,w) ((x)*(w)+M_PI-(w))
10634 #define line3(x,w) ((x)*(w))
10637 #define wparsam(x,r) (h0 + pow(0.15 + ((r)-0.0)/10.0 , -njet))
10639 #define thetasam(x,r,w,xp1,xp2) (line1(x,w)*(1.0-trans(x,xp1,xp2)) + line2(x,w)*trans(x,xp1,xp2))
10653 #define lineeq(x,w) ((x)*(w)+(0.5*M_PI)-(0.5*w))
10654 #define linepole(x,w) (line1(x,w))
10655 #define thetaL(x,wp,weq,xp1,xp2) ( linepole(x,wp)*(1.0-trans(x,xp1,xp2)) + lineeq(x,weq)*trans(x,xp1,xp2) )
10656 #define thetasam2(x,wp,weq,xp1,xp2) ( x<0.5 ? thetaL(x,wp,weq,xp1,xp2) : -thetaL(1.0-x,wp,weq,xp1,xp2)+M_PI )
10663 V[2] =
thetasam2(X[2], poleslope, eqslope, xpole, xeq);
10673 V[3]=2.0*M_PI*X[3];
10678 V[1] =
R0+exp(X[1]) ;
10679 V[2] = M_PI * X[2] + ((1. -
hslope) / 2.) * sin(2. * M_PI * X[2]);
10680 V[3]=2.0*M_PI*X[3];
10689 dualfprintf(fail_file,
"Should not be computing USERCOORDS analytically\n");
10691 dxdxp[3][3] = 2.0*M_PI;
10694 void set_points_user(
void)
10705 startx[1] = 0.3999985081775278946780799743777598329673;
10709 FTYPE endx1=21.40529883372801383045167738115556702610;
10727 trifprintf(
"ITERATIVE dx1: Rout=%21.15g R0=%21.15g npow=%21.15g cpow2=%21.15g cpow3=%21.15g npow2=%21.15g x1br=%21.15g rbr=%21.15g\n",
Rout,
R0,npow,cpow2,cpow3,npow2,x1br,rbr);
10729 FTYPE x1max0, x1max,dxmax;
10732 const int ITERMAX = 100;
10743 for( iter = 0; iter < ITERMAX; iter++ ) {
10747 if( fabs((x1max - x1max0)/x1max) < RELACC ) {
10753 dxmax= (pow( (log(
Rout-
R0) - npow*x1max0)/cpow2, 1./npow2 ) + x1br*1.0) - x1max0;
10758 FTYPE dVdx1=(npow + cpow2*npow2*cpow3*pow(cpow3*(x1max0-x1br*1.0),npow2-1.0)) * exp(npow*x1max0 + cpow2 * pow(cpow3*(x1max0-x1br*1.0),npow2));
10759 FTYPE V0 =
R0 + exp(npow*x1max0 + cpow2 * pow(cpow3*(x1max0-x1br*1.0),npow2));
10761 dxmax=(
Rout-V0)/dVdx1;
10763 dualfprintf(fail_file,
"dVdx1=%g V0=%g dxmax=%g x1max=%g x1max0=%g\n",dVdx1,V0,dxmax,x1max,x1max0);
10767 FTYPE dampingfactor=0.5;
10768 x1max = x1max0 + dampingfactor*dxmax;
10773 if( iter == ITERMAX ) {
10774 trifprintf(
"Error: iteration procedure for finding x1max has not converged: x1max = %g, dx1max/x1max = %g, iter = %d\n",
10775 x1max, (x1max-x1max0)/x1max, iter );
10779 trifprintf(
"x1max = %g (dx1max/x1max = %g, itno = %d)\n", x1max, (x1max-x1max0)/x1max, iter );
10804 dualfprintf(fail_file,
"ihoradjust=%21.15g totalsize[1]=%d Rhor=%21.15g R0=%21.15g npow=%21.15g Rout=%21.15g\n",ihoradjust,
totalsize[1],
Rhor,
R0,npow,
Rout);
10818 #define FRACN1 (0.1)
10819 #define ADJUSTFRACT (0.25)