From 5cc2bb52e6778abe20fc401913f28a2a5b06f23d Mon Sep 17 00:00:00 2001
From: mari2895 <mgogilashvili@fsu.edu>
Date: Thu, 7 Mar 2024 16:09:08 -0500
Subject: [PATCH] central density added in hystory output

---
 src/analysis/history.cpp           | 31 ++++++++++++++++++++++++++++++
 src/analysis/history.hpp           |  1 +
 src/fluid/fluid.cpp                |  6 ++++++
 src/radiation/cooling_function.cpp |  1 +
 4 files changed, 39 insertions(+)

diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp
index 5b044939..ace2031a 100644
--- a/src/analysis/history.cpp
+++ b/src/analysis/history.cpp
@@ -235,4 +235,35 @@ Real ReduceMagneticFluxPhi(MeshData<Real> *md) {
   return 0.5 * result; // 0.5 \int detg B^r dx2 dx3
 } // Phi
 
+Real ReduceCentralDensitySN(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 p = fluid_prim;
+  auto *pmb = rc->GetParentPointer();
+  using parthenon::MakePackDescriptor;
+  Mesh *pmesh = rc->GetMeshPointer();
+  auto &resolved_pkgs = pmesh->resolved_packages;
+  const int ndim = pmesh->ndim;
+
+  static auto desc = MakePackDescriptor<p::density>(resolved_pkgs.get());
+  auto v = desc.GetPack(rc);
+  const int nblocks = v.GetNBlocks();
+
+  parthenon::par_for(
+      DEFAULT_LOOP_PATTERN, "CentralDensity", 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) {
+        const Real x = coords.Xc<1>(k, j, i);
+        const Real y = coords.Xc<2>(k, j, i);
+        const Real z = coords.Xc<3>(k, j, i);
+        const Real sigma = analysis->Param<Real>("sigma");
+        const Real pi = 3.14;
+        const Real rhoc = v(b, p::density(), k, j, i) *
+                          std::exp(-(x * x + y * y + z * z) / sigma / sigma) / pi /
+                          sigma; // sigma > 0
+      });
+  return rhoc;
+} // rhoc
+
 } // namespace History
diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp
index 04a61d60..c154b3dd 100644
--- a/src/analysis/history.hpp
+++ b/src/analysis/history.hpp
@@ -41,6 +41,7 @@ Real ReduceMassAccretionRate(MeshData<Real> *md);
 Real ReduceJetEnergyFlux(MeshData<Real> *md);
 Real ReduceJetMomentumFlux(MeshData<Real> *md);
 Real ReduceMagneticFluxPhi(MeshData<Real> *md);
+Real ReduceCentralDensitySN(MeshData<Real> *md);
 
 template <typename Reducer_t>
 Real ReduceOneVar(MeshData<Real> *md, const std::string &varname, int idx = 0) {
diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp
index 2aade70d..340dd41b 100644
--- a/src/fluid/fluid.cpp
+++ b/src/fluid/fluid.cpp
@@ -321,6 +321,12 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceMass, "total baryon number"));
   hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceEn, "total conserved energy tau"));
 
+  auto CentralDensitySN = [](MeshData<Real> *md) {
+    Real rhoc = Hystory::ReduceCentralDensity(md);
+    return ReduceOneVar<Kokkos::Sum<Real>>(md, rhoc, 0);
+  };
+  hst_vars.emplace_back(HistoryOutputVar(HstSum, CentralDensitySN, "central density SN"));
+
   for (int d = 0; d < 3; ++d) {
     auto ReduceMom = [d](MeshData<Real> *md) {
       return History::ReduceOneVar<Kokkos::Sum<Real>>(md, fluid_cons::momentum::name(),
diff --git a/src/radiation/cooling_function.cpp b/src/radiation/cooling_function.cpp
index b6171ec6..03b30945 100644
--- a/src/radiation/cooling_function.cpp
+++ b/src/radiation/cooling_function.cpp
@@ -11,6 +11,7 @@
 // distribute copies to the public, perform publicly and display
 // publicly, and to permit others to do so.
 
+#include "analysis/history.hpp"
 #include "geometry/geometry.hpp"
 #include "light_bulb_constants.hpp"
 #include "phoebus_utils/variables.hpp"