Log edotcool

parthenon-hpc-lab · Nov 17, 2023 · a8a81e1 · a8a81e1
1 parent dfd6650
commit a8a81e1
Show file tree

Hide file tree

Showing 4 changed files with 57 additions and 14 deletions.
diff --git a/src/hydro/hydro_driver.cpp b/src/hydro/hydro_driver.cpp
@@ -230,6 +230,11 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) {
 #endif
   }
 
+  // Reset per cycle cooling logger
+  if (stage == 1 && hydro_pkg->AllParams().hasKey("cooling/total_deltaE_this_cycle")) {
+    hydro_pkg->UpdateParam("cooling/total_deltaE_this_cycle", 0.0);
+  }
+
   // First add split sources before the main time integration
   if (stage == 1) {
     TaskRegion &strang_init_region = tc.AddRegion(num_partitions);

diff --git a/src/hydro/srcterms/tabular_cooling.cpp b/src/hydro/srcterms/tabular_cooling.cpp
@@ -22,6 +22,7 @@
 // AthenaPK headers
 #include "../../units.hpp"
 #include "tabular_cooling.hpp"
+
 #include "utils/error_checking.hpp"
 
 namespace cooling {
@@ -264,6 +265,9 @@ TabularCooling::TabularCooling(ParameterInput *pin,
   cooling_table_obj_ = CoolingTableObj(log_lambdas_, log_temp_start_, log_temp_final_,
                                        d_log_temp_, n_temp_, mbar_over_kb,
                                        adiabatic_index, 1.0 - He_mass_fraction, units);
+
+  // Add mutable param to keep track of total energy lost
+  hydro_pkg->AddParam("cooling/total_deltaE_this_cycle", 0.0, true);
 }
 
 void TabularCooling::SrcTerm(MeshData<Real> *md, const Real dt) const {
@@ -313,12 +317,15 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData<Real> *md, const Real dt
   IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire);
   IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire);
 
-  par_for(
-      DEFAULT_LOOP_PATTERN, "TabularCooling::SubcyclingSplitSrcTerm", DevExecSpace(), 0,
-      cons_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 total_deltaE = NAN; // Total radiated energy (not a density)
+
+  md->GetBlockData(0)->GetBlockPointer()->par_reduce(
+      "TabularCooling::SubcyclingSplitSrcTerm", 0, cons_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 &ledot) {
         auto &cons = cons_pack(b);
         auto &prim = prim_pack(b);
+        const auto &coords = cons_pack.GetCoords(b);
         // Need to use `cons` here as prim may still contain state at t_0;
         const Real rho = cons(IDN, k, j, i);
         // TODO(pgrete) with potentially more EOS, a separate get_pressure (or similar)
@@ -470,11 +477,20 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData<Real> *md, const Real dt
         internal_e = (internal_e > internal_e_floor) ? internal_e : internal_e_floor;
 
         // Remove the cooling from the total energy density
-        cons(IEN, k, j, i) += rho * (internal_e - internal_e_initial);
+        const auto edotdensity = rho * (internal_e - internal_e_initial);
+        cons(IEN, k, j, i) += edotdensity;
+        ledot = edotdensity * coords.CellVolume(k, j, i);
         // Latter technically not required if no other tasks follows before
         // ConservedToPrim conversion, but keeping it for now (better safe than sorry).
         prim(IPR, k, j, i) = rho * internal_e * gm1;
-      });
+      },
+      Kokkos::Sum<Real>(total_deltaE));
+
+  // Updating current value as routine might be called multiple times (e.g. Strang split)
+  auto total_deltaE_this_cycle =
+      hydro_pkg->Param<Real>("cooling/total_deltaE_this_cycle");
+  hydro_pkg->UpdateParam("cooling/total_deltaE_this_cycle",
+                         total_deltaE_this_cycle + total_deltaE);
 }
 
 void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,
@@ -514,12 +530,15 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,
   const auto temp_final = std::pow(10.0, log_temp_final_);
   const auto lambda_final = lambda_final_;
 
-  par_for(
-      DEFAULT_LOOP_PATTERN, "TabularCooling::TownsendSrcTerm", DevExecSpace(), 0,
-      cons_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 total_deltaE = NAN; // Total radiated energy (not a density)
+
+  md->GetBlockData(0)->GetBlockPointer()->par_reduce(
+      "TabularCooling::TownsendSrcTerm", 0, cons_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 &ledot) {
         auto &cons = cons_pack(b);
         auto &prim = prim_pack(b);
+        const auto &coords = cons_pack.GetCoords(b);
         // Need to use `cons` here as prim may still contain state at t_0;
         const auto rho = cons(IDN, k, j, i);
         // TODO(pgrete) with potentially more EOS, a separate get_pressure (or similar)
@@ -587,11 +606,20 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,
         const auto internal_e_new = temp_new > temp_cool_floor
                                         ? temp_new / mbar_gm1_over_kb
                                         : temp_cool_floor / mbar_gm1_over_kb;
-        cons(IEN, k, j, i) += rho * (internal_e_new - internal_e);
+        const auto edotdensity = rho * (internal_e_new - internal_e);
+        cons(IEN, k, j, i) += edotdensity;
+        ledot = edotdensity * coords.CellVolume(k, j, i);
         // Latter technically not required if no other tasks follows before
         // ConservedToPrim conversion, but keeping it for now (better safe than sorry).
         prim(IPR, k, j, i) = rho * internal_e_new * gm1;
-      });
+      },
+      Kokkos::Sum<Real>(total_deltaE));
+
+  // Updating current value as routine might be called multiple times (e.g. Strang split)
+  auto total_deltaE_this_cycle =
+      hydro_pkg->Param<Real>("cooling/total_deltaE_this_cycle");
+  hydro_pkg->UpdateParam("cooling/total_deltaE_this_cycle",
+                         total_deltaE_this_cycle + total_deltaE);
 }
 
 Real TabularCooling::EstimateTimeStep(MeshData<Real> *md) const {

diff --git a/src/outputs/per_block.cpp b/src/outputs/per_block.cpp
@@ -410,7 +410,6 @@ void UserOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm,
     stats.emplace_back(Stats("vel_mag", "prim", {IV1, IV2, IV3}));
     stats.emplace_back(Stats("vel_mag_mw", "prim", {IV1, IV2, IV3}, Weight::Mass));
     stats.emplace_back(Stats("rho", "prim", 0));
-    stats.emplace_back(Stats("vx", "prim", 1));
 
     const std::vector<std::string> stat_types = {
         "min", "max", "absmin", "absmax", "mean", "rms", "stddev", "skew", "kurt"};

diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp
@@ -46,7 +46,7 @@ using utils::few_modes_ft::FewModesFT;
 
 // TODO(?) until we are able to process multiple variables in a single hst function call
 // we'll use this enum to identify the various vars.
-enum class HstQuan { Ms, Ma, pb };
+enum class HstQuan { Ms, Ma, pb, DeltaEcool };
 
 // Compute the local sum of either the sonic Mach number,
 // alfvenic Mach number, or plasma beta as specified by `hst_quan`.
@@ -55,6 +55,14 @@ Real TurbulenceHst(MeshData<Real> *md) {
   auto pmb = md->GetBlockData(0)->GetBlockPointer();
   auto hydro_pkg = pmb->packages.Get("Hydro");
   const auto gamma = hydro_pkg->Param<Real>("AdiabaticIndex");
+
+  if (hst_quan == HstQuan::DeltaEcool &&
+      hydro_pkg->AllParams().hasKey("cooling/total_deltaE_this_cycle")) {
+    return hydro_pkg->Param<Real>("cooling/total_deltaE_this_cycle");
+  } else {
+    return 0;
+  }
+
   const auto fluid = hydro_pkg->Param<Fluid>("fluid");
 
   const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
@@ -112,6 +120,9 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
 
   hst_vars.emplace_back(parthenon::HistoryOutputVar(parthenon::UserHistoryOperation::sum,
                                                     TurbulenceHst<HstQuan::Ms>, "Ms"));
+  hst_vars.emplace_back(parthenon::HistoryOutputVar(parthenon::UserHistoryOperation::sum,
+                                                    TurbulenceHst<HstQuan::DeltaEcool>,
+                                                    "DeltaEcool"));
   if (fluid == Fluid::glmmhd) {
     hst_vars.emplace_back(parthenon::HistoryOutputVar(
         parthenon::UserHistoryOperation::sum, TurbulenceHst<HstQuan::Ma>, "Ma"));