long delayed commit, way too many updates

PROBLEMS/MONOPOLE_2D/define.h

@@ -36,18 +36,18 @@
 //rmhd floors
 /************************************/
 
-#define B2RHORATIOMAXINIT 100
-#define B2UURATIOMAXINIT 100
+#define B2RHORATIOMAXINIT 100//100
+#define B2UURATIOMAXINIT 100//100
 
 #if defined(FORCEFREE)
 
+//#define ENFORCEENTROPY
 #define HYBRID_FORCEFREE
-#define HYBRID_FORCEFREE_SIGMACUT 50
+#define HYBRID_FORCEFREE_SIGMACUT 25
 #define FORCEFREE_SOLVE_PARALLEL
-//#define FORCEFREE_SOLVE_PARALLEL_OLD
+#//define FORCEFREE_SOLVE_PARALLEL_OLD
 
 //#define SKIPALLFLOORS // TODO seems critical for SOLVE_PARALLEL? 
-#define B2RHOFLOORFRAME DRIFTFRAME //ZAMOFRAME //DRIFTFRAME
 //#define CORRECT_POLARAXIS
 //#define NCCORRECTPOLAR 2
 
@@ -72,8 +72,8 @@
 #define B2UURATIOMIN 0.
 #define B2UURATIOMAX 2500.
 #define B2RHORATIOMIN 0.
-#define B2RHORATIOMAX 200.
-#define GAMMAMAXHD 50.
+#define B2RHORATIOMAX 500.
+#define GAMMAMAXHD 500.
 #endif
 
 /************************************/
@@ -145,7 +145,7 @@
 #define OUTVEL VEL4
 #define ALLSTEPSOUTPUT 0
 #define NSTEPSTOP 1.e10
-#define NOUTSTOP 102
+#define NOUTSTOP 110
 #define SILOOUTPUT 1
 #define PRIMOUTPUT 1
 //#define PRIMOUTPUTINMYCOORDS

finite.c

@@ -15,6 +15,7 @@ int
 reduce_order_check(ldouble *pm2,ldouble *pm1,ldouble *p0,ldouble *pp1,ldouble *pp2,int ix,int iy,int iz)
 {
   int reconstrpar=0;
+
 #ifdef REDUCEORDERTEMP
   ldouble t0,tp1,tm1,tmin;
   t0 =calc_PEQ_Tfromurho(p0[UU], p0[RHO]);
@@ -31,7 +32,22 @@ reduce_order_check(ldouble *pm2,ldouble *pm1,ldouble *p0,ldouble *pp1,ldouble *p
   if(xxBL[1]<REDUCEORDERRADIUS)
     reconstrpar=1;
 #endif
-
+
+  // ANDREW TODO -- different behavior on boundary vs interior of FF region? 
+#if defined(FORCEFREE) && defined(REDUCEORDERFF)
+  int ffflag;
+  ffflag = get_cflag(FFINVFLAG,ix,iy,iz);
+  if(ffflag>0)
+    reconstrpar=1;
+#endif
+
+#ifdef REDUCEORDERATBH
+  ldouble xxBL[4];
+  get_xx_arb(ix,iy,iz,xxBL,BLCOORDS);
+  if(xxBL[1]<rhorizonBL)
+    reconstrpar=1;
+#endif
+
   return reconstrpar;
 }
 
@@ -93,20 +109,19 @@ avg2point(ldouble *um2,ldouble *um1,ldouble *u0,ldouble *up1,ldouble *up2,
 {
   ldouble r0[NV],rm1[NV],rp1[NV];
 
-  if(param!=0) //overrule the standard reconstruction
+  int int_order_local = INT_ORDER;
+
+  if(param!=0) //reduce integration order
   {
-    int i;
-    if(param==1) //DONOR CELL
-    {
-      for(i=0;i<NV;i++)
-      {
-        ur[i]=u0[i];
-        ul[i]=u0[i];
-      }
-    }
+    //int_order_local = INT_ORDER-param;
+    int_order_local=0; // ANDREW TODO reduce directly to 0
   }  // if(param!=0)
+
+  // default to donor cell
+  if(int_order_local<0)
+    int_order_local = 0;
 
-  else if(INT_ORDER==0) //DONOR CELL
+  if(int_order_local==0) //DONOR CELL
   {
     int i;
 
@@ -117,7 +132,7 @@ avg2point(ldouble *um2,ldouble *um1,ldouble *u0,ldouble *up1,ldouble *up2,
     }
   }  // else if(INT_ORDER==0)
 
-  else if(INT_ORDER==1) //linear
+  else if(int_order_local==1) //linear
   {
     ldouble diffpar=theta;
     int i;
@@ -201,7 +216,7 @@ avg2point(ldouble *um2,ldouble *um1,ldouble *u0,ldouble *up1,ldouble *up2,
     } // for(i=0;i<NV;i++)
   }  // else if(INT_ORDER==1)
 
-  else if(INT_ORDER==2) //parabolic PPM
+  else if(int_order_local==2) //parabolic PPM
   {
     //The following is based on Colella & Woodward (J. Comp. Phys. 54, 174, 1984).
     //It uses five points: m2, m1, 0, p1, p2.
@@ -1876,7 +1891,7 @@ int
 set_grid(ldouble *mindx,ldouble *mindy, ldouble *mindz, ldouble *maxdtfac)
 {
   int i1,i2,ix,iy,iz;
-  ldouble mdx,mdy,mdz,dx,dy,dz,gloc[4][5],xx[4];
+  ldouble mdx,mdy,mdz,dx,dy,dz,gloc[4][5],xx[4],xxBL[4];
   mdx=mdy=mdz=-1;
   ldouble maxdt=-1;
 
@@ -1908,14 +1923,6 @@ set_grid(ldouble *mindx,ldouble *mindy, ldouble *mindz, ldouble *maxdtfac)
       if(i1>-NG) set_x(i1-1,2,.5*(get_xb(i1,2)+get_xb(i1-1,2)));
     }
 
-  //consistency check
-#if(MYCOORDS==BLCOORDS)
-  if(get_x(-1,0)<=rhorizonBL)
-    {
-      printf("ix %d > %f\n",-1,get_x(-1,0));
-      my_err("-1 cell inside horizon\n");
-    }
-#endif
 
   // what is the minimal cell size
   for(ix=ix1;ix<ix2;ix++)
@@ -2467,23 +2474,14 @@ int
 get_xx_arb(int ix,int iy,int iz,ldouble *xx,int COORDSOUT)
 {
 
-#ifdef PRECOMPUTE_MY2OUT // use precomputed coordinates if COORDS == OUTCOORDS
-  if(COORDSOUT == OUTCOORDS)
-  {
+#if defined(PRECOMPUTE_MY2OUT) && (COORDSOUT==OUTCOORDS) // use precomputed coordinates if COORDS == OUTCOORDS
     get_xxout(ix, iy, iz, xx); // time will be nonsense! seems ok everywhere this is used
-  }
-  else
-  {
-    ldouble xx0[4];
-    get_xx(ix,iy,iz,xx0);
-    coco_N(xx0,xx,MYCOORDS,COORDSOUT);
-  }
 #else
     ldouble xx0[4];
     get_xx(ix,iy,iz,xx0);
     coco_N(xx0,xx,MYCOORDS,COORDSOUT);
-#endif  
-  return 0;
+#endif
+    return 0;
 }
 
 
@@ -5187,8 +5185,7 @@ cell_fixup(int type)
 	  (get_cflag(RADIMPFIXUPFLAG,ix,iy,iz)!=0 && type==FIXUP_RADIMP)) && is_cell_active(ix,iy,iz)
 	)
 	{
-	  //ANDREW - I think maybe this should be here, to set the flag back to its default? 
-	  //this should not be here, should it?
+
 	  /*
 	  if(type==FIXUP_U2PMHD) set_cflag(HDFIXUPFLAG,ix,iy,iz,0); //try only once
 	  if(type==FIXUP_U2PRAD) set_cflag(RADFIXUPFLAG,ix,iy,iz,0); //try only once

frames.c

@@ -95,8 +95,9 @@ trans_pmhd_coco(ldouble *ppin, ldouble *ppout, int CO1,int CO2, ldouble *xxvec,
     }
 
 #ifdef FORCEFREE
-    pp2[UUFF]=pp1[UUFF];
-
+
+
+    // perpindicular
     ldouble uconff[4];
     uconff[0]=0.;
     uconff[1]=pp1[VXFF];
@@ -108,6 +109,13 @@ trans_pmhd_coco(ldouble *ppin, ldouble *ppout, int CO1,int CO2, ldouble *xxvec,
     pp2[VXFF]=uconff[1];
     pp2[VYFF]=uconff[2];
     pp2[VZFF]=uconff[3];
+
+    //parallel
+    //ANDREW TODO this is NOT right
+    //but it doesn't matter for now as long as FF prims are not used in evolution
+    //FF conserveds are all derived from normal prims
+    pp2[UUFF]=pp1[UUFF];
+
 #endif
 
 #endif

ko.c

@@ -96,7 +96,8 @@ main(int argc, char **argv)
 
   //precalculates metric, Christoffels, etc.
   calc_metric();
-
+  calc_cells_under_horiz();
+
   //save coordinate file
   #ifdef COORDOUTPUT
   fprint_coordfile(folder,"coord");
@@ -112,12 +113,6 @@ main(int argc, char **argv)
 #endif 
 
   //print scalings GU->CGS
-  #ifdef FORCEFREE
-  if(PROCID==0) printf("NVHD %d NVMHD %d NV %d\n",NVHD,NVMHD,NV);
-  if(PROCID==0) printf("%d %d %d %d %d %d %d %d %d %d %d %d %d\n",
-		       RHO,UU,VX,VY,VZ,ENTR,B1,B2,B3,UUFF,VXFF,VYFF,VZFF);
-  #endif
-
   if(PROCID==0) print_scalings();
 
   //**************
@@ -177,6 +172,11 @@ main(int argc, char **argv)
         printf("Sending initial data... ");
 	fflush(stdout);
       }
+
+#ifdef FORCEFREE
+      fill_ffprims(); // make force-free primitives consistent
+#endif
+
       //exchange initial state
       mpi_exchangedata();  
       calc_avgs_throughout();
@@ -284,6 +284,7 @@ main(int argc, char **argv)
   int iix,iiy,iiz;
   struct geometry geom;
   ldouble PERTURBATION;
+
 
   // reset Te so that ue/ugas is constant
 #ifdef RESETELECTRONTEMPERATURETOUEUGASRATIO

ko.h

@@ -177,6 +177,7 @@ ldouble min_dx,min_dy,min_dz;
 
 //precalculated metric parameters
 ldouble rhorizonBL,rISCOBL,rmboundBL,rphotonBL,etaNT;
+int cells_under_horizon;
 
 // jet coordinate specific
 #if (MYCOORDS==JETCOORDS)
@@ -676,6 +677,7 @@ int my_finger(ldouble t);
 ///////////////////////////////////////////////////////////////
 
 void initialize_constants();
+int calc_cells_under_horiz();
 int print_scalings();
 int set_initial_profile();
 void am_i_sane();
@@ -1003,15 +1005,16 @@ int calc_primitives(int ix,int iy,int iz,int type,int setflags);
 int u2p(ldouble *uu0, ldouble *pp,void *ggg,int corrected[3],int fixups[2],int type);
 int check_floors_mhd(ldouble *pp, int whichvel,void *ggg);
 
-int u2p_solver(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
+int u2p_solver_mhd(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
 int u2p_solver_nonrel(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
 int u2p_solver_Wp(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
 int u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
 int u2p_solver_Bonly(ldouble *uu, ldouble *pp, void *ggg);
-int u2p_solver_ff(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
+int u2p_solver_ff(ldouble *uu, ldouble *pp, void *ggg,int verbose);
 int fill_ffprims();
 int fill_ffprims_cell(ldouble *pp, void *gg);
-
+ldouble get_driftvel_velr(ldouble *pp,ldouble *velff,void* ggg);
+
 int count_entropy(int *n, int *n2);
 int copy_entropycount();
 int update_entropy();
@@ -1334,6 +1337,8 @@ ldouble pick_ViscousHeating(int ix,int iy,int iz);
 int test_gammagas();
 int test_calcgamma();
 
+int calc_faraday(ldouble F[4][4],int ix,int iy,int iz,int when);
+int calc_current(void* ggg,ldouble jcon[4],int *derdir);
 ///////////////////////////////////////////////////////////////
 // nonthermal.c ///////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////

misc.c

@@ -190,10 +190,55 @@ void initialize_constants()
   */
 
   #endif
+
 
   return;
 }
 
+//**********************************************************************
+/*! \fn int calc_cells_under_horiz()
+ \brief calculate the number of cells beneath the horizon for bhdisk problems
+*/
+//**********************************************************************
+
+int
+calc_cells_under_horiz()
+{
+    //consistency check -- how many cell centers are under horizon? 
+#ifdef BHDISK_PROBLEMTYPE
+  ldouble xx[4],xxBL[4];
+  int ix;
+  if(TOI==0 && TOJ==0 && TOK==0) // only do on 0th tile
+  {
+    cells_under_horizon=0;
+    for(ix=0;ix<NX;ix++)
+    {
+      get_xx(ix,0,0,xx); //BH problem types should have r coord independent of theta,phi
+      coco_N(xx,xxBL,MYCOORDS,BLCOORDS);
+      if(xxBL[1]>rhorizonBL)
+	break;
+      else
+	cells_under_horizon+=1;
+    }
+
+    if(cells_under_horizon<3)
+    {
+      printf("There are only %d cells under horizon at rh=%.2f! increase to at least 4\n",
+	     cells_under_horizon,rhorizonBL);
+      exit(-1);
+    }
+    if(cells_under_horizon==NX)
+    {
+      printf("All cells on inner tile are under horizon! \n");
+      exit(-1);
+    }
+
+    printf("There are %d cells under the horizon at rh=%.2f\n",
+	   cells_under_horizon,rhorizonBL);
+  }
+#endif
+  return 0;
+}
 
 //**********************************************************************
 /*! \fn int print_scalings()
@@ -358,7 +403,6 @@ if (NTZ % 2 != 0)
 #endif
 
 #ifdef FORCEFREE
-  printf("WARNING! FORCEFREE is untested and does not account for current sheets! use only with monopole!\n");
   #ifdef NONRELMHD
   printf("FORCEFREE does not work with NONRELMHD!\n");
   exit(-1);
@@ -539,6 +583,8 @@ if (NTZ % 2 != 0)
   if(PROCID==0) printf("RESTARTFROMMHD\n");
   if(PROCID==0) printf("urad/uu: %e ue/uu: %e\n", INITURADFRAC, INITUEFRAC);
 #endif
+
+
   return;
 }
 
@@ -1476,7 +1522,7 @@ print_primitives(ldouble *p)
   printf("B^2 = %.15e\n",p[B2]);
   printf("B^3 = %.15e\n",p[B3]);
 #ifdef FORCEFREE
-  printf("uu (ff) = %.15e\n", p[UUFF]);
+  printf("upar (ff) = %.15e\n", p[UUFF]);
   printf("u^1 (ff) = %.15e\n",p[VXFF]);
   printf("u^2 (ff) = %.15e\n",p[VYFF]);
   printf("u^3 (ff) = %.15e\n",p[VZFF]);
@@ -1525,6 +1571,7 @@ print_conserved(ldouble *u)
   printf("B^2 = %.15e\n",u[B2]);
   printf("B^3 = %.15e\n",u[B3]);
 #ifdef FORCEFREE
+  printf("mu b^0 (ff) = %.15e\n",u[UUFF]);
   printf("T^t_1 (ff) = %.15e\n",u[VXFF]);
   printf("T^t_2 (ff) = %.15e\n",u[VYFF]);
   printf("T^t_3 (ff) = %.15e\n",u[VZFF]);