calculating Mdot for SN post-processing

src/analysis/analysis.cpp

@@ -6,10 +6,17 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto analysis_pkg = std::make_shared<StateDescriptor>("analysis");
   Params &params = 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

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,100 @@ void ReduceLocalizationFunction(MeshData<Real> *md) {
 
 } // exp
 
+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;
+}
+
+// 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;
+  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<bool> *pdo_gain_reducer =
+                rad->MutableParam<parthenon::AllReduce<bool>>("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<Real>("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<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");
@@ -112,11 +115,20 @@ Real ReduceInGain(MeshData<Real> *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<Real>("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<Real>(result));
   return result;

src/progenitor/progenitordata.cpp

@@ -116,14 +116,25 @@ 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<class 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<class internal_variables::GcovHeat>(md, 0, 0) -
+           ReduceInGain<class internal_variables::GcovCool>(md, 0, 0);
+  };
+  auto Mdot400 = [](MeshData<Real> *md) {
+    Real rc = 0.04;
+    return History::CalculateMdot(md, rc, 0);
+  };
+
+  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, 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;