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;