fixed jetcoords -- still issues with cylindrification+rad

PROBLEMS/ADAFCRITRELEL/define.h

@@ -15,13 +15,13 @@
 /************************************/
 //radiation choices
 /************************************/
-//#define RADIATION
+#define RADIATION
 //#define SKIPRADSOURCE    //advective only
 //#define SKIPRADEVOLUTION //keeps initial values
 
 #ifdef RADIATION
 //#define BASICRADIMPLICIT
-//#define RADIMPSTOPWHENFAIL
+#define RADIMPSTOPWHENFAIL
 //#define RADIMP_START_WITH_BISECT
 //#define BALANCEENTROPYWITHRADIATION
 #define ALLOWRADCEILINGINIMPLICIT
@@ -81,6 +81,8 @@
 #define DONOTLIMITENTRINMIXING
 //#define MIXENTROPIES_CONST_PRESSURE
 
+
+#define DIVIDEVISCHEATBYDT
 #define PRINTVISCHEATINGTOSILO
 #define PRINTCOULOMBTOSILO
 #define CALCVISCHEATING
@@ -233,53 +235,26 @@
 #define BHSPIN 0.9375
 
 /************************************/
-//resolution
-/************************************/
-//total resolution
-#define TNX 64//320//288 //256 //16*16 //128 <->256
-#define TNY 64//256//224 //320 //14*9 //92 <->256
-#define TNZ 1//128//96
-
-//number of tiles
-#define NTX 32
-#define NTY 16//32
-#define NTZ 8
-
-/************************************/
-//coordinates
+//coordinates / resolution
 /************************************/
 #define myJETCOORDS
 //#define myMKS3COORDS
-//#define myMKS2COORDS
 #define METRICAXISYMMETRIC
 #define RMIN 1.
 #define RMAX 1.e5
-
 #define MKSR0 -1.35
 #define MKSH0 0.7
 #define MKSMY1 0.002
 #define MKSMY2 0.02
 #define MKSMP0 1.3
 
-#define HYPRBRK 5000
-#define FJET 0.4
-#define FDISK 0.5
-#define RUNI RMIN
-#define RCOLL_JET 1000
-#define RDECOLL_JET 2*RMIN
-#define RCOLL_DISK 5*RMIN
-#define RDECOLL_DISK 2*RMIN
-#define ALPHA_1 1
-#define ALPHA_2 0.375
-
-#ifdef myMKS2COORDS //modified Kerr-Shild
-#define METRICNUMERIC
-#define MYCOORDS MKS2COORDS
-#define MINX (log(RMIN-MKSR0))
-#define MAXX (log(RMAX-MKSR0))
-#define MINY (0.001)
-#define MAXY (1.-0.001)
-#endif
+//#ifdef myMKS2COORDS //modified Kerr-Shild
+//#define MYCOORDS MKS2COORDS
+//#define MINX (log(RMIN-MKSR0))
+//#define MAXX (log(RMAX-MKSR0))
+//#define MINY (0.001)
+//#define MAXY (1.-0.001)
+//#endif
 
 #ifdef myMKS3COORDS //modified Kerr-Shild further from axis
 #define METRICNUMERIC
@@ -290,58 +265,87 @@
 #define MAXY 1.
 #endif
 
+
 #ifdef myJETCOORDS //concentrate resolution in jet and disk zones
 #define MYCOORDS JETCOORDS
-//#define CYLINDRIFY
+
 #define METRICNUMERIC
-#define DERIVS_NOGSL
+#define DERIVS_NOGSL // use a faster numeric derivative for coordinate transformations
+#define PRECOMPUTE_MY2OUT // precompute transformation matrices from BL <--> JET
+
 #define MINX 0
 #define MAXX 1.
-#define MINY -(1.-1.e-8) //10^-9 away from the poles seems to be the last safe point
-#define MAXY 1.-1.e-8    //TODO -- can we fix this!??!?!?
-#define RCYL 20
-#define NCYL 1
+#define MINY -(1.-1.e-2) //10^-8 away from the poles seems to be the last safe point
+#define MAXY 1.-1.e-2  
+
+#define MKSR0 -1.35
+#define HYPRBRK 5000
+#define FJET 0.3
+#define FDISK 0.4
+#define RUNI RMIN
+#define RCOLL_JET 1000
+#define RDECOLL_JET 2*RMIN//2*RMIN
+#define RCOLL_DISK 20*RMIN//5*RMIN
+#define RDECOLL_DISK 2*RMIN//2*RMIN
+#define ALPHA_1 1
+#define ALPHA_2 0.375
+
+//#define CYLINDRIFY
+#define RCYL 5
+#define NCYL 1.//0.5
 #endif
 
+
 #define PHIWEDGE (2*M_PI)
 #define MINZ (-PHIWEDGE/2.)
 #define MAXZ (PHIWEDGE/2.)
 
-#define COPY_XBC
-//#define COPY_YBC
+//total resolution
+#define TNX 32//288 //256 //16*16 //128 <->256
+#define TNY 32//224 //320 //14*9 //92 <->256
+#define TNZ 1//128//96
+
+//number of tiles
+#define NTX 32
+#define NTY 16//32
+#define NTZ 8
+
 #define SPECIFIC_BC
 #define PERIODIC_ZBC
+#define COPY_XBC
+
 //#define PERIODIC_XBC
 //#define PERIODIC_YBC
 
 /************************************/
 //output
 /************************************/
 #define DTOUT1 10 //res
-#define DTOUT2 100 //avg
+#define DTOUT2 1000 //avg
 #define DTOUT3 100 //box,var
 
 #define OUTCOORDS BLCOORDS//KERRCOORDS                                                                    
 #define OUTVEL VEL4
 #define ALLSTEPSOUTPUT 0
-//#define RADOUTPUTINZAMO
+
 #define NSTEPSTOP 1.e20
-#define NOUTSTOP 5
+#define NOUTSTOP 1250
 
 //#define OUTPUTPERCORE
-#define COORDOUTPUT 1
-#define GRIDOUTPUT 1
+#define COORDOUTPUT 0
 #define SILOOUTPUT 1
 #define CGSOUTPUT
 #define OUTOUTPUT 0
-#define SIMOUTPUT 0//2
+#define SIMOUTPUT 2
 //#define SIMOUTPUT_GILAB2FF
 //#define SIMOUTPUTWITHINDTHETA .05
 #define GRTRANSSIMOUTPUT
+//#define GRTRANSSIMOUTPUT_2
+
 
-#define RADOUTPUT 0
+#define RADOUTPUT 1
 //#define RADOUTPUTWITHINDTHETA .05
-#define SCAOUTPUT 0
+#define SCAOUTPUT 1
 #define AVGOUTPUT 0
 #define NORELELAVGS
 #define THOUTPUT 0

PROBLEMS/JETCOORDS/bc.c

@@ -4,9 +4,11 @@
 
 /**********************/
 //geometries
-ldouble gdet_src,gdet_bc,Fx,Fy,Fz,ppback[NV];
+ldouble gdet_src,gdet_bc,Fx,Fy,Fz,ppback[NV],ppback2[NV],pp1[NV],pp2[NV];
 int iix,iiy,iiz,iv;  	  
 
+int iiytest=32;
+
 struct geometry geom;
 fill_geometry(ix,iy,iz,&geom);
 
@@ -48,18 +50,36 @@ if(BCtype==XBCHI) //outflow in magn, atm in rad., atm. in HD
     for(iv=0;iv<NV;iv++)
       {
 	pp[iv]=get_u(p,iv,iix,iiy,iiz);
+	ppback[iv]=pp[iv];
       }
 
     //!! begin rescale
 
     // AA -- example of how to use precomputed my2out and out2my for faster coco
     //first transform to BL
     #ifdef PRECOMPUTE_MY2OUT
-    trans_pmhd_coco_my2out(pp,pp,&geom,&geomBL);
+    trans_pmhd_coco_my2out(pp,pp, &geom, &geomBL);
     #else      
     trans_pmhd_coco(pp, pp, MYCOORDS,BLCOORDS, geom.xxvec,&geom,&geomBL);
     #endif
 
+    /*
+    if(iiy==iiytest)
+    {
+
+      printf("%d %d %d\n",iiy,geom.iy,geomBL.iy);
+      printf("\n");
+      trans_pmhd_coco_my2out(ppback, pp1, &geom, &geomBL);
+      trans_pmhd_coco(ppback, pp2, MYCOORDS,BLCOORDS, geom.xxvec,&geom,&geomBL);
+
+      for(iv=0;iv<NV;iv++) printf("%d %e | %e %e %e\n",iv, ppback[iv],pp1[iv],pp2[iv],pp[iv]);
+      printf("\n");
+
+      for(iv=0;iv<NV;iv++) ppback[iv]=pp[iv];
+      //exit(-1);
+    }
+    */
+
     struct geometry geombdBL;
     fill_geometry_arb(iix,iiy,iiz,&geombdBL,BLCOORDS);
     ldouble rghost = geomBL.xx;
@@ -86,26 +106,42 @@ if(BCtype==XBCHI) //outflow in magn, atm in rad., atm. in HD
     #endif
     //transform back after rescaling
 
+
+    /*
+    if(iiy==iiytest)
+    {
+      for(iv=0;iv<NV;iv++) ppback2[iv]=pp[iv];
+      
+      trans_pmhd_coco_out2my(ppback, pp1, &geomBL, &geom);
+      trans_pmhd_coco(ppback, pp2, BLCOORDS,MYCOORDS, geomBL.xxvec,&geomBL,&geom);
+
+      printf("scale1: %e scale2: %e\n",scale1,scale2); 
+      for(iv=0;iv<NV;iv++) printf("%d %e %e | %e %e\n",iv, ppback[iv],ppback2[iv],pp1[iv],pp2[iv]);
+      printf("\n");
+      exit(-1);
+    }
+    */
+
     #ifdef PRECOMPUTE_MY2OUT
     trans_pmhd_coco_out2my(pp,pp,&geomBL,&geom);
     #else      
     trans_pmhd_coco(pp, pp, BLCOORDS,MYCOORDS, geomBL.xxvec,&geomBL,&geom);
     #endif
-    
+
     //!! end rescale
 
-//checking for the gas inflow
+    //checking for the gas inflow
     ldouble ucon[4]={0.,pp[VX],pp[VY],pp[VZ]}; 
     conv_vels(ucon, ucon, VELPRIM, VEL4, geom.gg, geom.GG);
     if(MYCOORDS!=CYLCOORDS)
-      {
+    {
       #ifdef PRECOMPUTE_MY2OUT
       trans2_coco_my2out(ucon, ucon, geom.ix, geom.iy, geom.iz);
       #else      
       trans2_coco(geom.xxvec,ucon,ucon,MYCOORDS,BLCOORDS);
       #endif
-      }
-    if(ucon[1]<0.) //inflow, resseting to atmosphere
+    }
+    if(ucon[1]<0.) //inflow, resetting to atmosphere
       {
 	//atmosphere in rho,uint and velocities and zero magn. field
 	//set_hdatmosphere(pp,xxvec,gg,GG,4);