Enable floor reporting

Author: pgrete

src/eos/adiabatic_hydro.cpp

@@ -43,13 +43,29 @@ void AdiabaticHydroEOS::ConservedToPrimitive(MeshData<Real> *md) const {
 
   auto this_on_device = (*this);
 
-  parthenon::par_for(
+  std::int64_t floor_rho, floor_pres, floor_temp;
+  parthenon::par_reduce(
       DEFAULT_LOOP_PATTERN, "ConservedToPrimitive", parthenon::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) {
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i,
+                    std::int64_t &lfloor_rho, std::int64_t &lfloor_pres,
+                    std::int64_t &lfloor_temp) {
         const auto &cons = cons_pack(b);
         auto &prim = prim_pack(b);
 
-        return this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
-      });
+        auto floors_used =
+            this_on_device.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i);
+        if (floors_used & 1) lfloor_rho += 1;
+        if (floors_used & 2) lfloor_pres += 1;
+        if (floors_used & 4) lfloor_temp += 1;
+      },
+      floor_rho, floor_pres, floor_temp);
+  const auto floor_rho_pkg = pkg->Param<std::int64_t>("fixed_num_cells_floor_rho");
+  pkg->UpdateParam<std::int64_t>("fixed_num_cells_floor_rho", floor_rho_pkg + floor_rho);
+  const auto floor_pres_pkg = pkg->Param<std::int64_t>("fixed_num_cells_floor_pres");
+  pkg->UpdateParam<std::int64_t>("fixed_num_cells_floor_pres",
+                                 floor_pres_pkg + floor_pres);
+  const auto floor_temp_pkg = pkg->Param<std::int64_t>("fixed_num_cells_floor_temp");
+  pkg->UpdateParam<std::int64_t>("fixed_num_cells_floor_temp",
+                                 floor_temp_pkg + floor_temp);
 }

src/eos/adiabatic_hydro.hpp

@@ -49,9 +49,10 @@ class AdiabaticHydroEOS : public EquationOfState {
   // int& j, const int& i) \brief Fills an array of primitives given an array of
   // conserveds, potentially updating the conserved with floors
   template <typename View4D>
-  KOKKOS_INLINE_FUNCTION void ConsToPrim(View4D cons, View4D prim, const int &nhydro,
-                                         const int &nscalars, const int &k, const int &j,
-                                         const int &i) const {
+  KOKKOS_INLINE_FUNCTION int ConsToPrim(View4D cons, View4D prim, const int &nhydro,
+                                        const int &nscalars, const int &k, const int &j,
+                                        const int &i) const {
+    int floors_used = 0;
     Real gm1 = GetGamma() - 1.0;
     auto density_floor_ = GetDensityFloor();
     auto pressure_floor_ = GetPressureFloor();
@@ -72,13 +73,19 @@ class AdiabaticHydroEOS : public EquationOfState {
     Real &w_vz = prim(IV3, k, j, i);
     Real &w_p = prim(IPR, k, j, i);
 
+    PARTHENON_REQUIRE(u_d != 0.0,
+                      "Densities should never be exactly 0! This points to working with "
+                      "some default initialized and/or uninitialized data.");
     // Let's apply floors explicitly, i.e., by default floor will be disabled (<=0)
     // and the code will fail if a negative density is encountered.
     PARTHENON_REQUIRE(u_d > 0.0 || density_floor_ > 0.0,
                       "Got negative density. Consider enabling first-order flux "
                       "correction or setting a reasonble density floor.");
     // apply density floor, without changing momentum or energy
-    u_d = (u_d > density_floor_) ? u_d : density_floor_;
+    if (u_d < density_floor_) {
+      u_d = density_floor_;
+      floors_used = floors_used | 1; // set density floor flag
+    }
     w_d = u_d;
 
     Real di = 1.0 / u_d;
@@ -117,6 +124,7 @@ class AdiabaticHydroEOS : public EquationOfState {
       // apply pressure floor, correct total energy
       u_e = (pressure_floor_ / gm1) + e_k;
       w_p = pressure_floor_;
+      floors_used = floors_used | 2; // set pressure floor flag
     }
 
     // temperature (internal energy) based pressure floor
@@ -125,6 +133,7 @@ class AdiabaticHydroEOS : public EquationOfState {
       // apply temperature floor, correct total energy
       u_e = (u_d * e_floor_) + e_k;
       w_p = eff_pressure_floor;
+      floors_used = floors_used | 4; // set temperture floor flag
     }
 
     // temperature (internal energy) based pressure ceiling
@@ -139,6 +148,7 @@ class AdiabaticHydroEOS : public EquationOfState {
     for (auto n = nhydro; n < nhydro + nscalars; ++n) {
       prim(n, k, j, i) = cons(n, k, j, i) * di;
     }
+    return floors_used;
   }
 
  private:

src/hydro/hydro.cpp

@@ -522,6 +522,9 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
       auto mbar_over_kb = pkg->Param<Real>("mbar_over_kb");
       efloor = Tfloor / mbar_over_kb / (gamma - 1.0);
     }
+    pkg->AddParam<>("dfloor", dfloor);
+    pkg->AddParam<>("pfloor", pfloor);
+    pkg->AddParam<>("Tfloor", Tfloor);
 
     // By default disable ceilings by setting to infinity
     Real vceil =

src/hydro/hydro_driver.cpp

@@ -591,10 +591,25 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) {
           if (hydro_pkg->Param<bool>("first_order_flux_correct")) {
             msg << hydro_pkg->Param<std::int64_t>("fixed_num_cells_fofc") << " FOFC. ";
           }
+          if (hydro_pkg->Param<Real>("dfloor") > 0.0) {
+            msg << hydro_pkg->Param<std::int64_t>("fixed_num_cells_floor_rho")
+                << " dfloor. ";
+          }
+          if (hydro_pkg->Param<Real>("pfloor") > 0.0) {
+            msg << hydro_pkg->Param<std::int64_t>("fixed_num_cells_floor_pres")
+                << " pfloor. ";
+          }
+          if (hydro_pkg->Param<Real>("Tfloor") > 0.0) {
+            msg << hydro_pkg->Param<std::int64_t>("fixed_num_cells_floor_temp")
+                << " Tfloor. ";
+          }
           std::cout << msg.str() << "\n";
         }
 
         // reset counter for next stage
+        hydro_pkg->UpdateParam<std::int64_t>("fixed_num_cells_floor_rho", 0);
+        hydro_pkg->UpdateParam<std::int64_t>("fixed_num_cells_floor_pres", 0);
+        hydro_pkg->UpdateParam<std::int64_t>("fixed_num_cells_floor_temp", 0);
         hydro_pkg->UpdateParam<std::int64_t>("fixed_num_cells_fofc", 0);
         return TaskStatus::complete;
       },