diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index e4f07d5c..491e05d3 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -419,8 +419,11 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto units = pkg->Param("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 diff --git a/src/hydro/srcterms/tabular_cooling.cpp b/src/hydro/srcterms/tabular_cooling.cpp index 021967bd..bbdc6c30 100644 --- a/src/hydro/srcterms/tabular_cooling.cpp +++ b/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 *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 *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{"prim"}); @@ -336,13 +336,13 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData *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 *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 *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 *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 *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 *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 *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 *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{"prim"}); @@ -635,14 +635,14 @@ Real TabularCooling::EstimateTimeStep(MeshData *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; diff --git a/src/hydro/srcterms/tabular_cooling.hpp b/src/hydro/srcterms/tabular_cooling.hpp index 33af640b..dd0efeef 100644 --- a/src/hydro/srcterms/tabular_cooling.hpp +++ b/src/hydro/srcterms/tabular_cooling.hpp @@ -106,10 +106,10 @@ class TabularCooling { parthenon::ParArray1D 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) { diff --git a/src/pgen/cloud.cpp b/src/pgen/cloud.cpp index 1e271d97..ed39b2c2 100644 --- a/src/pgen/cloud.cpp +++ b/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("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); diff --git a/src/pgen/cluster/hydrostatic_equilibrium_sphere.cpp b/src/pgen/cluster/hydrostatic_equilibrium_sphere.cpp index 95b2c025..8df0250d 100644 --- a/src/pgen/cluster/hydrostatic_equilibrium_sphere.cpp +++ b/src/pgen/cluster/hydrostatic_equilibrium_sphere.cpp @@ -41,7 +41,7 @@ HydrostaticEquilibriumSphere:: : 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"); diff --git a/src/pgen/cluster/hydrostatic_equilibrium_sphere.hpp b/src/pgen/cluster/hydrostatic_equilibrium_sphere.hpp index 3d6db7ae..b3527cfa 100644 --- a/src/pgen/cluster/hydrostatic_equilibrium_sphere.hpp +++ b/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,7 +58,7 @@ 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; } @@ -66,13 +66,13 @@ class HydrostaticEquilibriumSphere { // 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; } diff --git a/src/units.hpp b/src/units.hpp index c3666f6b..77b0043a 100644 --- a/src/units.hpp +++ b/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(); } diff --git a/tst/regression/test_suites/cluster_hse/cluster_hse.py b/tst/regression/test_suites/cluster_hse/cluster_hse.py index 127a1100..b86f858a 100644 --- a/tst/regression/test_suites/cluster_hse/cluster_hse.py +++ b/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,14 +201,14 @@ 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) ) @@ -216,7 +216,7 @@ 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) diff --git a/tst/regression/test_suites/cluster_tabular_cooling/cluster_tabular_cooling.py b/tst/regression/test_suites/cluster_tabular_cooling/cluster_tabular_cooling.py index 4ea6e1eb..7c8d78b3 100644 --- a/tst/regression/test_suites/cluster_tabular_cooling/cluster_tabular_cooling.py +++ b/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