WIP: Test temperature based recon

docs/turbulence.md

@@ -0,0 +1,68 @@
+# Turbulence problem generator
+
+The problem generator has been refactored from its implementation in K-Athena.
+Until the documentation is fully adapted, please consult https://gitlab.com/pgrete/kathena/-/wikis/turbulence for general information.
+
+A sample input file is provided in [turbulence.in](../inputs/turbulence.in)
+
+## Blob injection
+
+In order to study the evolution of (cold) clouds in a turbulent environment,
+blobs can be injected to the simulation.
+The injection is a one off mechanism controlled by the following parameters
+
+```
+<problem/turbulence>
+# only one of the following three conditions can be set at a given time
+inject_once_at_time = -1.0
+inject_once_at_cycle = -1
+inject_once_on_restart = false # should not be set in the input file, but only via the command line upon restart
+
+inject_n_blobs = -1 # number of blob to inject
+
+# then for the given number of blobs follow parameters need to be given (starting to count with 0)
+inject_blob_radius_0 = ... # float, in code length units, no default value
+inject_blob_loc_0 = ...,...,... # location vector of three comma-separated floats, in code length units, no default value
+inject_blob_chi_0 = ... # float, dimensionless, no default value, density ratio to existing value
+
+inject_blob_radius_1 = ...
+...
+```
+
+In practice, this will result in blobs being seeded at a given time, cycle, or upon restart
+by adjusting the density within the blob's radius by a factor of $\chi$ and
+at the same time adjusting the temperature by a factor of $1/\chi$ so that the
+blob remain in pressure equilibrium.
+
+While this is an action that is performed once, it can be repeated upon restart (or a later
+time) by resetting the variables.
+
+A current restriction is that the blobs cannot be seeded across a domain boundary (i.e.,
+the periodicity of the box is not taken into account).
+
+## Rescaling **not recommended*
+
+*The rescaling described in the following is generally not recommended, as it result in a
+state that is not naturally reached.
+Moreover, given the artificial nature of a hard reset, some time after the rescaling is
+required for the system to readjust.
+
+For non-isothermal simulations, the plasma will eventually heat up over time due to dissipation.
+One possibility to remove that extra heat (or add heat), is to rescale the temperature in the simulation.
+This can be done via the following parameters:
+
+```
+<problem/turbulence>
+# only one of the following three conditions can be set at a given time
+rescale_once_at_time = -1.0
+rescale_once_at_cycle = -1
+rescale_once_on_restart = false # should not be set in the input file, but only via the command line upon restart
+
+rescale_to_rms_Ms = -1.0
+```
+
+As the parameters suggest, rescaling is a one off action (though it can be repeated when
+the parameters are set again for a subsequent restart).
+The density and velocity field are not changed, only the specific internal energy is
+adjusted so that the volume-weighted root mean squared sonic Mach number matches
+the target value.

src/CMakeLists.txt

@@ -16,6 +16,7 @@ add_executable(
         hydro/srcterms/gravitational_field.hpp
         hydro/srcterms/tabular_cooling.hpp
         hydro/srcterms/tabular_cooling.cpp
+        #outputs/per_block.cpp
         refinement/gradient.cpp
         refinement/other.cpp
         utils/few_modes_ft.cpp

src/hydro/hydro.cpp

@@ -587,6 +587,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
                       "refinement/maxdensity_refine_above");
     pkg->AddParam<Real>("refinement/maxdensity_deref_below", deref_below);
     pkg->AddParam<Real>("refinement/maxdensity_refine_above", refine_above);
+  } else if (refine_str == "always") {
+    pkg->CheckRefinementBlock = refinement::other::Always;
   } else if (refine_str == "user") {
     pkg->CheckRefinementBlock = Hydro::ProblemCheckRefinementBlock;
   }
@@ -775,8 +777,8 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
   int il, iu, jl, ju, kl, ku;
   jl = jb.s, ju = jb.e, kl = kb.s, ku = kb.e;
   // TODO(pgrete): are these looop limits are likely too large for 2nd order
-  if (pmb->block_size.nx2 > 1) {
-    if (pmb->block_size.nx3 == 1) // 2D
+  if (pmb->block_size.nx(X2DIR) > 1) {
+    if (pmb->block_size.nx(X3DIR) == 1) // 2D
       jl = jb.s - 1, ju = jb.e + 1, kl = kb.s, ku = kb.e;
     else // 3D
       jl = jb.s - 1, ju = jb.e + 1, kl = kb.s - 1, ku = kb.e + 1;
@@ -848,7 +850,7 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
         parthenon::ScratchPad2D<Real>::shmem_size(num_scratch_vars, nx1) * 3;
     // set the loop limits
     il = ib.s - 1, iu = ib.e + 1, kl = kb.s, ku = kb.e;
-    if (pmb->block_size.nx3 == 1) // 2D
+    if (pmb->block_size.nx(X3DIR) == 1) // 2D
       kl = kb.s, ku = kb.e;
     else // 3D
       kl = kb.s - 1, ku = kb.e + 1;

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);

src/hydro/prolongation/custom_ops.hpp

@@ -81,34 +81,40 @@ struct ProlongateCellMinModMultiD {
 
     Real dx1fm = 0;
     Real dx1fp = 0;
+    [[maybe_unused]] Real gx1m = 0, gx1p = 0;
     Real gx1c = 0;
     if constexpr (INCLUDE_X1) {
       Real dx1m, dx1p;
       GetGridSpacings<1, el>(coords, coarse_coords, cib, ib, i, fi, &dx1m, &dx1p, &dx1fm,
                              &dx1fp);
-      gx1c = GradMinMod(fc, coarse(element_idx, l, m, n, k, j, i - 1),
-                        coarse(element_idx, l, m, n, k, j, i + 1), dx1m, dx1p);
+      gx1c =
+          GradMinMod(fc, coarse(element_idx, l, m, n, k, j, i - 1),
+                     coarse(element_idx, l, m, n, k, j, i + 1), dx1m, dx1p, gx1m, gx1p);
     }
 
     Real dx2fm = 0;
     Real dx2fp = 0;
+    [[maybe_unused]] Real gx2m = 0, gx2p = 0;
     Real gx2c = 0;
     if constexpr (INCLUDE_X2) {
       Real dx2m, dx2p;
       GetGridSpacings<2, el>(coords, coarse_coords, cjb, jb, j, fj, &dx2m, &dx2p, &dx2fm,
                              &dx2fp);
-      gx2c = GradMinMod(fc, coarse(element_idx, l, m, n, k, j - 1, i),
-                        coarse(element_idx, l, m, n, k, j + 1, i), dx2m, dx2p);
+      gx2c =
+          GradMinMod(fc, coarse(element_idx, l, m, n, k, j - 1, i),
+                     coarse(element_idx, l, m, n, k, j + 1, i), dx2m, dx2p, gx2m, gx2p);
     }
     Real dx3fm = 0;
     Real dx3fp = 0;
+    [[maybe_unused]] Real gx3m = 0, gx3p = 0;
     Real gx3c = 0;
     if constexpr (INCLUDE_X3) {
       Real dx3m, dx3p;
       GetGridSpacings<3, el>(coords, coarse_coords, ckb, kb, k, fk, &dx3m, &dx3p, &dx3fm,
                              &dx3fp);
-      gx3c = GradMinMod(fc, coarse(element_idx, l, m, n, k - 1, j, i),
-                        coarse(element_idx, l, m, n, k + 1, j, i), dx3m, dx3p);
+      gx3c =
+          GradMinMod(fc, coarse(element_idx, l, m, n, k - 1, j, i),
+                     coarse(element_idx, l, m, n, k + 1, j, i), dx3m, dx3p, gx3m, dx3p);
     }
 
     // Max. expected total difference. (dx#fm/p are positive by construction)

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 {