Print the wavespeed that sets the CFL timestep

src/hydro/hydro.cpp

@@ -7,6 +7,9 @@
 #include <algorithm>
 #include <limits>
 #include <memory>
+#ifdef MPI_PARALLEL
+#include <mpi.h>
+#endif
 #include <string>
 #include <vector>
 
@@ -26,9 +29,11 @@
 #include "../recon/wenoz_simple.hpp"
 #include "../refinement/refinement.hpp"
 #include "../units.hpp"
+#include "basic_types.hpp"
 #include "defs.hpp"
 #include "diffusion/diffusion.hpp"
 #include "glmmhd/glmmhd.hpp"
+#include "globals.hpp"
 #include "hydro.hpp"
 #include "interface/params.hpp"
 #include "outputs/outputs.hpp"
@@ -186,6 +191,9 @@ TaskStatus AddSplitSourcesStrang(MeshData<Real> *md, const SimTime &tm) {
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto pkg = std::make_shared<StateDescriptor>("Hydro");
 
+  CellPrimValues cell_values{};
+  pkg->AddParam<CellPrimValues>("cfl_cell_properties", cell_values, true);
+
   Real cfl = pin->GetOrAddReal("parthenon/time", "cfl", 0.3);
   pkg->AddParam<>("cfl", cfl);
 
@@ -751,7 +759,10 @@ Real EstimateHyperbolicTimestep(MeshData<Real> *md) {
   IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior);
   IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior);
 
-  Real min_dt_hyperbolic = std::numeric_limits<Real>::max();
+  // reduce ValPropPair<Real, CellPrimValues> so we can obtain the properties of
+  // the cell that is limiting the timestep
+  ValPropPair<Real, CellPrimValues> min_dt_hyperbolic{std::numeric_limits<Real>::max(),
+                                                      {}};
 
   const auto ndim_ = prim_pack.GetNdim();
   Kokkos::parallel_reduce(
@@ -760,7 +771,8 @@ Real EstimateHyperbolicTimestep(MeshData<Real> *md) {
           DevExecSpace(), {0, kb.s, jb.s, ib.s},
           {prim_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1},
           {1, 1, 1, ib.e + 1 - ib.s}),
-      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &min_dt) {
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i,
+                    ValPropPair<Real, CellPrimValues> &min_dt) {
         const auto &prim = prim_pack(b);
         const auto &coords = prim_pack.GetCoords(b);
         // Need to reference variables here so that they are properly caught by
@@ -774,37 +786,114 @@ Real EstimateHyperbolicTimestep(MeshData<Real> *md) {
         w[IV2] = prim(IV2, k, j, i);
         w[IV3] = prim(IV3, k, j, i);
         w[IPR] = prim(IPR, k, j, i);
-        Real lambda_max_x, lambda_max_y, lambda_max_z;
+
+        const Real B1 = prim(IB1, k, j, i);
+        const Real B2 = prim(IB2, k, j, i);
+        const Real B3 = prim(IB3, k, j, i);
+
+        Real lambda_max_x{};
+        Real lambda_max_y{};
+        Real lambda_max_z{};
+
         if constexpr (fluid == Fluid::euler) {
           lambda_max_x = eos.SoundSpeed(w);
           lambda_max_y = lambda_max_x;
           lambda_max_z = lambda_max_x;
 
         } else if constexpr (fluid == Fluid::glmmhd) {
-          lambda_max_x = eos.FastMagnetosonicSpeed(
-              w[IDN], w[IPR], prim(IB1, k, j, i), prim(IB2, k, j, i), prim(IB3, k, j, i));
+          lambda_max_x = eos.FastMagnetosonicSpeed(w[IDN], w[IPR], B1, B2, B3);
           if (ndim > 1) {
-            lambda_max_y =
-                eos.FastMagnetosonicSpeed(w[IDN], w[IPR], prim(IB2, k, j, i),
-                                          prim(IB3, k, j, i), prim(IB1, k, j, i));
+            lambda_max_y = eos.FastMagnetosonicSpeed(w[IDN], w[IPR], B2, B3, B1);
           }
           if (ndim > 2) {
-            lambda_max_z =
-                eos.FastMagnetosonicSpeed(w[IDN], w[IPR], prim(IB3, k, j, i),
-                                          prim(IB1, k, j, i), prim(IB2, k, j, i));
+            lambda_max_z = eos.FastMagnetosonicSpeed(w[IDN], w[IPR], B3, B1, B2);
           }
         } else {
           PARTHENON_FAIL("Unknown fluid in EstimateTimestep");
         }
-        min_dt = fmin(min_dt, coords.Dxc<1>(k, j, i) / (fabs(w[IV1]) + lambda_max_x));
+        min_dt.value =
+            fmin(min_dt.value, coords.Dxc<1>(k, j, i) / (fabs(w[IV1]) + lambda_max_x));
         if (ndim > 1) {
-          min_dt = fmin(min_dt, coords.Dxc<2>(k, j, i) / (fabs(w[IV2]) + lambda_max_y));
+          min_dt.value =
+              fmin(min_dt.value, coords.Dxc<2>(k, j, i) / (fabs(w[IV2]) + lambda_max_y));
         }
         if (ndim > 2) {
-          min_dt = fmin(min_dt, coords.Dxc<3>(k, j, i) / (fabs(w[IV3]) + lambda_max_z));
+          min_dt.value =
+              fmin(min_dt.value, coords.Dxc<3>(k, j, i) / (fabs(w[IV3]) + lambda_max_z));
         }
+
+        CellPrimValues this_cell{w[IDN], w[IV1], w[IV2], w[IV3], w[IPR], B1, B2, B3};
+        min_dt.index = this_cell;
       },
-      Kokkos::Min<Real>(min_dt_hyperbolic));
+      Kokkos::Min<ValPropPair<Real, CellPrimValues>>(min_dt_hyperbolic));
+
+  // Locally, min_dt_hyperbolic.value now contains the min timestep and
+  // min_dt_hyperbolic.index contains the primitive vars for the cell that set the min
+  // timestep on this local partition.
+
+#ifdef MPI_PARALLEL
+  MPI_Datatype mpi_valproppair_types[2] = {MPI_PARTHENON_REAL, MPI_BYTE};
+  MPI_Aint mpi_valproppair_disps[2] = {
+      offsetof(valprop_reduce_type, value),
+      offsetof(valprop_reduce_type, index),
+  };
+  int mpi_valproppair_lens[2] = {1, sizeof(CellPrimValues)};
+
+  // create MPI datatype
+  MPI_Datatype mpi_valproppair;
+  MPI_Type_create_struct(2, mpi_valproppair_lens, mpi_valproppair_disps,
+                         mpi_valproppair_types, &mpi_valproppair);
+  MPI_Type_commit(&mpi_valproppair);
+
+  // create MPI reduction op
+  MPI_Op mpi_minloc_valproppair;
+  MPI_Op_create(ValPropPairMPIReducer, 1, &mpi_minloc_valproppair);
+
+  // do MPI reduction
+  MPI_Allreduce(MPI_IN_PLACE, &min_dt_hyperbolic, 1, mpi_valproppair,
+                mpi_minloc_valproppair, MPI_COMM_WORLD);
+#endif // MPI_PARALLEL
+
+  if (parthenon::Globals::my_rank == 0) {
+    // print cell properties
+    CellPrimValues &props = min_dt_hyperbolic.index;
+
+    Real w[(NHYDRO)];
+    w[IDN] = props.rho;
+    w[IV1] = props.v1;
+    w[IV2] = props.v2;
+    w[IV3] = props.v3;
+    w[IPR] = props.P;
+    const Real B1 = props.B1;
+    const Real B2 = props.B2;
+    const Real B3 = props.B3;
+    const Real B_sq = SQR(B1) + SQR(B2) + SQR(B3);
+
+    // print sound speed
+    const parthenon::Real cs = eos_.SoundSpeed(w);
+    const parthenon::Real v_A = std::sqrt(B_sq / props.rho);
+    const parthenon::Real v_abs = std::sqrt(SQR(props.v1) + SQR(props.v2) + SQR(props.v3));
+
+    // print timestep
+    std::cout << "\n----> dt = " << cfl_hyp * min_dt_hyperbolic.value;
+    std::cout << " (c_s = " << cs << ", v_A = " << v_A << ", |v| = " << v_abs << ")\n";
+
+    // find which is largest
+    std::vector<parthenon::Real> v{cs, v_A, v_abs};
+    std::vector<parthenon::Real>::iterator result = std::max_element(v.begin(), v.end());
+    const int max_idx = std::distance(v.begin(), result);
+
+    if (max_idx == 0) {
+      std::cout << "\tTimestep limited by sound speed.\n";
+    } else if (max_idx == 1) {
+      std::cout << "\tTimestep limited by Alfven speed.\n";
+    } else if (max_idx == 2) {
+      std::cout << "\tTimestep limited by fluid velocity.\n";
+    }
+  }
+
+  // save the result
+  hydro_pkg->UpdateParam("cfl_cell_properties", min_dt_hyperbolic.index);
 
   // TODO(pgrete) THIS WORKAROUND IS NOT THREAD SAFE (though this will only become
   // relevant once parthenon uses host-multithreading in the driver).
@@ -814,11 +903,11 @@ Real EstimateHyperbolicTimestep(MeshData<Real> *md) {
   // processes.
   if constexpr (fluid == Fluid::glmmhd) {
     auto dt_hyp_pkg = hydro_pkg->Param<Real>("dt_hyp");
-    if (cfl_hyp * min_dt_hyperbolic < dt_hyp_pkg) {
-      hydro_pkg->UpdateParam("dt_hyp", cfl_hyp * min_dt_hyperbolic);
+    if (cfl_hyp * min_dt_hyperbolic.value < dt_hyp_pkg) {
+      hydro_pkg->UpdateParam("dt_hyp", cfl_hyp * min_dt_hyperbolic.value);
     }
   }
-  return cfl_hyp * min_dt_hyperbolic;
+  return cfl_hyp * min_dt_hyperbolic.value;
 }
 
 // provide the routine that estimates a stable timestep for this package

src/hydro/hydro.hpp

@@ -75,6 +75,61 @@ constexpr size_t GetNVars<Fluid::glmmhd>() {
   return 9; // above plus B_x, B_y, B_z, psi
 }
 
+struct CellPrimValues {
+  Real rho{};
+  Real v1{};
+  Real v2{};
+  Real v3{};
+  Real P{};
+  Real B1{};
+  Real B2{};
+  Real B3{};
+};
+
+template <typename TV, typename TI>
+struct ValPropPair {
+  TV value;
+  TI index;
+
+  static constexpr ValPropPair<TV, TI> max() {
+    return ValPropPair<TV, TI>{std::numeric_limits<TV>::max(), TI()};
+  }
+
+  friend constexpr bool operator<(ValPropPair<TV, TI> const &a,
+                                  ValPropPair<TV, TI> const &b) {
+    return a.value < b.value;
+  }
+
+  friend constexpr bool operator>(ValPropPair<TV, TI> const &a,
+                                  ValPropPair<TV, TI> const &b) {
+    return a.value > b.value;
+  }
+};
+
+typedef ValPropPair<Real, CellPrimValues> valprop_reduce_type;
+
+#ifdef MPI_PARALLEL
+inline void ValPropPairMPIReducer(void *in, void *inout, int *len, MPI_Datatype *type) {
+  valprop_reduce_type *invals = static_cast<valprop_reduce_type *>(in);
+  valprop_reduce_type *inoutvals = static_cast<valprop_reduce_type *>(inout);
+
+  for (int i = 0; i < *len; i++) {
+    if (invals[i] < inoutvals[i]) {
+      inoutvals[i] = invals[i];
+    }
+  }
+};
+#endif // MPI_PARALLEL
+
 } // namespace Hydro
 
+namespace Kokkos { // reduction identity must be defined in Kokkos namespace
+template <typename TV, typename TI>
+struct reduction_identity<Hydro::ValPropPair<TV, TI>> {
+  KOKKOS_FORCEINLINE_FUNCTION static Hydro::ValPropPair<TV, TI> min() {
+    return Hydro::ValPropPair<TV, TI>::max(); // confusingly, this is correct
+  }
+};
+} // namespace Kokkos
+
 #endif // HYDRO_HYDRO_HPP_