diff --git a/PROBLEMS/JETCOORDS/init.c b/PROBLEMS/JETCOORDS/init.c index 124e0cf..8da03ae 100644 --- a/PROBLEMS/JETCOORDS/init.c +++ b/PROBLEMS/JETCOORDS/init.c @@ -178,7 +178,7 @@ else // we are inside the donut //midplane at r_mag init_dsandvels_fishbone_moncrief(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); //midplane at rchop - init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); + init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_chop, &ell); //vector potential follows contours of UU ldouble uchop = u_av - u_av_chop; //vpot->zero on contour of radius r=rchop ldouble uchopmid = u_av_mid - u_av_chop; //vpot->zero away from midplane @@ -213,7 +213,7 @@ else // we are inside the donut //midplane at r_mag init_dsandvels_fishbone_moncrief(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); //midplane at rchop - init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); + init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_chop, &ell); //vector potential follows contours of UU ldouble uchop = u_av - u_av_chop; //vpot->zero on contour of radius r=rchop ldouble uchopmid = u_av_mid - u_av_chop; //vpot->zero away from midplane diff --git a/PROBLEMS/MONOPOLE_2D/define.h b/PROBLEMS/MONOPOLE_2D/define.h index f1adbde..16bf356 100755 --- a/PROBLEMS/MONOPOLE_2D/define.h +++ b/PROBLEMS/MONOPOLE_2D/define.h @@ -20,14 +20,14 @@ //#define U2PCONV 1.e-14 #define FORCEFREE -//#define NOLOGINS + /************************************/ //reconstruction / Courant /************************************/ -#define INT_ORDER 1 //TODO why is int_order 1 more stable??? +#define INT_ORDER 2 //TODO why is int_order 1 more stable??? #define TIMESTEPPING RK2IMEX -#define TSTEPLIM .9 +#define TSTEPLIM .8 #define FLUXLIMITER 1 #define FLUXMETHOD LAXF_FLUX #define MINMOD_THETA 1.5 @@ -35,22 +35,29 @@ /************************************/ //rmhd floors /************************************/ +#define SPLIT_MONOPOLE -#define B2RHORATIOMAXINIT 100//100 -#define B2UURATIOMAXINIT 100//100 +#define RHOATMMIN 1.e-20 +#define UINTATMMIN 1.e-20 + +#define B2RHORATIOMAXINIT 500//100//100 +#define B2UURATIOMAXINIT 500//100//100 #if defined(FORCEFREE) +#define NOLOGINS //#define ENFORCEENTROPY -#define HYBRID_FORCEFREE -#define HYBRID_FORCEFREE_SIGMACUT 25 +//#define HYBRID_FORCEFREE +//#define HYBRID_FORCEFREE_SIGMACUT 25 #define FORCEFREE_SOLVE_PARALLEL -#//define FORCEFREE_SOLVE_PARALLEL_OLD +#define FORCEFREE_PARALLEL_COLD //#define SKIPALLFLOORS // TODO seems critical for SOLVE_PARALLEL? -//#define CORRECT_POLARAXIS -//#define NCCORRECTPOLAR 2 +#define CORRECT_POLARAXIS +#define NCCORRECTPOLAR 2 + +#define B2RHOFLOORFRAME FFFRAME #define UURHORATIOMIN 0. #define UURHORATIOMAX 50. #define B2UURATIOMIN 0. @@ -58,22 +65,23 @@ #define B2RHORATIOMIN 0. #define B2RHORATIOMAX 1.e100 -#define GAMMAMAXFF 500. //lower than GAMMAMAXHD? -#define GAMMAMAXHD 500. //why can't this be pushed higher on the monopole? +#define GAMMAMAXFF 10//100. //lower than GAMMAMAXHD? +#define GAMMAMAXHD 100//100. //why can't this be pushed higher on the monopole? #else #define CORRECT_POLARAXIS #define NCCORRECTPOLAR 2 -#define B2RHOFLOORFRAME ZAMOFRAME +#define B2RHOFLOORFRAME DRIFTFRAME//ZAMOFRAME #define UURHORATIOMIN 0. #define UURHORATIOMAX 50. #define B2UURATIOMIN 0. -#define B2UURATIOMAX 2500. +#define B2UURATIOMAX B2UURATIOMAXINIT #define B2RHORATIOMIN 0. -#define B2RHORATIOMAX 500. -#define GAMMAMAXHD 500. +#define B2RHORATIOMAX B2RHORATIOMAXINIT + +#define GAMMAMAXHD 100. #endif /************************************/ @@ -92,15 +100,16 @@ #define RMAX 200. #define MKS1R0 MKSR0 +#define TNX 128 +#define TNY 256//128 +#define TNZ 1 + #ifdef myMKS1COORDS //modified Kerr-Shild #define MYCOORDS MKS1COORDS #define MINX (log(RMIN - MKSR0)) #define MAXX (log(RMAX - MKSR0)) -#define MINY 0.001*M_PI -#define MAXY M_PI - 0.001*M_PI -#define TNX 128 -#define TNY 64//256 -#define TNZ 1 +#define MINY 0.005*M_PI +#define MAXY M_PI - 0.005*M_PI #elif defined(myMKS2COORDS) //modified Kerr-Shild #define MYCOORDS MKS2COORDS @@ -111,9 +120,6 @@ #define MAXX (log(RMAX - MKSR0)) #define MINY (0.001) #define MAXY (1.-0.001) -#define TNX 128 -#define TNY 128 -#define TNZ 1 #else //Schwarzschild #define MYCOORDS SCHWCOORDS @@ -121,9 +127,6 @@ #define MAXX (25.3) #define MINY (0.01*Pi/2.) #define MAXY (Pi-0.01*Pi/2.) -#define TNX 64 -#define TNY 32 -#define TNZ 1 #endif #define PHIWEDGE (M_PI/2.) @@ -158,5 +161,3 @@ #define GAMMA (4./3.) #define DTOUT1 1 #define DTOUT2 1.e0 -#define RHOATMMIN 1.e-20 -#define UINTATMMIN 1.e-20 diff --git a/PROBLEMS/RADSURVEY/init.c b/PROBLEMS/RADSURVEY/init.c index 7afebe4..fd8b9de 100644 --- a/PROBLEMS/RADSURVEY/init.c +++ b/PROBLEMS/RADSURVEY/init.c @@ -208,12 +208,12 @@ else // we are inside the donut #ifdef LIMOTORUS init_dsandvels_limotorus(r_mag, th_mag, BHSPIN, &rho, &u_av, &ell); // current loc init_dsandvels_limotorus(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane - init_dsandvels_limotorus(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane at rchop + init_dsandvels_limotorus(rchop, M_PI/2., BHSPIN, &rho, &u_av_chop, &ell); // midplane at rchop rin = LT_RIN; #else init_dsandvels_fishbone_moncrief(r_mag, th_mag, BHSPIN, &rho, &u_av, &ell); // current loc init_dsandvels_fishbone_moncrief(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane - init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane at rchop + init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_chop, &ell); // midplane at rchop rin = FM_rin; #endif @@ -249,12 +249,12 @@ else // we are inside the donut #ifdef LIMOTORUS init_dsandvels_limotorus(r_mag, th_mag, BHSPIN, &rho, &u_av, &ell); // current loc init_dsandvels_limotorus(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane - init_dsandvels_limotorus(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane at rchop + init_dsandvels_limotorus(rchop, M_PI/2., BHSPIN, &rho, &u_av_chop, &ell); // midplane at rchop rin = LT_RIN; #else init_dsandvels_fishbone_moncrief(r_mag, th_mag, BHSPIN, &rho, &u_av, &ell); // current loc init_dsandvels_fishbone_moncrief(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane - init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane at rchop + init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_chop, &ell); // midplane at rchop rin = FM_rin; #endif diff --git a/magn.c b/magn.c index e348db8..1462bb3 100644 --- a/magn.c +++ b/magn.c @@ -631,6 +631,16 @@ calc_BfromA_core() B[2]=(dA[3][1] - dA[1][3])/geom.gdet; B[3]=(dA[1][2] - dA[2][1])/geom.gdet; +#if PROBLEM==132 && defined(SPLIT_MONOPOLE) + ldouble xxBL[4]; + get_xx_arb(ix,iy,iz,xxBL,BLCOORDS); + if(xxBL[2]>0.5*M_PI) + { + B[1]*=-1; + B[2]*=-2; + B[3]*=-3; + } +#endif set_u(pvecpot,1,ix,iy,iz,B[1]); set_u(pvecpot,2,ix,iy,iz,B[2]); set_u(pvecpot,3,ix,iy,iz,B[3]); diff --git a/physics.c b/physics.c index 77efade..f672773 100644 --- a/physics.c +++ b/physics.c @@ -1220,6 +1220,7 @@ int f_flux_prime(ldouble *pp, int idim, int ix, int iy, int iz,ldouble *ff,int l { if(isnan(T[ii][jj])) { + #if !(defined(FORCEFREE) && !defined(HYBRID_FORCEFREE)) printf("%d > nan tmunu: %d %d %e at %d %d %d\n",PROCID,ii,jj,T[ii][jj],ix+TOI,iy+TOJ,iz+TOK); //printf("%d > nan tmunu: %e %e %e %e\n",PROCID,gamma,pre,w,eta); //print_4vector(ucon); @@ -1227,6 +1228,7 @@ int f_flux_prime(ldouble *pp, int idim, int ix, int iy, int iz,ldouble *ff,int l //print_Nvector(pp,NV); my_err("nan in flux_prime\n"); exit(-1); + #endif } } diff --git a/problem.h b/problem.h index 7087025..7654bea 100644 --- a/problem.h +++ b/problem.h @@ -152,8 +152,8 @@ //if you are using a problem older than 100 with a different kappa defined in kappa.c //make sure your kappa.c ends with (kappa=(.....)) NOT (return kappa)!! -//#define PROBLEM 132 -#define PROBLEM 144 +#define PROBLEM 132 +//#define PROBLEM 144 #if(PROBLEM==144)