diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index d88e3d65..8f578dd0 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,76 @@ void ReduceLocalizationFunction(MeshData *md) { } // exp +// 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; + 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 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)", + 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) { + 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)) { + + 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..55dfdae3 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"); @@ -105,6 +108,9 @@ Real ReduceInGain(MeshData *md, int idx = 0) { Real result = 0.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(), @@ -113,15 +119,43 @@ Real ReduceInGain(MeshData *md, int idx = 0) { 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); + ((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); 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; } +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 92043be8..39de2da3 100644 --- a/src/progenitor/progenitordata.cpp +++ b/src/progenitor/progenitordata.cpp @@ -1,14 +1,15 @@ #include #include +#include "analysis/analysis.hpp" #include "analysis/history.hpp" #include "ascii_reader.hpp" #include "geometry/geometry.hpp" #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) { @@ -83,6 +84,19 @@ 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 net_heat_threshold = pin->GetOrAddReal("progenitor", "net_heat_threshold", + 1e-8); // corresponds to r < 80 km + UnitConversions units(pin); + Real LengthCGSToCode = units.GetLengthCGSToCode(); + auto mdot_radii = pin->GetOrAddVector("progenitor", "mdot_radii", + std::vector{400}); // default 400km + // Add Params params.Add("mass_density", mass_density); params.Add("temp", temp); @@ -108,6 +122,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("net_heat_threshold", net_heat_threshold); + params.Add("mdot_radii", mdot_radii); + // Reductions auto HstSum = parthenon::UserHistoryOperation::sum; auto HstMax = parthenon::UserHistoryOperation::max; @@ -116,14 +135,28 @@ 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); + }; + for (auto rc : mdot_radii) { + 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)) + "km")); + } + + 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, Mdot_gain, "Mdot gain")); params.Add(parthenon::hist_param_key, hst_vars); return progenitor_pkg;