added force free velocities as additional primitives

Makefile

@@ -2,7 +2,7 @@
 
 ifneq ($(SERIAL),1)
 CC=mpicc 
-CFLAGS=-O3 -DMPI
+CFLAGS=-O3 -DMPI -fcommon
 
 else
 //CC=gcc
@@ -13,7 +13,7 @@ else
 -fsanitize=address -fno-omit-frame-pointer
 
 CC=/usr/bin/h5cc
-CFLAGS = -O3 -Wno-unused-result -I/usr/lib/gcc/x86_64-linux-gnu/5.4.0/include -I/usr/include/hdf5/serial -Wunused-function -fopenmp -fcommon
+CFLAGS = -O3 -Wno-unused-result -I/usr/lib/gcc/x86_64-linux-gnu/5.4.0/include -I/usr/include/hdf5/serial -Wunused-function -w -fopenmp -fcommon
 
 endif
 

PROBLEMS/MONOPOLE_2D/bc.c

@@ -2,6 +2,8 @@
 //int calc_bc(int ix,int iy,int iz,ldouble t,
 //	ldouble *uu,ldouble *pp,int ifinit,int BCtype)
 
+#define BCCOORDS KSCOORDS
+
 /**********************/
 //geometries
 ldouble gdet_src,gdet_bc,Fx,Fy,Fz,ppback[NV];
@@ -11,16 +13,16 @@ 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,BCCOORDS);
 
-ldouble gg[4][5],GG[4][5],ggsrc[4][5],eup[4][4],elo[4][4];
+ldouble gg[4][5],GG[4][5];
 pick_g(ix,iy,iz,gg);
 pick_G(ix,iy,iz,GG);
 
 /**********************/
 
 //radius
-if(BCtype==XBCHI) //outflow in magn, atm in rad., atm. in HD
+  if(BCtype==XBCHI) //outflow in magn, atm in rad., atm. in HD
   {
     iix=NX-1;
     iiy=iy;
@@ -32,54 +34,62 @@ if(BCtype==XBCHI) //outflow in magn, atm in rad., atm. in HD
 	pp[iv]=get_u(p,iv,iix,iiy,iiz);
       }
 
+    //ANDREW TODO what to do about outer BC
+    /*
     //!! begin rescale
     //first transform to BL
-    trans_pmhd_coco(pp, pp, MYCOORDS,BLCOORDS, geom.xxvec,&geom,&geomBL);
-
+    trans_pmhd_coco(pp, pp, MYCOORDS,BCCOORDS, geom.xxvec,&geom,&geomBL);
+
+    // boundary cell and scale factors
     struct geometry geombdBL;
-    fill_geometry_arb(iix,iiy,iiz,&geombdBL,BLCOORDS);
+    fill_geometry_arb(iix,iiy,iiz,&geombdBL,BCCOORDS);
     ldouble rghost = geomBL.xx;
     ldouble rbound = geombdBL.xx;
     ldouble deltar = rbound - rghost;
-    ldouble scale1 =  rbound*rbound/rghost/rghost;
-    ldouble scale2 = rbound/rghost;
-
-    pp[RHO]*=scale1;
-    pp[UU] *=scale1;
+    ldouble scale1 =  rbound/rghost; //r^-1
+    ldouble scale2 =  rbound*rbound/rghost/rghost; //r^-2
+    ldouble scale3 =  rbound*rbound*rbound/rghost/rghost/rghost; //r^-3 -- inital atm
+
+    //  bsq = br*br + r^2bth^2 + r^2*sin^2(th)*bph^2
+      
+    pp[RHO]*=scale3;//scale2;
+    pp[UU] *=scale3;//scale2;
     #ifdef MAGNFIELD
-    pp[B1] *=scale1;
-    pp[B2] *=scale2;
-    pp[B3] *=scale2;
+    pp[B1] *=scale2;
+    pp[B2] *=scale1;
+    pp[B3] *=scale1;
     #endif
-    //pp[VX] *=1.;
+    pp[VX] *=1.;
     pp[VY] *=scale1;
     pp[VZ] *=scale1;
     #ifdef RADIATION
-    pp[EE0] *=scale1;
-    //pp[FX0] *=1.;
-    pp[FY0] *=scale1;
-    pp[FZ0] *=scale1;
+    pp[EE0] *=scale2;
+    pp[FX0] *=1.;
+    pp[FY0] *=scale2;
+    pp[FZ0] *=scale2;
     #endif
     //transform back after rescaling
-    trans_pmhd_coco(pp, pp,BLCOORDS, MYCOORDS, geomBL.xxvec,&geomBL, &geom);
-    //!! end rescale
+    trans_pmhd_coco(pp, pp,BCCOORDS, MYCOORDS, geomBL.xxvec,&geomBL, &geom);
 
-    //checking for the gas inflow
+    //!! end rescale
+    /*
+    //check for gas inflow
     ldouble ucon[4]={0.,pp[VX],pp[VY],pp[VZ]};    
     conv_vels(ucon,ucon,VELPRIM,VEL4,geom.gg,geom.GG);
-    trans2_coco(geom.xxvec,ucon,ucon,MYCOORDS,BLCOORDS);
-    if(ucon[1]<0.) //inflow, resseting to atmosphere
+    trans2_coco(geom.xxvec,ucon,ucon,MYCOORDS,BCCOORDS);
+    if(ucon[1]<0.) //inflow, reset to zero velocity
       {
 	//atmosphere in rho,uint and velocities and zero magn. field
 	//set_hdatmosphere(pp,xxvec,gg,GG,4);
 	ucon[1]=0.;
-	trans2_coco(geomBL.xxvec,ucon,ucon,BLCOORDS,MYCOORDS);
+	trans2_coco(geomBL.xxvec,ucon,ucon,BCCOORDS,MYCOORDS);
 	conv_vels(ucon,ucon,VEL4,VELPRIM,geom.gg,geom.GG);
 	pp[VX]=ucon[1];
 	pp[VY]=ucon[2];
 	pp[VZ]=ucon[3];//atmosphere in rho,uint and velocities and zero magn. field
       }
-
+    */
+
     p2u(pp,uu,&geom);
     return 0;  
   }
@@ -91,7 +101,7 @@ if(BCtype==XBCHI) //outflow in magn, atm in rad., atm. in HD
 
       for(iv=0;iv<NV;iv++)
        {
-	 pp[iv]=get_u(p,iv,0,iiy,iiz);
+	 pp[iv]=get_u(p,iv,iix,iiy,iiz);
        }
 
      p2u(pp,uu,&geom);

PROBLEMS/MONOPOLE_2D/define.h

@@ -2,83 +2,113 @@
 //restart
 /************************************/
 #define RESTART
-#define RESTARTNUM 44
+#define RESTARTNUM -1
 
-/************************************/
-//radiation choices
-/************************************/
 //#define RADIATION
-//#define SKIPRADSOURCE
+
 
 /************************************/
 //magnetic choices
 /************************************/
-//#define MAGNFIELD
-//#define VECPOTGIVEN
-//#define MONOPOLE_FIELD_CORNERS // generates monopole field potential directly on corners
-                               // no cell averaging
+#define MAGNFIELD
+#define GDETIN 1    //must be 1 for MAGNFIELD
+#define VECPOTGIVEN
+//#define INIT_MAGN_CORNERS //initialize magnetic field on corners/edges (which?)
+
+//#define FORCEFREE
 
-//#define MAXBETA 0.01 //close to the target pmag/pgas
+//#define NOLOGINS 
 
 /************************************/
 //reconstruction / Courant
 /************************************/
-#define INT_ORDER 1
+#define INT_ORDER 1 //TODO why is int_order 1 more stable???
 #define TIMESTEPPING RK2IMEX
-#define TSTEPLIM .8
-#define FLUXLIMITER 0
+#define TSTEPLIM .9
+#define FLUXLIMITER 1
 #define FLUXMETHOD LAXF_FLUX
 #define MINMOD_THETA 1.5
 
-/************************************/
-//viscosity choices
-/************************************/
-#define HDVISCOSITY NOVISCOSITY
-#define RADVISCOSITY NOVISCOSITY
-
 /************************************/
 //rmhd floors
 /************************************/
-#define UURHORATIOMIN 1.e-15
+
+#define B2RHORATIOMAXINIT 500 
+#define B2UURATIOMAXINIT 500
+//#define SIGMAWCONSTINIT 1.e4
+
+#if defined(FORCEFREE)
+
+
+//DIFTFRAME not compatible with FORCEFREE yet
+#define B2RHOFLOORFRAME FFFRAME // ZAMOFRAME 
+
+//#define CORRECT_POLARAXIS
+//#define NCCORRECTPOLAR 2
+
+#define UURHORATIOMIN 0.
+#define UURHORATIOMAX 50.
+#define B2UURATIOMIN 0.
+#define B2UURATIOMAX 1.e10
+#define B2RHORATIOMIN 0.
+#define B2RHORATIOMAX 1.e10
+
+#define GAMMAMAXFF 1000. //lower than GAMMAMAXHD? 
+#define GAMMAMAXHD 1000. //why can't this be pushed higher on the monopole? 
+
+#define ALLOWENTROPYU2P 0
+
+#else
+#define CORRECT_POLARAXIS
+#define NCCORRECTPOLAR 2
+
+#define B2RHOFLOORFRAME ZAMOFRAME 
+#define UURHORATIOMIN 0.
 #define UURHORATIOMAX 50.
-//#define EERHORATIOMIN 1.e-20
-//#define EERHORATIOMAX 1.e20
-//#define EEUURATIOMIN 1.e-20
-//#define EEUURATIOMAX 1.e20
 #define B2UURATIOMIN 0.
 #define B2UURATIOMAX 2500.
 #define B2RHORATIOMIN 0.
-#define B2RHORATIOMAX 50.
-#define GAMMAMAXRAD 50.
+#define B2RHORATIOMAX 500.
+#define GAMMAMAXHD 50.
+#endif
 
 /************************************/
 //blackhole
 /************************************/
 #define MASS 10.
-//#define BHSPIN 0.
-#define BHSPIN 0.9375
+#define BHSPIN 0.5 // 0.9375
 #define RHOR (1.+sqrt(1. - BHSPIN*BHSPIN))
+
 /************************************/
 //coordinates / resolution
 /************************************/
-#define myMKS2COORDS
-#define MKSR0 -2
-#define MKSH0 0.8
+#define myMKS1COORDS
 #define METRICAXISYMMETRIC
 #define RMIN 0.7*RHOR //1.8<->6 ANDREW
-#define RMAX 1000.
+#define RMAX 100.
 #define MKS1R0 MKSR0
 
-#ifdef myMKS2COORDS //modified Kerr-Shild
+#ifdef myMKS1COORDS //modified Kerr-Shild
+#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 TNX 128//256
+#define TNY 128//256
+#define TNZ 1
+
+#elif defined(myMKS2COORDS) //modified Kerr-Shild
 #define MYCOORDS MKS2COORDS
+#define MKSR0 -4.*RHOR
+#define MKSH0 0.8
+
 #define MINX (log(RMIN - MKSR0))
 #define MAXX (log(RMAX - MKSR0))
-//#define MINY (0.1*Pi/2.)
-//#define MAXY (Pi-0.1*Pi/2.)
 #define MINY (0.001)
 #define MAXY (1.-0.001)
-#define TNX 512
-#define TNY 40
+#define TNX 128
+#define TNY 128
 #define TNZ 1
 
 #else //Schwarzschild
@@ -95,42 +125,34 @@
 #define PHIWEDGE (M_PI/2.)
 #define MINZ (-PHIWEDGE/2.)
 #define MAXZ (PHIWEDGE/2.)
-//#define MINZ -1.
-//#define MAXZ 1.
 
 #define SPECIFIC_BC
-#define PEROIDIC_ZBC
+#define PERIODIC_ZBC
 
 //number of tiles
 #define NTX 4//48
-#define NTY 1//24
+#define NTY 2//24
 #define NTZ 1
 
 /************************************/
 //output
 /************************************/
-#define OUTCOORDS KERRCOORDS    
+#define OUTCOORDS KSCOORDS //MKS1COORDS
 #define OUTVEL VEL4
 #define ALLSTEPSOUTPUT 0
 #define NSTEPSTOP 1.e10
-#define NOUTSTOP 5000
+#define NOUTSTOP 100
 #define SILOOUTPUT 1
+#define PRIMOUTPUT 1
+//#define PRIMOUTPUTINMYCOORDS
 #define SCAOUTPUT 0
 #define SILO2D_XZPLANE
 
 /************************************/
 //common physics / atmosphere
 /************************************/
 #define GAMMA (4./3.)
-#define NODONUT 0
-#define INFLOWING 0
-#define ELL 3.85
-#define URIN (0.)
-#define KKK 7127.
-#define UTPOT .965
-#define DTOUT1 1.e0
+#define DTOUT1 1
 #define DTOUT2 1.e0
 #define RHOATMMIN  1.e-20
 #define UINTATMMIN 1.e-20
-//#define UINTATMMIN  (calc_PEQ_ufromTrho(1.e10,RHOATMMIN))
-//#define ERADATMMIN  (calc_LTE_EfromT(3.e6)/10)