fixed bug with u_av_chop in multiple init.c

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

PROBLEMS/MONOPOLE_2D/define.h

@@ -20,60 +20,68 @@
 //#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
 
 /************************************/
 //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.
 #define B2UURATIOMAX 1.e100
 #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,19 +120,13 @@
 #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
 #define MINX (1.5*r_horizon_BL(BHSPIN))
 #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

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
 

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]);

physics.c

@@ -1220,13 +1220,15 @@ 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);
 	    //print_metric(geom.gg);
 	    //print_Nvector(pp,NV);
 	    my_err("nan in flux_prime\n");
 	    exit(-1);
+	    #endif
 	}
     }
 

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)