Merge pull request #48 from lanl/dempsey/rad_cgs

Convert radiation problems to cgs

.github/workflows/ci.yml

@@ -81,14 +81,14 @@ jobs:
             --log_file=ci_cpu_log.txt
       - name: Upload logs
         if: always()
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: logs
           path: tst/testing/logs
           retention-days: 3
       - name: Upload figures
         if: always()
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: figs
           path: tst/testing/figs

.github/workflows/nightly.yml

@@ -56,14 +56,14 @@ jobs:
             --log_file=ci_cpu_log.txt
       - name: Upload CPU test log
         if: always()
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: ci_cpu_log.txt
           path: tst/ci_cpu_log.txt
           retention-days: 3
       - name: Upload figures
         if: always()
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: figs
           path: tst/figs

inputs/radiation/rad_shock.in

@@ -14,6 +14,12 @@
 <artemis>
 problem = shock          # name of the pgen
 coordinates = cartesian  # coordinate system
+physical_units = cgs     # not scale-free
+unit_conversion = base   # [L/T/M/Temp] = 1
+time = 1.0e-10           # time unit
+length = 0.01575         # length unit
+mass = 3.906984375e-6    # mass unit (implying dunit=1.0)
+temperature = 2.18e6     # temperature unit
 
 <parthenon/job>
 problem_id = shock  # problem ID: basename of output filenames
@@ -23,65 +29,72 @@ variables = gas.prim.density,  &
             gas.prim.velocity, &
             gas.prim.sie,      &
             field.jaybenne.energy_tally
-file_type = hdf5   # HDF5 data dump
-dt        = 0.01   # time increment between outputs
+file_type = hdf5  # HDF5 data dump
+dt        = 1.0   # time increment between outputs
 
 <parthenon/time>
-nlim       = -1    # cycle limit
-tlim       = 0.05  # time limit
-integrator = rk2   # time integration algorithm
+nlim       = -1   # cycle limit
+tlim       = 6.0  # time limit
+integrator = rk2  # time integration algorithm
 
 <parthenon/mesh>
-nx1    = 512       # Number of zones in X1-direction
-x1min  = -0.01     # minimum value of X1
-x1max  = 0.01      # maximum value of X1
+nx1    = 128       # Number of zones in X1-direction
+x1min  = 0.0       # minimum value of X1
+x1max  = 1.0       # maximum value of X1
 ix1_bc = ic        # inner-X1 boundary flag
 ox1_bc = ic        # outer-X1 boundary flag
 
 nx2    = 1         # Number of zones in X2-direction
-x2min  = -0.01     # minimum value of X2
-x2max  = 0.01      # maximum value of X2
+x2min  = 0.0       # minimum value of X2
+x2max  = 1.0       # maximum value of X2
 ix2_bc = periodic  # inner-X2 boundary flag
 ox2_bc = periodic  # outer-X2 boundary flag
 
 nx3    = 1         # Number of zones in X3-direction
-x3min  = -0.01     # minimum value of X3
-x3max  = 0.01      # maximum value of X3
+x3min  = 0.0       # minimum value of X3
+x3max  = 1.0       # maximum value of X3
 ix3_bc = periodic  # inner-X3 boundary flag
 ox3_bc = periodic  # outer-X3 boundary flag
 
 <parthenon/swarm>
-ix1_bc = jaybenne_reflecting
-ox1_bc = jaybenne_reflecting
-ix2_bc = periodic
-ox2_bc = periodic
-ix3_bc = periodic
-ox3_bc = periodic
+ix1_bc = jaybenne_reflecting  # inner-X1 boundary flag for swarms
+ox1_bc = jaybenne_reflecting  # outer-X1 boundary flag for swarms
+ix2_bc = periodic  # inner-X2 boundary flag for swarms
+ox2_bc = periodic  # outer-X2 boundary flag for swarms
+ix3_bc = periodic  # inner-X3 boundary flag for swarms
+ox3_bc = periodic  # outer-X3 boundary flag for swarms
 
 <parthenon/meshblock>
-nx1 = 64  # Number of cells in each MeshBlock, X1-dir
-nx2 = 1   # Number of cells in each MeshBlock, X2-dir
-nx3 = 1   # Number of cells in each MeshBlock, X3-dir
+nx1 = 128  # Number of cells in each MeshBlock, X1-dir
+nx2 = 1    # Number of cells in each MeshBlock, X2-dir
+nx3 = 1    # Number of cells in each MeshBlock, X3-dir
 
 <physics>
-gas = true
-radiation = true
+gas = true        # enable gas hydrodynamics
+radiation = true  # enable IMC radiation
 
 <gas>
-gamma = 1.666666  # adiabatic index
-cv = 1.5          # specific heat
-cfl = 0.8
-reconstruct = plm
-riemann = hllc
-dfloor = 1.0e-10
-siefloor = 1.0e-10
+gamma = 1.666666    # adiabatic index
+mu = 1.008          # mean molecular weight
+cfl = 0.8           # CFL number
+reconstruct = plm   # reconstruction algorithm
+riemann = hllc      # Riemann solver
 
 <gas/opacity/absorption>
-opacity_model = shocktube_a  # absorption opacity model
-coef_kappa_a = 577.35  # kappa_a,0 constant
-rho_exp = -1.0         # dens exponent for opacity powerlaw
-temp_exp = 0.0         # temp exponent for opacity powerlaw
+opacity_model = powerlaw  # absorption opacity model
+coef_kappa_a = 577.0      # kappa_a,0 constant (must be passed in cgs)
+rho_exp = -1.0            # dens exponent for opacity powerlaw
+temp_exp = 0.0            # temp exponent for opacity powerlaw
 
 <jaybenne>
-num_particles = 10000
-use_ddmc = false
+num_particles = 100000  # particle resolution
+use_ddmc = false        # use DDMC?
+
+<problem>
+rhol = 5.69      # density (left)
+rhor = 17.1      # density (right)
+vxl = 0.329524   # x-velocity (left)
+vxr = 0.109841   # x-velocity (right)
+tl = 1.0         # gas and radiation temp (left)
+tr = 3.66055     # gas and radiation temp (right)
+xdisc = 0.828571 # discontinuity position

inputs/radiation/thermalization.in

@@ -13,6 +13,8 @@
 
 <artemis>
 problem = thermalization  # name of the pgen
+physical_units = cgs
+unit_conversion = base
 coordinates = cartesian   # coordinate system
 
 <parthenon/job>
@@ -24,12 +26,12 @@ variables = gas.prim.density,  &
             gas.prim.sie,      &
             field.jaybenne.energy_tally
 file_type = hdf5   # HDF5 data dump
-dt        = 0.01   # time increment between outputs
+dt        = 1.0e-9 # time increment between outputs
 
 <parthenon/time>
-nlim       = -1    # cycle limit
-tlim       = 5.0   # time limit
-integrator = rk2   # time integration algorithm
+nlim       = -1      # cycle limit
+tlim       = 6.0e-8  # time limit
+integrator = rk1     # time integration algorithm
 
 <parthenon/mesh>
 nx1    = 4         # Number of zones in X1-direction
@@ -60,22 +62,23 @@ gas = true
 radiation = true
 
 <gas>
-gamma = 2.0  # adiabatic index
-cv = 8.0     # specific heat
+gamma = 1.4  # adiabatic index
+mu = 28.96   # mean molecular weight
 cfl = 0.8
 reconstruct = plm
 riemann = hllc
 dfloor = 1.0e-10
 siefloor = 1.0e-10
 
 <gas/opacity/absorption>
-opacity_model = thermalization
-kappa_a = 1.0
+opacity_model = constant
+kappa_a = 2.0
 
 # <gas/opacity/scattering>
 # scattering_model = constant
 # kappa_s = 10.0
 
 <jaybenne>
-dt = 0.01
+dt = 1.0e-10
 num_particles = 1000
+use_ddmc = true

src/gas/gas.cpp

@@ -103,18 +103,22 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
   if (eos_name == "ideal") {
     const Real gamma = pin->GetOrAddReal("gas", "gamma", 1.66666666667);
     auto cv = Null<Real>();
+    auto mu = Null<Real>();
     if (pin->DoesParameterExist("gas", "cv")) {
-      PARTHENON_REQUIRE(!pin->DoesParameterExist("gas", "mmw"),
-                        "Cannot specify both cv and mmw");
+      PARTHENON_REQUIRE(!pin->DoesParameterExist("gas", "mu"),
+                        "Cannot specify both cv and mu");
       cv = pin->GetReal("gas", "cv");
       PARTHENON_REQUIRE(cv > 0, "Only positive cv allowed!");
+      mu = constants.GetKBCode() / ((gamma - 1.) * constants.GetAMUCode() * cv);
     } else {
-      const Real mu = pin->GetOrAddReal("gas", "mu", 1.);
+      mu = pin->GetOrAddReal("gas", "mu", 1.);
       PARTHENON_REQUIRE(mu > 0, "Only positive mean molecular weight allowed!");
       cv = constants.GetKBCode() / ((gamma - 1.) * constants.GetAMUCode() * mu);
     }
+    params.Add("mu", mu);
+    params.Add("cv", cv);
     EOS eos_host = singularity::UnitSystem<singularity::IdealGas>(
-        singularity::IdealGas(gamma - 1., cv),
+        singularity::IdealGas(gamma - 1., cv * units.GetSpecificHeatCodeToPhysical()),
         singularity::eos_units_init::LengthTimeUnitsInit(), units.GetTimeCodeToPhysical(),
         units.GetMassCodeToPhysical(), units.GetLengthCodeToPhysical(),
         units.GetTemperatureCodeToPhysical());
@@ -139,15 +143,13 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
   } else if (opacity_model_name == "constant") {
     const Real kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "kappa_a", 0.0);
     opacity = NonCGSUnits<Gray>(Gray(kappa_a), time, mass, length, temp);
-  } else if (opacity_model_name == "shocktube_a") {
+  } else if (opacity_model_name == "powerlaw") {
     const Real coef_kappa_a =
         pin->GetOrAddReal("gas/opacity/absorption", "coef_kappa_a", 0.0);
     const Real rho_exp = pin->GetOrAddReal("gas/opacity/absorption", "rho_exp", 0.0);
     const Real temp_exp = pin->GetOrAddReal("gas/opacity/absorption", "temp_exp", 0.0);
-    opacity = ArtemisUtils::ShocktubeAOpacity(coef_kappa_a, rho_exp, temp_exp);
-  } else if (opacity_model_name == "thermalization") {
-    const Real kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "kappa_a", 0.0);
-    opacity = ArtemisUtils::ThermalizationOpacity(kappa_a);
+    opacity = NonCGSUnits<PowerLaw>(PowerLaw(coef_kappa_a, rho_exp, temp_exp), time, mass,
+                                    length, temp);
   } else {
     PARTHENON_FAIL("Opacity model not recognized!");
   }

src/pgen/shock.hpp

@@ -14,6 +14,15 @@
 #define PGEN_SHOCK_HPP_
 //! \file shock.hpp
 //! \brief
+//!
+//! This is the Mach=3 problem from Lowrie & Edwards (2008).
+//! The specific values are taken from the Fornax and Quokka code papers
+//!
+//!  mu = mH, gamma = 5/3, rho*kappa = 577 /cm
+//!  left state:         |  right state:
+//!      T = 2.18e6 K    |   T = 7.98e6 K
+//!    rho = 5.69 g/cc   | rho = 17.1 g/cc
+//!     vx = 5.19e7 cm/s |  vx = 1.73e7 cm/s
 
 // artemis headers
 #include "artemis.hpp"
@@ -32,7 +41,6 @@ namespace shock {
 struct ShockParams {
   Real rhol, vxl, tl;
   Real rhor, vxr, tr;
-  Real cv;
   Real xdisc;
 };
 
@@ -44,14 +52,13 @@ inline void InitShockParams(MeshBlock *pmb, ParameterInput *pin) {
   Params &params = artemis_pkg->AllParams();
   if (!(params.hasKey("shock_params"))) {
     ShockParams shock_params;
-    shock_params.rhol = pin->GetOrAddReal("problem", "rhol", 1.0);
-    shock_params.vxl = pin->GetOrAddReal("problem", "vxl", 2.0);
-    shock_params.tl = pin->GetOrAddReal("problem", "tl", 0.6);
-    shock_params.rhor = pin->GetOrAddReal("problem", "rhor", 2.285714);
-    shock_params.vxr = pin->GetOrAddReal("problem", "vxr", 0.875000);
-    shock_params.tr = pin->GetOrAddReal("problem", "tr", 1.246875);
+    shock_params.rhol = pin->GetOrAddReal("problem", "rhol", 5.69);
+    shock_params.vxl = pin->GetOrAddReal("problem", "vxl", 5.19e7);
+    shock_params.tl = pin->GetOrAddReal("problem", "tl", 2.18e6);
+    shock_params.rhor = pin->GetOrAddReal("problem", "rhor", 17.1);
+    shock_params.vxr = pin->GetOrAddReal("problem", "vxr", 1.73e7);
+    shock_params.tr = pin->GetOrAddReal("problem", "tr", 7.98e6);
     shock_params.xdisc = pin->GetOrAddReal("problem", "xdisc", 0.0005);
-    shock_params.cv = pin->GetOrAddReal("gas", "cv", 1.5);
     params.Add("shock_params", shock_params);
   }
 }
@@ -98,12 +105,13 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
         const bool upwind = (xi[0] <= shkp.xdisc);
         const Real rho = upwind ? shkp.rhol : shkp.rhor;
         const Real vx = upwind ? shkp.vxl : shkp.vxr;
-        const Real sie = upwind ? shkp.cv * shkp.tl : shkp.cv * shkp.tr;
+        const Real T = upwind ? shkp.tl : shkp.tr;
         v(0, gas::prim::density(0), k, j, i) = rho;
         v(0, gas::prim::velocity(0), k, j, i) = vx;
         v(0, gas::prim::velocity(1), k, j, i) = 0.0;
         v(0, gas::prim::velocity(2), k, j, i) = 0.0;
-        v(0, gas::prim::sie(0), k, j, i) = sie;
+        v(0, gas::prim::sie(0), k, j, i) =
+            eos_d.InternalEnergyFromDensityTemperature(rho, T);
       });
 
   if (do_radiation) jaybenne::InitializeRadiation(md.get(), true);
@@ -119,6 +127,7 @@ inline void ShockInnerX1(std::shared_ptr<MeshBlockData<Real>> &mbd, bool coarse)
 
   auto artemis_pkg = pmb->packages.Get("artemis");
   auto shkp = artemis_pkg->Param<ShockParams>("shock_params");
+  auto eos_d = pmb->packages.Get("gas")->Param<EOS>("eos_d");
   const auto nb = IndexRange{0, 0};
 
   static auto descriptors =
@@ -134,7 +143,8 @@ inline void ShockInnerX1(std::shared_ptr<MeshBlockData<Real>> &mbd, bool coarse)
           v(0, gas::prim::velocity(0), k, j, i) = shkp.vxl;
           v(0, gas::prim::velocity(1), k, j, i) = 0.0;
           v(0, gas::prim::velocity(2), k, j, i) = 0.0;
-          v(0, gas::prim::sie(0), k, j, i) = shkp.cv * shkp.tl;
+          v(0, gas::prim::sie(0), k, j, i) =
+              eos_d.InternalEnergyFromDensityTemperature(shkp.rhol, shkp.tl);
         });
   }
 
@@ -151,6 +161,7 @@ inline void ShockOuterX1(std::shared_ptr<MeshBlockData<Real>> &mbd, bool coarse)
 
   auto artemis_pkg = pmb->packages.Get("artemis");
   auto shkp = artemis_pkg->Param<ShockParams>("shock_params");
+  auto eos_d = pmb->packages.Get("gas")->Param<EOS>("eos_d");
   const auto nb = IndexRange{0, 0};
 
   static auto descriptors =
@@ -166,7 +177,8 @@ inline void ShockOuterX1(std::shared_ptr<MeshBlockData<Real>> &mbd, bool coarse)
           v(0, gas::prim::velocity(0), k, j, i) = shkp.vxr;
           v(0, gas::prim::velocity(1), k, j, i) = 0.0;
           v(0, gas::prim::velocity(2), k, j, i) = 0.0;
-          v(0, gas::prim::sie(0), k, j, i) = shkp.cv * shkp.tr;
+          v(0, gas::prim::sie(0), k, j, i) =
+              eos_d.InternalEnergyFromDensityTemperature(shkp.rhor, shkp.tr);
         });
   }