From 346737dcad51639c4499caedd249feeed40dd5ca Mon Sep 17 00:00:00 2001
From: mari2895 <mgogilashvili@fsu.edu>
Date: Mon, 4 Mar 2024 13:56:43 -0500
Subject: [PATCH] entropy added in outputs

---
 src/fluid/fluid.cpp             |  1 +
 src/phoebus_driver.cpp          | 34 +++++++++++++++++++++++++++++++++
 src/phoebus_utils/variables.hpp |  1 +
 3 files changed, 36 insertions(+)

diff --git a/src/fluid/fluid.cpp b/src/fluid/fluid.cpp
index 82aa5826..cb4fe21f 100644
--- a/src/fluid/fluid.cpp
+++ b/src/fluid/fluid.cpp
@@ -204,6 +204,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   }
   physics->AddField(p::pressure::name(), mprim_scalar);
   physics->AddField(p::temperature::name(), mprim_scalar);
+  physics->AddField(p::entropy::name(), mprim_scalar);
   physics->AddField(p::gamma1::name(), mprim_scalar);
   if (ye) {
     physics->AddField(p::ye::name(), mprim_scalar);
diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp
index 432f9655..8e88f818 100644
--- a/src/phoebus_driver.cpp
+++ b/src/phoebus_driver.cpp
@@ -36,6 +36,7 @@
 #include "phoebus_driver.hpp"
 #include "phoebus_package.hpp"
 #include "phoebus_utils/robust.hpp"
+#include "phoebus_utils/variables.hpp"
 #include "progenitor/progenitordata.hpp"
 #include "radiation/radiation.hpp"
 #include "tov/tov.hpp"
@@ -1019,6 +1020,39 @@ TaskListStatus PhoebusDriver::MonteCarloStep() {
 void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) {
   auto tracer_pkg = pmb->packages.Get("tracers");
   bool do_tracers = tracer_pkg->Param<bool>("active");
+  bool is_stellarcollapse = (pin->GetString("eos", "type") == "StellarCollapse");
+
+  if (is_stellarcollapse) {
+    namespace p = fluid_prim;
+    auto rc = pmb->meshblock_data.Get().get();
+    PackIndexMap imap;
+    auto v = rc->PackVariables(
+        {p::density::name(), p::ye::name(), p::temperature::name(), p::entropy::name()},
+        imap);
+
+    const int irho = imap[fluid_prim::density::name()].first;
+    const int iye = imap[fluid_prim::ye::name()].first;
+    const int itemp = imap[fluid_prim::temperature::name()].first;
+    const int is = imap[fluid_prim::entropy::name()].first;
+
+    IndexRange ib = rc->GetBoundsI(IndexDomain::interior);
+    IndexRange jb = rc->GetBoundsJ(IndexDomain::interior);
+    IndexRange kb = rc->GetBoundsK(IndexDomain::interior);
+
+    auto eos = pmb->packages.Get("eos")->Param<Microphysics::EOS::EOS>("d.EOS");
+    singularity::StellarCollapse eos_sc =
+        eos.GetUnmodifiedObject().get<singularity::StellarCollapse>();
+    parthenon::par_for(
+        DEFAULT_LOOP_PATTERN, "CoolingFunctionCalculateFourForce", DevExecSpace(), kb.s,
+        kb.e, jb.s, jb.e, ib.s, ib.e,
+        KOKKOS_LAMBDA(const int k, const int j, const int i) {
+          Real lambda[2];
+          lambda[0] = v(iye, k, j, i);
+          const Real s = eos_sc.EntropyFromDensityTemperature(v(irho, k, j, i),
+                                                              v(itemp, k, j, i), lambda);
+          v(is, k, j, i) = s;
+        });
+  }
   if (do_tracers) {
     auto &mbd = pmb->meshblock_data.Get();
     tracers::FillTracers(mbd.get());
diff --git a/src/phoebus_utils/variables.hpp b/src/phoebus_utils/variables.hpp
index d510af9e..a6dde283 100644
--- a/src/phoebus_utils/variables.hpp
+++ b/src/phoebus_utils/variables.hpp
@@ -49,6 +49,7 @@ VARIABLE(p, bfield);
 VARIABLE(p, ye);
 VARIABLE_NONS(pressure);
 VARIABLE_NONS(temperature);
+VARIABLE_NONS(entropy);
 VARIABLE_NONS(gamma1);
 } // namespace fluid_prim