outputs for SN postprocessing: max density, M_gain, and Qdot_gain

src/analysis/history.cpp

@@ -12,6 +12,7 @@
 // publicly, and to permit others to do so.
 
 #include "history.hpp"
+#include "analysis/analysis.hpp"
 #include "geometry/geometry.hpp"
 #include "geometry/geometry_utils.hpp"
 #include "history_utils.hpp"
@@ -235,4 +236,46 @@ Real ReduceMagneticFluxPhi(MeshData<Real> *md) {
   return 0.5 * result; // 0.5 \int detg B^r dx2 dx3
 } // Phi
 
+// SN analysis
+
+void ReduceLocalizationFunction(MeshData<Real> *md) {
+  const auto ib = md->GetBoundsI(IndexDomain::interior);
+  const auto jb = md->GetBoundsJ(IndexDomain::interior);
+  const auto kb = md->GetBoundsK(IndexDomain::interior);
+  namespace diag = diagnostic_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<diag::localization_function>(resolved_pkgs.get());
+  auto v = desc.GetPack(md);
+  const int nblocks = v.GetNBlocks();
+  auto geom = Geometry::GetCoordinateSystem(md);
+
+  parthenon::par_for(
+      parthenon::LoopPatternMDRange(), "Central Density for 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) {
+        auto &coords = v.GetCoordinates(b);
+        auto analysis = pmb->packages.Get("analysis").get();
+        const Real x[3] = {coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
+                           coords.Xc<3>(k, j, i)};
+        const Real sigma = analysis->Param<Real>("sigma");
+        Real gam[3][3];
+        Real r2 = 0;
+        geom.Metric(CellLocation::Cent, 0, k, j, i, gam);
+        for (int n = 0; n < 3; ++n) {
+          for (int m = 0; m < 3; ++m) {
+            r2 += gam[n][m] * x[n] * x[m];
+          }
+        }
+        v(b, diag::localization_function(), k, j, i) = std::exp(-r2 / sigma / sigma);
+      });
+
+} // exp (this function returns normalization function that is used for localizing
+  // quantities at the center, or at some particular case. For example SN diagnostics
+  // oftec computes quantities at 400 km.)
+
 } // namespace History

src/analysis/history.hpp

@@ -17,12 +17,14 @@
 #include <string>
 #include <vector>
 
+#include "geometry/geometry.hpp"
+#include "geometry/geometry_utils.hpp"
+#include "phoebus_utils/variables.hpp"
 #include <kokkos_abstraction.hpp>
+#include <parthenon/driver.hpp>
 #include <parthenon/package.hpp>
 #include <utils/error_checking.hpp>
 
-#include "phoebus_utils/variables.hpp"
-
 using namespace parthenon::package::prelude;
 
 /*
@@ -41,6 +43,7 @@ Real ReduceMassAccretionRate(MeshData<Real> *md);
 Real ReduceJetEnergyFlux(MeshData<Real> *md);
 Real ReduceJetMomentumFlux(MeshData<Real> *md);
 Real ReduceMagneticFluxPhi(MeshData<Real> *md);
+void ReduceLocalizationFunction(MeshData<Real> *md);
 
 template <typename Reducer_t>
 Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
@@ -58,21 +61,71 @@ Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
   // We choose to apply volume weighting when using the sum reduction.
   // This assumes that for "Sum" variables, the variable is densitized, meaning
   // it already contains a factor of the measure: sqrt(det(gamma))
+  Real result = 0.0;
+  Reducer_t reducer(result);
+
   const bool volume_weighting =
-      std::is_same<Reducer_t, Kokkos::Sum<Real, HostExecSpace>>::value;
+      std::is_same<Reducer_t, Kokkos::Sum<Real, HostExecSpace>>::value ||
+      std::is_same<Reducer_t, Kokkos::Sum<Real, DevExecSpace>>::value ||
+      std::is_same<Reducer_t, Kokkos::Sum<Real, Kokkos::DefaultExecutionSpace>>::value;
 
+  parthenon::par_reduce(
+      parthenon::LoopPatternMDRange(), "Phoebus History for " + varname, DevExecSpace(),
+      0, pack.GetDim(5) - 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) {
+        // join is a Kokkos construct
+        // that automatically does the
+        // reduction operation locally
+        const auto &coords = pack.GetCoords(b);
+        const Real vol = volume_weighting ? coords.CellVolume(k, j, i) : 1.0;
+        reducer.join(lresult, pack(b, ivar.first + idx, k, j, i) * vol);
+      },
+      reducer);
+  return result;
+}
+
+template <typename Reducer_t>
+Real ReduceInGain(MeshData<Real> *md, const std::string &varname, int idx = 0) {
+  const auto ib = md->GetBoundsI(IndexDomain::interior);
+  const auto jb = md->GetBoundsJ(IndexDomain::interior);
+  const auto kb = md->GetBoundsK(IndexDomain::interior);
+
+  PackIndexMap imap;
+  namespace iv = internal_variables;
+  std::vector<std::string> vars = {varname, iv::GcovHeat::name(), iv::GcovCool::name()};
+  const auto pack = md->PackVariables(vars, imap);
+  const auto ivar = imap[varname];
+  const auto iheat = imap[iv::GcovHeat::name()].first;
+  const auto icool = imap[iv::GcovCool::name()].first;
+  auto *pmb = md->GetParentPointer();
+  PARTHENON_REQUIRE(ivar.first >= 0, "Var must exist");
+  PARTHENON_REQUIRE(ivar.second >= ivar.first + idx, "Var must exist");
+
+  // We choose to apply volume weighting when using the sum reduction.
+  // This assumes that for "Sum" variables, the variable is densitized, meaning
+  // it already contains a factor of the measure: sqrt(det(gamma))
   Real result = 0.0;
   Reducer_t reducer(result);
+
+  const bool volume_weighting =
+      std::is_same<Reducer_t, Kokkos::Sum<Real, HostExecSpace>>::value ||
+      std::is_same<Reducer_t, Kokkos::Sum<Real, DevExecSpace>>::value ||
+      std::is_same<Reducer_t, Kokkos::Sum<Real, Kokkos::DefaultExecutionSpace>>::value;
+  auto geom = Geometry::GetCoordinateSystem(md);
+
   parthenon::par_reduce(
       parthenon::LoopPatternMDRange(), "Phoebus History for " + varname, DevExecSpace(),
       0, pack.GetDim(5) - 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) {
         // join is a Kokkos construct
         // that automatically does the
         // reduction operation locally
+        Real gdet = geom.DetGamma(CellLocation::Cent, 0, k, j, i);
+        bool is_netheat = (pack(b, iheat, k, j, i) - pack(b, icool, k, j, i) > 0);
         const auto &coords = pack.GetCoords(b);
         const Real vol = volume_weighting ? coords.CellVolume(k, j, i) : 1.0;
-        reducer.join(lresult, pack(b, ivar.first + idx, k, j, i) * vol);
+        reducer.join(lresult,
+                     pack(b, ivar.first + idx, k, j, i) * gdet * vol * is_netheat);
       },
       reducer);
   return result;

src/fluid/fluid.cpp

@@ -207,6 +207,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   physics->AddField(p::entropy::name(), mprim_scalar);
   physics->AddField(p::cs::name(), mprim_scalar);
   physics->AddField(diag::ratio_divv_cs::name(), mprim_scalar);
+  physics->AddField(diag::localization_function::name(), mprim_scalar);
   physics->AddField(diag::entropy_z_0::name(), mprim_scalar);
   physics->AddField(p::gamma1::name(), mprim_scalar);
   if (ye) {
@@ -309,6 +310,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   // Reductions
   // By default compute integrated value of scalar conserved vars
   auto HstSum = parthenon::UserHistoryOperation::sum;
+  auto HstMax = parthenon::UserHistoryOperation::max;
+  using History::ReduceInGain;
   using History::ReduceOneVar;
   using parthenon::HistoryOutputVar;
   parthenon::HstVar_list hst_vars = {};
@@ -319,6 +322,23 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto ReduceEn = [](MeshData<Real> *md) {
     return ReduceOneVar<Kokkos::Sum<Real>>(md, fluid_cons::energy::name(), 0);
   };
+  auto MaxDensitySN = [](MeshData<Real> *md) {
+    return ReduceOneVar<Kokkos::Max<Real>>(md, fluid_prim::density::name(), 0);
+  };
+  auto Mgain = [](MeshData<Real> *md) {
+    return ReduceInGain<Kokkos::Sum<Real, HostExecSpace>>(md, fluid_prim::density::name(),
+                                                          0);
+  };
+  auto Qgain = [](MeshData<Real> *md) {
+    return ReduceInGain<Kokkos::Sum<Real, HostExecSpace>>(
+               md, internal_variables::GcovHeat::name(), 0) -
+           ReduceInGain<Kokkos::Sum<Real, HostExecSpace>>(
+               md, internal_variables::GcovCool::name(), 0);
+  };
+
+  hst_vars.emplace_back(HistoryOutputVar(HstMax, MaxDensitySN, "maximum density SN"));
+  hst_vars.emplace_back(HistoryOutputVar(HstSum, Mgain, "Mgain"));
+  hst_vars.emplace_back(HistoryOutputVar(HstSum, Qgain, "total net heat"));
   hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceMass, "total baryon number"));
   hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceEn, "total conserved energy tau"));
 

src/phoebus_utils/variables.hpp

@@ -122,6 +122,7 @@ namespace diagnostic_variables {
 VARIABLE_NONS(divb);
 VARIABLE_NONS(ratio_divv_cs);
 VARIABLE_NONS(entropy_z_0);
+VARIABLE_NONS(localization_function);
 VARIABLE_CUSTOM(divf, flux_divergence);
 VARIABLE_NONS(src_terms);
 VARIABLE_CUSTOM(r_divf, r.flux_divergence);