From 8f04e18f75f4e655eff8e289759d06af87f4fbd9 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Sat, 3 Jun 2023 00:21:42 +1000 Subject: [PATCH] add derived fields to cluster pgen (#63) * add derived fields * add Mach number * use log10 radius instead * Added temperature as derived field * fix previous commit * Tabular Cooling DeDt refactor * Added cooling_time derived field * cpp-py-formatter * Update src/pgen/cluster.cpp Co-authored-by: forrestglines * Update src/pgen/cluster.cpp Co-authored-by: Philipp Grete * Update src/pgen/cluster.cpp Co-authored-by: Philipp Grete * Update src/pgen/cluster.cpp Co-authored-by: Philipp Grete * Update src/pgen/cluster.cpp Co-authored-by: Philipp Grete * Update src/pgen/cluster.cpp Co-authored-by: Philipp Grete * fix indices * Refactor cooling_table_obj, add mach_alfven * Removed extra 0, indexers * cpp-py-formatter * Explicility set edot==0 -> t_cool=NAN, and M_A * Fix calc of alfven mach in cluster output Piggyback removal of labels for acceleration fields in turb generator as the (new) default ones are effectively identical. * WIP: Rolled back tabular cooling * Compute cooling time with old DeDt * cpp-py-formatter * Reapplied cooling table refactor, fixed bug * cpp-py-formatter --------- Co-authored-by: Forrest Glines Co-authored-by: Forrest Glines Co-authored-by: par-hermes Co-authored-by: Philipp Grete Co-authored-by: Forrest Glines --- src/hydro/hydro.cpp | 2 +- src/hydro/srcterms/tabular_cooling.cpp | 112 +++++++---------- src/hydro/srcterms/tabular_cooling.hpp | 159 +++++++++++++++---------- src/main.cpp | 1 + src/pgen/cluster.cpp | 132 ++++++++++++++++++++ src/pgen/pgen.hpp | 1 + src/pgen/turbulence.cpp | 7 +- 7 files changed, 278 insertions(+), 136 deletions(-) diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index 08736ede..40cdd27b 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -503,7 +503,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { pkg->AddParam<>("enable_cooling", cooling); if (cooling == Cooling::tabular) { - TabularCooling tabular_cooling(pin); + TabularCooling tabular_cooling(pin, pkg); pkg->AddParam<>("tabular_cooling", tabular_cooling); } diff --git a/src/hydro/srcterms/tabular_cooling.cpp b/src/hydro/srcterms/tabular_cooling.cpp index f8981cb2..9164f4cc 100644 --- a/src/hydro/srcterms/tabular_cooling.cpp +++ b/src/hydro/srcterms/tabular_cooling.cpp @@ -27,17 +27,9 @@ namespace cooling { using namespace parthenon; -TabularCooling::TabularCooling(ParameterInput *pin) { - Units units(pin); - - const Real He_mass_fraction = pin->GetReal("hydro", "He_mass_fraction"); - const Real H_mass_fraction = 1.0 - He_mass_fraction; - const Real mu = 1 / (He_mass_fraction * 3. / 4. + (1 - He_mass_fraction) * 2); - - gm1_ = pin->GetReal("hydro", "gamma") - 1.0; - - mu_mh_gm1_by_k_B_ = mu * units.mh() * gm1_ / units.k_boltzmann(); - X_by_mh_ = H_mass_fraction / units.mh(); +TabularCooling::TabularCooling(ParameterInput *pin, + std::shared_ptr hydro_pkg) { + auto units = hydro_pkg->Param("units"); const std::string table_filename = pin->GetString("cooling", "table_filename"); @@ -197,7 +189,10 @@ TabularCooling::TabularCooling(ParameterInput *pin) { lambda_final_ = std::pow(10.0, log_lambdas[n_temp_ - 1]); // Setup log_lambdas_ used in Dedt() - if ((integrator_ != CoolIntegrator::townsend) || (cooling_time_cfl_ > 0.0)) { + { + // log_lambdas is used if the integrator isn't Townsend, if the cooling CFL + // is set, or if cooling time is a extra derived field. Since we don't have + // a good way to check the last condition we always initialize log_lambdas_ log_lambdas_ = ParArray1D("log_lambdas_", n_temp_); // Read log_lambdas in host_log_lambdas, changing to code units along the way @@ -260,6 +255,15 @@ TabularCooling::TabularCooling(ParameterInput *pin) { Kokkos::deep_copy(townsend_alpha_k_, host_townsend_alpha_k); Kokkos::deep_copy(townsend_Y_k_, host_townsend_Y_k); } + + // Create a lightweight object for computing cooling rates within kernels + const auto mbar_over_kb = hydro_pkg->Param("mbar_over_kb"); + const auto adiabatic_index = hydro_pkg->Param("AdiabaticIndex"); + const auto He_mass_fraction = hydro_pkg->Param("He_mass_fraction"); + + cooling_table_obj_ = CoolingTableObj(log_lambdas_, log_temp_start_, log_temp_final_, + d_log_temp_, n_temp_, mbar_over_kb, + adiabatic_index, 1.0 - He_mass_fraction, units); } void TabularCooling::SrcTerm(MeshData *md, const Real dt) const { @@ -283,16 +287,10 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData *md, const Real dt const bool mhd_enabled = hydro_pkg->Param("fluid") == Fluid::glmmhd; // Grab member variables for compiler - // Everything needed by DeDt - 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_; - const unsigned int n_temp = n_temp_; - const auto log_lambdas = log_lambdas_; - - const Real gm1 = gm1_; + const CoolingTableObj cooling_table_obj = cooling_table_obj_; + const auto gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + const auto mbar_gm1_over_kb = hydro_pkg->Param("mbar_over_kb") * gm1; + const unsigned int max_iter = max_iter_; const Real min_sub_dt = dt / max_iter; @@ -304,7 +302,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_mh_gm1_by_k_B; // specific internal en. + const Real internal_e_floor = temp_floor / mbar_gm1_over_kb; // specific internal en. // Grab some necessary variables const auto &prim_pack = md->PackVariables(std::vector{"prim"}); @@ -337,14 +335,11 @@ 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_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_mh_gm1_by_k_B, n_h2_by_rho, log_temp_start, log_temp_final, - d_log_temp, n_temp, log_lambdas, valid); + return cooling_table_obj.DeDt(e, rho, valid); }; Real sub_t = 0; // current subcycle time @@ -489,17 +484,19 @@ 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_mh_gm1_by_k_B = mu_mh_gm1_by_k_B_; - const auto X_by_mh = X_by_mh_; + + const auto units = hydro_pkg->Param("units"); + const auto gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + const auto mbar_gm1_over_kb = hydro_pkg->Param("mbar_over_kb") * gm1; + const Real X_by_mh2 = + std::pow((1 - hydro_pkg->Param("He_mass_fraction")) / units.mh(), 2); const auto lambdas = lambdas_; const auto temps = temps_; const auto alpha_k = townsend_alpha_k_; const auto Y_k = townsend_Y_k_; - const auto gm1 = gm1_; - - const auto internal_e_floor = T_floor_ / mu_mh_gm1_by_k_B; + const auto internal_e_floor = T_floor_ / mbar_gm1_over_kb; const auto temp_cool_floor = std::pow(10.0, log_temp_start_); // low end of cool table // Grab some necessary variables @@ -548,13 +545,13 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData *md, return; } - auto temp = mu_mh_gm1_by_k_B * internal_e; + auto temp = mbar_gm1_over_kb * 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_mh * X_by_mh; + const Real n_h2_by_rho = rho * X_by_mh2; // Get the index of the right temperature bin // TODO(?) this could be optimized for using a binary search @@ -571,7 +568,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_mh_gm1_by_k_B * n_h2_by_rho; + tef + lambda_final * dt / temp_final * mbar_gm1_over_kb * 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 @@ -588,8 +585,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_mh_gm1_by_k_B - : temp_cool_floor / mu_mh_gm1_by_k_B; + ? temp_new / mbar_gm1_over_kb + : temp_cool_floor / mbar_gm1_over_kb; 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). @@ -606,23 +603,18 @@ Real TabularCooling::EstimateTimeStep(MeshData *md) const { return std::numeric_limits::infinity(); } - // Everything needed by DeDt - 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_; - const unsigned int n_temp = n_temp_; - const auto log_lambdas = log_lambdas_; + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); - const Real gm1 = gm1_; + const CoolingTableObj cooling_table_obj = cooling_table_obj_; + const auto gm1 = (hydro_pkg->Param("AdiabaticIndex") - 1.0); + const auto mbar_gm1_over_kb = hydro_pkg->Param("mbar_over_kb") * gm1; // Determine the cooling floor, whichever is higher of the cooling table floor // or fluid solver floor 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_mh_gm1_by_k_B; // specific internal en. + const Real internal_e_floor = temp_floor / mbar_gm1_over_kb; // specific internal en. // Grab some necessary variables const auto &prim_pack = md->PackVariables(std::vector{"prim"}); @@ -644,15 +636,10 @@ 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_mh * X_by_mh; const Real internal_e = pres / (rho * gm1); - bool dedt_valid = true; - - const Real de_dt = - 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); + const Real de_dt = cooling_table_obj.DeDt(internal_e, rho); // Compute cooling time // If de_dt is zero (temperature is smaller than lower end of cooling table) or @@ -683,15 +670,8 @@ void TabularCooling::TestCoolingTable(ParameterInput *pin) const { // Grab member variables for compiler // Everything needed by DeDt - 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_; - const unsigned int n_temp = n_temp_; - const auto log_lambdas = log_lambdas_; - - const Real gm1 = gm1_; + const CoolingTableObj cooling_table_obj = cooling_table_obj_; + const auto gm1 = pin->GetReal("hydro", "gamma") - 1.0; // Make some device arrays to store the test data ParArray2D d_rho("d_rho", n_rho, n_pres), d_pres("d_pres", n_rho, n_pres), @@ -706,17 +686,11 @@ 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_mh * X_by_mh; - const Real internal_e = pres / (rho * gm1); d_internal_e(j, i) = internal_e; - bool dedt_valid = true; - - const Real de_dt = - 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); + const Real de_dt = cooling_table_obj.DeDt(internal_e, rho); d_de_dt(j, i) = de_dt; }); diff --git a/src/hydro/srcterms/tabular_cooling.hpp b/src/hydro/srcterms/tabular_cooling.hpp index 85ea87b7..c1ffa8b2 100644 --- a/src/hydro/srcterms/tabular_cooling.hpp +++ b/src/hydro/srcterms/tabular_cooling.hpp @@ -25,6 +25,7 @@ // AthenaPK headers #include "../../main.hpp" +#include "../../units.hpp" #ifdef MPI_PARALLEL #include @@ -86,6 +87,98 @@ struct RK45Stepper { enum class CoolIntegrator { undefined, rk12, rk45, townsend }; +class CoolingTableObj { + /************************************************************ + * Cooling Table Object, for interpolating a cooling rate out of a cooling + * table. Currently assumes evenly space log_temperatures in cooling table + * + * Lightweight object intended for inlined computation within kernels + ************************************************************/ + private: + // Log cooling rate/ne^3 + parthenon::ParArray1D log_lambdas_; + + // Spacing of cooling table + // TODO: assumes evenly spaced cooling table + parthenon::Real log_temp_start_, log_temp_final_, d_log_temp_; + unsigned int n_temp_; + + // Mean molecular mass * ( adiabatic_index -1) / boltzmann_constant + parthenon::Real mbar_gm1_over_k_B_; + + // (Hydrogen mass fraction / hydrogen atomic mass)^2 + parthenon::Real x_H_over_m_h2_; + + public: + CoolingTableObj() + : log_lambdas_(), log_temp_start_(NAN), log_temp_final_(NAN), d_log_temp_(NAN), + n_temp_(0), mbar_gm1_over_k_B_(NAN), x_H_over_m_h2_(NAN) {} + CoolingTableObj(const parthenon::ParArray1D log_lambdas, + const parthenon::Real log_temp_start, + const parthenon::Real log_temp_final, const parthenon::Real d_log_temp, + const unsigned int n_temp, const parthenon::Real mbar_over_kb, + const parthenon::Real adiabatic_index, const parthenon::Real x_H, + const Units units) + : log_lambdas_(log_lambdas), log_temp_start_(log_temp_start), + log_temp_final_(log_temp_final), d_log_temp_(d_log_temp), n_temp_(n_temp), + mbar_gm1_over_k_B_(mbar_over_kb * (adiabatic_index - 1)), + x_H_over_m_h2_(SQR(x_H / units.mh())) {} + + // Interpolate a cooling rate from the table + // from internal energy density and density + KOKKOS_INLINE_FUNCTION parthenon::Real + DeDt(const parthenon::Real &e, const parthenon::Real &rho, bool &is_valid) const { + using namespace parthenon; + + if (e < 0 || std::isnan(e)) { + is_valid = false; + return 0; + } + + const Real temp = mbar_gm1_over_k_B_ * e; + const Real log_temp = log10(temp); + Real log_lambda; + if (log_temp < log_temp_start_) { + return 0; + } else if (log_temp > log_temp_final_) { + // Above table + // Return de/dt + // TODO(forrestglines):Currently free-free cooling is used for + // temperatures above the table. This behavior could be generalized via + // templates + log_lambda = 0.5 * log_temp - 0.5 * log_temp_final_ + log_lambdas_(n_temp_ - 1); + } else { + // Inside table, interpolate assuming log spaced temperatures + + // Determine where temp is in the table + const unsigned int i_temp = + static_cast((log_temp - log_temp_start_) / d_log_temp_); + const Real log_temp_i = log_temp_start_ + d_log_temp_ * i_temp; + + // log_temp should be between log_temps[i_temp] and log_temps[i_temp+1] + PARTHENON_REQUIRE(log_temp >= log_temp_i && log_temp <= log_temp_i + d_log_temp_, + "FATAL ERROR in [CoolingTable::DeDt]: Failed to find log_temp"); + + const Real log_lambda_i = log_lambdas_(i_temp); + const Real log_lambda_ip1 = log_lambdas_(i_temp + 1); + + // Linearly interpolate lambda at log_temp + log_lambda = log_lambda_i + (log_temp - log_temp_i) * + (log_lambda_ip1 - log_lambda_i) / d_log_temp_; + } + // Return de/dt + const Real lambda = pow(10., log_lambda); + const Real de_dt = -lambda * x_H_over_m_h2_ * rho; + return de_dt; + } + + KOKKOS_INLINE_FUNCTION parthenon::Real DeDt(const parthenon::Real &e, + const parthenon::Real &rho) const { + bool is_valid = true; + return DeDt(e, rho, is_valid); + } +}; + class TabularCooling { private: // Defines uniformly spaced log temperature range of the table @@ -105,14 +198,6 @@ class TabularCooling { // Townsend cooling power law indices parthenon::ParArray1D townsend_alpha_k_; - // Some constants - // 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_; - CoolIntegrator integrator_; // Temperature floor (assumed in Kelvin and only used in cooling function) @@ -137,60 +222,11 @@ class TabularCooling { // Used for roundoff as subcycle approaches end of timestep static constexpr parthenon::Real KEpsilon_ = 1e-12; - // Interpolate a cooling rate from the table - static KOKKOS_INLINE_FUNCTION parthenon::Real - 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, - const parthenon::ParArray1D &log_lambdas, bool &valid) { - using namespace parthenon; - - if (e < 0 || std::isnan(e)) { - valid = false; - return 0; - } - - 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) { - return 0; - } else if (log_temp > log_temp_final) { - // Above table - // Return de/dt - // TODO(forrestglines):Currently free-free cooling is used for - // temperatures above the table. This behavior could be generalized via - // templates - log_lambda = 0.5 * log_temp - 0.5 * log_temp_final + log_lambdas(n_temp - 1); - } else { - // Inside table - - // Determine where temp is in the table - const unsigned int i_temp = - static_cast((log_temp - log_temp_start) / d_log_temp); - const Real log_temp_i = log_temp_start + d_log_temp * i_temp; - - // log_temp should be between log_temps[i_temp] and log_temps[i_temp+1] - if (log_temp < log_temp_i || log_temp > log_temp_i + d_log_temp) { - // FIXME exception - } - - const Real log_lambda_i = log_lambdas(i_temp); - const Real log_lambda_ip1 = log_lambdas(i_temp + 1); - - // Linearly interpolate lambda at log_temp - log_lambda = log_lambda_i + - (log_temp - log_temp_i) * (log_lambda_ip1 - log_lambda_i) / d_log_temp; - } - // Return de/dt - const Real lambda = pow(10., log_lambda); - const Real de_dt = -lambda * n_h2_by_rho; - return de_dt; - } + CoolingTableObj cooling_table_obj_; public: - TabularCooling(parthenon::ParameterInput *pin); + TabularCooling(parthenon::ParameterInput *pin, + std::shared_ptr hydro_pkg); void SrcTerm(parthenon::MeshData *md, const parthenon::Real dt) const; @@ -206,6 +242,9 @@ class TabularCooling { parthenon::Real EstimateTimeStep(parthenon::MeshData *md) const; + // Get a lightweight object for computing cooling rate from the cooling table + const CoolingTableObj GetCoolingTableObj() const { return cooling_table_obj_; } + void TestCoolingTable(parthenon::ParameterInput *pin) const; }; diff --git a/src/main.cpp b/src/main.cpp index 8cb090a3..1d9ef467 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -86,6 +86,7 @@ int main(int argc, char *argv[]) { Hydro::ProblemSourceFirstOrder = rand_blast::RandomBlasts; } else if (problem == "cluster") { pman.app_input->ProblemGenerator = cluster::ProblemGenerator; + pman.app_input->MeshBlockUserWorkBeforeOutput = cluster::UserWorkBeforeOutput; Hydro::ProblemInitPackageData = cluster::ProblemInitPackageData; Hydro::ProblemSourceUnsplit = cluster::ClusterSrcTerm; Hydro::ProblemEstimateTimestep = cluster::ClusterEstimateTimestep; diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp index 757f203a..69385ea0 100644 --- a/src/pgen/cluster.cpp +++ b/src/pgen/cluster.cpp @@ -33,6 +33,7 @@ // AthenaPK headers #include "../hydro/hydro.hpp" #include "../hydro/srcterms/gravitational_field.hpp" +#include "../hydro/srcterms/tabular_cooling.hpp" #include "../main.hpp" // Cluster headers @@ -214,6 +215,35 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd ************************************************************/ SNIAFeedback snia_feedback(pin, hydro_pkg); + + /************************************************************ + * Add derived fields + * NOTE: these must be filled in UserWorkBeforeOutput + ************************************************************/ + + auto m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector({1})); + + // log10 of cell-centered radius + hydro_pkg->AddField("log10_cell_radius", m); + // entropy + hydro_pkg->AddField("entropy", m); + // sonic Mach number v/c_s + hydro_pkg->AddField("mach_sonic", m); + // temperature + hydro_pkg->AddField("temperature", m); + + if (hydro_pkg->Param("enable_cooling") == Cooling::tabular) { + // cooling time + hydro_pkg->AddField("cooling_time", m); + } + + if (hydro_pkg->Param("fluid") == Fluid::glmmhd) { + // alfven Mach number v_A/c_s + hydro_pkg->AddField("mach_alfven", m); + + // plasma beta + hydro_pkg->AddField("plasma_beta", m); + } } //======================================================================================== @@ -405,4 +435,106 @@ void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin) { } // END if(hydro_pkg->Param("fluid") == Fluid::glmmhd) } +void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { + // get hydro + auto pkg = pmb->packages.Get("Hydro"); + const Real gam = pin->GetReal("hydro", "gamma"); + const Real gm1 = (gam - 1.0); + + // get prim vars + auto &data = pmb->meshblock_data.Get(); + auto const &prim = data->Get("prim").data; + + // get derived fields + auto &log10_radius = data->Get("log10_cell_radius").data; + auto &entropy = data->Get("entropy").data; + auto &mach_sonic = data->Get("mach_sonic").data; + auto &temperature = data->Get("temperature").data; + + // for computing temperature from primitives + auto mbar_over_kb = pkg->Param("mbar_over_kb"); + + // fill derived vars (*including ghost cells*) + auto &coords = pmb->coords; + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); + + pmb->par_for( + "Cluster::UserWorkBeforeOutput", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + // get gas properties + const Real rho = prim(IDN, k, j, i); + const Real v1 = prim(IV1, k, j, i); + const Real v2 = prim(IV2, k, j, i); + const Real v3 = prim(IV3, k, j, i); + const Real P = prim(IPR, k, j, i); + + // compute radius + const Real x = coords.Xc<1>(i); + const Real y = coords.Xc<2>(j); + const Real z = coords.Xc<3>(k); + const Real r2 = SQR(x) + SQR(y) + SQR(z); + log10_radius(k, j, i) = 0.5 * std::log10(r2); + + // compute entropy + const Real K = P / std::pow(rho, gam); + entropy(k, j, i) = K; + + const Real v_mag = std::sqrt(SQR(v1) + SQR(v2) + SQR(v3)); + const Real c_s = std::sqrt(gam * P / rho); // ideal gas EOS + const Real M_s = v_mag / c_s; + mach_sonic(k, j, i) = M_s; + + // compute temperature + temperature(k, j, i) = mbar_over_kb * P / rho; + }); + if (pkg->Param("enable_cooling") == Cooling::tabular) { + auto &cooling_time = data->Get("cooling_time").data; + + // get cooling function + const cooling::TabularCooling &tabular_cooling = + pkg->Param("tabular_cooling"); + const auto cooling_table_obj = tabular_cooling.GetCoolingTableObj(); + + pmb->par_for( + "Cluster::UserWorkBeforeOutput::CoolingTime", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + // get gas properties + const Real rho = prim(IDN, k, j, i); + const Real P = prim(IPR, k, j, i); + + // compute cooling time + const Real eint = P / (rho * gm1); + const Real edot = cooling_table_obj.DeDt(eint, rho); + cooling_time(k, j, i) = (edot != 0) ? -eint / edot : NAN; + }); + } + + if (pkg->Param("fluid") == Fluid::glmmhd) { + auto &plasma_beta = data->Get("plasma_beta").data; + auto &mach_alfven = data->Get("mach_alfven").data; + + pmb->par_for( + "Cluster::UserWorkBeforeOutput::MHD", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + // get gas properties + const Real rho = prim(IDN, k, j, i); + const Real P = prim(IPR, k, j, i); + const Real Bx = prim(IB1, k, j, i); + const Real By = prim(IB2, k, j, i); + const Real Bz = prim(IB3, k, j, i); + const Real B2 = (SQR(Bx) + SQR(By) + SQR(Bz)); + + // compute Alfven mach number + const Real v_A = std::sqrt(B2 / rho); + const Real c_s = std::sqrt(gam * P / rho); // ideal gas EOS + mach_alfven(k, j, i) = mach_sonic(k, j, i) * c_s / v_A; + + // compute plasma beta + plasma_beta(k, j, i) = (B2 != 0) ? P / (0.5 * B2) : NAN; + }); + } +} + } // namespace cluster diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index d8eba455..84e36b00 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -99,6 +99,7 @@ using namespace parthenon::driver::prelude; void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg); void InitUserMeshData(ParameterInput *pin); void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin); +void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin); void ClusterSrcTerm(MeshData *md, const parthenon::SimTime &tm, const Real beta_dt); parthenon::Real ClusterEstimateTimestep(MeshData *md); } // namespace cluster diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index d74fac64..91bb935a 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -120,18 +120,13 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg pkg->UpdateParam(parthenon::hist_param_key, hst_vars); // Step 2. Add appropriate fields required by this pgen - std::vector acc_labels(3); - acc_labels[0] = "Acceleration1"; - acc_labels[1] = "Acceleration2"; - acc_labels[2] = "Acceleration3"; - // Using OneCopy here to save memory. We typically don't need to update/evolve the // acceleration field for various stages in a cycle as the "model" error of the // turbulence driver is larger than the numerical one any way. This may need to be // changed if an "as close as possible" comparison between methods/codes is the goal and // not turbulence from a physical point of view. Metadata m({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, - std::vector({3}), acc_labels); + std::vector({3})); pkg->AddField("acc", m); auto num_modes =