From e341ca54880c070eaa84d0025b222f411f1ca2ea Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Mon, 1 Mar 2021 12:51:59 -0700
Subject: [PATCH 01/14] Add code for a radiation pressure dominated disk

---
 core/decs.h              |  7 +++++--
 core/defs.h              |  2 +-
 core/eos.c               |  6 ++++--
 core/eos_gamma.c         |  6 +++++-
 core/fixup.c             |  4 ++--
 core/io.c                |  6 +++---
 core/main.c              |  6 ++++++
 core/phys.c              |  2 +-
 core/utop.c              |  2 +-
 prob/torus_cbc/build.py  | 36 +++++++++++++++++++++++++++++-------
 prob/torus_cbc/problem.c | 16 +++++++++++-----
 11 files changed, 68 insertions(+), 25 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index e55e35d..e471ec5 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -196,9 +196,11 @@
 
 // EOS
 #define EOS_TYPE_GAMMA (0)
+#define EOS_TYPE_GAMMA_GASPRESS (0)
+#define EOS_TYPE_GAMMA_RADPRESS (0)
 #define EOS_TYPE_POLYTROPE (1)
 #define EOS_TYPE_TABLE (2)
-#if EOS == EOS_TYPE_GAMMA
+#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS
 #define EOS_NUM_EXTRA (0)
 #define POLYTROPE_FALLBACK (0)
 #elif EPS == EOS_TYPE_POLYTROPE
@@ -787,7 +789,8 @@ double EOS_Theta_unit();
 #if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
 double EOS_Gamma_pressure_rho0_u(double rho, double u);
 double EOS_Gamma_pressure_rho0_w(double rho, double w);
-double EOS_Gamma_entropy_rho0_u(double rho, double u);
+double EOS_Gamma_Gaspress_entropy_rho0_u(double rho, double u);
+double EOS_Gamma_Radpress_entropy_rho0_u(double rho, double u);
 double EOS_Gamma_enthalpy_rho0_u(double rho, double u);
 double EOS_Gamma_adiabatic_constant_rho0_u(double rho, double u);
 double EOS_Gamma_sound_speed_rho0_u(double rho, double u);
diff --git a/core/defs.h b/core/defs.h
index a26d977..f9fc69c 100644
--- a/core/defs.h
+++ b/core/defs.h
@@ -98,7 +98,7 @@ char opac_file[STRLEN];
 #endif // neutrinos
 #endif // radiation
 
-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK
 double gam;
 #endif
 #if EOS == EOS_TYPE_POLYTROPE
diff --git a/core/eos.c b/core/eos.c
index 2279072..358bb0e 100644
--- a/core/eos.c
+++ b/core/eos.c
@@ -159,8 +159,10 @@ double EOS_enthalpy_rho0_u(double rho, double u, const double *extra) {
 
 double EOS_entropy_rho0_u(double rho, double u, const double *extra) {
   double ent;
-#if EOS == EOS_TYPE_GAMMA
-  ent = EOS_Gamma_entropy_rho0_u(rho, u);
+#if EOS == EOS_TYPE_GAMMA_GASPRESS
+  ent = EOS_Gamma_Gaspress_entropy_rho0_u(rho, u);
+#elif EOS == EOS_TYPE_GAMMA_RADPRESS
+  ent = EOS_Gamma_Radpress_entropy_rho0_u(rho, u);
 #elif EOS == EOS_TYPE_POLYTROPE
   ent   = EOS_Poly_entropy_rho0_u(rho, u, poly_K, poly_gam);
 #elif EOS == EOS_TYPE_TABLE
diff --git a/core/eos_gamma.c b/core/eos_gamma.c
index 225748e..bf2abb0 100644
--- a/core/eos_gamma.c
+++ b/core/eos_gamma.c
@@ -17,10 +17,14 @@ double EOS_Gamma_pressure_rho0_w(double rho, double w) {
   return ((w - rho) * (gam - 1.) / gam);
 }
 
-double EOS_Gamma_entropy_rho0_u(double rho, double u) {
+double EOS_Gamma_Gaspress_entropy_rho0_u(double rho, double u) {  
   return (gam - 1.0) * EOS_Gamma_adiabatic_constant_rho0_u(rho, u);
 }
 
+double EOS_Gamma_Radpress_entropy_rho0_u(double rho, double u) {
+  return pow((64./3) * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
+}
+
 double EOS_Gamma_enthalpy_rho0_u(double rho, double u) { return rho + u * gam; }
 
 double EOS_Gamma_adiabatic_constant_rho0_u(double rho, double u) {
diff --git a/core/fixup.c b/core/fixup.c
index 86ad1bd..d027ed4 100644
--- a/core/fixup.c
+++ b/core/fixup.c
@@ -229,9 +229,9 @@ void fixup1zone(
   fixup_passive(i, j, k, pv, pv_prefloor);
 #endif
 
-#if ELECTRONS && EOS == EOS_TYPE_GAMMA
+#if ELECTRONS && EOS == EOS_TYPE_GAMMA_GASPRESS  
   // Reset entropy after floors
-  pv[KTOT] = EOS_Gamma_entropy_rho0_u(pv[RHO], pv[UU]);
+  pv[KTOT] = EOS_Gamma_Gaspress_entropy_rho0_u(pv[RHO], pv[UU]);
 
   // Set KTOTMAX to 3 by controlling u, to avoid anomalous cooling from funnel
   // wall
diff --git a/core/io.c b/core/io.c
index 58143f7..ed3a82a 100644
--- a/core/io.c
+++ b/core/io.c
@@ -730,7 +730,7 @@ void dump() {
 #endif // RADIATION
 #endif // METRIC
 
-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  
   WRITE_HDR(gam, TYPE_DBL);
 #endif // GAMMA
 #if EOS == EOS_TYPE_POLYTROPE
@@ -1099,7 +1099,7 @@ void restart_write(int restart_type) {
 #endif // METRIC
 
 // EOS
-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  
   WRITE_HDR(gam, TYPE_DBL);
 #endif // GAMMA EOS
 #if EOS == EOS_TYPE_POLYTROPE
@@ -1381,7 +1381,7 @@ void restart_read(char *fname) {
 
   READ_HDR(cour, TYPE_DBL);
 
-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  
   READ_HDR(gam, TYPE_DBL);
 #endif // EOS
 
diff --git a/core/main.c b/core/main.c
index fd72f64..fedc7c1 100644
--- a/core/main.c
+++ b/core/main.c
@@ -60,6 +60,12 @@ int main(int argc, char *argv[]) {
     #if EOS == EOS_TYPE_GAMMA
     fprintf(stdout, "          *  IDEAL GAS EOS                                           *\n");
     #endif
+    #if EOS == EOS_TYPE_GAMMA_GASPRESS
+    fprintf(stdout, "          *  GAS-PRESSURE DOMINATED EOS                    *\n");
+    #endif
+    #if EOS == EOS_TYPE_GAMMA_RADPRESS
+    fprintf(stdout, "          *  RADIATION-PRESSURE DOMINATED EOS              *\n");
+    #endif
     #if RADIATION == RADTYPE_LIGHT
     fprintf(stdout, "          *  RADIATION: PHOTON TRANSPORT                             *\n");
     #elif RADIATION == RADTYPE_NEUTRINOS
diff --git a/core/phys.c b/core/phys.c
index cd2ce98..fbbcc4d 100644
--- a/core/phys.c
+++ b/core/phys.c
@@ -211,7 +211,7 @@ void mhd_vchar(double *Pr, struct of_state *q, struct of_geom *geom, int js,
   bsq = dot(q->bcon, q->bcov);
   rho = fabs(Pr[RHO] + SMALL);
   u   = Pr[UU];
-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  
   u = fabs(u + SMALL);
 #endif
   ef  = EOS_enthalpy_rho0_u(rho, u, extra);
diff --git a/core/utop.c b/core/utop.c
index 95e25cc..3c946b9 100644
--- a/core/utop.c
+++ b/core/utop.c
@@ -161,7 +161,7 @@ int Utoprim(double U[NVAR], struct of_geom *geom, double prim[NVAR]) {
   // Return without updating non-B primitives
   if (rho0 < 0)
     return (5);
-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  
   if (u < 0)
     return (5);
 #endif
diff --git a/prob/torus_cbc/build.py b/prob/torus_cbc/build.py
index 84ad77e..b3943ef 100644
--- a/prob/torus_cbc/build.py
+++ b/prob/torus_cbc/build.py
@@ -18,9 +18,10 @@
 from math import log10,ceil
 PROB = 'torus_cbc'
 
-DO_GAMMA = '-gamma' in sys.argv # No table
+DO_GAMMA_RADPRESS = '-gammaradpress' in sys.argv # No table
+DO_GAMMA_GASPRESS = '-gammagaspress' in sys.argv # No table
 GAMTABLE = '-gamtable' in sys.argv # fake table
-RELTABLE = '-reltable' in sys.argv or (not DO_GAMMA or GAMTABLE)
+RELTABLE = '-reltable' in sys.argv or (not DO_GAMMA_RADPRESS or DO_GAMMA_GASPRESS or GAMTABLE)
 INITIAL = '-initial' in sys.argv # initial data only
 NOB = '-nob' in sys.argv # or INITIAL # no B fields
 TOROIDALB = "-toroidalb" in sys.argv
@@ -51,6 +52,7 @@
 KILL = '-kill' in sys.argv # kill all packets
 
 ENT_FROM_CLI  = '-ent' in sys.argv
+KAPPA_FROM_CLI = '-kappa' in sys.argv
 RHO_FROM_CLI  = '-rho' in sys.argv
 RIN_FROM_CLI  = '-rin' in sys.argv
 RMAX_FROM_CLI = '-rmax' in sys.argv
@@ -146,12 +148,31 @@
 if USE_TABLE:
     bhl.report_var('Disk Ye',YE)
 
+if USE_TABLE:
+    EOS_TYPE = "EOS_TYPE_TABLE"
+elif DO_GAMMA_RADPRESS:
+    EOS_TYPE = "EOS_TYPE_GAMMA_RADPRESS"
+elif DO_GAMMA_GASPRESS:
+    EOS_TYPE = "EOS_TYPE_GAMMA_GASPRESS"
+
 if ENT_FROM_CLI:
     ENTROPY = float(sys.argv[sys.argv.index('-ent') + 1])
 else:
     ENTROPY = 4
+        if EOS_TYPE == "EOS_TYPE_GAMMA_RADPRESS":
+        print("Using the default initial entropy = 4.", end=' ')
+        print("Entropy is a required parameter in setting up radiation-pressure dominated disks.")
+
 bhl.report_var('ENTROPY',ENTROPY)
 
+if KAPPA_FROM_CLI:
+    KAPPA_EOS = float(sys.argv[sys.argv.index('-kappa') + 1])
+else:
+    KAPPA_EOS = 1.e-3
+    if EOS_TYPE == "EOS_TYPE_GAMMA_GASPRESS":
+        print("Using the default initial kappa = 1.e-3.", end=' ')
+        print("kappa is a required parameter in setting up gas-pressure dominated disks.")
+
 if CLASSIC: # classic harm disk
     # Rmax and rho fixed
     bhl.report_var('DISK_TYPE','CLASSIC')
@@ -211,8 +232,7 @@
 # use selection instead
 Rout_rad = ENTROPY*ceil(Rmax) # not safe to use 3x
 Rout_vis = Rout
-GAMMA = 13./9.
-KAPPA = 1.e-3
+GAMMA = 4./3.
 L_unit = cgs['GNEWT']*cgs['MSOLAR']*MBH/(cgs['CL']**2)
 M_UNIT = RHO_unit*(L_unit**3)
 
@@ -314,7 +334,6 @@
 if not NEUTRINOS:
     NPH_PER_PROC = 0.0
 
-EOS = "EOS_TYPE_TABLE" if USE_TABLE else "EOS_TYPE_GAMMA"
 NVP0 = 3
 NVAR_PASSIVE = NVP0 if USE_TABLE else 0
 #GAMMA_FALLBACK = True
@@ -389,7 +408,7 @@
 bhl.config.set_cparm('QUADRANT_SYMMETRY', QUAD)
 
 # EOS
-bhl.config.set_cparm("EOS", EOS)
+bhl.config.set_cparm("EOS_TYPE", EOS_TYPE)
 bhl.config.set_cparm('NVAR_PASSIVE', NVAR_PASSIVE)
 bhl.config.set_cparm('GAMMA_FALLBACK', GAMMA_FALLBACK)
 
@@ -466,7 +485,10 @@
         bhl.config.set_rparm('lrho_guess','double',LRHO_GUESS)
 if USE_GAMMA or GAMMA_FALLBACK:
     bhl.config.set_rparm('gam', 'double', default = GAMMA)
-    bhl.config.set_rparm("kappa", "double", default = KAPPA)
+        if EOS_TYPE == "EOS_TYPE_GAMMA_RADPRESS":
+        bhl.config.set_rparm('entropy','double',ENTROPY)
+    elif EOS_TYPE == "EOS_TYPE_GAMMA_GASPRESS":
+        bhl.config.set_rparm("kappa_eos", "double", default = KAPPA_EOS)
 
 # Opacities
 if FORTRAN or HDF:
diff --git a/prob/torus_cbc/problem.c b/prob/torus_cbc/problem.c
index ce50867..f2d1126 100644
--- a/prob/torus_cbc/problem.c
+++ b/prob/torus_cbc/problem.c
@@ -16,8 +16,8 @@ static char   bfield_type[STRLEN];
 static int    renormalize_densities;
 static double rin;
 static double rmax;
-static double beta;
-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+static double beta; // desired minimum ratio of gas to magnetic pressure
+#if EOS == EOS_TYPE_GAMMA_GASPRESS || GAMMA_FALLBACK
 static double kappa_eos;
 #endif
 #if EOS == EOS_TYPE_TABLE
@@ -44,7 +44,7 @@ void set_problem_params() {
   set_param("rmax", &rmax);
   set_param("beta", &beta);
   set_param("renorm_dens", &renormalize_densities);
-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA_GASPRESS || GAMMA_FALLBACK
   set_param("kappa_eos", &kappa_eos);
 #endif
 #if EOS == EOS_TYPE_TABLE
@@ -221,7 +221,11 @@ void init_prob() {
       disk_cell[i][j][k] = 1;
 
       hm1 = exp(lnh[i][j][k]) - 1.;
-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA_RADPRESS
+      fprintf(stdout, "entropy: %f\n", entropy);
+      rho = (64./3) * (pow(hm1, 3)/pow(entropy, 4)) * (pow((gam - 1.)/gam), 3);
+      u = hm1*rho/gam;
+#elif EOS == EOS_TYPE_GAMMA_GASPRESS || GAMMA_FALLBACK
       rho = pow(hm1 * (gam - 1.) / (kappa_eos * gam), 1. / (gam - 1.));
       u   = kappa_eos * pow(rho, gam) / (gam - 1.);
 #elif EOS == EOS_TYPE_TABLE
@@ -316,6 +320,8 @@ void init_prob() {
       PsaveLocal[i][j][k][UU] /= rhomax;
     }
 
+    // Only let cells inside the torus (radius greater than rin and non-minimal 
+    //     // enthalpy) contribute to pressmax
     if (r > rin && lnh[i][j][k] >= 0.) {
 #if EOS == EOS_TYPE_TABLE
       EOS_SC_fill(PsaveLocal[i][j][k], extra[i][j][k]);
@@ -376,7 +382,7 @@ void init_prob() {
         mtot, mtot * M_unit, mtot * M_unit / MSUN);
   }
 // debug
-#if EOS == EOS_TYPE_TABLE
+#if EOS == EOS_TYPE_TABLE || EOS_TYPE_GAMMA_RADPRESS
   if (mpi_io_proc()) {
     fprintf(stdout, "Calculating max entropy:\n");
   }

From 2b612f214e536d66a53f748c5f3a413861ff62ab Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Mon, 1 Mar 2021 12:58:20 -0700
Subject: [PATCH 02/14] Fix a commented line

---
 prob/torus_cbc/problem.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/prob/torus_cbc/problem.c b/prob/torus_cbc/problem.c
index f2d1126..b3002f5 100644
--- a/prob/torus_cbc/problem.c
+++ b/prob/torus_cbc/problem.c
@@ -321,7 +321,7 @@ void init_prob() {
     }
 
     // Only let cells inside the torus (radius greater than rin and non-minimal 
-    //     // enthalpy) contribute to pressmax
+    // enthalpy) contribute to pressmax
     if (r > rin && lnh[i][j][k] >= 0.) {
 #if EOS == EOS_TYPE_TABLE
       EOS_SC_fill(PsaveLocal[i][j][k], extra[i][j][k]);

From be323b413c29aaf196feab8352103efd46fe04ae Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Mon, 1 Mar 2021 20:33:21 -0700
Subject: [PATCH 03/14] Keep EOS_TYPE_GAMMA and break that into GAMMA_EOS =
 RADPRESS and GASPRESS

---
 core/decs.h              |  8 +++----
 core/defs.h              |  2 +-
 prob/torus_cbc/build.py  | 45 +++++++++++++++++++++++++++-------------
 prob/torus_cbc/problem.c | 18 ++++++++--------
 4 files changed, 45 insertions(+), 28 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index e471ec5..7a2e9b5 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -196,14 +196,14 @@
 
 // EOS
 #define EOS_TYPE_GAMMA (0)
-#define EOS_TYPE_GAMMA_GASPRESS (0)
-#define EOS_TYPE_GAMMA_RADPRESS (0)
 #define EOS_TYPE_POLYTROPE (1)
 #define EOS_TYPE_TABLE (2)
-#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS
+#if EOS == EOS_TYPE_GAMMA
 #define EOS_NUM_EXTRA (0)
 #define POLYTROPE_FALLBACK (0)
-#elif EPS == EOS_TYPE_POLYTROPE
+#define GASPRESS (1)
+#define RADPRESS (2)
+#elif EOS == EOS_TYPE_POLYTROPE
 #define EOS_NUM_EXTRA (0)
 #define POLYTROPE_FALLBACK (0)
 #elif EOS == EOS_TYPE_TABLE
diff --git a/core/defs.h b/core/defs.h
index f9fc69c..a26d977 100644
--- a/core/defs.h
+++ b/core/defs.h
@@ -98,7 +98,7 @@ char opac_file[STRLEN];
 #endif // neutrinos
 #endif // radiation
 
-#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
 double gam;
 #endif
 #if EOS == EOS_TYPE_POLYTROPE
diff --git a/prob/torus_cbc/build.py b/prob/torus_cbc/build.py
index b3943ef..1075061 100644
--- a/prob/torus_cbc/build.py
+++ b/prob/torus_cbc/build.py
@@ -18,10 +18,11 @@
 from math import log10,ceil
 PROB = 'torus_cbc'
 
-DO_GAMMA_RADPRESS = '-gammaradpress' in sys.argv # No table
-DO_GAMMA_GASPRESS = '-gammagaspress' in sys.argv # No table
+DO_GAMMA = '-gamma' in sys.argv # No table
+GAMMA_GASPRESS = '-gammagaspress' in sys.argv # No table
+GAMMA_RADPRESS = '-gammaradpress' in sys.argv # No table
 GAMTABLE = '-gamtable' in sys.argv # fake table
-RELTABLE = '-reltable' in sys.argv or (not DO_GAMMA_RADPRESS or DO_GAMMA_GASPRESS or GAMTABLE)
+RELTABLE = '-reltable' in sys.argv or (not DO_GAMMA or GAMTABLE)
 INITIAL = '-initial' in sys.argv # initial data only
 NOB = '-nob' in sys.argv # or INITIAL # no B fields
 TOROIDALB = "-toroidalb" in sys.argv
@@ -85,6 +86,9 @@
 if any([N1TOT_FROM_CLI, N2TOT_FROM_CLI, N3TOT_FROM_CLI]) and not N1N2N3TOT_FROM_CLI:
     raise ValueError("Must turn on the -n1n2n3tot flag if inputting -n1tot, -n2tot, -n3tot")
 
+if any([GAMMA_GASPRESS, GAMMA_RADPRESS]) and not DO_GAMMA:
+    raise ValueError("Must turn on the -gamma flag if inputting -gammagaspress or -gammaradpress")
+
 if NOB:
     BFIELD = "none"
 elif TOROIDALB:
@@ -150,16 +154,28 @@
 
 if USE_TABLE:
     EOS_TYPE = "EOS_TYPE_TABLE"
-elif DO_GAMMA_RADPRESS:
-    EOS_TYPE = "EOS_TYPE_GAMMA_RADPRESS"
-elif DO_GAMMA_GASPRESS:
-    EOS_TYPE = "EOS_TYPE_GAMMA_GASPRESS"
-
+if DO_GAMMA:
+    EOS_TYPE = "EOS_TYPE_GAMMA"
+    if GAMMA_RADPRESS:
+        EOS_GAMMA = "GAMMA_RADPRESS"
+        GAMMA = 4./3.
+    elif GAMMA_GASPRESS:
+        EOS_GAMMA = "GAMMA_GASPRESS"
+        GAMMA = 5./3.
+    else:
+        print("Using the default gas-pressure dominated ideal gas EoS. 
+              You can choose gas-pressure dominated (add flag -gammagaspress),
+              radiation-pressure dominated (add flag -gammaradpress).")
+        EOS_GAMMA = "GAMMA_GASPRESS"
+        GAMMA = 4./3.
+else:
+    raise ValueError("Bad EOS chosen")
+        
 if ENT_FROM_CLI:
     ENTROPY = float(sys.argv[sys.argv.index('-ent') + 1])
 else:
     ENTROPY = 4
-        if EOS_TYPE == "EOS_TYPE_GAMMA_RADPRESS":
+        if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "GAMMA_RADPRESS":
         print("Using the default initial entropy = 4.", end=' ')
         print("Entropy is a required parameter in setting up radiation-pressure dominated disks.")
 
@@ -169,7 +185,7 @@
     KAPPA_EOS = float(sys.argv[sys.argv.index('-kappa') + 1])
 else:
     KAPPA_EOS = 1.e-3
-    if EOS_TYPE == "EOS_TYPE_GAMMA_GASPRESS":
+    if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "GAMMA_GASPRESS":
         print("Using the default initial kappa = 1.e-3.", end=' ')
         print("kappa is a required parameter in setting up gas-pressure dominated disks.")
 
@@ -232,7 +248,6 @@
 # use selection instead
 Rout_rad = ENTROPY*ceil(Rmax) # not safe to use 3x
 Rout_vis = Rout
-GAMMA = 4./3.
 L_unit = cgs['GNEWT']*cgs['MSOLAR']*MBH/(cgs['CL']**2)
 M_UNIT = RHO_unit*(L_unit**3)
 
@@ -408,7 +423,9 @@
 bhl.config.set_cparm('QUADRANT_SYMMETRY', QUAD)
 
 # EOS
-bhl.config.set_cparm("EOS_TYPE", EOS_TYPE)
+bhl.config.set_cparm("EOS", EOS_TYPE)
+if EOS_TYPE == 'EOS_TYPE_GAMMA':
+    bhl.config.set_cparm("EOS_GAMMA", EOS_GAMMA)
 bhl.config.set_cparm('NVAR_PASSIVE', NVAR_PASSIVE)
 bhl.config.set_cparm('GAMMA_FALLBACK', GAMMA_FALLBACK)
 
@@ -485,9 +502,9 @@
         bhl.config.set_rparm('lrho_guess','double',LRHO_GUESS)
 if USE_GAMMA or GAMMA_FALLBACK:
     bhl.config.set_rparm('gam', 'double', default = GAMMA)
-        if EOS_TYPE == "EOS_TYPE_GAMMA_RADPRESS":
+    if EOS_TYPE == "EOS_TYPE_GAMMA" and GAMMA_EOS == "RADPRESS":
         bhl.config.set_rparm('entropy','double',ENTROPY)
-    elif EOS_TYPE == "EOS_TYPE_GAMMA_GASPRESS":
+    elif EOS_TYPE == "EOS_TYPE_GAMMA" and GAMMA_EOS == "GASPRESS":
         bhl.config.set_rparm("kappa_eos", "double", default = KAPPA_EOS)
 
 # Opacities
diff --git a/prob/torus_cbc/problem.c b/prob/torus_cbc/problem.c
index b3002f5..26848ab 100644
--- a/prob/torus_cbc/problem.c
+++ b/prob/torus_cbc/problem.c
@@ -17,7 +17,7 @@ static int    renormalize_densities;
 static double rin;
 static double rmax;
 static double beta; // desired minimum ratio of gas to magnetic pressure
-#if EOS == EOS_TYPE_GAMMA_GASPRESS || GAMMA_FALLBACK
+#if (EOS == EOS_TYPE_GAMMA && GAMMA_EOS == GASPRESS) || GAMMA_FALLBACK
 static double kappa_eos;
 #endif
 #if EOS == EOS_TYPE_TABLE
@@ -44,10 +44,10 @@ void set_problem_params() {
   set_param("rmax", &rmax);
   set_param("beta", &beta);
   set_param("renorm_dens", &renormalize_densities);
-#if EOS == EOS_TYPE_GAMMA_GASPRESS || GAMMA_FALLBACK
+#if (EOS_TYPE == EOS_TYPE_GAMMA && GAMMA_EOS == GASPRESS) || GAMMA_FALLBACK
   set_param("kappa_eos", &kappa_eos);
 #endif
-#if EOS == EOS_TYPE_TABLE
+#if EOS_TYPE == EOS_TYPE_TABLE
   set_param("const_ye", &const_ye);
 #if !GAMMA_FALLBACK
   set_param("entropy", &entropy);
@@ -68,7 +68,7 @@ void init_prob() {
   double SSin, hm1;
 
 // Diagnostics for entropy
-#if EOS == EOS_TYPE_TABLE
+#if EOS_TYPE == EOS_TYPE_TABLE
   double ent, entmax;
 #endif
 
@@ -82,7 +82,7 @@ void init_prob() {
   double rhomin, umin;
   rhomin = RHOMINLIMIT;
   umin   = UUMINLIMIT;
-#if EOS == EOS_TYPE_TABLE
+#if EOS_TYPE == EOS_TYPE_TABLE
   umin = EOS_SC_get_minu(rhomin, const_ye, umin);
 #endif
 
@@ -94,7 +94,7 @@ void init_prob() {
   PASSTYPE(PASSIVE_START + 1) = PASSTYPE_NUMBER;
 #endif
 
-#if EOS == EOS_TYPE_TABLE
+#if EOS_TYPE == EOS_TYPE_TABLE
   double ye, ye_atm;
 #if !GAMMA_FALLBACK
   // guesses. Should be function local.
@@ -221,11 +221,11 @@ void init_prob() {
       disk_cell[i][j][k] = 1;
 
       hm1 = exp(lnh[i][j][k]) - 1.;
-#if EOS == EOS_TYPE_GAMMA_RADPRESS
+#if EOS == EOS_TYPE_GAMMA && GAMMA_EOS == RADPRESS
       fprintf(stdout, "entropy: %f\n", entropy);
       rho = (64./3) * (pow(hm1, 3)/pow(entropy, 4)) * (pow((gam - 1.)/gam), 3);
       u = hm1*rho/gam;
-#elif EOS == EOS_TYPE_GAMMA_GASPRESS || GAMMA_FALLBACK
+#elif (EOS == EOS_TYPE_GAMMA && GAMMA_EOS == GASPRESS) || GAMMA_FALLBACK
       rho = pow(hm1 * (gam - 1.) / (kappa_eos * gam), 1. / (gam - 1.));
       u   = kappa_eos * pow(rho, gam) / (gam - 1.);
 #elif EOS == EOS_TYPE_TABLE
@@ -382,7 +382,7 @@ void init_prob() {
         mtot, mtot * M_unit, mtot * M_unit / MSUN);
   }
 // debug
-#if EOS == EOS_TYPE_TABLE || EOS_TYPE_GAMMA_RADPRESS
+#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && GAMMA_EOS == RADPRESS)
   if (mpi_io_proc()) {
     fprintf(stdout, "Calculating max entropy:\n");
   }

From 79aec3b0a095585a15aee81a3dea9367829b4d37 Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Mon, 1 Mar 2021 20:38:37 -0700
Subject: [PATCH 04/14] Fix the gamma option names consistent to be consistent

---
 prob/torus_cbc/build.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/prob/torus_cbc/build.py b/prob/torus_cbc/build.py
index 1075061..05d6644 100644
--- a/prob/torus_cbc/build.py
+++ b/prob/torus_cbc/build.py
@@ -157,16 +157,16 @@
 if DO_GAMMA:
     EOS_TYPE = "EOS_TYPE_GAMMA"
     if GAMMA_RADPRESS:
-        EOS_GAMMA = "GAMMA_RADPRESS"
+        EOS_GAMMA = "RADPRESS"
         GAMMA = 4./3.
     elif GAMMA_GASPRESS:
-        EOS_GAMMA = "GAMMA_GASPRESS"
+        EOS_GAMMA = "GASPRESS"
         GAMMA = 5./3.
     else:
         print("Using the default gas-pressure dominated ideal gas EoS. 
               You can choose gas-pressure dominated (add flag -gammagaspress),
               radiation-pressure dominated (add flag -gammaradpress).")
-        EOS_GAMMA = "GAMMA_GASPRESS"
+        EOS_GAMMA = "GASPRESS"
         GAMMA = 4./3.
 else:
     raise ValueError("Bad EOS chosen")
@@ -175,7 +175,7 @@
     ENTROPY = float(sys.argv[sys.argv.index('-ent') + 1])
 else:
     ENTROPY = 4
-        if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "GAMMA_RADPRESS":
+        if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "RADPRESS":
         print("Using the default initial entropy = 4.", end=' ')
         print("Entropy is a required parameter in setting up radiation-pressure dominated disks.")
 
@@ -185,7 +185,7 @@
     KAPPA_EOS = float(sys.argv[sys.argv.index('-kappa') + 1])
 else:
     KAPPA_EOS = 1.e-3
-    if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "GAMMA_GASPRESS":
+    if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "GASPRESS":
         print("Using the default initial kappa = 1.e-3.", end=' ')
         print("kappa is a required parameter in setting up gas-pressure dominated disks.")
 

From 461f50518251e869f8bbd135ef7f388f6f945613 Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Mon, 1 Mar 2021 20:45:49 -0700
Subject: [PATCH 05/14] GAMMA_EOS --> EOS_GAMMA in all occurences

---
 prob/torus_cbc/build.py  |  4 ++--
 prob/torus_cbc/problem.c | 10 +++++-----
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/prob/torus_cbc/build.py b/prob/torus_cbc/build.py
index 05d6644..57bdffc 100644
--- a/prob/torus_cbc/build.py
+++ b/prob/torus_cbc/build.py
@@ -502,9 +502,9 @@
         bhl.config.set_rparm('lrho_guess','double',LRHO_GUESS)
 if USE_GAMMA or GAMMA_FALLBACK:
     bhl.config.set_rparm('gam', 'double', default = GAMMA)
-    if EOS_TYPE == "EOS_TYPE_GAMMA" and GAMMA_EOS == "RADPRESS":
+    if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "RADPRESS":
         bhl.config.set_rparm('entropy','double',ENTROPY)
-    elif EOS_TYPE == "EOS_TYPE_GAMMA" and GAMMA_EOS == "GASPRESS":
+    elif EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "GASPRESS":
         bhl.config.set_rparm("kappa_eos", "double", default = KAPPA_EOS)
 
 # Opacities
diff --git a/prob/torus_cbc/problem.c b/prob/torus_cbc/problem.c
index 26848ab..18d96ca 100644
--- a/prob/torus_cbc/problem.c
+++ b/prob/torus_cbc/problem.c
@@ -17,7 +17,7 @@ static int    renormalize_densities;
 static double rin;
 static double rmax;
 static double beta; // desired minimum ratio of gas to magnetic pressure
-#if (EOS == EOS_TYPE_GAMMA && GAMMA_EOS == GASPRESS) || GAMMA_FALLBACK
+#if (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == GASPRESS) || GAMMA_FALLBACK
 static double kappa_eos;
 #endif
 #if EOS == EOS_TYPE_TABLE
@@ -44,7 +44,7 @@ void set_problem_params() {
   set_param("rmax", &rmax);
   set_param("beta", &beta);
   set_param("renorm_dens", &renormalize_densities);
-#if (EOS_TYPE == EOS_TYPE_GAMMA && GAMMA_EOS == GASPRESS) || GAMMA_FALLBACK
+#if (EOS_TYPE == EOS_TYPE_GAMMA && EOS_GAMMA == GASPRESS) || GAMMA_FALLBACK
   set_param("kappa_eos", &kappa_eos);
 #endif
 #if EOS_TYPE == EOS_TYPE_TABLE
@@ -221,11 +221,11 @@ void init_prob() {
       disk_cell[i][j][k] = 1;
 
       hm1 = exp(lnh[i][j][k]) - 1.;
-#if EOS == EOS_TYPE_GAMMA && GAMMA_EOS == RADPRESS
+#if EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS
       fprintf(stdout, "entropy: %f\n", entropy);
       rho = (64./3) * (pow(hm1, 3)/pow(entropy, 4)) * (pow((gam - 1.)/gam), 3);
       u = hm1*rho/gam;
-#elif (EOS == EOS_TYPE_GAMMA && GAMMA_EOS == GASPRESS) || GAMMA_FALLBACK
+#elif (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == GASPRESS) || GAMMA_FALLBACK
       rho = pow(hm1 * (gam - 1.) / (kappa_eos * gam), 1. / (gam - 1.));
       u   = kappa_eos * pow(rho, gam) / (gam - 1.);
 #elif EOS == EOS_TYPE_TABLE
@@ -382,7 +382,7 @@ void init_prob() {
         mtot, mtot * M_unit, mtot * M_unit / MSUN);
   }
 // debug
-#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && GAMMA_EOS == RADPRESS)
+#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
   if (mpi_io_proc()) {
     fprintf(stdout, "Calculating max entropy:\n");
   }

From 24aee2a8ef92e1795425e0b69091b1b9ec9c8672 Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Mon, 1 Mar 2021 20:49:34 -0700
Subject: [PATCH 06/14] EOS_TYPE --> EOS everywhere

---
 prob/torus_cbc/problem.c | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/prob/torus_cbc/problem.c b/prob/torus_cbc/problem.c
index 18d96ca..9105438 100644
--- a/prob/torus_cbc/problem.c
+++ b/prob/torus_cbc/problem.c
@@ -44,10 +44,10 @@ void set_problem_params() {
   set_param("rmax", &rmax);
   set_param("beta", &beta);
   set_param("renorm_dens", &renormalize_densities);
-#if (EOS_TYPE == EOS_TYPE_GAMMA && EOS_GAMMA == GASPRESS) || GAMMA_FALLBACK
+#if (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == GASPRESS) || GAMMA_FALLBACK
   set_param("kappa_eos", &kappa_eos);
 #endif
-#if EOS_TYPE == EOS_TYPE_TABLE
+#if EOS == EOS_TYPE_TABLE
   set_param("const_ye", &const_ye);
 #if !GAMMA_FALLBACK
   set_param("entropy", &entropy);
@@ -68,7 +68,7 @@ void init_prob() {
   double SSin, hm1;
 
 // Diagnostics for entropy
-#if EOS_TYPE == EOS_TYPE_TABLE
+#if EOS == EOS_TYPE_TABLE
   double ent, entmax;
 #endif
 
@@ -82,7 +82,7 @@ void init_prob() {
   double rhomin, umin;
   rhomin = RHOMINLIMIT;
   umin   = UUMINLIMIT;
-#if EOS_TYPE == EOS_TYPE_TABLE
+#if EOS == EOS_TYPE_TABLE
   umin = EOS_SC_get_minu(rhomin, const_ye, umin);
 #endif
 
@@ -94,7 +94,7 @@ void init_prob() {
   PASSTYPE(PASSIVE_START + 1) = PASSTYPE_NUMBER;
 #endif
 
-#if EOS_TYPE == EOS_TYPE_TABLE
+#if EOS == EOS_TYPE_TABLE
   double ye, ye_atm;
 #if !GAMMA_FALLBACK
   // guesses. Should be function local.

From d136659b7766dfb8659b94c2597c3298beb96be6 Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Tue, 9 Mar 2021 00:59:36 -0700
Subject: [PATCH 07/14] Update the remaining code with the two updated
 EOS_TYPE_GAMMA options

---
 core/eos.c              | 10 ++++++----
 core/fixup.c            |  6 +++++-
 core/io.c               |  6 +++---
 core/main.c             | 11 +++++------
 core/phys.c             |  2 +-
 core/utop.c             |  2 +-
 prob/torus_cbc/build.py | 16 ++++++++--------
 7 files changed, 29 insertions(+), 24 deletions(-)

diff --git a/core/eos.c b/core/eos.c
index 358bb0e..a35924e 100644
--- a/core/eos.c
+++ b/core/eos.c
@@ -159,10 +159,12 @@ double EOS_enthalpy_rho0_u(double rho, double u, const double *extra) {
 
 double EOS_entropy_rho0_u(double rho, double u, const double *extra) {
   double ent;
-#if EOS == EOS_TYPE_GAMMA_GASPRESS
-  ent = EOS_Gamma_Gaspress_entropy_rho0_u(rho, u);
-#elif EOS == EOS_TYPE_GAMMA_RADPRESS
-  ent = EOS_Gamma_Radpress_entropy_rho0_u(rho, u);
+#if EOS == EOS_TYPE_GAMMA
+  #if EOS_GAMMA == GASPRESS
+    ent = EOS_Gamma_Gaspress_entropy_rho0_u(rho, u);
+  #elif EOS_GAMMA == RADPRESS
+    ent = EOS_Gamma_Radpress_entropy_rho0_u(rho, u);
+  #endif
 #elif EOS == EOS_TYPE_POLYTROPE
   ent   = EOS_Poly_entropy_rho0_u(rho, u, poly_K, poly_gam);
 #elif EOS == EOS_TYPE_TABLE
diff --git a/core/fixup.c b/core/fixup.c
index d027ed4..2791a9f 100644
--- a/core/fixup.c
+++ b/core/fixup.c
@@ -229,7 +229,7 @@ void fixup1zone(
   fixup_passive(i, j, k, pv, pv_prefloor);
 #endif
 
-#if ELECTRONS && EOS == EOS_TYPE_GAMMA_GASPRESS  
+#if ELECTRONS && EOS == EOS_TYPE_GAMMA && GAMMA_EOS == GASPRESS
   // Reset entropy after floors
   pv[KTOT] = EOS_Gamma_Gaspress_entropy_rho0_u(pv[RHO], pv[UU]);
 
@@ -240,6 +240,10 @@ void fixup1zone(
     pv[UU]   = KTOTMAX * pow(pv[RHO], gam) / (gam - 1.);
     pv[KTOT] = KTOTMAX;
   }
+
+#elif ELECTRONS && EOS == EOS_TYPE_GAMMA && GAMMA_EOS == RADPRESS
+  // Reset entropy after floors
+  pv[KTOT] = EOS_Gamma_Radpress_entropy_rho0_u(pv[RHO], pv[UU]);
 #endif // ELECTRONS
 
   // Limit gamma with respect to normal observer
diff --git a/core/io.c b/core/io.c
index ed3a82a..fb8f2d4 100644
--- a/core/io.c
+++ b/core/io.c
@@ -730,7 +730,7 @@ void dump() {
 #endif // RADIATION
 #endif // METRIC
 
-#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  
+#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK  
   WRITE_HDR(gam, TYPE_DBL);
 #endif // GAMMA
 #if EOS == EOS_TYPE_POLYTROPE
@@ -1099,7 +1099,7 @@ void restart_write(int restart_type) {
 #endif // METRIC
 
 // EOS
-#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  
+#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK  
   WRITE_HDR(gam, TYPE_DBL);
 #endif // GAMMA EOS
 #if EOS == EOS_TYPE_POLYTROPE
@@ -1381,7 +1381,7 @@ void restart_read(char *fname) {
 
   READ_HDR(cour, TYPE_DBL);
 
-#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  
+#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK  
   READ_HDR(gam, TYPE_DBL);
 #endif // EOS
 
diff --git a/core/main.c b/core/main.c
index fedc7c1..0773642 100644
--- a/core/main.c
+++ b/core/main.c
@@ -59,12 +59,11 @@ int main(int argc, char *argv[]) {
     #endif
     #if EOS == EOS_TYPE_GAMMA
     fprintf(stdout, "          *  IDEAL GAS EOS                                           *\n");
-    #endif
-    #if EOS == EOS_TYPE_GAMMA_GASPRESS
-    fprintf(stdout, "          *  GAS-PRESSURE DOMINATED EOS                    *\n");
-    #endif
-    #if EOS == EOS_TYPE_GAMMA_RADPRESS
-    fprintf(stdout, "          *  RADIATION-PRESSURE DOMINATED EOS              *\n");
+      #if EOS_GAMMA == GASPRESS
+        fprintf(stdout, "          *  GAS-PRESSURE DOMINATED IDEAL GAS EOS                    *\n");
+      #elif EOS_GAMMA == RADPRESS
+        fprintf(stdout, "          *  RADIATION-PRESSURE DOMINATED IDEAL GAS EOS              *\n");
+      #endif
     #endif
     #if RADIATION == RADTYPE_LIGHT
     fprintf(stdout, "          *  RADIATION: PHOTON TRANSPORT                             *\n");
diff --git a/core/phys.c b/core/phys.c
index fbbcc4d..b181d15 100644
--- a/core/phys.c
+++ b/core/phys.c
@@ -211,7 +211,7 @@ void mhd_vchar(double *Pr, struct of_state *q, struct of_geom *geom, int js,
   bsq = dot(q->bcon, q->bcov);
   rho = fabs(Pr[RHO] + SMALL);
   u   = Pr[UU];
-#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  
+#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK  
   u = fabs(u + SMALL);
 #endif
   ef  = EOS_enthalpy_rho0_u(rho, u, extra);
diff --git a/core/utop.c b/core/utop.c
index 3c946b9..a3e9224 100644
--- a/core/utop.c
+++ b/core/utop.c
@@ -161,7 +161,7 @@ int Utoprim(double U[NVAR], struct of_geom *geom, double prim[NVAR]) {
   // Return without updating non-B primitives
   if (rho0 < 0)
     return (5);
-#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  
+#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK  
   if (u < 0)
     return (5);
 #endif
diff --git a/prob/torus_cbc/build.py b/prob/torus_cbc/build.py
index 57bdffc..2cbc96f 100644
--- a/prob/torus_cbc/build.py
+++ b/prob/torus_cbc/build.py
@@ -163,9 +163,9 @@
         EOS_GAMMA = "GASPRESS"
         GAMMA = 5./3.
     else:
-        print("Using the default gas-pressure dominated ideal gas EoS. 
-              You can choose gas-pressure dominated (add flag -gammagaspress),
-              radiation-pressure dominated (add flag -gammaradpress).")
+        print('Using the default gas-pressure dominated ideal gas EoS.',
+              'You can choose gas-pressure dominated (add flag -gammagaspress)',
+              'radiation-pressure dominated (add flag -gammaradpress).')
         EOS_GAMMA = "GASPRESS"
         GAMMA = 4./3.
 else:
@@ -175,9 +175,9 @@
     ENTROPY = float(sys.argv[sys.argv.index('-ent') + 1])
 else:
     ENTROPY = 4
-        if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "RADPRESS":
-        print("Using the default initial entropy = 4.", end=' ')
-        print("Entropy is a required parameter in setting up radiation-pressure dominated disks.")
+    if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "RADPRESS":
+        print('Using the default initial entropy = 4.',
+              'Entropy is a required parameter in setting up radiation-pressure dominated disks.')
 
 bhl.report_var('ENTROPY',ENTROPY)
 
@@ -186,8 +186,8 @@
 else:
     KAPPA_EOS = 1.e-3
     if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "GASPRESS":
-        print("Using the default initial kappa = 1.e-3.", end=' ')
-        print("kappa is a required parameter in setting up gas-pressure dominated disks.")
+        print('Using the default initial kappa = 1.e-3.',
+              'kappa is a required parameter in setting up gas-pressure dominated disks.')
 
 if CLASSIC: # classic harm disk
     # Rmax and rho fixed

From d51238ceb54874f34d976c9a554cf6408f21d534 Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Tue, 9 Mar 2021 14:26:50 -0700
Subject: [PATCH 08/14] Trying to fix ent uninitialized err from unittest

---
 core/eos.c | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/core/eos.c b/core/eos.c
index a35924e..519c0f8 100644
--- a/core/eos.c
+++ b/core/eos.c
@@ -160,11 +160,11 @@ double EOS_enthalpy_rho0_u(double rho, double u, const double *extra) {
 double EOS_entropy_rho0_u(double rho, double u, const double *extra) {
   double ent;
 #if EOS == EOS_TYPE_GAMMA
-  #if EOS_GAMMA == GASPRESS
-    ent = EOS_Gamma_Gaspress_entropy_rho0_u(rho, u);
-  #elif EOS_GAMMA == RADPRESS
-    ent = EOS_Gamma_Radpress_entropy_rho0_u(rho, u);
-  #endif
+#if EOS_GAMMA == GASPRESS
+  ent = EOS_Gamma_Gaspress_entropy_rho0_u(rho, u);
+#elif EOS_GAMMA == RADPRESS
+  ent = EOS_Gamma_Radpress_entropy_rho0_u(rho, u);
+#endif
 #elif EOS == EOS_TYPE_POLYTROPE
   ent   = EOS_Poly_entropy_rho0_u(rho, u, poly_K, poly_gam);
 #elif EOS == EOS_TYPE_TABLE

From f42c723875ec582c627bf83042ef94cb2bf9f421 Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Tue, 9 Mar 2021 15:01:20 -0700
Subject: [PATCH 09/14] Add a default value for EOS_GAMMA

---
 script/config.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/script/config.py b/script/config.py
index cf6a8c5..ec5041d 100644
--- a/script/config.py
+++ b/script/config.py
@@ -349,6 +349,10 @@ def build(PROBLEM, PATHS):
      and CPARMS['EOS'] != 'EOS_TYPE_GAMMA':
     raise ValueError("ELECTRONS only compatible with Gamma law EOS.\n"
                      +"Please set EOS = EOS_TYPE_GAMMA.\n")
+  if CPARMS['EOS'] == 'EOS_TYPE_GAMMA':
+    if not util.parm_is_active(CPARMS, 'EOS_GAMMA'):
+      set_cparm('EOS_GAMMA', "GASPRESS")
+      print_config("Gamma law entropy type", CPARMS["EOS_GAMMA"])
 
   if CPARMS['EOS'] == 'EOS_TYPE_TABLE' and CPARMS['NVAR_PASSIVE'] < 2:
     raise ValueError("Tabulated EOS requires at least two passive scalars\n"

From a76212a12deceda1120a702fd223473816483825 Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Wed, 10 Mar 2021 19:58:59 -0700
Subject: [PATCH 10/14] Multiplied in rad constant in rho eqn and added rad
 dominated option in EOS_Gamma_temp

---
 core/decs.h              |  2 +-
 core/defs.h              |  2 +-
 core/io.c                |  2 ++
 core/util.c              |  2 ++
 prob/torus_cbc/problem.c | 15 ++++++++++-----
 5 files changed, 16 insertions(+), 7 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index 7a2e9b5..644a3c8 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -429,7 +429,7 @@ extern double M_unit;
 extern double Reh;
 extern double Risco;
 #if NEED_UNITS
-extern double mbh, Mbh, L_unit, T_unit, M_unit, RHO_unit, U_unit, B_unit;
+extern double mbh, Mbh, L_unit, T_unit, M_unit, RHO_unit, U_unit, B_unit, TEMP_unit;
 #endif
 #if EOS == EOS_TYPE_TABLE
 extern double TEMP_unit;
diff --git a/core/defs.h b/core/defs.h
index a26d977..e077ab8 100644
--- a/core/defs.h
+++ b/core/defs.h
@@ -114,7 +114,7 @@ double Reh;
 double Risco;
 
 #if NEED_UNITS
-double mbh, Mbh, L_unit, T_unit, M_unit, RHO_unit, U_unit, B_unit;
+double mbh, Mbh, L_unit, T_unit, M_unit, RHO_unit, U_unit, B_unit, TEMP_unit;
 #endif
 
 #if EOS == EOS_TYPE_TABLE
diff --git a/core/io.c b/core/io.c
index fb8f2d4..cc47c59 100644
--- a/core/io.c
+++ b/core/io.c
@@ -753,6 +753,7 @@ void dump() {
   WRITE_HDR(T_unit, TYPE_DBL);
   WRITE_HDR(U_unit, TYPE_DBL);
   WRITE_HDR(B_unit, TYPE_DBL);
+  WRITE_HDR(TEMP_unit, TYPE_DBL);
 #if EOS == EOS_TYPE_TABLE
   WRITE_HDR(TEMP_unit, TYPE_DBL);
 #endif // EOS_TYPE_TABLE
@@ -1117,6 +1118,7 @@ void restart_write(int restart_type) {
   WRITE_HDR(T_unit, TYPE_DBL);
   WRITE_HDR(U_unit, TYPE_DBL);
   WRITE_HDR(B_unit, TYPE_DBL);
+  WRITE_HDR(TEMP_unit, TYPE_DBL);
 #if EOS == EOS_TYPE_TABLE
   WRITE_HDR(TEMP_unit, TYPE_DBL);
 #endif // EOS_TYPE_TABLE
diff --git a/core/util.c b/core/util.c
index 73c0af5..3a2bd68 100644
--- a/core/util.c
+++ b/core/util.c
@@ -145,6 +145,8 @@ void set_units() {
   B_unit   = CL * sqrt(4. * M_PI * RHO_unit);
 #if EOS == EOS_TYPE_TABLE
   TEMP_unit = MEV;
+#elif EOS == EOS_TYPE_GAMMA
+  TEMP_unit = GNEWT * KBOL / pow(CL, 4);
 #endif
 
 #if RADIATION
diff --git a/prob/torus_cbc/problem.c b/prob/torus_cbc/problem.c
index 9105438..fa2b151 100644
--- a/prob/torus_cbc/problem.c
+++ b/prob/torus_cbc/problem.c
@@ -23,10 +23,12 @@ static double kappa_eos;
 #if EOS == EOS_TYPE_TABLE
 static double const_ye; // if > 0, set ye this way
 #if !GAMMA_FALLBACK
-static double entropy;
 static double lrho_guess;
 #endif
 #endif
+#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
+static double entropy;
+#endif
 #if RADIATION && TRACERS
 static int ntracers;
 #endif
@@ -68,7 +70,7 @@ void init_prob() {
   double SSin, hm1;
 
 // Diagnostics for entropy
-#if EOS == EOS_TYPE_TABLE
+#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
   double ent, entmax;
 #endif
 
@@ -222,8 +224,9 @@ void init_prob() {
 
       hm1 = exp(lnh[i][j][k]) - 1.;
 #if EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS
-      fprintf(stdout, "entropy: %f\n", entropy);
-      rho = (64./3) * (pow(hm1, 3)/pow(entropy, 4)) * (pow((gam - 1.)/gam), 3);
+      a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4)/M_unit;
+      fprintf(stdout, "radiation density constant: %f\n", a);
+      rho = (64./3) * a * (pow(hm1, 3)/pow(entropy, 4)) * pow((gam - 1.)/gam, 3);
       u = hm1*rho/gam;
 #elif (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == GASPRESS) || GAMMA_FALLBACK
       rho = pow(hm1 * (gam - 1.) / (kappa_eos * gam), 1. / (gam - 1.));
@@ -391,7 +394,9 @@ void init_prob() {
     coord(i, j, k, CENT, X);
     bl_coord(X, &r, &th);
     if (r > rin && lnh[i][j][k] >= 0.) {
-      EOS_SC_fill(P[i][j][k], extra[i][j][k]);
+      #if EOS == EOS_TYPE_TABLE
+        EOS_SC_fill(P[i][j][k], extra[i][j][k]);
+      #endif
       ent = EOS_entropy_rho0_u(P[i][j][k][RHO], P[i][j][k][UU], extra[i][j][k]);
       if (ent > entmax)
         entmax = ent;

From 33909467bbfe5d23b69eb226b94f0e3246c54a23 Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Wed, 10 Mar 2021 20:01:11 -0700
Subject: [PATCH 11/14] Forgot to add eos_gamma in the prev commit

---
 core/eos_gamma.c | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/core/eos_gamma.c b/core/eos_gamma.c
index bf2abb0..1c0623f 100644
--- a/core/eos_gamma.c
+++ b/core/eos_gamma.c
@@ -22,7 +22,8 @@ double EOS_Gamma_Gaspress_entropy_rho0_u(double rho, double u) {
 }
 
 double EOS_Gamma_Radpress_entropy_rho0_u(double rho, double u) {
-  return pow((64./3) * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
+  double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4)/M_unit;
+  return pow((64./3) * a * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
 }
 
 double EOS_Gamma_enthalpy_rho0_u(double rho, double u) { return rho + u * gam; }
@@ -66,8 +67,15 @@ double EOS_Gamma_sound_speed_rho0_u(double rho, double u) {
 double EOS_Gamma_u_press(double press) { return press / (gam - 1.); }
 
 double EOS_Gamma_temp(double rho, double u) {
-  // return TEMP_unit*(gam - 1.)*u/rho;
-  return (gam - 1.) * u / rho;
+  #if EOS_GAMMA == RADPRESS
+    double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4)/M_unit;
+    double press = (gam - 1.) * u;
+    double temp = pow((3 * press)/a, 1./4);
+  #else
+    // return TEMP_unit*(gam - 1.)*u/rho;
+    double temp = (gam - 1.) * u / rho;
+  #endif
+  return temp;
 }
 
 #if RADIATION

From 06164ffa60c58b7cfaab5ed5a333ea18ca98bd3d Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Wed, 17 Mar 2021 22:11:09 -0600
Subject: [PATCH 12/14] edits to make entropy be declared in problem.c and some
 debugging prints

---
 core/decs.h              |  3 +++
 core/defs.h              |  3 +++
 core/eos_gamma.c         | 10 ++++++----
 core/input.c             |  4 ++++
 prob/torus_cbc/problem.c | 25 ++++++++++++++++++-------
 5 files changed, 34 insertions(+), 11 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index 644a3c8..4c267d8 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -419,6 +419,9 @@ extern double a;
 #if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
 extern double gam;
 #endif
+#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
+extern double entropy;
+#endif
 #if EOS == EOS_TYPE_POLYTROPE
 extern double poly_K, poly_gam;
 #endif
diff --git a/core/defs.h b/core/defs.h
index e077ab8..fa5d00e 100644
--- a/core/defs.h
+++ b/core/defs.h
@@ -101,6 +101,9 @@ char opac_file[STRLEN];
 #if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
 double gam;
 #endif
+#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
+double entropy;
+#endif
 #if EOS == EOS_TYPE_POLYTROPE
 double poly_K, poly_gam;
 #endif
diff --git a/core/eos_gamma.c b/core/eos_gamma.c
index 1c0623f..be66c6e 100644
--- a/core/eos_gamma.c
+++ b/core/eos_gamma.c
@@ -22,8 +22,9 @@ double EOS_Gamma_Gaspress_entropy_rho0_u(double rho, double u) {
 }
 
 double EOS_Gamma_Radpress_entropy_rho0_u(double rho, double u) {
-  double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4)/M_unit;
-  return pow((64./3) * a * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
+  double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
+  //return pow((64./3) * a * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
+  return pow((64./3) * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
 }
 
 double EOS_Gamma_enthalpy_rho0_u(double rho, double u) { return rho + u * gam; }
@@ -68,9 +69,10 @@ double EOS_Gamma_u_press(double press) { return press / (gam - 1.); }
 
 double EOS_Gamma_temp(double rho, double u) {
   #if EOS_GAMMA == RADPRESS
-    double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4)/M_unit;
+    double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
     double press = (gam - 1.) * u;
-    double temp = pow((3 * press)/a, 1./4);
+    //double temp = pow((3 * press)/a, 1./4);
+    double temp = pow((3 * press), 1./4);
   #else
     // return TEMP_unit*(gam - 1.)*u/rho;
     double temp = (gam - 1.) * u / rho;
diff --git a/core/input.c b/core/input.c
index 7e3989d..1c039f6 100644
--- a/core/input.c
+++ b/core/input.c
@@ -95,6 +95,10 @@ void set_core_params() {
 #if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
   set_param("gam", &gam);
 #endif
+#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
+  set_param("entropy", &entropy);
+#endif
+
 #if EOS == EOS_TYPE_GAMMA
   sprintf(eos, "GAMMA");
 #elif EOS == EOS_TYPE_POLYTROPE
diff --git a/prob/torus_cbc/problem.c b/prob/torus_cbc/problem.c
index fa2b151..f11f70c 100644
--- a/prob/torus_cbc/problem.c
+++ b/prob/torus_cbc/problem.c
@@ -26,9 +26,9 @@ static double const_ye; // if > 0, set ye this way
 static double lrho_guess;
 #endif
 #endif
-#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
-static double entropy;
-#endif
+//#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
+//static double entropy;
+//#endif
 #if RADIATION && TRACERS
 static int ntracers;
 #endif
@@ -52,7 +52,7 @@ void set_problem_params() {
 #if EOS == EOS_TYPE_TABLE
   set_param("const_ye", &const_ye);
 #if !GAMMA_FALLBACK
-  set_param("entropy", &entropy);
+  //set_param("entropy", &entropy);
   set_param("lrho_guess", &lrho_guess);
 #endif
 #endif
@@ -180,6 +180,7 @@ void init_prob() {
       lnh[i][j][k] = 1.;
     }
 
+//fprintf(stdout, "entropy: %f\n", entropy);
 #if EOS == EOS_TYPE_TABLE
     ye     = const_ye; // may need to change this eventually
     ye_atm = 0.5;
@@ -224,9 +225,19 @@ void init_prob() {
 
       hm1 = exp(lnh[i][j][k]) - 1.;
 #if EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS
-      a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4)/M_unit;
-      fprintf(stdout, "radiation density constant: %f\n", a);
-      rho = (64./3) * a * (pow(hm1, 3)/pow(entropy, 4)) * pow((gam - 1.)/gam, 3);
+      fprintf(stdout, "AR: %e\n", AR);
+      fprintf(stdout, "T_unit: %e\n", T_unit);
+      fprintf(stdout, "L_unit: %e\n", L_unit);
+      fprintf(stdout, "TEMP_unit: %e\n", TEMP_unit);
+      fprintf(stdout, "M_unit: %e\n", M_unit);
+      double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
+      fprintf(stdout, "radiation density constant: %e\n", a);
+      fprintf(stdout, "hm1: %e\n", hm1);
+      fprintf(stdout, "entropy: %f\n", entropy);
+      fprintf(stdout, "gam: %f\n", gam);
+      //rho = (64./3) * a * (pow(hm1, 3)/pow(entropy, 4)) * pow((gam - 1.)/gam, 3);
+      rho = (64./3) * (pow(hm1, 3)/pow(entropy, 4)) * pow((gam - 1.)/gam, 3);
+      fprintf(stdout, "rho: %e\n", rho);
       u = hm1*rho/gam;
 #elif (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == GASPRESS) || GAMMA_FALLBACK
       rho = pow(hm1 * (gam - 1.) / (kappa_eos * gam), 1. / (gam - 1.));

From 87b77cb72fd0f86c840d49e8780b95a6c959180a Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Tue, 4 May 2021 01:41:50 -0600
Subject: [PATCH 13/14] switching back to correctly initializing entropy in
 problem.c instead of globally + experimenting with entropy CU

---
 core/decs.h              |  6 +++---
 core/defs.h              |  6 +++---
 core/eos_gamma.c         |  8 ++++----
 core/input.c             |  6 +++---
 core/util.c              |  3 ++-
 prob/torus_cbc/problem.c | 29 +++++++++++++++++++++--------
 6 files changed, 36 insertions(+), 22 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index 4c267d8..ea38928 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -419,9 +419,9 @@ extern double a;
 #if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
 extern double gam;
 #endif
-#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
-extern double entropy;
-#endif
+//#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
+//extern double entropy;
+//#endif
 #if EOS == EOS_TYPE_POLYTROPE
 extern double poly_K, poly_gam;
 #endif
diff --git a/core/defs.h b/core/defs.h
index fa5d00e..43f1a83 100644
--- a/core/defs.h
+++ b/core/defs.h
@@ -101,9 +101,9 @@ char opac_file[STRLEN];
 #if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
 double gam;
 #endif
-#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
-double entropy;
-#endif
+//#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
+//double entropy;
+//#endif
 #if EOS == EOS_TYPE_POLYTROPE
 double poly_K, poly_gam;
 #endif
diff --git a/core/eos_gamma.c b/core/eos_gamma.c
index be66c6e..a9c8807 100644
--- a/core/eos_gamma.c
+++ b/core/eos_gamma.c
@@ -23,8 +23,8 @@ double EOS_Gamma_Gaspress_entropy_rho0_u(double rho, double u) {
 
 double EOS_Gamma_Radpress_entropy_rho0_u(double rho, double u) {
   double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
-  //return pow((64./3) * a * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
-  return pow((64./3) * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
+  return pow((64./3) * a * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
+  //return pow((64./3) * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
 }
 
 double EOS_Gamma_enthalpy_rho0_u(double rho, double u) { return rho + u * gam; }
@@ -71,8 +71,8 @@ double EOS_Gamma_temp(double rho, double u) {
   #if EOS_GAMMA == RADPRESS
     double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
     double press = (gam - 1.) * u;
-    //double temp = pow((3 * press)/a, 1./4);
-    double temp = pow((3 * press), 1./4);
+    double temp = pow((3 * press)/a, 1./4);
+    //double temp = pow((3 * press), 1./4);
   #else
     // return TEMP_unit*(gam - 1.)*u/rho;
     double temp = (gam - 1.) * u / rho;
diff --git a/core/input.c b/core/input.c
index 1c039f6..fe7f868 100644
--- a/core/input.c
+++ b/core/input.c
@@ -95,9 +95,9 @@ void set_core_params() {
 #if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
   set_param("gam", &gam);
 #endif
-#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
-  set_param("entropy", &entropy);
-#endif
+//#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
+//  set_param("entropy", &entropy);
+//#endif
 
 #if EOS == EOS_TYPE_GAMMA
   sprintf(eos, "GAMMA");
diff --git a/core/util.c b/core/util.c
index 3a2bd68..b8c6ffb 100644
--- a/core/util.c
+++ b/core/util.c
@@ -146,7 +146,8 @@ void set_units() {
 #if EOS == EOS_TYPE_TABLE
   TEMP_unit = MEV;
 #elif EOS == EOS_TYPE_GAMMA
-  TEMP_unit = GNEWT * KBOL / pow(CL, 4);
+  //TEMP_unit = GNEWT * KBOL / pow(CL, 4);
+  TEMP_unit = (1. / KBOL) * (M_unit * pow(L_unit, 2.) / pow(T_unit, 2.)); 
 #endif
 
 #if RADIATION
diff --git a/prob/torus_cbc/problem.c b/prob/torus_cbc/problem.c
index f11f70c..31f4270 100644
--- a/prob/torus_cbc/problem.c
+++ b/prob/torus_cbc/problem.c
@@ -26,9 +26,9 @@ static double const_ye; // if > 0, set ye this way
 static double lrho_guess;
 #endif
 #endif
-//#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
-//static double entropy;
-//#endif
+#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
+static double entropy;
+#endif
 #if RADIATION && TRACERS
 static int ntracers;
 #endif
@@ -49,10 +49,13 @@ void set_problem_params() {
 #if (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == GASPRESS) || GAMMA_FALLBACK
   set_param("kappa_eos", &kappa_eos);
 #endif
+#if (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
+  set_param("entropy", &entropy);
+#endif
 #if EOS == EOS_TYPE_TABLE
   set_param("const_ye", &const_ye);
 #if !GAMMA_FALLBACK
-  //set_param("entropy", &entropy);
+  set_param("entropy", &entropy);
   set_param("lrho_guess", &lrho_guess);
 #endif
 #endif
@@ -74,6 +77,10 @@ void init_prob() {
   double ent, entmax;
 #endif
 
+#if (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
+  double entropy_CU;
+#endif
+
   // Magnetic field
   double rho_av, rhomax, umax, pressmax, hm1max, bsq_ij, bsq_max, q;
 
@@ -225,22 +232,27 @@ void init_prob() {
 
       hm1 = exp(lnh[i][j][k]) - 1.;
 #if EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS
+      a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
       fprintf(stdout, "AR: %e\n", AR);
       fprintf(stdout, "T_unit: %e\n", T_unit);
       fprintf(stdout, "L_unit: %e\n", L_unit);
       fprintf(stdout, "TEMP_unit: %e\n", TEMP_unit);
       fprintf(stdout, "M_unit: %e\n", M_unit);
-      double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
       fprintf(stdout, "radiation density constant: %e\n", a);
       fprintf(stdout, "hm1: %e\n", hm1);
-      fprintf(stdout, "entropy: %f\n", entropy);
+      fprintf(stdout, "entropy_CU: %f\n", entropy_CU);
       fprintf(stdout, "gam: %f\n", gam);
-      //rho = (64./3) * a * (pow(hm1, 3)/pow(entropy, 4)) * pow((gam - 1.)/gam, 3);
-      rho = (64./3) * (pow(hm1, 3)/pow(entropy, 4)) * pow((gam - 1.)/gam, 3);
+      entropy_CU = (entropy * KBOL / MP) * pow(T_unit,2) * TEMP_unit * pow(L_unit, -2) * pow(M_unit, -1);
+      rho = (64./3) * a * (pow(hm1, 3)/pow(entropy_CU, 4)) * pow((gam - 1.)/gam, 3);
+      //rho = (64./3) * (pow(hm1, 3)/pow(entropy_CU, 4)) * pow((gam - 1.)/gam, 3);
       fprintf(stdout, "rho: %e\n", rho);
       u = hm1*rho/gam;
 #elif (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == GASPRESS) || GAMMA_FALLBACK
       rho = pow(hm1 * (gam - 1.) / (kappa_eos * gam), 1. / (gam - 1.));
+      fprintf(stdout, "rho: %e\n", rho);
+      fprintf(stdout, "hm1: %e\n", hm1);
+      fprintf(stdout, "kappa_eos: %e\n", kappa_eos);
+      fprintf(stdout, "gam: %f\n", gam);
       u   = kappa_eos * pow(rho, gam) / (gam - 1.);
 #elif EOS == EOS_TYPE_TABLE
       hm1 += hm1_min;
@@ -342,6 +354,7 @@ void init_prob() {
 #endif
       press = EOS_pressure_rho0_u(
           PsaveLocal[i][j][k][RHO], PsaveLocal[i][j][k][UU], extra[i][j][k]);
+      fprintf(stdout, "press: %e\n", press);
       if (press > pressmax)
         pressmax = press;
     }

From 6a331713bde4c33157dc0d18bafdd6523454e339 Mon Sep 17 00:00:00 2001
From: Soumi De <soumide@lanl.gov>
Date: Tue, 25 May 2021 02:49:38 -0600
Subject: [PATCH 14/14] Calculate in CGS then divide by code units

---
 core/eos_gamma.c         | 15 +++++++++------
 prob/torus_cbc/build.py  |  1 +
 prob/torus_cbc/problem.c | 24 +++++++++---------------
 3 files changed, 19 insertions(+), 21 deletions(-)

diff --git a/core/eos_gamma.c b/core/eos_gamma.c
index a9c8807..e5dd0c9 100644
--- a/core/eos_gamma.c
+++ b/core/eos_gamma.c
@@ -22,9 +22,12 @@ double EOS_Gamma_Gaspress_entropy_rho0_u(double rho, double u) {
 }
 
 double EOS_Gamma_Radpress_entropy_rho0_u(double rho, double u) {
-  double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
-  return pow((64./3) * a * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
-  //return pow((64./3) * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
+  //double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
+  double u_cgs = u * U_unit;
+  double rho_cgs = rho * RHO_unit;
+  double entropy_cgs = pow((64./3) * AR * (pow(u_cgs, 3)*pow(gam, 3)/pow(rho_cgs, 4)) * pow((gam -1)/gam, 3), 1./4);
+  //return pow((64./3) * a * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
+  return entropy_cgs / (KBOL / MP);
 }
 
 double EOS_Gamma_enthalpy_rho0_u(double rho, double u) { return rho + u * gam; }
@@ -69,9 +72,9 @@ double EOS_Gamma_u_press(double press) { return press / (gam - 1.); }
 
 double EOS_Gamma_temp(double rho, double u) {
   #if EOS_GAMMA == RADPRESS
-    double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
-    double press = (gam - 1.) * u;
-    double temp = pow((3 * press)/a, 1./4);
+    double press_cgs = (gam - 1.) * u * U_unit;
+    double temp_cgs = pow((3 * press_cgs)/AR, 1./4);
+    double temp = temp_cgs * KBOL / MEV;
     //double temp = pow((3 * press), 1./4);
   #else
     // return TEMP_unit*(gam - 1.)*u/rho;
diff --git a/prob/torus_cbc/build.py b/prob/torus_cbc/build.py
index 2cbc96f..9a20655 100644
--- a/prob/torus_cbc/build.py
+++ b/prob/torus_cbc/build.py
@@ -424,6 +424,7 @@
 
 # EOS
 bhl.config.set_cparm("EOS", EOS_TYPE)
+bhl.config.set_cparm('OUTPUT_EOSVARS', True)
 if EOS_TYPE == 'EOS_TYPE_GAMMA':
     bhl.config.set_cparm("EOS_GAMMA", EOS_GAMMA)
 bhl.config.set_cparm('NVAR_PASSIVE', NVAR_PASSIVE)
diff --git a/prob/torus_cbc/problem.c b/prob/torus_cbc/problem.c
index 31f4270..2dd15ec 100644
--- a/prob/torus_cbc/problem.c
+++ b/prob/torus_cbc/problem.c
@@ -78,7 +78,7 @@ void init_prob() {
 #endif
 
 #if (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
-  double entropy_CU;
+  double entropy_cgs, hm1_cgs;
 #endif
 
   // Magnetic field
@@ -232,21 +232,15 @@ void init_prob() {
 
       hm1 = exp(lnh[i][j][k]) - 1.;
 #if EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS
-      a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
-      fprintf(stdout, "AR: %e\n", AR);
-      fprintf(stdout, "T_unit: %e\n", T_unit);
-      fprintf(stdout, "L_unit: %e\n", L_unit);
-      fprintf(stdout, "TEMP_unit: %e\n", TEMP_unit);
-      fprintf(stdout, "M_unit: %e\n", M_unit);
-      fprintf(stdout, "radiation density constant: %e\n", a);
-      fprintf(stdout, "hm1: %e\n", hm1);
-      fprintf(stdout, "entropy_CU: %f\n", entropy_CU);
-      fprintf(stdout, "gam: %f\n", gam);
-      entropy_CU = (entropy * KBOL / MP) * pow(T_unit,2) * TEMP_unit * pow(L_unit, -2) * pow(M_unit, -1);
-      rho = (64./3) * a * (pow(hm1, 3)/pow(entropy_CU, 4)) * pow((gam - 1.)/gam, 3);
-      //rho = (64./3) * (pow(hm1, 3)/pow(entropy_CU, 4)) * pow((gam - 1.)/gam, 3);
+      //a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
+      hm1_cgs = hm1 * pow(CL, 2);
+      fprintf(stdout, "hm1_cgs: %e\n", hm1_cgs);
+      entropy_cgs = entropy * KBOL / MP;
+      fprintf(stdout, "entropy_cgs: %e\n", entropy_cgs);
+      rho = ((64./3) * AR * (pow(hm1_cgs, 3)/pow(entropy_cgs, 4)) * pow((gam - 1.)/gam, 3)) / RHO_unit;
       fprintf(stdout, "rho: %e\n", rho);
-      u = hm1*rho/gam;
+      u = (hm1 * rho / gam);
+      fprintf(stdout, "u: %e\n", u);
 #elif (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == GASPRESS) || GAMMA_FALLBACK
       rho = pow(hm1 * (gam - 1.) / (kappa_eos * gam), 1. / (gam - 1.));
       fprintf(stdout, "rho: %e\n", rho);