From 7f3bc17e54529a89f1827d43ca7d41422baff876 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 20 Mar 2024 19:34:36 -0400 Subject: [PATCH] few more details chainged to address suggestions --- src/analysis/history.cpp | 7 ++++--- src/analysis/history.hpp | 6 ++++-- src/progenitor/progenitordata.cpp | 16 +++++++++------- 3 files changed, 17 insertions(+), 12 deletions(-) diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index 4d35fcff..8f578dd0 100644 --- a/src/analysis/history.cpp +++ b/src/analysis/history.cpp @@ -309,6 +309,7 @@ Real CalculateMdot(MeshData *md, Real rc, bool gain) { const bool do_gain = pdo_gain_reducer->val; auto progenitor = pmb->packages.Get("progenitor").get(); const Real inside_pns_threshold = progenitor->Param("inside_pns_threshold"); + const Real net_heat_threshold = progenitor->Param("net_heat_threshold"); parthenon::par_reduce( parthenon::LoopPatternMDRange(), "Calculates mass accretion rate (SN)", @@ -322,9 +323,9 @@ Real CalculateMdot(MeshData *md, Real rc, bool gain) { Real r = std::sqrt(C[1] * C[1] + C[2] * C[2] + C[3] * C[3]); if (r <= rc) { if (gain) { - bool is_netheat = (v(b, internal_variables::GcovHeat(), k, j, i) - - v(b, internal_variables::GcovCool(), k, j, i) > - 1.e-8); // checks that in the gain region + bool is_netheat = ((v(b, internal_variables::GcovHeat(), k, j, i) - + v(b, internal_variables::GcovCool(), k, j, i)) > + net_heat_threshold); // checks that in the gain region bool is_inside_pns = (r < inside_pns_threshold); // checks that inside PNS if (do_gain && (is_inside_pns || is_netheat)) { diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index 98901a8d..55dfdae3 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -110,6 +110,7 @@ Real ReduceInGain(MeshData *md, bool is_conserved, int idx = 0) { auto geom = Geometry::GetCoordinateSystem(md); auto progenitor = pmb->packages.Get("progenitor").get(); const Real outside_pns_threshold = progenitor->Param("outside_pns_threshold"); + const Real net_heat_threshold = progenitor->Param("net_heat_threshold"); parthenon::par_reduce( parthenon::LoopPatternMDRange(), @@ -117,8 +118,9 @@ Real ReduceInGain(MeshData *md, bool is_conserved, int idx = 0) { kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lresult) { Real gdet = geom.DetGamma(CellLocation::Cent, 0, k, j, i); - bool is_netheat = (v(b, iv::GcovHeat(), k, j, i) - v(b, iv::GcovCool(), k, j, i) > - 1.e-8); // checks that in the gain region + bool is_netheat = + ((v(b, iv::GcovHeat(), k, j, i) - v(b, iv::GcovCool(), k, j, i)) > + net_heat_threshold); // checks that in the gain region bool is_outside_pns = (v(b, fluid_prim::entropy(), k, j, i) > outside_pns_threshold); // checks that outside PNS const auto &coords = v.GetCoordinates(b); diff --git a/src/progenitor/progenitordata.cpp b/src/progenitor/progenitordata.cpp index 79d91f79..535a1424 100644 --- a/src/progenitor/progenitordata.cpp +++ b/src/progenitor/progenitordata.cpp @@ -90,11 +90,12 @@ std::shared_ptr Initialize(ParameterInput *pin) { 2.42e-5); // corresponds to entropy > 3 kb/baryon Real inside_pns_threshold = pin->GetOrAddReal("progenitor", "inside_pns_threshold", 0.008); // corresponds to r < 80 km + Real net_heat_threshold = pin->GetOrAddReal("progenitor", "inside_pns_threshold", + 1e-8); // corresponds to r < 80 km UnitConversions units(pin); Real LengthCGSToCode = units.GetLengthCGSToCode(); - Real rc = 400e5 * LengthCGSToCode; auto mdot_radii = pin->GetOrAddVector("progenitor", "mdot_radii", - std::vector{rc}); // default 400km + std::vector{400}); // default 400km // Add Params params.Add("mass_density", mass_density); @@ -123,6 +124,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { params.Add("outside_pns_threshold", outside_pns_threshold); params.Add("inside_pns_threshold", inside_pns_threshold); + params.Add("net_heat_threshold", net_heat_threshold); params.Add("mdot_radii", mdot_radii); // Reductions @@ -140,12 +142,12 @@ std::shared_ptr Initialize(ParameterInput *pin) { ReduceInGain(md, 0, 0); }; for (auto rc : mdot_radii) { - auto Mdot = [rc, LengthCGSToCode](MeshData *md) { - return History::CalculateMdot(md, rc, 0); + auto rc_code = rc * 1e5 * LengthCGSToCode; + auto Mdot = [rc_code](MeshData *md) { + return History::CalculateMdot(md, rc_code, 0); }; - hst_vars.emplace_back(HistoryOutputVar( - HstSum, Mdot, - "Mdot at r = " + std::to_string(int(rc / LengthCGSToCode / 1.e5)) + "km")); + hst_vars.emplace_back( + HistoryOutputVar(HstSum, Mdot, "Mdot at r = " + std::to_string(int(rc)) + "km")); } Real x1max = pin->GetReal("parthenon/mesh", "x1max");