From 4f761b2e8c47b90bbd13f15a1f72605c1da0f689 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Tue, 19 Mar 2024 18:23:08 -0400 Subject: [PATCH 1/7] calculating Mdot for SN post-processing --- src/analysis/analysis.cpp | 7 +++ src/analysis/history.cpp | 97 +++++++++++++++++++++++++++++++ src/analysis/history.hpp | 22 +++++-- src/progenitor/progenitordata.cpp | 17 +++++- 4 files changed, 135 insertions(+), 8 deletions(-) diff --git a/src/analysis/analysis.cpp b/src/analysis/analysis.cpp index 54837fb3..e7adde50 100644 --- a/src/analysis/analysis.cpp +++ b/src/analysis/analysis.cpp @@ -6,10 +6,17 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto analysis_pkg = std::make_shared("analysis"); Params ¶ms = analysis_pkg->AllParams(); Real sigma = pin->GetOrAddReal("analysis", "sigma", 2e-3); + Real outside_pns_threshold = + pin->GetOrAddReal("analysis", "outside_pns_threshold", + 2.42e-5); // corresponds to entropy > 3 kb/baryon + Real inside_pns_threshold = pin->GetOrAddReal("analysis", "inside_pns_threshold", + 0.008); // corresponds to r < 80 km if (sigma < 0) { PARTHENON_THROW("sigma must be greater than 0"); } params.Add("sigma", sigma); + params.Add("outside_pns_threshold", outside_pns_threshold); + params.Add("inside_pns_threshold", inside_pns_threshold); return analysis_pkg; } // initialize diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index d88e3d65..aa0f10ac 100644 --- a/src/analysis/history.cpp +++ b/src/analysis/history.cpp @@ -18,6 +18,7 @@ #include "history_utils.hpp" #include "phoebus_utils/relativity_utils.hpp" #include +#include namespace History { @@ -279,4 +280,100 @@ void ReduceLocalizationFunction(MeshData *md) { } // exp +using namespace parthenon::package::prelude; +template +KOKKOS_INLINE_FUNCTION Real ComputeDivInPillbox(int ndim, int b, int k, int j, int i, + const parthenon::Coordinates_t &coords, + const F &f) { + Real div_mass_flux_integral; + div_mass_flux_integral = + -(f(b, 1, k, j, i + 1) - f(b, 1, k, j, i)) * coords.FaceArea(k, j, i); + + if (ndim >= 2) { + div_mass_flux_integral += + -(f(b, 2, k, j + 1, i) - f(b, 2, k, j, i)) * coords.FaceArea(k, j, i); + } + if (ndim >= 3) { + div_mass_flux_integral += + -(f(b, 3, k + 1, j, i) - f(b, 3, k, j, i)) * coords.FaceArea(k, j, i); + } + return div_mass_flux_integral; +} + +// This function calculates mass accretion rate which I defined as +// Mdot=Int(dV*d/dx^i{detg*rho*U^i}) where detg is the determinant of four metric, U is +// four-velocity, and dV=d^3x + +Real CalculateMdot(MeshData *md, Real rc, bool gain) { + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + namespace c = fluid_cons; + using parthenon::MakePackDescriptor; + auto *pmb = md->GetParentPointer(); + Mesh *pmesh = md->GetMeshPointer(); + auto &resolved_pkgs = pmesh->resolved_packages; + const int ndim = pmesh->ndim; + using parthenon::PDOpt; + static auto desc = MakePackDescriptor( + resolved_pkgs.get(), {}, {PDOpt::WithFluxes}); + auto v = desc.GetPack(md); + + const int nblocks = v.GetNBlocks(); + auto geom = Geometry::GetCoordinateSystem(md); + Real result = 0.0; + const bool is_monopole_cart = + (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleCart)); + const bool is_monopole_sph = + (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleSph)); + + parthenon::par_reduce( + parthenon::LoopPatternMDRange(), "Calculates mass accretion rate (SN)", + DevExecSpace(), 0, nblocks - 1, 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) { + const auto &coords = v.GetCoordinates(b); + const Real vol = coords.CellVolume(k, j, i); + + Real C[NDFULL]; + geom.Coords(CellLocation::Cent, b, k, j, i, C); + Real r = std::sqrt(C[1] * C[1] + C[2] * C[2] + C[3] * C[3]); + if (r <= rc) { + if (gain) { + auto rad = pmb->packages.Get("radiation").get(); + const parthenon::AllReduce *pdo_gain_reducer = + rad->MutableParam>("do_gain_reducer"); + const bool do_gain = pdo_gain_reducer->val; + 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 + auto analysis = pmb->packages.Get("analysis").get(); + const Real inside_pns_threshold = + analysis->Param("inside_pns_threshold"); + bool is_inside_pns = (r < inside_pns_threshold); // checks that inside PNS + if (do_gain && (is_inside_pns || is_netheat)) { + + lresult += ComputeDivInPillbox( + ndim, b, k, j, i, coords, [=](int b, int dir, int k, int j, int i) { + return v.flux(b, dir, c::density(), k, j, i); + }); + } else { + lresult += 0; + } + + } else { + lresult += ComputeDivInPillbox(ndim, b, k, j, i, coords, + [=](int b, int dir, int k, int j, int i) { + return v.flux(b, dir, c::density(), k, j, i); + }); + } + + } else { + lresult += 0.0; + } + }, + Kokkos::Sum(result)); + return result; +} // mdot + } // namespace History diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index 8dbfb164..a6cd5e2a 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -44,6 +44,7 @@ Real ReduceJetEnergyFlux(MeshData *md); Real ReduceJetMomentumFlux(MeshData *md); Real ReduceMagneticFluxPhi(MeshData *md); void ReduceLocalizationFunction(MeshData *md); +Real CalculateMdot(MeshData *md, Real rc, bool gain); template Real ReduceOneVar(MeshData *md, const std::string &varname, int idx = 0) { @@ -85,18 +86,20 @@ Real ReduceOneVar(MeshData *md, const std::string &varname, int idx = 0) { } template -Real ReduceInGain(MeshData *md, int idx = 0) { +Real ReduceInGain(MeshData *md, bool is_conserved, int idx = 0) { const auto ib = md->GetBoundsI(IndexDomain::interior); const auto jb = md->GetBoundsJ(IndexDomain::interior); const auto kb = md->GetBoundsK(IndexDomain::interior); namespace iv = internal_variables; using parthenon::MakePackDescriptor; + auto *pmb = md->GetParentPointer(); Mesh *pmesh = md->GetMeshPointer(); auto &resolved_pkgs = pmesh->resolved_packages; const int ndim = pmesh->ndim; static auto desc = - MakePackDescriptor(resolved_pkgs.get()); + MakePackDescriptor( + resolved_pkgs.get()); auto v = desc.GetPack(md); PARTHENON_REQUIRE_THROWS(v.ContainsHost(0, iv::GcovHeat(), iv::GcovCool()), "Must be doing SN simulation"); @@ -112,11 +115,20 @@ Real ReduceInGain(MeshData *md, 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) > 0); + 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 + auto analysis = pmb->packages.Get("analysis").get(); + const Real outside_pns_threshold = analysis->Param("outside_pns_threshold"); + 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); const Real vol = coords.CellVolume(k, j, i); - lresult += gdet * vol * is_netheat * v(b, Varname(idx), k, j, i); + if (is_conserved) { + lresult += vol * is_netheat * is_outside_pns * v(b, Varname(idx), k, j, i); + } else { + lresult += + gdet * vol * is_netheat * is_outside_pns * v(b, Varname(idx), k, j, i); + } }, Kokkos::Sum(result)); return result; diff --git a/src/progenitor/progenitordata.cpp b/src/progenitor/progenitordata.cpp index 92043be8..92d6641e 100644 --- a/src/progenitor/progenitordata.cpp +++ b/src/progenitor/progenitordata.cpp @@ -116,14 +116,25 @@ std::shared_ptr Initialize(ParameterInput *pin) { using parthenon::HistoryOutputVar; parthenon::HstVar_list hst_vars = {}; auto Mgain = [](MeshData *md) { - return ReduceInGain(md, 0); + return ReduceInGain(md, 1, 0); }; auto Qgain = [](MeshData *md) { - return ReduceInGain(md, 0) - - ReduceInGain(md, 0); + return ReduceInGain(md, 0, 0) - + ReduceInGain(md, 0, 0); + }; + auto Mdot400 = [](MeshData *md) { + Real rc = 0.04; + return History::CalculateMdot(md, rc, 0); + }; + + Real x1max = pin->GetReal("parthenon/mesh", "x1max"); + auto Mdot_gain = [x1max](MeshData *md) { + return History::CalculateMdot(md, x1max, 1); }; hst_vars.emplace_back(HistoryOutputVar(HstSum, Mgain, "Mgain")); hst_vars.emplace_back(HistoryOutputVar(HstSum, Qgain, "total net heat")); + hst_vars.emplace_back(HistoryOutputVar(HstSum, Mdot400, "Mdot at 400km")); + hst_vars.emplace_back(HistoryOutputVar(HstSum, Mdot_gain, "Mdot gain")); params.Add(parthenon::hist_param_key, hst_vars); return progenitor_pkg; From 2b604be4f87299c007e03e4df8af05dd8c5dead5 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 20 Mar 2024 12:29:16 -0400 Subject: [PATCH 2/7] comments applied --- src/analysis/analysis.cpp | 3 ++ src/analysis/history.cpp | 47 ++++++++----------------------- src/analysis/history.hpp | 24 ++++++++++++++-- src/progenitor/progenitordata.cpp | 13 ++++----- 4 files changed, 42 insertions(+), 45 deletions(-) diff --git a/src/analysis/analysis.cpp b/src/analysis/analysis.cpp index e7adde50..5ae92d8a 100644 --- a/src/analysis/analysis.cpp +++ b/src/analysis/analysis.cpp @@ -11,12 +11,15 @@ std::shared_ptr Initialize(ParameterInput *pin) { 2.42e-5); // corresponds to entropy > 3 kb/baryon Real inside_pns_threshold = pin->GetOrAddReal("analysis", "inside_pns_threshold", 0.008); // corresponds to r < 80 km + Real radius_cutoff_mdot = + pin->GetOrAddReal("analysis", "radius_cutoff_mdot", 0.04); // default 400km if (sigma < 0) { PARTHENON_THROW("sigma must be greater than 0"); } params.Add("sigma", sigma); params.Add("outside_pns_threshold", outside_pns_threshold); params.Add("inside_pns_threshold", inside_pns_threshold); + params.Add("radius_cutoff_mdot", radius_cutoff_mdot); return analysis_pkg; } // initialize diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index aa0f10ac..8fe1a098 100644 --- a/src/analysis/history.cpp +++ b/src/analysis/history.cpp @@ -280,26 +280,6 @@ void ReduceLocalizationFunction(MeshData *md) { } // exp -using namespace parthenon::package::prelude; -template -KOKKOS_INLINE_FUNCTION Real ComputeDivInPillbox(int ndim, int b, int k, int j, int i, - const parthenon::Coordinates_t &coords, - const F &f) { - Real div_mass_flux_integral; - div_mass_flux_integral = - -(f(b, 1, k, j, i + 1) - f(b, 1, k, j, i)) * coords.FaceArea(k, j, i); - - if (ndim >= 2) { - div_mass_flux_integral += - -(f(b, 2, k, j + 1, i) - f(b, 2, k, j, i)) * coords.FaceArea(k, j, i); - } - if (ndim >= 3) { - div_mass_flux_integral += - -(f(b, 3, k + 1, j, i) - f(b, 3, k, j, i)) * coords.FaceArea(k, j, i); - } - return div_mass_flux_integral; -} - // This function calculates mass accretion rate which I defined as // Mdot=Int(dV*d/dx^i{detg*rho*U^i}) where detg is the determinant of four metric, U is // four-velocity, and dV=d^3x @@ -323,10 +303,12 @@ Real CalculateMdot(MeshData *md, Real rc, bool gain) { const int nblocks = v.GetNBlocks(); auto geom = Geometry::GetCoordinateSystem(md); Real result = 0.0; - const bool is_monopole_cart = - (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleCart)); - const bool is_monopole_sph = - (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleSph)); + auto rad = pmb->packages.Get("radiation").get(); + const parthenon::AllReduce *pdo_gain_reducer = + rad->MutableParam>("do_gain_reducer"); + const bool do_gain = pdo_gain_reducer->val; + auto analysis = pmb->packages.Get("analysis").get(); + const Real inside_pns_threshold = analysis->Param("inside_pns_threshold"); parthenon::par_reduce( parthenon::LoopPatternMDRange(), "Calculates mass accretion rate (SN)", @@ -340,20 +322,13 @@ 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) { - auto rad = pmb->packages.Get("radiation").get(); - const parthenon::AllReduce *pdo_gain_reducer = - rad->MutableParam>("do_gain_reducer"); - const bool do_gain = pdo_gain_reducer->val; 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 - auto analysis = pmb->packages.Get("analysis").get(); - const Real inside_pns_threshold = - analysis->Param("inside_pns_threshold"); bool is_inside_pns = (r < inside_pns_threshold); // checks that inside PNS if (do_gain && (is_inside_pns || is_netheat)) { - lresult += ComputeDivInPillbox( + lresult += -ComputeDivInPillbox( ndim, b, k, j, i, coords, [=](int b, int dir, int k, int j, int i) { return v.flux(b, dir, c::density(), k, j, i); }); @@ -362,10 +337,10 @@ Real CalculateMdot(MeshData *md, Real rc, bool gain) { } } else { - lresult += ComputeDivInPillbox(ndim, b, k, j, i, coords, - [=](int b, int dir, int k, int j, int i) { - return v.flux(b, dir, c::density(), k, j, i); - }); + lresult += -ComputeDivInPillbox( + ndim, b, k, j, i, coords, [=](int b, int dir, int k, int j, int i) { + return v.flux(b, dir, c::density(), k, j, i); + }); } } else { diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index a6cd5e2a..bc5cfcfe 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -108,6 +108,8 @@ Real ReduceInGain(MeshData *md, bool is_conserved, int idx = 0) { Real result = 0.0; auto geom = Geometry::GetCoordinateSystem(md); + auto analysis = pmb->packages.Get("analysis").get(); + const Real outside_pns_threshold = analysis->Param("outside_pns_threshold"); parthenon::par_reduce( parthenon::LoopPatternMDRange(), @@ -117,8 +119,6 @@ Real ReduceInGain(MeshData *md, bool is_conserved, int idx = 0) { 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 - auto analysis = pmb->packages.Get("analysis").get(); - const Real outside_pns_threshold = analysis->Param("outside_pns_threshold"); 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); @@ -134,6 +134,26 @@ Real ReduceInGain(MeshData *md, bool is_conserved, int idx = 0) { return result; } +using namespace parthenon::package::prelude; +template +KOKKOS_INLINE_FUNCTION Real ComputeDivInPillbox(int ndim, int b, int k, int j, int i, + const parthenon::Coordinates_t &coords, + const F &f) { + Real div_mass_flux_integral; + div_mass_flux_integral = + (f(b, 1, k, j, i + 1) - f(b, 1, k, j, i)) * coords.FaceArea(k, j, i); + + if (ndim >= 2) { + div_mass_flux_integral += + (f(b, 2, k, j + 1, i) - f(b, 2, k, j, i)) * coords.FaceArea(k, j, i); + } + if (ndim >= 3) { + div_mass_flux_integral += + (f(b, 3, k + 1, j, i) - f(b, 3, k, j, i)) * coords.FaceArea(k, j, i); + } + return div_mass_flux_integral; +} + } // namespace History #endif // ANALYSIS_HISTORY_HPP_ diff --git a/src/progenitor/progenitordata.cpp b/src/progenitor/progenitordata.cpp index 92d6641e..5c67fbdc 100644 --- a/src/progenitor/progenitordata.cpp +++ b/src/progenitor/progenitordata.cpp @@ -1,6 +1,7 @@ #include #include +#include "analysis/analysis.hpp" #include "analysis/history.hpp" #include "ascii_reader.hpp" #include "geometry/geometry.hpp" @@ -109,6 +110,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { params.Add("Srr_adm_dev", Srr_adm_dev); // Reductions + const Real rc = pin->GetReal("analysis", "radius_cutoff_mdot"); auto HstSum = parthenon::UserHistoryOperation::sum; auto HstMax = parthenon::UserHistoryOperation::max; using History::ReduceInGain; @@ -116,16 +118,13 @@ std::shared_ptr Initialize(ParameterInput *pin) { using parthenon::HistoryOutputVar; parthenon::HstVar_list hst_vars = {}; auto Mgain = [](MeshData *md) { - return ReduceInGain(md, 1, 0); + return ReduceInGain(md, 1, 0); }; auto Qgain = [](MeshData *md) { - return ReduceInGain(md, 0, 0) - - ReduceInGain(md, 0, 0); - }; - auto Mdot400 = [](MeshData *md) { - Real rc = 0.04; - return History::CalculateMdot(md, rc, 0); + return ReduceInGain(md, 0, 0) - + ReduceInGain(md, 0, 0); }; + auto Mdot400 = [rc](MeshData *md) { return History::CalculateMdot(md, rc, 0); }; Real x1max = pin->GetReal("parthenon/mesh", "x1max"); auto Mdot_gain = [x1max](MeshData *md) { From d4073c5c37a15786b2075b88374a2000d6ffcc6a Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 20 Mar 2024 16:20:15 -0400 Subject: [PATCH 3/7] moved thresholds from analysis pkg to progenitor pkg --- src/analysis/analysis.cpp | 10 ---------- src/analysis/history.cpp | 4 ++-- src/analysis/history.hpp | 4 ++-- src/progenitor/progenitordata.cpp | 18 ++++++++++++++++-- 4 files changed, 20 insertions(+), 16 deletions(-) diff --git a/src/analysis/analysis.cpp b/src/analysis/analysis.cpp index 5ae92d8a..54837fb3 100644 --- a/src/analysis/analysis.cpp +++ b/src/analysis/analysis.cpp @@ -6,20 +6,10 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto analysis_pkg = std::make_shared("analysis"); Params ¶ms = analysis_pkg->AllParams(); Real sigma = pin->GetOrAddReal("analysis", "sigma", 2e-3); - Real outside_pns_threshold = - pin->GetOrAddReal("analysis", "outside_pns_threshold", - 2.42e-5); // corresponds to entropy > 3 kb/baryon - Real inside_pns_threshold = pin->GetOrAddReal("analysis", "inside_pns_threshold", - 0.008); // corresponds to r < 80 km - Real radius_cutoff_mdot = - pin->GetOrAddReal("analysis", "radius_cutoff_mdot", 0.04); // default 400km if (sigma < 0) { PARTHENON_THROW("sigma must be greater than 0"); } params.Add("sigma", sigma); - params.Add("outside_pns_threshold", outside_pns_threshold); - params.Add("inside_pns_threshold", inside_pns_threshold); - params.Add("radius_cutoff_mdot", radius_cutoff_mdot); return analysis_pkg; } // initialize diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index 8fe1a098..4d35fcff 100644 --- a/src/analysis/history.cpp +++ b/src/analysis/history.cpp @@ -307,8 +307,8 @@ Real CalculateMdot(MeshData *md, Real rc, bool gain) { const parthenon::AllReduce *pdo_gain_reducer = rad->MutableParam>("do_gain_reducer"); const bool do_gain = pdo_gain_reducer->val; - auto analysis = pmb->packages.Get("analysis").get(); - const Real inside_pns_threshold = analysis->Param("inside_pns_threshold"); + auto progenitor = pmb->packages.Get("progenitor").get(); + const Real inside_pns_threshold = progenitor->Param("inside_pns_threshold"); parthenon::par_reduce( parthenon::LoopPatternMDRange(), "Calculates mass accretion rate (SN)", diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index bc5cfcfe..98901a8d 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -108,8 +108,8 @@ Real ReduceInGain(MeshData *md, bool is_conserved, int idx = 0) { Real result = 0.0; auto geom = Geometry::GetCoordinateSystem(md); - auto analysis = pmb->packages.Get("analysis").get(); - const Real outside_pns_threshold = analysis->Param("outside_pns_threshold"); + auto progenitor = pmb->packages.Get("progenitor").get(); + const Real outside_pns_threshold = progenitor->Param("outside_pns_threshold"); parthenon::par_reduce( parthenon::LoopPatternMDRange(), diff --git a/src/progenitor/progenitordata.cpp b/src/progenitor/progenitordata.cpp index 5c67fbdc..9cb55316 100644 --- a/src/progenitor/progenitordata.cpp +++ b/src/progenitor/progenitordata.cpp @@ -84,6 +84,15 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto S_adm_dev = S_adm.getOnDevice(); auto Srr_adm_dev = Srr_adm.getOnDevice(); + // Post-processing params + Real outside_pns_threshold = + pin->GetOrAddReal("progenitor", "outside_pns_threshold", + 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 radius_cutoff_mdot = + pin->GetOrAddReal("progenitor", "radius_cutoff_mdot", 0.04); // default 400km + // Add Params params.Add("mass_density", mass_density); params.Add("temp", temp); @@ -109,8 +118,11 @@ std::shared_ptr Initialize(ParameterInput *pin) { params.Add("S_adm_dev", S_adm_dev); params.Add("Srr_adm_dev", Srr_adm_dev); + params.Add("outside_pns_threshold", outside_pns_threshold); + params.Add("inside_pns_threshold", inside_pns_threshold); + params.Add("radius_cutoff_mdot", radius_cutoff_mdot); + // Reductions - const Real rc = pin->GetReal("analysis", "radius_cutoff_mdot"); auto HstSum = parthenon::UserHistoryOperation::sum; auto HstMax = parthenon::UserHistoryOperation::max; using History::ReduceInGain; @@ -124,7 +136,9 @@ std::shared_ptr Initialize(ParameterInput *pin) { return ReduceInGain(md, 0, 0) - ReduceInGain(md, 0, 0); }; - auto Mdot400 = [rc](MeshData *md) { return History::CalculateMdot(md, rc, 0); }; + auto Mdot400 = [radius_cutoff_mdot](MeshData *md) { + return History::CalculateMdot(md, radius_cutoff_mdot, 0); + }; Real x1max = pin->GetReal("parthenon/mesh", "x1max"); auto Mdot_gain = [x1max](MeshData *md) { From 994e2ccbba75087d0e10f371fa10196c80062cb7 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 20 Mar 2024 17:17:44 -0400 Subject: [PATCH 4/7] Mdot can now take a set or thresholds --- src/progenitor/progenitordata.cpp | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/src/progenitor/progenitordata.cpp b/src/progenitor/progenitordata.cpp index 9cb55316..79d91f79 100644 --- a/src/progenitor/progenitordata.cpp +++ b/src/progenitor/progenitordata.cpp @@ -8,8 +8,8 @@ #include "microphysics/eos_phoebus/eos_phoebus.hpp" #include "monopole_gr/monopole_gr.hpp" #include "pgen/pgen.hpp" +#include "phoebus_utils/unit_conversions.hpp" #include "progenitordata.hpp" - namespace Progenitor { std::shared_ptr Initialize(ParameterInput *pin) { @@ -90,8 +90,11 @@ 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 radius_cutoff_mdot = - pin->GetOrAddReal("progenitor", "radius_cutoff_mdot", 0.04); // default 400km + UnitConversions units(pin); + Real LengthCGSToCode = units.GetLengthCGSToCode(); + Real rc = 400e5 * LengthCGSToCode; + auto mdot_radii = pin->GetOrAddVector("progenitor", "mdot_radii", + std::vector{rc}); // default 400km // Add Params params.Add("mass_density", mass_density); @@ -120,7 +123,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("radius_cutoff_mdot", radius_cutoff_mdot); + params.Add("mdot_radii", mdot_radii); // Reductions auto HstSum = parthenon::UserHistoryOperation::sum; @@ -136,9 +139,14 @@ std::shared_ptr Initialize(ParameterInput *pin) { return ReduceInGain(md, 0, 0) - ReduceInGain(md, 0, 0); }; - auto Mdot400 = [radius_cutoff_mdot](MeshData *md) { - return History::CalculateMdot(md, radius_cutoff_mdot, 0); - }; + for (auto rc : mdot_radii) { + auto Mdot = [rc, LengthCGSToCode](MeshData *md) { + return History::CalculateMdot(md, rc, 0); + }; + hst_vars.emplace_back(HistoryOutputVar( + HstSum, Mdot, + "Mdot at r = " + std::to_string(int(rc / LengthCGSToCode / 1.e5)) + "km")); + } Real x1max = pin->GetReal("parthenon/mesh", "x1max"); auto Mdot_gain = [x1max](MeshData *md) { @@ -146,7 +154,6 @@ std::shared_ptr Initialize(ParameterInput *pin) { }; hst_vars.emplace_back(HistoryOutputVar(HstSum, Mgain, "Mgain")); hst_vars.emplace_back(HistoryOutputVar(HstSum, Qgain, "total net heat")); - hst_vars.emplace_back(HistoryOutputVar(HstSum, Mdot400, "Mdot at 400km")); hst_vars.emplace_back(HistoryOutputVar(HstSum, Mdot_gain, "Mdot gain")); params.Add(parthenon::hist_param_key, hst_vars); From 7f3bc17e54529a89f1827d43ca7d41422baff876 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 20 Mar 2024 19:34:36 -0400 Subject: [PATCH 5/7] 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"); From 4658fa09ebfe72b7600a132ebebaa42f577a8c80 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 20 Mar 2024 20:56:28 -0400 Subject: [PATCH 6/7] typo --- src/analysis/history.hpp | 4 ++++ src/progenitor/progenitordata.cpp | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index 55dfdae3..44491e59 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -123,6 +123,10 @@ Real ReduceInGain(MeshData *md, bool is_conserved, int idx = 0) { 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 + std::cout << "is_outside_pns=" << is_outside_pns + << "outside_pns_threshold=" << outside_pns_threshold << std::endl; + std::cout << "is_netheat=" << is_netheat + << "net_heat_threshold=" << net_heat_threshold << std::endl; const auto &coords = v.GetCoordinates(b); const Real vol = coords.CellVolume(k, j, i); if (is_conserved) { diff --git a/src/progenitor/progenitordata.cpp b/src/progenitor/progenitordata.cpp index 535a1424..39de2da3 100644 --- a/src/progenitor/progenitordata.cpp +++ b/src/progenitor/progenitordata.cpp @@ -90,7 +90,7 @@ 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", + Real net_heat_threshold = pin->GetOrAddReal("progenitor", "net_heat_threshold", 1e-8); // corresponds to r < 80 km UnitConversions units(pin); Real LengthCGSToCode = units.GetLengthCGSToCode(); From c6f72c55fce4b3a5f3142599c45588608d0fa4da Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 20 Mar 2024 20:57:33 -0400 Subject: [PATCH 7/7] print statements deleted --- src/analysis/history.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index 44491e59..55dfdae3 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -123,10 +123,6 @@ Real ReduceInGain(MeshData *md, bool is_conserved, int idx = 0) { 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 - std::cout << "is_outside_pns=" << is_outside_pns - << "outside_pns_threshold=" << outside_pns_threshold << std::endl; - std::cout << "is_netheat=" << is_netheat - << "net_heat_threshold=" << net_heat_threshold << std::endl; const auto &coords = v.GetCoordinates(b); const Real vol = coords.CellVolume(k, j, i); if (is_conserved) {