way too many updates

Makefile

@@ -20,7 +20,7 @@ endif
 LIBS=-lm -lgsl -lgslcblas -lsiloh5 -lfftw3 -lrt -lhdf5_serial
 
 RM=/bin/rm
-OBJS = mpi.o u2prad.o magn.o silo.o postproc.o fileop.o misc.o physics.o finite.o problem.o metric.o relele.o rad.o opacities.o u2p.o frames.o p2u.o nonthermal.o
+OBJS = mpi.o u2prad.o magn.o silo.o postproc.o fileop.o misc.o physics.o finite.o problem.o metric.o relele.o rad.o opacities.o u2p.o u2p_ff.o frames.o p2u.o nonthermal.o
 
 all: ko ana avg outavg phisli thsli phiavg regrid dumps2hdf5
 

PROBLEMS/MONOPOLE_2D/bc.c

@@ -112,16 +112,19 @@ pick_G(ix,iy,iz,GG);
 //in 3D polar cells overwritten with #define CORRECT_POLARAXIS_3D
 if(BCtype==YBCLO) //upper spin axis 
   {      
-
-    iiy=-iy-1;
+    iiy=-1*iy-1; 
     iiz=iz;
     iix=ix;
     gdet_src=get_g(g,3,4,iix,iiy,iiz);  
     gdet_bc=get_g(g,3,4,ix,iy,iz);  
     for(iv=0;iv<NV;iv++)
       {
 	//v_theta
+#ifdef FORCEFREE
+        if(iv==VY || iv==B2 || iv==VYFF || iv==FY0)
+#else
 	if(iv==VY || iv==B2 || iv==FY0)
+#endif
 	  pp[iv]=-get_u(p,iv,iix,iiy,iiz);
 	else
 	  pp[iv]=get_u(p,iv,iix,iiy,iiz);
@@ -142,7 +145,11 @@ if(BCtype==YBCHI) //lower spin axis
 
     for(iv=0;iv<NV;iv++)
       {
+#ifdef FORCEFREE
+        if(iv==VY || iv==B2 || iv==VYFF || iv==FY0)
+#else
 	if(iv==VY || iv==B2 || iv==FY0)
+#endif
 	  pp[iv]=-get_u(p,iv,iix,iiy,iiz);
 	else
 	  pp[iv]=get_u(p,iv,iix,iiy,iiz);

PROBLEMS/MONOPOLE_2D/define.h

@@ -3,6 +3,7 @@
 /************************************/
 #define RESTART
 #define RESTARTNUM -1
+#define BHDISK_PROBLEMTYPE
 
 //#define RADIATION
 
@@ -13,8 +14,11 @@
 #define MAGNFIELD
 #define GDETIN 1    //must be 1 for MAGNFIELD
 #define VECPOTGIVEN
+#define MONOPOLE_FIELD_CORNERS
 //#define INIT_MAGN_CORNERS
 
+//#define U2PCONV 1.e-14
+
 #define FORCEFREE
 //#define NOLOGINS 
 
@@ -32,35 +36,37 @@
 //rmhd floors
 /************************************/
 
-#define B2RHORATIOMAXINIT 500 
-#define B2UURATIOMAXINIT 500
-//#define SIGMAWCONSTINIT 1.e4
+#define B2RHORATIOMAXINIT 200
+#define B2UURATIOMAXINIT 200
 
 #if defined(FORCEFREE)
 
-//DIFTFRAME not compatible with FORCEFREE yet
-#define B2RHOFLOORFRAME FFFRAME // ZAMOFRAME 
+#define HYBRID_FORCEFREE
+#define HYBRID_FORCEFREE_SIGMACUT 50
+#define FORCEFREE_SOLVE_PARALLEL
+//#define FORCEFREE_SOLVE_PARALLEL_OLD
 
+//#define SKIPALLFLOORS // TODO seems critical for SOLVE_PARALLEL? 
+#define B2RHOFLOORFRAME ZAMOFRAME //DRIFTFRAME
 //#define CORRECT_POLARAXIS
 //#define NCCORRECTPOLAR 2
 
 #define UURHORATIOMIN 0.
-#define UURHORATIOMAX 50.
+#define UURHORATIOMAX 50. 
 #define B2UURATIOMIN 0.
 #define B2UURATIOMAX 1.e10
 #define B2RHORATIOMIN 0.
 #define B2RHORATIOMAX 1.e10
 
-#define GAMMAMAXFF 1000. //lower than GAMMAMAXHD? 
+#define GAMMAMAXFF 1000.  //lower than GAMMAMAXHD? 
 #define GAMMAMAXHD 10000. //why can't this be pushed higher on the monopole? 
 
-#define ALLOWENTROPYU2P 0
-
 #else
+
 #define CORRECT_POLARAXIS
 #define NCCORRECTPOLAR 2
 
-#define B2RHOFLOORFRAME DRIFTFRAME //ZAMOFRAME 
+#define B2RHOFLOORFRAME ZAMOFRAME 
 #define UURHORATIOMIN 0.
 #define UURHORATIOMAX 50.
 #define B2UURATIOMIN 0.
@@ -90,8 +96,8 @@
 #define MYCOORDS MKS1COORDS
 #define MINX (log(RMIN - MKSR0))
 #define MAXX (log(RMAX - MKSR0))
-#define MINY 0.01*M_PI/2.
-#define MAXY M_PI - 0.01*M_PI/2.
+#define MINY 0.001*M_PI
+#define MAXY M_PI - 0.001*M_PI
 #define TNX 128
 #define TNY 64//256
 #define TNZ 1
@@ -139,7 +145,7 @@
 #define OUTVEL VEL4
 #define ALLSTEPSOUTPUT 0
 #define NSTEPSTOP 1.e10
-#define NOUTSTOP 100
+#define NOUTSTOP 102
 #define SILOOUTPUT 1
 #define PRIMOUTPUT 1
 //#define PRIMOUTPUTINMYCOORDS

PROBLEMS/MONOPOLE_2D/init.c

@@ -6,7 +6,7 @@ struct geometry geom;
 fill_geometry(ix,iy,iz,&geom);
 
 struct geometry geomBL;
-fill_geometry_arb(ix,iy,iz,&geomBL,BLCOORDS);
+fill_geometry_arb(ix,iy,iz,&geomBL,KSCOORDS);
 
 ldouble r = geomBL.xx;
 ldouble th = geomBL.yy;
@@ -17,11 +17,13 @@ uint = UINTATMMIN + (r/10./rhor)/pow(r,4)/B2UURATIOMAXINIT;
 
 pp[0]=rho;
 pp[1]=uint;
-pp[2]=0.; //zero initial velocity in co-rotating frame
+pp[2]=0.; 
 pp[3]=0.;
 pp[4]=0.;
 
+// These will be made consistent with B direction after init
 #ifdef FORCEFREE
+pp[UUFF]=uint;
 pp[VXFF]=0.;
 pp[VYFF]=0.;
 pp[VZFF]=0.;
@@ -35,6 +37,7 @@ pp[B1]=pp[B2]=pp[B3]=0.;
 #endif
 
 //transform primitives from BL to MYCOORDS
+trans_pmhd_coco(pp, pp, KSCOORDS,MYCOORDS, geomBL.xxvec,&geomBL,&geom);
 //trans_pall_coco(pp, pp, BLCOORDS, MYCOORDS,xxvecBL,&geomBL,&geom);
 
 #ifdef MAGNFIELD //calculate vector potential

choices.h

@@ -218,7 +218,7 @@
 //number of magneto-hydro variables
 #ifdef MAGNFIELD
 #ifdef FORCEFREE
-#define NVMHD (NVHD+6)
+#define NVMHD (NVHD+7)
 #else
 #define NVMHD (NVHD+3)
 #endif 

frames.c

@@ -51,8 +51,8 @@ trans_pmhd_coco(ldouble *ppin, ldouble *ppout, int CO1,int CO2, ldouble *xxvec,
   }
   else
   {
-    pp2[0]=pp1[0];
-    pp2[1]=pp1[1];
+    pp2[RHO]=pp1[RHO];
+    pp2[UU]=pp1[UU];
 
     //bcon in CO1
     ldouble ucon[4], ucov[4],uconback[4];    
@@ -80,10 +80,10 @@ trans_pmhd_coco(ldouble *ppin, ldouble *ppout, int CO1,int CO2, ldouble *xxvec,
     //to VELPRIM
     conv_vels_ut(ucon,ucon,VEL4,VELPRIM,geom2->gg,geom2->GG);
 
-    pp2[2]=ucon[1];
-    pp2[3]=ucon[2];
-    pp2[4]=ucon[3];
-    
+    pp2[VX]=ucon[1];
+    pp2[VY]=ucon[2];
+    pp2[VZ]=ucon[3];
+
 #ifdef MAGNFIELD
 
     //coming back to primitive B^i
@@ -93,7 +93,26 @@ trans_pmhd_coco(ldouble *ppin, ldouble *ppout, int CO1,int CO2, ldouble *xxvec,
     {
       pp2[B1+i] = Bcon[1+i];
     }
+
+#ifdef FORCEFREE
+    pp2[UUFF]=pp1[UUFF];
+
+    ldouble uconff[4];
+    uconff[0]=0.;
+    uconff[1]=pp1[VXFF];
+    uconff[2]=pp1[VYFF];
+    uconff[3]=pp1[VZFF];
+    conv_vels(uconff,uconff,VELPRIM,VEL4,geom1->gg,geom1->GG);
+    trans2_coco(xxvec,uconff,uconff,CO1,CO2);
+    conv_vels(uconff,uconff,VEL4,VELPRIM,geom2->gg,geom2->GG);
+    pp2[VXFF]=uconff[1];
+    pp2[VYFF]=uconff[2];
+    pp2[VZFF]=uconff[3];
+#endif
+
 #endif
+
+
   }
 
   for(i=0;i<NVMHD;i++)
@@ -1651,6 +1670,21 @@ trans_pmhd_coco_precompute(ldouble *ppin, ldouble *ppout, void* ggg1,void* ggg2,
     {
       pp2[B1+i] = Bcon[1+i];
     }
+
+#ifdef FORCEFREE
+    ldouble uconff[4];
+    uconff[0]=0.;
+    uconff[1]=pp1[VXFF];
+    uconff[2]=pp1[VYFF];
+    uconff[3]=pp1[VZFF];
+    conv_vels(uconff,uconff,VELPRIM,VEL4,geom1->gg,geom1->GG);
+    trans2_coco_precompute(uconff,uconff,geom1->ix,geom1->iy,geom1->iz,which);
+    conv_vels(uconff,uconff,VEL4,VELPRIM,geom2->gg,geom2->GG);
+    pp2[VXFF]=uconff[1];
+    pp2[VYFF]=uconff[2];
+    pp2[VZFF]=uconff[3];
+#endif
+
 #endif
   }
 

ko.c

@@ -114,10 +114,11 @@ main(int argc, char **argv)
   //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\n",RHO,UU,VX,VY,VZ,ENTR,B1,B2,B3,VXFF,VYFF,VZFF);
+  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();
-  //exit(1);
 
   //**************
   //tests
@@ -245,6 +246,9 @@ main(int argc, char **argv)
       }
       calc_BfromA(p,1);
 
+#ifdef FORCEFREE
+      update_ffprims(); // make force-free primitives consistent
+#endif
       //exchange magn. field calculated in domain
       mpi_exchangedata();
       calc_avgs_throughout();

ko.h

@@ -897,23 +897,6 @@ int trans_prad_coco(ldouble *ppin, ldouble *ppout, int CO1,int CO2, ldouble *xxv
 
 //deprecated
 int prad_ff2lab(ldouble *pp1, ldouble *pp2, void* ggg);
-//int prad_lab2ff(ldouble *pp1, ldouble *pp2, void *ggg);
-//int prad_on2lab(ldouble *pp1, ldouble *pp2, void* ggg);
-//int prad_lab2on(ldouble *pp1, ldouble *pp2, void *ggg);
-//int prad_ff2zamo(ldouble *pp1, ldouble *pp2, ldouble gg[][5], ldouble GG[][5], ldouble eup[][4]);
-//int prad_zamo2ff(ldouble *pp1, ldouble *pp2, ldouble gg[][5], ldouble GG[][5], ldouble eup[][4]);
-//int boost2_zamo2ff(ldouble* A1,ldouble* A2,ldouble *pp,ldouble gg[][5],ldouble GG[][5],ldouble eup[][4]);
-//int boost2_ff2zamo(ldouble A1[4],ldouble A2[4],ldouble *pp,ldouble gg[][5],ldouble GG[][5],ldouble eup[][4]);
-//int boost22_zamo2ff(ldouble T1[][4],ldouble T2[][4],ldouble *pp,ldouble gg[][5],ldouble GG[][5],ldouble eup[][4]);
-//int boost22_ff2zamo(ldouble T1[][4],ldouble T2[][4],ldouble *pp,ldouble gg[][5],ldouble GG[][5],ldouble eup[][4]);
-//int trans22_zamo2lab(ldouble T1[][4],ldouble T2[][4],ldouble elo[][4]);
-//int trans22_lab2zamo(ldouble T1[][4],ldouble T2[][4],ldouble eup[][4]);
-//int trans2_lab2zamo(ldouble *u1,ldouble *u2,ldouble eup[][4]);
-//int trans2_zamo2lab(ldouble *u1,ldouble *u2,ldouble elo[][4]);
-//int trans22_on2cc(ldouble T1[][4],ldouble T2[][4],ldouble tlo[][4]);
-//int trans22_cc2on(ldouble T1[][4],ldouble T2[][4],ldouble tup[][4]);
-//int trans2_cc2on(ldouble *u1,ldouble *u2,ldouble tup[][4]);
-//int trans2_on2cc(ldouble *u1,ldouble *u2,ldouble tlo[][4]);
 
 int calc_Lorentz_lab2ff(ldouble *pp,ldouble gg[][5],ldouble GG[][5],ldouble L[][4]);
 int calc_Lorentz_lab2ff_4vel(ldouble *pp, ldouble gg[][5], ldouble GG[][5], ldouble L[][4], ldouble ucon[4], ldouble ucov[4]);
@@ -1024,8 +1007,10 @@ int u2p_solver(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_ff(ldouble *uu, ldouble *pp, void *ggg,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 update_ffprims();
+int update_ffprims_cell(ldouble *pp, ldouble *uu, void *gg);
 
 int count_entropy(int *n, int *n2);
 int copy_entropycount();

misc.c

@@ -1476,6 +1476,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("u^1 (ff) = %.15e\n",p[VXFF]);
   printf("u^2 (ff) = %.15e\n",p[VYFF]);
   printf("u^3 (ff) = %.15e\n",p[VZFF]);

mnemonics.h

@@ -35,9 +35,10 @@
 #define B3 (NVHD+2)
 
 #ifdef FORCEFREE
-#define VXFF (NVHD+3)
-#define VYFF (NVHD+4)
-#define VZFF (NVHD+5)
+#define UUFF (NVHD+3)
+#define VXFF (NVHD+4)
+#define VYFF (NVHD+5)
+#define VZFF (NVHD+6)
 #endif
 
 #define EE (NVMHD)