Remove deltaEcool this cycle due to issues with handling substages

parthenon-hpc-lab · Aug 7, 2024 · f9bc877 · f9bc877
1 parent 47681b3
commit f9bc877
Show file tree

Hide file tree

Showing 3 changed files with 4 additions and 36 deletions.
diff --git a/src/hydro/hydro_driver.cpp b/src/hydro/hydro_driver.cpp
@@ -538,11 +538,6 @@ 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) {
     const auto &diffint = hydro_pkg->Param<DiffInt>("diffint");

diff --git a/src/hydro/srcterms/tabular_cooling.cpp b/src/hydro/srcterms/tabular_cooling.cpp
@@ -265,9 +265,6 @@ 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 {
@@ -319,10 +316,10 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData<Real> *md, const Real dt
 
   Real total_deltaE = NAN; // Total radiated energy (not a density)
 
-  md->GetBlockData(0)->GetBlockPointer()->par_reduce(
+  md->GetBlockData(0)->GetBlockPointer()->par_for(
       "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) {
+      KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) {
         auto &cons = cons_pack(b);
         auto &prim = prim_pack(b);
         const auto &coords = cons_pack.GetCoords(b);
@@ -479,18 +476,10 @@ void TabularCooling::SubcyclingFixedIntSrcTerm(MeshData<Real> *md, const Real dt
         // Remove the cooling from the total energy density
         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,
@@ -614,12 +603,6 @@ void TabularCooling::TownsendSrcTerm(parthenon::MeshData<parthenon::Real> *md,
         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/pgen/turbulence.cpp b/src/pgen/turbulence.cpp
@@ -47,7 +47,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, DeltaEcool, temperature };
+enum class HstQuan { Ms, Ma, pb, temperature };
 
 // Compute the local sum of either the sonic Mach number,
 // alfvenic Mach number, or plasma beta as specified by `hst_quan`.
@@ -57,11 +57,6 @@ Real TurbulenceHst(MeshData<Real> *md) {
   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");
-  }
-
   const auto fluid = hydro_pkg->Param<Fluid>("fluid");
 
   const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
@@ -126,11 +121,6 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
 
   hst_vars.emplace_back(parthenon::HistoryOutputVar(parthenon::UserHistoryOperation::sum,
                                                     TurbulenceHst<HstQuan::Ms>, "Ms"));
-  if (pkg->Param<Cooling>("enable_cooling") == Cooling::tabular) {
-    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"));