changes for madcc and jetcoords

choices.h

@@ -515,7 +515,7 @@
 
 //number of profiles
 #ifndef NRADPROFILES
-#define NRADPROFILES (66-1)
+#define NRADPROFILES (68-1)
 #endif
 
 #ifndef NTHPROFILES
@@ -715,6 +715,7 @@
 #ifdef RADIATION
 #ifndef NO_COMPTONIZATION
 #define COMPTONIZATION
+#define COMPTONIZATIONFLAG
 #endif
 #endif
 

finite.c

@@ -5419,7 +5419,7 @@ smooth_polaraxis()
   //spherical like coords
   if (MYCOORDS==SCHWCOORDS || MYCOORDS==KSCOORDS || MYCOORDS==KERRCOORDS || MYCOORDS==SPHCOORDS ||
       MYCOORDS==MKS1COORDS || MYCOORDS==MKS2COORDS || MYCOORDS==MKS3COORDS || MYCOORDS==TKS3COORDS ||
-      MYCOORDS==MSPH1COORDS)
+      MYCOORDS==MSPH1COORDS || MYCOORDS==JETCOORDS)
     {
       int ix;
 //#pragma omp parallel for private(ix)
@@ -5652,7 +5652,7 @@ correct_polaraxis()
       //spherical like coords
       if (MYCOORDS==SCHWCOORDS || MYCOORDS==KSCOORDS || MYCOORDS==KERRCOORDS || MYCOORDS==SPHCOORDS ||
 	  MYCOORDS==MKS1COORDS || MYCOORDS==MKS2COORDS || MYCOORDS==MKS3COORDS || MYCOORDS==TKS3COORDS ||
-	  MYCOORDS==MSPH1COORDS)
+	  MYCOORDS==MSPH1COORDS || MYCOORDS==JETCOORDS)
 	{
 	  #pragma omp parallel for private(ic,ix,iy,iz,iv,iysrc) schedule (static)
 	  for(ix=0;ix<NX;ix++)
@@ -5931,7 +5931,7 @@ correct_polaraxis_3d()
       //spherical like coords
       if (MYCOORDS==SCHWCOORDS || MYCOORDS==KSCOORDS || MYCOORDS==KERRCOORDS || MYCOORDS==SPHCOORDS ||
 	  MYCOORDS==MKS1COORDS || MYCOORDS==MKS2COORDS|| MYCOORDS==MKS3COORDS|| MYCOORDS==TKS3COORDS ||
-	  MYCOORDS==MSPH1COORDS)
+	  MYCOORDS==MSPH1COORDS || MYCOORDS==JETCOORDS)
 	{
           int ix;
 //#pragma omp parallel for private(ix)

metric.c

@@ -4478,7 +4478,7 @@ calc_gttpert_arb(double *xx, int COORDS)
       gpert[0][0]=gpert[1][1]=2.*r/Sigma;
       gpert[2][2]=gbase[2][2]-1.0;
       gpert[3][3]=gbase[3][3]-1.0;
-      //PIN
+
       //transformation matrix
       if(COORDS==MKS1COORDS) dxdx_MKS12KS(xx,dxdxp);
       if(COORDS==MKS2COORDS) dxdx_MKS22KS(xx,dxdxp);
@@ -4524,7 +4524,6 @@ calc_gttpert_arb(double *xx, int COORDS)
   return gttpert;
 }
 
-
 /******************************************************/
 //helper functions for jet coordinates and cylindrification
 /******************************************************/
@@ -4815,6 +4814,7 @@ int set_cyl_params()
 
 
    //test
+   /*
    printf("RCYL %e rmidcyl %e x2cyl %e \n",RCYL,rmidcyl,x2cyl);
    printf("THETACYL: %e THETAAX %e\n",thetaCYL,thetaAX);
    
@@ -4831,8 +4831,9 @@ int set_cyl_params()
    printf("thcyl %.14e\n",cylindrify(rtest,x2test,&tpar));
    if(sinth0(rtest,x2test,&tpar)>1)
      {printf("nan in cylindrify for r=%e, x2test=%e!\n",rtest,x2test); exit(-1);}
-
-  return 0;
+   */
+
+   return 0;
 }
 
 

misc.c

@@ -358,8 +358,7 @@ if (NTZ % 2 != 0)
 #endif
 
 #ifdef RADIATION  
-#if !defined(NO_COMPTONIZATION) && !defined(COMPTONIZATION)
-#define COMPTONIZATION
+#if defined(COMPTONIZATIONFLAG)
   if (PROCID == 0)
   {
     printf("\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");

physics.c

@@ -1654,7 +1654,7 @@ calc_TfromS3rho(ldouble S3,ldouble rho, int type,int ix,int iy,int iz)
     }
 
   if(!isfinite(T))
-    T=BIG;
+    T = BIG;
 
   return T;
 }
@@ -1775,29 +1775,37 @@ entri_from_entre_energy_cons(ldouble *pp, int ix, int iy, int iz)
 ldouble calc_PEQ_Teifrompp(ldouble* pp, ldouble* Te, ldouble* Ti,int ix, int iy, int iz)
 {
   ldouble Tgas=calc_PEQ_Tfromurho(pp[UU],pp[RHO],ix,iy,iz);
-
+
+  ldouble Tiloc,Teloc;
 #ifndef EVOLVEELECTRONS
-  *Ti=*Te=Tgas;
+  Tiloc=Teloc=Tgas;
 
 #ifdef FIXEDTETIRATIO
   ldouble factor=FIXEDTETIRATIO;
-  *Te=  (factor*MU_E*MU_I*Tgas)/(MU_GAS * (MU_E + factor * MU_I));
-  *Ti = *Te/factor;
+  Teloc=  (factor*MU_E*MU_I*Tgas)/(MU_GAS * (MU_E + factor * MU_I));
+  Tiloc = Teloc/factor;
 #endif
 
 #else //EVOLVEELECTRONS
 
   ldouble neth = calc_thermal_ne(pp);
   ldouble rhoeth = MU_E * M_PROTON * neth;
-  ldouble Teloc=calc_TfromSerho(pp[ENTRE],rhoeth,ELECTRONS,ix,iy,iz);
-  ldouble Tiloc=calc_TfromSerho(pp[ENTRI],pp[RHO],IONS,ix,iy,iz);
+  Teloc=calc_TfromSerho(pp[ENTRE],rhoeth,ELECTRONS,ix,iy,iz);
+  Tiloc=calc_TfromSerho(pp[ENTRI],pp[RHO],IONS,ix,iy,iz);
 
+#endif //EVOLVEELECTRONS
+
+  //!AC -- modified for  issues with gammagas in  corners
+  //!AC -- this will effectively make all corner gammagas 4/3
+  // I don't think that matters though...
+  if(!isfinite(Teloc)) Teloc=BIG;
+  if(!isfinite(Tiloc)) Tiloc=BIG;
+
   if(Teloc<TEMPEMINIMAL) Teloc=TEMPEMINIMAL;
   if(Tiloc<TEMPIMINIMAL) Tiloc=TEMPIMINIMAL;
-  
+
   *Te=Teloc;
   *Ti=Tiloc;
-#endif //EVOLVEELECTRONS
 
   return Tgas;
 }
@@ -1953,7 +1961,6 @@ ldouble pick_gammagas(int ix,int iy,int iz)
 //********************************************************************************
 // Calculate gamma_int for total gas
 //********************************************************************************
-
 ldouble
 calc_gammagas(ldouble* pp,int ix,int iy,int iz)
 {
@@ -1970,9 +1977,26 @@ calc_gammagas(ldouble* pp,int ix,int iy,int iz)
   gamma=calc_gammaintfromTei(Te,Ti);
   #endif
 
+  //!AC
+  /*
+  ldouble Tgas=calc_PEQ_Tfromurho(pp[UU],pp[RHO],ix,iy,iz);
+  if(!isfinite(Te) || !isfinite(Ti) || !isfinite(Tgas))
+    {
+      int gix, giy, giz;
+      mpi_local2globalidx(ix, iy, iz, &gix, &giy, &giz);
+      printf("Te/Ti not finite >>> %d: > %d %d %d > %e %e %e %e %e %e\n",PROCID,gix,giy,giz,Te,Ti,Tgas,pp[RHO],pp[UU],gamma); 
+      fflush(stdout);
+      //getch();
+      //exit(-1);
+    }
+  */
+
   if(!isfinite(gamma))
     {
-      printf("gammagas not finite >>> %d: > %d %d %d > %e %e %e %e %e\n",PROCID,ix,iy,iz,Te,Ti,calc_relel_ne(pp), calc_relel_uint(pp),gamma); 
+      ldouble Tgas=calc_PEQ_Tfromurho(pp[UU],pp[RHO],ix,iy,iz);
+      int gix, giy, giz;
+      mpi_local2globalidx(ix, iy, iz, &gix, &giy, &giz);
+      printf("gammagas not finite >>> %d: > %d %d %d > %e %e %e %e %e %e\n",PROCID,gix,giy,giz,Te,Ti,Tgas,pp[RHO],pp[UU],gamma); 
       fflush(stdout);
       //getch();
       //exit(-1);
@@ -2347,7 +2371,7 @@ heat_electronions_with_state(ldouble dtin)
 	  p_index=RELEL_HEAT_INDEX;
 	  #endif
 
-          #endif
+          #endif //RELEL_HEAT_RECONNECTION
 
 	  //Determine gamma_inj_min & gamma_inj_max
 	  #ifndef RELEL_HEAT_FIX_LIMITS

postproc.c

@@ -521,9 +521,10 @@ int calc_radialprofiles(ldouble profiles[][NX])
         //if(muBe>0.05 && (xxBL[2]<M_PI/4. || xxBL[2]>3.*M_PI/4.))
         //  isjet=1;
 	//new criterion based on beta*gamma
-	if(betagamma2>1.)
+	if(betagamma2>1. && xxBL[2]<M_PI/3.) // top jet only
 	  isjet=1;
-        else isjet=0;
+        else 
+          isjet=0;
 
         ldouble pregas = GAMMAM1*uint;
         ldouble ptot = pregas;
@@ -631,7 +632,7 @@ int calc_radialprofiles(ldouble profiles[][NX])
 
 	//density-weighted MRI wavelength using code comparison paper defn
 	//ANDREW?? -- should it be bsq/2?
-	ldouble lambdaMRI = 2*M_PI * bcon[2] / sqrt(bsq + rho + uint + pregas) / Omega;
+	ldouble lambdaMRI = 2*M_PI * fabs(bcon[2]) / sqrt(bsq + rho + uint + pregas) / fabs(Omega);
 	rholambda134[ix] += rho*lambdaMRI* dx[1]*dx[2]*geomBL.gdet; 
 
         //total mhd energy flux (14)
@@ -914,7 +915,7 @@ int calc_radialprofiles(ldouble profiles[][NX])
     profiles[11][ix]=calc_photloc(ix);
 
     // special quantities for problem 134/139
-    if (PROBLEM == 134 || PROBLEM ==  139)
+    if (PROBLEM == 134 || PROBLEM == 139)
     {
       profiles[7][ix] = rho134[ix] / normalize[ix];
       profiles[8][ix] = pgas134[ix] / normalize[ix];

problem.h

@@ -146,7 +146,7 @@
 //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 87
+#define PROBLEM 139
 
 
 #if(PROBLEM==139)

rad.c

@@ -737,8 +737,9 @@ solve_implicit_lab_4dprim(ldouble *uu00,ldouble *pp00,void *ggg,
 	  print_state_implicit_lab_4dprim(iter-1,xxx,f1,err,N); 
 	  printf("\n kappaCGS: %e ", kappaGU2CGS(calc_kappa(pp,geom,&opac)/pp[RHO]) );
 	  printf("rhoCGS: %e ", rhoGU2CGS(pp[RHO]));
-	  ldouble Ti0,Te0;
           ldouble Tgas0=state0.Tgas;
+	  ldouble Te0=state0.Te;
+	  ldouble Ti0=state0.Ti;
 	  ldouble Trad0=state0.Trad;
 	  printf("T: %e Te: %e ", Tgas0,Te0);
 	  printf("Trad: %e ", Trad0);
@@ -1644,7 +1645,12 @@ f_implicit_lab_4dprim_with_state(ldouble *uu, ldouble *pp, void *sss,
 
     //ANDREW - in the entropy equation, should there be a chemical potential term due to change in volume?
     //is this equation without relel consistent with equation in terms of entropy per particle? 
-    ldouble Qe = CC + EA + CC_relel + cool_relel_dq - mu*cool_relel_dn;
+    ldouble Qe;
+    Qe = CC + EA + CC_relel + cool_relel_dq - mu*cool_relel_dn;
+    #ifdef SKIP_ELECTRON_COOL
+    Qe=0.;
+    #endif
+
     f[ientre]=Te*(pp[ENTRE] - pp0[ENTRE]) - dtau * Qe;
     if(fabs(f[ientre])>SMALL)
       err[ientre]=fabs(f[ientre])/(fabs(pp[ENTRE]*Te)+fabs(pp0[ENTRE]*Te)+fabs(dtau*Qe));
@@ -2755,7 +2761,7 @@ calc_all_Gi_with_state(ldouble *pp, void *sss, void* ggg,
   fac=step_function(xxBL[1]-1.2*rhorizonBL,.1*rhorizonBL);
   if(xxBL[1]<rhorizonBL) fac=0.;
 #endif
-  
+
   for(i = 0; i < 4; i++)
   {
     Gith_lab[i] += fac * Gic_lab[i];
@@ -2931,14 +2937,16 @@ calc_Compt_Gi_with_state(ldouble *pp, void *sss, void* ggg, ldouble *Gic, ldoubl
   ldouble kappaes = state->kappaes;
   ldouble ThatradBB = state->TradBB;
   ldouble Thatrad = state->Trad;
-
+
+  /*
   ldouble urfcon[4], uffcov[4];
   for (i = 0; i < 4; i++)
   {
     urfcon[i] = state->urfcon[i];
     uffcov[i] = state->ucov[i];
   }
-
+  */
+
   //ANDREW correction factor for the nonthermal electrons present in kappa_es
   ldouble relel_corr = 1.0;
 #ifdef RELELECTRONS

relele.c

@@ -269,6 +269,11 @@ conv_vels_core(ldouble *u1,ldouble *u2conout,int which1,int which2,ldouble gg[][
     if(u2con[0]<0)
       u2con[0] = fabs(u2con[0]);
 
+    //ANDREW's version
+    //u2con[1]=(u1[1]-alpgam*GG[0][1])/u2con[0];
+    //u2con[2]=(u1[2]-alpgam*GG[0][2])/u2con[0];
+    //u2con[3]=(u1[3]-alpgam*GG[0][3])/u2con[0];
+
     //Identical to Olek's
     u2con[1]=u1[1]/u2con[0] + GG[0][1]/GG[0][0];
     u2con[2]=u1[2]/u2con[0] + GG[0][2]/GG[0][0];
@@ -511,6 +516,43 @@ calc_normalobs_relvel(ldouble GG[][5], ldouble *ncon)
 int
 set_hdatmosphere(ldouble *pp,ldouble *xx,ldouble gg[][5],ldouble GG[][5],int atmtype)
 {
+  if(atmtype==-1) // for MAD code comparison project
+    {
+      // normal observer -- zero in BL
+      ldouble ucon[4];
+      ldouble xx2[4];
+      ldouble GGBL[4][5];
+
+      // BL coords
+      coco_N(xx,xx2,MYCOORDS,BLCOORDS);
+      calc_G_arb(xx2,GGBL,BLCOORDS);
+
+      // normal observer in BL = stationary observer
+      calc_normalobs_4vel(GGBL,ucon);
+
+      // to MYCOORDS
+      trans2_coco(xx2,ucon,ucon,BLCOORDS,MYCOORDS);
+
+      // to VELPRIM
+      conv_vels(ucon,ucon,VEL4,VELPRIM,gg,GG);
+
+      pp[2]=ucon[1];
+      pp[3]=ucon[2];
+      pp[4]=ucon[3];
+
+      // density & pressure
+      ldouble r=xx2[1];
+      ldouble rout=1.; //RHOATMMIN etc. given at rout=2
+
+      pp[0] = RHOATMMIN*pow(r/rout,-1.5);
+      pp[1] = UINTATMMIN*pow(r/rout,-2.5);
+
+      #ifdef MAGNFIELD
+      pp[B1]=pp[B2]=pp[B3]=0.;
+      #endif
+
+      return 0;
+    }
   if(atmtype==0) //normal observer plur rho \propto r^-1.5
     {
       //normal observer

u2p.c

@@ -1848,7 +1848,7 @@ u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose)
     return -102;
   }
 
-  if(!isfinite(W) || !isfinite(W))
+  if(!isfinite(W))
   {
     if(verbose) printf("nan/inf W in u2p_solver with Etype: %d\n",Etype);
     return -103;