Maintain sanity (or madness)

inputs/torus.pin

@@ -25,9 +25,9 @@ variables = p.density,       &
             c.energy,        &
             pressure,      	 &
             g.c.coord,       &
-            g.f1.coord,       &
-            g.f2.coord,       &
-            g.f3.coord,       &
+            g.f1.coord,      &
+            g.f2.coord,      &
+            g.f3.coord,      &
             g.n.coord,       &
             g.c.detgam,      &
             g.c.bcon,        &
@@ -40,6 +40,10 @@ variables = p.density,       &
 file_type   = hdf5      # Tabular data dump
 dt          = 5.0       # time increment between outputs
 
+<parthenon/output2>
+file_type   = hst       # History data dump
+dt          = 1.0       # Time increment between outputs
+
 <parthenon/time>
 nlim        = -1        # cycle limit
 tlim        = 2000.0    # time limit
@@ -136,6 +140,10 @@ ceiling_type = ConstantGamSie
 sie0_ceiling = 1.0
 gam0_ceiling = 10.0
 enable_ox1_fmks_inflow_check = true
+enable_phi_enforcement = false
+enforced_phi = 1.
+enforced_phi_timescale = 500.
+enforced_phi_cadence = 1.
 
 <torus>
 target_beta = 100.

scripts/python/phoedf.py

@@ -75,7 +75,7 @@ def flatten_indices(mu, nu):
     self.gcov = np.zeros([self.NumBlocks, 4, 4, self.Nx3, self.Nx2, self.Nx1])
     for mu in range(4):
       for nu in range(4):
-        self.gcov[:,mu,nu,:,:,:] = self.flatgcov[:,flatten_indices(mu,nu),:,:,:]
+        self.gcov[:,mu,nu,:,:,:] = self.flatgcov[:,flatten_indices(mu,nu),:,4:-4,4:-4]
     del(self.flatgcov)
     self.gcon = np.zeros([self.NumBlocks, 4, 4, self.Nx3, self.Nx2, self.Nx1])
     for b in range(self.NumBlocks):

src/fixup/fixup.cpp

@@ -18,12 +18,14 @@
 #include <bvals/bvals_interfaces.hpp>
 #include <defs.hpp>
 
+#include "analysis/history.hpp"
 #include "fluid/con2prim_robust.hpp"
 #include "fluid/fluid.hpp"
 #include "fluid/prim2con.hpp"
 #include "geometry/geometry.hpp"
 #include "microphysics/eos_phoebus/eos_phoebus.hpp"
 #include "phoebus_utils/programming_utils.hpp"
+#include "phoebus_utils/reduction.hpp"
 #include "phoebus_utils/relativity_utils.hpp"
 #include "phoebus_utils/robust.hpp"
 #include "phoebus_utils/variables.hpp"
@@ -71,6 +73,23 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   bool enable_mhd_floors = pin->GetOrAddBoolean("fixup", "enable_mhd_floors", false);
   bool enable_rad_floors = pin->GetOrAddBoolean("fixup", "enable_rad_floors", false);
   bool enable_rad_ceilings = pin->GetOrAddBoolean("fixup", "enable_rad_ceilings", false);
+  bool enable_phi_enforcement =
+      pin->GetOrAddBoolean("fixup", "enable_phi_enforcement", false);
+  params.Add("enable_phi_enforcement", enable_phi_enforcement);
+  Real enforced_phi = pin->GetOrAddReal("fixup", "enforced_phi", 0.);
+  params.Add("enforced_phi", enforced_phi);
+  Real enforced_phi_timescale =
+      pin->GetOrAddReal("fixup", "enforced_phi_timescale", 1.e3);
+  params.Add("enforced_phi_timescale", enforced_phi_timescale);
+  Real enforced_phi_cadence = pin->GetOrAddReal("fixup", "enforced_phi_cadence", 10.);
+  params.Add("enforced_phi_cadence", enforced_phi_cadence);
+  Real enforced_phi_start_time =
+      pin->GetOrAddReal("fixup", "enforced_phi_start_time", 175.);
+  params.Add("enforced_phi_start_time", enforced_phi_start_time);
+  if (enable_phi_enforcement) {
+    PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::FMKS),
+                      "Phi enforcement only supported for BH geometry!");
+  }
 
   if (enable_mhd_floors && !enable_mhd) {
     enable_mhd_floors = false;
@@ -645,6 +664,156 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {
   return TaskStatus::complete;
 }
 
+static Real dphi_dt = 0.;
+static Real next_dphi_dt_update_time = 0.;
+static Real phi_factor = 0.;
+
+TaskStatus EndOfStepModify(MeshData<Real> *md, const Real t, const Real dt,
+                           bool last_stage) {
+  auto *pm = md->GetParentPointer();
+  StateDescriptor *fix_pkg = pm->packages.Get("fixup").get();
+
+  const bool enable_phi_enforcement = fix_pkg->Param<bool>("enable_phi_enforcement");
+
+  namespace p = fluid_prim;
+  namespace c = fluid_cons;
+  const std::vector<std::string> vars({p::bfield, c::bfield});
+  PackIndexMap imap;
+  auto pack = md->PackVariables(vars, imap);
+
+  const int pblo = imap[p::bfield].first;
+  const int cblo = imap[c::bfield].first;
+
+  const auto ib = md->GetBoundsI(IndexDomain::interior);
+  const auto jb = md->GetBoundsJ(IndexDomain::interior);
+  const auto kb = md->GetBoundsK(IndexDomain::interior);
+
+  auto geom = Geometry::GetCoordinateSystem(md);
+  auto gpkg = pm->packages.Get("geometry");
+  bool derefine_poles = gpkg->Param<bool>("derefine_poles");
+  Real h = gpkg->Param<Real>("h");
+  Real xt = gpkg->Param<Real>("xt");
+  Real alpha = gpkg->Param<Real>("alpha");
+  Real x0 = gpkg->Param<Real>("x0");
+  Real smooth = gpkg->Param<Real>("smooth");
+  auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth);
+
+  if (enable_phi_enforcement) {
+    const Real enforced_phi = fix_pkg->Param<Real>("enforced_phi");
+    const Real enforced_phi_timescale = fix_pkg->Param<Real>("enforced_phi_timescale");
+    const Real enforced_phi_cadence = fix_pkg->Param<Real>("enforced_phi_cadence");
+    const Real enforced_phi_start_time = fix_pkg->Param<Real>("enforced_phi_start_time");
+    if (t > enforced_phi_start_time) {
+
+      // This only needs to be done at initialization
+      ParArrayND<Real> A("vector potential", pack.GetDim(5), jb.e + 2, ib.e + 2);
+      // Update B field
+      parthenon::par_for(
+          parthenon::LoopPatternMDRange(), "Phoebus::Fixup::EndOfStepModify::EvaluateQ",
+          DevExecSpace(), 0, pack.GetDim(5) - 1, jb.s, jb.e + 1, ib.s, ib.e + 1,
+          KOKKOS_LAMBDA(const int b, const int j, const int i) {
+            const auto &coords = pack.GetCoords(b);
+            const Real x1 = coords.Xf<1>(0, j, i);
+            const Real x2 = coords.Xf<2>(0, j, i);
+            const Real r = tr.bl_radius(x1);
+            const Real th = tr.bl_theta(x1, x2);
+
+            const Real x = r * sin(th);
+            const Real z = r * cos(th);
+            const Real a_hyp = 1.2;
+            const Real b_hyp = 3. * a_hyp;
+            const Real x_hyp = a_hyp * sqrt(1. + pow(z / b_hyp, 2.));
+            const Real q = (pow(x, 2) - pow(x_hyp, 2)) / pow(x_hyp, 2.);
+            if (x < x_hyp) {
+              A(b, j, i) = q;
+            }
+          });
+
+      // Recalculate dphi_dt occasionally
+      if (last_stage && t >= next_dphi_dt_update_time) {
+        next_dphi_dt_update_time += enforced_phi_cadence;
+
+        // Record mdot
+        const Real mdot_proc = History::ReduceMassAccretionRate(md);
+        const Real mdot = reduction::Sum(mdot_proc);
+
+        // Measure phi
+        const Real phi_proc = History::ReduceMagneticFluxPhi(md);
+        const Real phi = reduction::Sum(phi_proc) / std::sqrt(mdot);
+
+        // Calculate amount of phi to add this timestep
+        const Real dphi = (enforced_phi - phi) * dt / enforced_phi_timescale;
+        dphi_dt = dphi / dt;
+
+        parthenon::par_for(
+            parthenon::LoopPatternMDRange(),
+            "Phoebus::Fixup::EndOfStepModify::EvaluateBField", DevExecSpace(), 0,
+            pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e - 1, ib.s, ib.e - 1,
+            KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+              const auto &coords = pack.GetCoords(b);
+              const Real gammadet = geom.DetGamma(CellLocation::Cent, b, k, j, i);
+              pack(b, pblo, k, j, i) +=
+                  (-(A(b, j, i) - A(b, j + 1, i) + A(b, j, i + 1) - A(b, j + 1, i + 1)) /
+                   (2.0 * coords.CellWidthFA(X2DIR, k, j, i) * gammadet));
+              pack(b, pblo + 1, k, j, i) +=
+                  ((A(b, j, i) + A(b, j + 1, i) - A(b, j, i + 1) - A(b, j + 1, i + 1)) /
+                   (2.0 * coords.CellWidthFA(X1DIR, k, j, i) * gammadet));
+
+              // Used by ReduceMagneticFluxPhi
+              pack(b, cblo, k, j, i) = pack(b, pblo, k, j, i) * gammadet;
+              pack(b, cblo + 1, k, j, i) = pack(b, pblo + 1, k, j, i) * gammadet;
+            });
+
+        // Measure change in phi
+        const Real fiducial_phi_proc = History::ReduceMagneticFluxPhi(md);
+        const Real fiducial_phi = reduction::Sum(fiducial_phi_proc) / std::sqrt(mdot);
+        phi_factor = dphi / (fiducial_phi - phi) / dt;
+
+        // Remove the fiducial contribution
+        parthenon::par_for(
+            parthenon::LoopPatternMDRange(),
+            "Phoebus::Fixup::EndOfStepModify::EvaluateBField", DevExecSpace(), 0,
+            pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e - 1, ib.s, ib.e - 1,
+            KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+              const auto &coords = pack.GetCoords(b);
+              const Real gammadet = geom.DetGamma(CellLocation::Cent, b, k, j, i);
+
+              pack(b, pblo, k, j, i) -=
+                  (-(A(b, j, i) - A(b, j + 1, i) + A(b, j, i + 1) - A(b, j + 1, i + 1)) /
+                   (2.0 * coords.CellWidthFA(X2DIR, k, j, i) * gammadet));
+              pack(b, pblo + 1, k, j, i) -=
+                  ((A(b, j, i) + A(b, j + 1, i) - A(b, j, i + 1) - A(b, j + 1, i + 1)) /
+                   (2.0 * coords.CellWidthFA(X1DIR, k, j, i) * gammadet));
+
+              pack(b, cblo, k, j, i) = pack(b, pblo, k, j, i) * gammadet;
+              pack(b, cblo + 1, k, j, i) = pack(b, pblo + 1, k, j, i) * gammadet;
+            });
+      }
+
+      // Properly update B field
+      parthenon::par_for(
+          parthenon::LoopPatternMDRange(),
+          "Phoebus::Fixup::EndOfStepModify::EvaluateBField", DevExecSpace(), 0,
+          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) {
+            const auto &coords = pack.GetCoords(b);
+            const Real gammadet = geom.DetGamma(CellLocation::Cent, b, k, j, i);
+
+            pack(b, pblo, k, j, i) +=
+                phi_factor * dt *
+                (-(A(b, j, i) - A(b, j + 1, i) + A(b, j, i + 1) - A(b, j + 1, i + 1)) /
+                 (2.0 * coords.CellWidthFA(X2DIR, k, j, i) * gammadet));
+            pack(b, pblo + 1, k, j, i) +=
+                phi_factor * dt *
+                ((A(b, j, i) + A(b, j + 1, i) - A(b, j, i + 1) - A(b, j + 1, i + 1)) /
+                 (2.0 * coords.CellWidthFA(X1DIR, k, j, i) * gammadet));
+          });
+    } // enforced_phi_start_time
+  }   // enable_phi_enforcement
+
+  return TaskStatus::complete;
+}
+
 template <typename T>
 TaskStatus ApplyFloors(T *rc) {
   auto *pm = rc->GetParentPointer().get();

src/fixup/fixup.hpp

@@ -38,6 +38,8 @@ template <typename T>
 TaskStatus ConservedToPrimitiveFixup(T *rc);
 template <typename T>
 TaskStatus SourceFixup(T *rc);
+TaskStatus EndOfStepModify(MeshData<Real> *, const Real t, const Real dt,
+                           const bool last_stage);
 
 static struct ConstantRhoSieFloor {
 } constant_rho_sie_floor_tag;

src/phoebus_driver.cpp

@@ -180,6 +180,7 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) {
   BlockList_t &blocks = pmesh->block_list;
 
   const Real beta = integrator->beta[stage - 1];
+  const Real t = tm.time;
   const Real dt = integrator->dt;
   const auto &stage_name = integrator->stage_name;
 
@@ -351,6 +352,28 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) {
     }
   }
 
+  // Extra per-step user work
+  TaskRegion &sync_region_5 = tc.AddRegion(num_partitions);
+  PARTHENON_REQUIRE(num_partitions == 1,
+                    "Reductions don't work for multiple partitions?");
+  for (int i = 0; i < num_partitions; i++) {
+    auto &tl = sync_region_5[i];
+
+    auto &md = pmesh->mesh_data.GetOrAdd(stage_name[stage], i);
+
+    auto end_mods = tl.AddTask(none, fixup::EndOfStepModify, md.get(), t, dt,
+                               stage == integrator->nstages);
+  }
+  // This is a bad pattern. Having a per-mesh data p2c would help.
+  TaskRegion &async_region_4 = tc.AddRegion(num_independent_task_lists);
+  for (int i = 0; i < blocks.size(); i++) {
+    auto &pmb = blocks[i];
+    auto &tl = async_region_4[i];
+    auto &sc = pmb->meshblock_data.Get(stage_name[stage]);
+
+    auto p2c = tl.AddTask(none, fluid::PrimitiveToConserved, sc.get());
+  }
+
   TaskRegion &sync_region_2 = tc.AddRegion(num_partitions);
   for (int ib = 0; ib < num_partitions; ib++) {
     auto &base = pmesh->mesh_data.GetOrAdd("base", ib);

src/phoebus_utils/reduction.hpp

@@ -1,4 +1,4 @@
-// © 2022. Triad National Security, LLC. All rights reserved.  This
+// © 2022-2023. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract
 // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
 // is operated by Triad National Security, LLC for the U.S.
@@ -32,10 +32,17 @@ Real inline Min(const Real &x) {
   return xmin;
 }
 
+Real inline Sum(const Real &x) {
+  Real xsum;
+  MPI_Allreduce(&x, &xsum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  return xsum;
+}
+
 #else
 
 Real inline Max(const Real &x) { return x; }
 Real inline Min(const Real &x) { return x; }
+Real inline Sum(const Real &x) { return x; }
 
 #endif // MPI_PARALLEL