calculating Mass accretion rate for SN post-processing

src/analysis/history.cpp

@@ -18,6 +18,7 @@
 #include "history_utils.hpp"
 #include "phoebus_utils/relativity_utils.hpp"
 #include <interface/sparse_pack.hpp>
+#include <parthenon/package.hpp>
 
 namespace History {
 
@@ -279,4 +280,75 @@ void ReduceLocalizationFunction(MeshData<Real> *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<Real> *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<c::density, internal_variables::GcovHeat,
+                                        internal_variables::GcovCool>(
+      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<bool> *pdo_gain_reducer =
+      rad->MutableParam<parthenon::AllReduce<bool>>("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<Real>("inside_pns_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) >
+                               1.e-8); // 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<Real>(result));
+  return result;
+} // mdot
+
 } // namespace History

src/analysis/history.hpp

@@ -44,6 +44,7 @@ Real ReduceJetEnergyFlux(MeshData<Real> *md);
 Real ReduceJetMomentumFlux(MeshData<Real> *md);
 Real ReduceMagneticFluxPhi(MeshData<Real> *md);
 void ReduceLocalizationFunction(MeshData<Real> *md);
+Real CalculateMdot(MeshData<Real> *md, Real rc, bool gain);
 
 template <typename Reducer_t>
 Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
@@ -85,18 +86,20 @@ Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
 }
 
 template <typename Varname>
-Real ReduceInGain(MeshData<Real> *md, int idx = 0) {
+Real ReduceInGain(MeshData<Real> *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<Varname, iv::GcovHeat, iv::GcovCool>(resolved_pkgs.get());
+      MakePackDescriptor<Varname, fluid_prim::entropy, iv::GcovHeat, iv::GcovCool>(
+          resolved_pkgs.get());
   auto v = desc.GetPack(md);
   PARTHENON_REQUIRE_THROWS(v.ContainsHost(0, iv::GcovHeat(), iv::GcovCool()),
                            "Must be doing SN simulation");
@@ -105,23 +108,52 @@ Real ReduceInGain(MeshData<Real> *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<Real>("outside_pns_threshold");
 
   parthenon::par_reduce(
       parthenon::LoopPatternMDRange(),
       "Phoebus History for integrals in gain region (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) {
         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
+        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<Real>(result));
   return result;
 }
 
+using namespace parthenon::package::prelude;
+template <typename F>
+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<X1DIR>(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<X2DIR>(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<X3DIR>(k, j, i);
+  }
+  return div_mass_flux_integral;
+}
+
 } // namespace History
 
 #endif // ANALYSIS_HISTORY_HPP_

src/progenitor/progenitordata.cpp

@@ -1,14 +1,15 @@
 #include <memory>
 #include <vector>
 
+#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<StateDescriptor> Initialize(ParameterInput *pin) {
@@ -83,6 +84,18 @@ std::shared_ptr<StateDescriptor> 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
+  UnitConversions units(pin);
+  Real LengthCGSToCode = units.GetLengthCGSToCode();
+  Real rc = 400e5 * LengthCGSToCode;
+  auto mdot_radii = pin->GetOrAddVector("progenitor", "mdot_radii",
+                                        std::vector<Real>{rc}); // default 400km
+
   // Add Params
   params.Add("mass_density", mass_density);
   params.Add("temp", temp);
@@ -108,6 +121,10 @@ std::shared_ptr<StateDescriptor> 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("mdot_radii", mdot_radii);
+
   // Reductions
   auto HstSum = parthenon::UserHistoryOperation::sum;
   auto HstMax = parthenon::UserHistoryOperation::max;
@@ -116,14 +133,28 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   using parthenon::HistoryOutputVar;
   parthenon::HstVar_list hst_vars = {};
   auto Mgain = [](MeshData<Real> *md) {
-    return ReduceInGain<class fluid_prim::density>(md, 0);
+    return ReduceInGain<fluid_cons::density>(md, 1, 0);
   };
   auto Qgain = [](MeshData<Real> *md) {
-    return ReduceInGain<class internal_variables::GcovHeat>(md, 0) -
-           ReduceInGain<class internal_variables::GcovCool>(md, 0);
+    return ReduceInGain<internal_variables::GcovHeat>(md, 0, 0) -
+           ReduceInGain<internal_variables::GcovCool>(md, 0, 0);
+  };
+  for (auto rc : mdot_radii) {
+    auto Mdot = [rc, LengthCGSToCode](MeshData<Real> *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<Real> *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;