Add floor reporting for MHD

Author: pgrete

src/eos/adiabatic_glmmhd.cpp

@@ -43,14 +43,30 @@ void AdiabaticGLMMHDEOS::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);
         // auto &nu = entropy_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_glmmhd.hpp

@@ -63,9 +63,10 @@ class AdiabaticGLMMHDEOS : 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;
     auto gam = GetGamma();
     auto gm1 = gam - 1.0;
     auto density_floor_ = GetDensityFloor();
@@ -95,13 +96,19 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
     Real &w_Bz = prim(IB3, k, j, i);
     Real &w_psi = prim(IPS, 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;
@@ -146,6 +153,7 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
       // apply pressure floor, correct total energy
       u_e = (pressure_floor_ / gm1) + e_k + e_B;
       w_p = pressure_floor_;
+      floors_used = floors_used | 2; // set pressure floor flag
     }
 
     // temperature (internal energy) based pressure floor
@@ -154,6 +162,7 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
       // apply temperature floor, correct total energy
       u_e = (u_d * e_floor_) + e_k + e_B;
       w_p = eff_pressure_floor;
+      floors_used = floors_used | 4; // set temperture floor flag
     }
 
     // temperature (internal energy) based pressure ceiling
@@ -168,6 +177,7 @@ class AdiabaticGLMMHDEOS : 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/eos/adiabatic_hydro.hpp

@@ -187,7 +187,6 @@ class AdiabaticHydroEOS : public EquationOfState {
     // apply density floor, without changing momentum or energy
     u_d = w_d;
 
-    Real di = 1.0 / u_d;
     u_m1 = w_d * w_vx;
     u_m2 = w_d * w_vy;
     u_m3 = w_d * w_vz;