Merge branch 'lhllc' of github.com:parthenon-hpc-lab/athenapk into lhllc

src/hydro/hydro.cpp

@@ -419,8 +419,11 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
       auto units = pkg->Param<Units>("units");
       const auto He_mass_fraction = pin->GetReal("hydro", "He_mass_fraction");
       const auto mu = 1 / (He_mass_fraction * 3. / 4. + (1 - He_mass_fraction) * 2);
-      pkg->AddParam<>("mbar_over_kb",
-                      mu * units.atomic_mass_unit() / units.k_boltzmann());
+      pkg->AddParam<>("mu", mu);
+      pkg->AddParam<>("He_mass_fraction", He_mass_fraction);
+      // Following convention in the astro community, we're using mh as unit for the mean
+      // molecular weight
+      pkg->AddParam<>("mbar_over_kb", mu * units.mh() / units.k_boltzmann());
     }
 
     // By default disable floors by setting a negative value

src/hydro/srcterms/tabular_cooling.cpp

@@ -35,8 +35,8 @@ TabularCooling::TabularCooling(ParameterInput *pin) {
 
   gm1_ = pin->GetReal("hydro", "gamma") - 1.0;
 
-  mu_m_u_gm1_by_k_B_ = mu * units.atomic_mass_unit() * gm1_ / units.k_boltzmann();
-  X_by_m_u_ = H_mass_fraction / units.atomic_mass_unit();
+  mu_mh_gm1_by_k_B_ = mu * units.mh() * gm1_ / units.k_boltzmann();
+  X_by_mh_ = H_mass_fraction / units.mh();
 
   const std::string table_filename = pin->GetString("cooling", "table_filename");
 
@@ -283,8 +283,8 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData<Real> *md, const Real dt
   // Grab member variables for compiler
 
   // Everything needed by DeDt
-  const Real mu_m_u_gm1_by_k_B = mu_m_u_gm1_by_k_B_;
-  const Real X_by_m_u = X_by_m_u_;
+  const Real mu_mh_gm1_by_k_B = mu_mh_gm1_by_k_B_;
+  const Real X_by_mh = X_by_mh_;
   const Real log_temp_start = log_temp_start_;
   const Real log_temp_final = log_temp_final_;
   const Real d_log_temp = d_log_temp_;
@@ -303,7 +303,7 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData<Real> *md, const Real dt
   const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table
   const Real temp_floor = (T_floor_ > temp_cool_floor) ? T_floor_ : temp_cool_floor;
 
-  const Real internal_e_floor = temp_floor / mu_m_u_gm1_by_k_B; // specific internal en.
+  const Real internal_e_floor = temp_floor / mu_mh_gm1_by_k_B; // specific internal en.
 
   // Grab some necessary variables
   const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
@@ -336,13 +336,13 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData<Real> *md, const Real dt
         internal_e /= rho;
         const Real internal_e_initial = internal_e;
 
-        const Real n_h2_by_rho = rho * X_by_m_u * X_by_m_u;
+        const Real n_h2_by_rho = rho * X_by_mh * X_by_mh;
 
         bool dedt_valid = true;
 
         // Wrap DeDt into a functor for the RKStepper
         auto DeDt_wrapper = [&](const Real t, const Real e, bool &valid) {
-          return DeDt(e, mu_m_u_gm1_by_k_B, n_h2_by_rho, log_temp_start, log_temp_final,
+          return DeDt(e, mu_mh_gm1_by_k_B, n_h2_by_rho, log_temp_start, log_temp_final,
                       d_log_temp, n_temp, log_lambdas, valid);
         };
 
@@ -483,8 +483,8 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,
 
   // Grab member variables for compiler
   const auto dt = dt_; // HACK capturing parameters still broken with Cuda 11.6 ...
-  const auto mu_m_u_gm1_by_k_B = mu_m_u_gm1_by_k_B_;
-  const auto X_by_m_u = X_by_m_u_;
+  const auto mu_mh_gm1_by_k_B = mu_mh_gm1_by_k_B_;
+  const auto X_by_mh = X_by_mh_;
 
   const auto lambdas = lambdas_;
   const auto temps = temps_;
@@ -493,7 +493,7 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,
 
   const auto gm1 = gm1_;
 
-  const auto internal_e_floor = T_floor_ / mu_m_u_gm1_by_k_B;
+  const auto internal_e_floor = T_floor_ / mu_mh_gm1_by_k_B;
   const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table
 
   // Grab some necessary variables
@@ -542,13 +542,13 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,
           return;
         }
 
-        auto temp = mu_m_u_gm1_by_k_B * internal_e;
+        auto temp = mu_mh_gm1_by_k_B * internal_e;
         // Temperature is above floor (see conditional above) but below cooling table:
         // -> no cooling
         if (temp < temp_cool_floor) {
           return;
         }
-        const Real n_h2_by_rho = rho * X_by_m_u * X_by_m_u;
+        const Real n_h2_by_rho = rho * X_by_mh * X_by_mh;
 
         // Get the index of the right temperature bin
         // TODO(?) this could be optimized for using a binary search
@@ -565,7 +565,7 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,
 
         // Compute the adjusted TEF for new timestep (Eqn. 26) (term in brackets)
         const auto tef_adj =
-            tef + lambda_final * dt / temp_final * mu_m_u_gm1_by_k_B * n_h2_by_rho;
+            tef + lambda_final * dt / temp_final * mu_mh_gm1_by_k_B * n_h2_by_rho;
 
         // TEF is a strictly decreasing function and new_tef > tef
         // Check if the new TEF falls into a lower bin, i.e., find the right bin for A7
@@ -582,8 +582,8 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,
                      1.0 / (1.0 - alpha_k(idx)));
         // Set new temp (at the lowest to the lower end of the cooling table)
         const auto internal_e_new = temp_new > temp_cool_floor
-                                        ? temp_new / mu_m_u_gm1_by_k_B
-                                        : temp_cool_floor / mu_m_u_gm1_by_k_B;
+                                        ? temp_new / mu_mh_gm1_by_k_B
+                                        : temp_cool_floor / mu_mh_gm1_by_k_B;
         cons(IEN, k, j, i) += rho * (internal_e_new - internal_e);
         // Latter technically not required if no other tasks follows before
         // ConservedToPrim conversion, but keeping it for now (better safe than sorry).
@@ -598,8 +598,8 @@ Real TabularCooling::EstimateTimeStep(MeshData<Real> *md) const {
   // Grab member variables for compiler
 
   // Everything needed by DeDt
-  const Real mu_m_u_gm1_by_k_B = mu_m_u_gm1_by_k_B_;
-  const Real X_by_m_u = X_by_m_u_;
+  const Real mu_mh_gm1_by_k_B = mu_mh_gm1_by_k_B_;
+  const Real X_by_mh = X_by_mh_;
   const Real log_temp_start = log_temp_start_;
   const Real log_temp_final = log_temp_final_;
   const Real d_log_temp = d_log_temp_;
@@ -613,7 +613,7 @@ Real TabularCooling::EstimateTimeStep(MeshData<Real> *md) const {
   const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table
   const Real temp_floor = (T_floor_ > temp_cool_floor) ? T_floor_ : temp_cool_floor;
 
-  const Real internal_e_floor = temp_floor / mu_m_u_gm1_by_k_B; // specific internal en.
+  const Real internal_e_floor = temp_floor / mu_mh_gm1_by_k_B; // specific internal en.
 
   // Grab some necessary variables
   const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
@@ -635,14 +635,14 @@ Real TabularCooling::EstimateTimeStep(MeshData<Real> *md) const {
 
         const Real rho = prim(IDN, k, j, i);
         const Real pres = prim(IPR, k, j, i);
-        const Real n_h2_by_rho = rho * X_by_m_u * X_by_m_u;
+        const Real n_h2_by_rho = rho * X_by_mh * X_by_mh;
 
         const Real internal_e = pres / (rho * gm1);
 
         bool dedt_valid = true;
 
         const Real de_dt =
-            DeDt(internal_e, mu_m_u_gm1_by_k_B, n_h2_by_rho, log_temp_start,
+            DeDt(internal_e, mu_mh_gm1_by_k_B, n_h2_by_rho, log_temp_start,
                  log_temp_final, d_log_temp, n_temp, log_lambdas, dedt_valid);
 
         // Compute cooling time
@@ -674,8 +674,8 @@ void TabularCooling::TestCoolingTable(ParameterInput *pin) const {
   // Grab member variables for compiler
 
   // Everything needed by DeDt
-  const auto mu_m_u_gm1_by_k_B = mu_m_u_gm1_by_k_B_;
-  const auto X_by_m_u = X_by_m_u_;
+  const auto mu_mh_gm1_by_k_B = mu_mh_gm1_by_k_B_;
+  const auto X_by_mh = X_by_mh_;
   const auto log_temp_start = log_temp_start_;
   const auto log_temp_final = log_temp_final_;
   const auto d_log_temp = d_log_temp_;
@@ -697,7 +697,7 @@ void TabularCooling::TestCoolingTable(ParameterInput *pin) const {
         d_rho(j, i) = rho;
         d_pres(j, i) = pres;
 
-        const Real n_h2_by_rho = rho * X_by_m_u * X_by_m_u;
+        const Real n_h2_by_rho = rho * X_by_mh * X_by_mh;
 
         const Real internal_e = pres / (rho * gm1);
 
@@ -706,7 +706,7 @@ void TabularCooling::TestCoolingTable(ParameterInput *pin) const {
         bool dedt_valid = true;
 
         const Real de_dt =
-            DeDt(internal_e, mu_m_u_gm1_by_k_B, n_h2_by_rho, log_temp_start,
+            DeDt(internal_e, mu_mh_gm1_by_k_B, n_h2_by_rho, log_temp_start,
                  log_temp_final, d_log_temp, n_temp, log_lambdas, dedt_valid);
 
         d_de_dt(j, i) = de_dt;

src/hydro/srcterms/tabular_cooling.hpp

@@ -106,10 +106,10 @@ class TabularCooling {
   parthenon::ParArray1D<parthenon::Real> townsend_alpha_k_;
 
   // Some constants
-  // mean_molecular_mass*atomic_mass_unit*(adiabatic_index-1)/k_B
-  parthenon::Real mu_m_u_gm1_by_k_B_;
-  // H_mass_fraction/atomic_mass_unit
-  parthenon::Real X_by_m_u_;
+  // mean_molecular_mass*mh*(adiabatic_index-1)/k_B
+  parthenon::Real mu_mh_gm1_by_k_B_;
+  // H_mass_fraction/mh
+  parthenon::Real X_by_mh_;
   // adiabatic_index -1
   parthenon::Real gm1_;
 
@@ -132,7 +132,7 @@ class TabularCooling {
 
   // Interpolate a cooling rate from the table
   static KOKKOS_INLINE_FUNCTION parthenon::Real
-  DeDt(const parthenon::Real &e, const parthenon::Real &mu_m_u_gm1_by_k_B,
+  DeDt(const parthenon::Real &e, const parthenon::Real &mu_mh_gm1_by_k_B,
        const parthenon::Real &n_h2_by_rho, const parthenon::Real &log_temp_start,
        const parthenon::Real &log_temp_final, const parthenon::Real &d_log_temp,
        const unsigned int n_temp,
@@ -144,7 +144,7 @@ class TabularCooling {
       return 0;
     }
 
-    const Real temp = mu_m_u_gm1_by_k_B * e;
+    const Real temp = mu_mh_gm1_by_k_B * e;
     const Real log_temp = log10(temp);
     Real log_lambda;
     if (log_temp < log_temp_start) {

src/pgen/cloud.cpp

@@ -45,12 +45,8 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin) {
 
   auto gamma = pin->GetReal("hydro", "gamma");
   auto gm1 = (gamma - 1.0);
-  // TODO(pgrete): reuse code from tabular_cooling
-  const auto He_mass_fraction = pin->GetReal("hydro", "He_mass_fraction");
-  const auto H_mass_fraction = 1.0 - He_mass_fraction;
-  const auto mu = 1 / (He_mass_fraction * 3. / 4. + (1 - He_mass_fraction) * 2);
-  const auto mu_m_u_gm1_by_k_B_ =
-      mu * units.atomic_mass_unit() * gm1 / units.k_boltzmann();
+  const auto &pkg = mesh->packages.Get("Hydro");
+  const auto mbar_over_kb = pkg->Param<Real>("mbar_over_kb");
 
   r_cloud = pin->GetReal("problem/cloud", "r0_cgs") / units.code_length_cgs();
   rho_cloud = pin->GetReal("problem/cloud", "rho_cloud_cgs") / units.code_density_cgs();
@@ -59,15 +55,15 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin) {
   auto v_wind = pin->GetReal("problem/cloud", "v_wind_cgs") /
                 (units.code_length_cgs() / units.code_time_cgs());
 
-  // mu_m_u_gm1_by_k_B is already in code units
-  rhoe_wind = T_wind * rho_wind / mu_m_u_gm1_by_k_B_;
+  // mu_mh_gm1_by_k_B is already in code units
+  rhoe_wind = T_wind * rho_wind / mbar_over_kb / gm1;
   const auto c_s_wind = std::sqrt(gamma * gm1 * rhoe_wind / rho_wind);
   const auto chi_0 = rho_cloud / rho_wind;               // cloud to wind density ratio
   const auto t_cc = r_cloud * std::sqrt(chi_0) / v_wind; // cloud crushting time (code)
   const auto pressure =
       gm1 * rhoe_wind; // one value for entire domain given initial pressure equil.
 
-  const auto T_cloud = pressure / gm1 / rho_cloud * mu_m_u_gm1_by_k_B_;
+  const auto T_cloud = pressure / rho_cloud * mbar_over_kb;
 
   auto plasma_beta = pin->GetOrAddReal("problem/cloud", "plasma_beta", -1.0);
 

src/pgen/cluster/hydrostatic_equilibrium_sphere.cpp

@@ -41,7 +41,7 @@ HydrostaticEquilibriumSphere<GravitationalField, EntropyProfile>::
     : gravitational_field_(gravitational_field), entropy_profile_(entropy_profile) {
   Units units(pin);
 
-  atomic_mass_unit_ = units.atomic_mass_unit();
+  mh_ = units.mh();
   k_boltzmann_ = units.k_boltzmann();
 
   const Real He_mass_fraction = pin->GetReal("hydro", "He_mass_fraction");

src/pgen/cluster/hydrostatic_equilibrium_sphere.hpp

@@ -36,7 +36,7 @@ class HydrostaticEquilibriumSphere {
   const EntropyProfile entropy_profile_;
 
   // Physical constants
-  parthenon::Real atomic_mass_unit_, k_boltzmann_;
+  parthenon::Real mh_, k_boltzmann_;
 
   // Density to fix baryons at a radius (change to temperature?)
   parthenon::Real R_fix_, rho_fix_;
@@ -58,21 +58,21 @@ class HydrostaticEquilibriumSphere {
   // of entropy
   parthenon::Real P_from_rho_K(const parthenon::Real rho, const parthenon::Real K) const {
     const parthenon::Real P =
-        K * pow(mu_ / mu_e_, 2. / 3.) * pow(rho / (mu_ * atomic_mass_unit_), 5. / 3.);
+        K * pow(mu_ / mu_e_, 2. / 3.) * pow(rho / (mu_ * mh_), 5. / 3.);
     return P;
   }
 
   // Get density from pressure and entropy, using ideal gas law and definition
   // of entropy
   parthenon::Real rho_from_P_K(const parthenon::Real P, const parthenon::Real K) const {
     const parthenon::Real rho =
-        pow(P / K, 3. / 5.) * mu_ * atomic_mass_unit_ / pow(mu_ / mu_e_, 2. / 5);
+        pow(P / K, 3. / 5.) * mu_ * mh_ / pow(mu_ / mu_e_, 2. / 5);
     return rho;
   }
 
   // Get total number density from density
   parthenon::Real n_from_rho(const parthenon::Real rho) const {
-    const parthenon::Real n = rho / (mu_ * atomic_mass_unit_);
+    const parthenon::Real n = rho / (mu_ * mh_);
     return n;
   }
 

src/units.hpp

@@ -1,6 +1,6 @@
 
 // Athena-Parthenon - a performance portable block structured AMR MHD code
-// Copyright (c) 2020, Athena Parthenon Collaboration. All rights reserved.
+// Copyright (c) 2020-2023, Athena Parthenon Collaboration. All rights reserved.
 // Licensed under the 3-Clause License (the "LICENSE")
 
 #ifndef UNITS_HPP_
@@ -34,6 +34,8 @@ class Units {
   static constexpr parthenon::Real erg_cgs = 1;                            // erg
   static constexpr parthenon::Real gauss_cgs = 1;                          // gauss
   static constexpr parthenon::Real microgauss_cgs = 1e-6;                  // gauss
+  static constexpr parthenon::Real mh_cgs =
+      1.007947 * atomic_mass_unit_cgs; // g (matching the definition used in yt)
 
   // PHYSICAL CONSTANTS
   static constexpr parthenon::Real k_boltzmann_cgs = 1.3806488e-16; // erg/k
@@ -122,6 +124,7 @@ class Units {
   parthenon::Real atomic_mass_unit() const {
     return atomic_mass_unit_cgs / code_mass_cgs();
   }
+  parthenon::Real mh() const { return mh_cgs / code_mass_cgs(); }
   parthenon::Real erg() const { return erg_cgs / code_energy_cgs(); }
   parthenon::Real gauss() const { return gauss_cgs / code_magnetic_cgs(); }
   parthenon::Real microgauss() const { return microgauss_cgs / code_magnetic_cgs(); }

tst/regression/test_suites/cluster_hse/cluster_hse.py

@@ -46,7 +46,7 @@ def __init__(self):
         # Setup constants
         self.k_b = unyt.kb_cgs
         self.G = unyt.G_cgs
-        self.m_u = unyt.amu
+        self.mh = unyt.mh
 
         self.adiabatic_index = 5.0 / 3.0
 
@@ -201,22 +201,22 @@ def P_from_rho_K(rho, K):
             return (
                 K
                 * (self.mu / self.mu_e) ** (2.0 / 3.0)
-                * (rho / (self.mu * self.m_u)) ** (5.0 / 3.0)
+                * (rho / (self.mu * self.mh)) ** (5.0 / 3.0)
             )
 
         def rho_from_P_K(P, K):
             return (
                 (P / K) ** (3.0 / 5.0)
                 * self.mu
-                * self.m_u
+                * self.mh
                 / (self.mu / self.mu_e) ** (2.0 / 5.0)
             )
 
         def K_from_r(r):
             return self.K_0 + self.K_100 * (r / self.R_K) ** self.alpha_K
 
         def n_from_rho(rho):
-            return rho / (self.mu * self.m_u)
+            return rho / (self.mu * self.mh)
 
         def ne_from_rho(rho):
             return self.mu / self.mu_e * n_from_rho(rho)

tst/regression/test_suites/cluster_tabular_cooling/cluster_tabular_cooling.py

@@ -215,10 +215,10 @@ def Analyse(self, parameters):
         initial_internal_e = (
             self.uniform_gas_pres / (self.uniform_gas_rho * (self.adiabatic_index - 1))
         ).in_units("erg/g")
-        n_h = (self.uniform_gas_rho * H_mass_fraction / unyt.amu).in_units("1/cm**3")
+        n_h = (self.uniform_gas_rho * H_mass_fraction / unyt.mh).in_units("1/cm**3")
 
         # Compute temperature
-        initial_temp = initial_internal_e * (mu * unyt.amu * gm1 / unyt.kb)
+        initial_temp = initial_internal_e * (mu * unyt.mh * gm1 / unyt.kb)
 
         # Unit helpers
         Ta = unyt.unyt_quantity(1, "K")
@@ -230,7 +230,7 @@ def Analyse(self, parameters):
         )
         cooling_b = self.log_lambda1 - self.log_temp1 * cooling_m
         B = unyt.unyt_quantity(10**cooling_b, "erg*cm**3/s")
-        X = mu * unyt.amu * (self.adiabatic_index - 1) / (unyt.kb * Ta)
+        X = mu * unyt.mh * (self.adiabatic_index - 1) / (unyt.kb * Ta)
         Y = B * n_h**2 / self.uniform_gas_rho
 
         # Integrate for the final internal energy