diff --git a/inputs/torus.pin b/inputs/torus.pin index 93555091..82dd178f 100644 --- a/inputs/torus.pin +++ b/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 + +file_type = hst # History data dump +dt = 1.0 # Time increment between outputs + 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. target_beta = 100. diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index d7dbd54b..2d627f5c 100644 --- a/scripts/python/phoedf.py +++ b/scripts/python/phoedf.py @@ -25,6 +25,7 @@ from phoebus_eos import * from phoebus_opacities import * + # ---------------------------------------------------------------------------- # # Phoebus-specific derived class from Parthenon's phdf datafile reader class phoedf(phdf.phdf): @@ -62,8 +63,7 @@ def __init__(self, filename, geomfile=None): if self.RadiationActive: self.NumSpecies = self.Params["radiation/num_species"] - # TODO(BRR) store output data only once with get/load functions - + # Construct geometry from gcov if geomfile is None: self.flatgcov = self.Get("g.c.gcov", flatten=False) @@ -71,6 +71,7 @@ def flatten_indices(mu, nu): ind = [[0, 1, 3, 6], [1, 2, 4, 7], [3, 4, 5, 8], [6, 7, 8, 9]] return ind[mu][nu] + # Not necessary to provide Minkowski geometry if self.flatgcov is None: if self.Params["geometry/geometry_name"] == "Minkowski": self.flatgcov = np.zeros( @@ -80,11 +81,37 @@ def flatten_indices(mu, nu): self.flatgcov[:, flatten_indices(1, 1), :, :, :] = 1.0 self.flatgcov[:, flatten_indices(2, 2), :, :, :] = 1.0 self.flatgcov[:, flatten_indices(3, 3), :, :, :] = 1.0 + else: + assert ( + False + ), "Geometry not provided in dump file but also not Minkowski!" self.gcov = np.zeros([self.NumBlocks, 4, 4, self.Nx3, self.Nx2, self.Nx1]) + + # Account for included ghost zones in flattened geometry data + if self.flatgcov.shape[-1] != self.Nx1: + istart = self.NGhost + iend = -self.NGhost + else: + istart = None + iend = None + if self.flatgcov.shape[-2] != self.Nx2: + jstart = self.NGhost + jend = -self.NGhost + else: + jstart = None + jend = None + if self.flatgcov.shape[-3] != self.Nx3: + kstart = self.NGhost + kend = -self.NGhost + else: + kstart = None + kend = None + + # Unroll gcov for convenience for mu in range(4): for nu in range(4): self.gcov[:, mu, nu, :, :, :] = self.flatgcov[ - :, flatten_indices(mu, nu), :, :, : + :, flatten_indices(mu, nu), kstart:kend, jstart:jend, istart:iend ] del self.flatgcov self.gcon = np.zeros([self.NumBlocks, 4, 4, self.Nx3, self.Nx2, self.Nx1]) diff --git a/scripts/python/plot_torus.py b/scripts/python/plot_torus.py index 4ceb9b69..44ba73b7 100644 --- a/scripts/python/plot_torus.py +++ b/scripts/python/plot_torus.py @@ -117,7 +117,6 @@ def plot_frame(ifname, fname, savefig, geomfile=None, rlim=40, coords="cartesian ax.set_title(r"$\Gamma$") if rad_active: - ax = axes[0, 2] Pg = dfile.GetPg() Pm = np.clip(dfile.GetPm(), 1.0e-20, 1.0e20) @@ -272,7 +271,6 @@ def plot_frame(ifname, fname, savefig, geomfile=None, rlim=40, coords="cartesian if __name__ == "__main__": - parser = argparse.ArgumentParser(description="Plot neutrino thermalization") parser.add_argument( "-gf", diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 718387a0..d1e558d5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -20,6 +20,7 @@ add_executable(phoebus fixup/fixup_c2p.cpp fixup/fixup_radc2p.cpp fixup/fixup_src.cpp + fixup/fixup_netfield.cpp fluid/con2prim.hpp fluid/con2prim_robust.hpp diff --git a/src/analysis/history_utils.hpp b/src/analysis/history_utils.hpp index 4537d133..5eba068b 100644 --- a/src/analysis/history_utils.hpp +++ b/src/analysis/history_utils.hpp @@ -37,16 +37,16 @@ KOKKOS_INLINE_FUNCTION Real CalcMassFlux(Pack &pack, Geometry &geom, const int p const int b, const int k, const int j, const int i) { - Real gdet = geom.DetGamma(CellLocation::Cent, k, j, i); - Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); + Real gdet = geom.DetGamma(CellLocation::Cent, b, k, j, i); + Real lapse = geom.Lapse(CellLocation::Cent, b, k, j, i); Real shift[3]; - geom.ContravariantShift(CellLocation::Cent, k, j, i, shift); + geom.ContravariantShift(CellLocation::Cent, b, k, j, i, shift); const Real vel[] = {pack(b, pvel_lo, k, j, i), pack(b, pvel_lo + 1, k, j, i), pack(b, pvel_hi, k, j, i)}; Real gcov4[4][4]; - geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov4); + geom.SpacetimeMetric(CellLocation::Cent, b, k, j, i, gcov4); const Real W = phoebus::GetLorentzFactor(vel, gcov4); const Real ucon = vel[0] - shift[0] * W / lapse; @@ -60,7 +60,7 @@ KOKKOS_INLINE_FUNCTION Real CalcMagnetization(Pack &pack, Geometry &geom, const int prho, const int b, const int k, const int j, const int i) { Real gam[3][3]; - geom.Metric(CellLocation::Cent, k, j, i, gam); + geom.Metric(CellLocation::Cent, b, k, j, i, gam); const Real Bp[] = {pack(b, pb_lo, k, j, i), pack(b, pb_lo + 1, k, j, i), pack(b, pb_hi, k, j, i)}; const Real vcon[] = {pack(b, pvel_lo, k, j, i), pack(b, pvel_lo + 1, k, j, i), @@ -80,11 +80,11 @@ KOKKOS_INLINE_FUNCTION Real CalcEMEnergyFlux(Pack &pack, Geometry &geom, const int b, const int k, const int j, const int i) { Real gam[3][3]; - geom.Metric(CellLocation::Cent, k, j, i, gam); - Real gdet = geom.DetGamma(CellLocation::Cent, k, j, i); - Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); + geom.Metric(CellLocation::Cent, b, k, j, i, gam); + Real gdet = geom.DetGamma(CellLocation::Cent, b, k, j, i); + Real lapse = geom.Lapse(CellLocation::Cent, b, k, j, i); Real shift[3]; - geom.ContravariantShift(CellLocation::Cent, k, j, i, shift); + geom.ContravariantShift(CellLocation::Cent, b, k, j, i, shift); const Real vcon[] = {pack(b, pvel_lo, k, j, i), pack(b, pvel_lo + 1, k, j, i), pack(b, pvel_hi, k, j, i)}; @@ -118,11 +118,11 @@ KOKKOS_INLINE_FUNCTION Real CalcEMMomentumFlux(Pack &pack, Geometry &geom, const int b, const int k, const int j, const int i) { Real gam[3][3]; - geom.Metric(CellLocation::Cent, k, j, i, gam); - Real gdet = geom.DetGamma(CellLocation::Cent, k, j, i); - Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); + geom.Metric(CellLocation::Cent, b, k, j, i, gam); + Real gdet = geom.DetGamma(CellLocation::Cent, b, k, j, i); + Real lapse = geom.Lapse(CellLocation::Cent, b, k, j, i); Real shift[3]; - geom.ContravariantShift(CellLocation::Cent, k, j, i, shift); + geom.ContravariantShift(CellLocation::Cent, b, k, j, i, shift); const Real vcon[] = {pack(b, pvel_lo, k, j, i), pack(b, pvel_lo + 1, k, j, i), pack(b, pvel_hi, k, j, i)}; @@ -158,7 +158,7 @@ template KOKKOS_INLINE_FUNCTION Real CalcMagneticFluxPhi(Pack &pack, Geometry &geom, const int cb_lo, const int b, const int k, const int j, const int i) { - Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); + Real lapse = geom.Lapse(CellLocation::Cent, b, k, j, i); // \int B^r sqrt(-gdet) := (detgam B^r) * lapse return std::abs(pack(b, cb_lo, k, j, i)) * lapse; diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index b4380da4..1ec75a90 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -18,12 +18,14 @@ #include #include +#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,30 @@ std::shared_ptr 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); + if (enable_phi_enforcement) { + PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::FMKS), + "Phi enforcement only supported for BH geometry!"); + Real enforced_phi = pin->GetReal("fixup", "enforced_phi"); + 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.); + PARTHENON_REQUIRE(enforced_phi_start_time >= 0., + "Phi enforcement start time must be non-negative!"); + params.Add("enforced_phi_start_time", enforced_phi_start_time); + + // Mutable parameters for following sanity maintenance + params.Add("dphi_dt", 0., true); + params.Add("next_dphi_dt_update_time", enforced_phi_start_time, true); + params.Add("phi_factor", 0., true); + } if (enable_mhd_floors && !enable_mhd) { enable_mhd_floors = false; diff --git a/src/fixup/fixup.hpp b/src/fixup/fixup.hpp index 6060454a..c063ebce 100644 --- a/src/fixup/fixup.hpp +++ b/src/fixup/fixup.hpp @@ -1,4 +1,4 @@ -// © 2021-2022. Triad National Security, LLC. All rights reserved. +// © 2021-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. @@ -17,9 +17,11 @@ #include #include +#include #include #include using namespace parthenon::package::prelude; +using namespace parthenon::driver::prelude; namespace fixup { @@ -38,6 +40,20 @@ template TaskStatus ConservedToPrimitiveFixup(T *rc); template TaskStatus SourceFixup(T *rc); +TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, const int stage, + std::vector *sums); +TaskStatus ModifyNetField(MeshData *, const Real t, const Real dt, const int stage, + const bool fiducial, const Real fiducial_factor = 1.); +TaskStatus EvaluateNetFieldScaleControls(MeshData *, const Real t, const Real dt, + const int stage, const Real Mdot, + const Real Phi); +TaskStatus NetFieldStartReduce(MeshData *md, const Real t, const int stage, + AllReduce> *net_field_totals); +TaskStatus NetFieldCheckReduce(MeshData *md, const Real t, const int stage, + AllReduce> *net_field_totals); +TaskStatus UpdateNetFieldScaleControls(MeshData *md, const Real t, const Real dt, + const int stage, std::vector *vals0, + std::vector *vals1); static struct ConstantRhoSieFloor { } constant_rho_sie_floor_tag; diff --git a/src/fixup/fixup_c2p.cpp b/src/fixup/fixup_c2p.cpp index deeead29..5a110a07 100644 --- a/src/fixup/fixup_c2p.cpp +++ b/src/fixup/fixup_c2p.cpp @@ -1,4 +1,4 @@ -// © 2022. Triad National Security, LLC. All rights reserved. +// © 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. diff --git a/src/fixup/fixup_netfield.cpp b/src/fixup/fixup_netfield.cpp new file mode 100644 index 00000000..c253145f --- /dev/null +++ b/src/fixup/fixup_netfield.cpp @@ -0,0 +1,254 @@ +// © 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. +// Department of Energy/National Nuclear Security Administration. All +// rights in the program are reserved by Triad National Security, LLC, +// and the U.S. Department of Energy/National Nuclear Security +// Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, +// distribute copies to the public, perform publicly and display +// publicly, and to permit others to do so. + +#include + +#include "fixup.hpp" + +#include +#include +#include + +#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" + +using Microphysics::EOS::EOS; +using robust::ratio; + +namespace fixup { + +TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, const int stage, + std::vector *sums) { + auto *pm = md->GetParentPointer(); + StateDescriptor *fix_pkg = pm->packages.Get("fixup").get(); + + const bool enable_phi_enforcement = fix_pkg->Param("enable_phi_enforcement"); + + if (stage == 1) { + if (enable_phi_enforcement) { + const Real enforced_phi_start_time = + fix_pkg->Param("enforced_phi_start_time"); + const Real next_dphi_dt_update_time = + fix_pkg->Param("next_dphi_dt_update_time"); + if (t >= next_dphi_dt_update_time && t >= enforced_phi_start_time) { + const Real Mdot = History::ReduceMassAccretionRate(md); + const Real Phi = History::ReduceMagneticFluxPhi(md); + + (*sums)[0] += Mdot; + (*sums)[1] += Phi; + } + } + } + + return TaskStatus::complete; +} + +TaskStatus NetFieldStartReduce(MeshData *md, const Real t, const int stage, + AllReduce> *net_field_totals) { + auto *pm = md->GetParentPointer(); + StateDescriptor *fix_pkg = pm->packages.Get("fixup").get(); + + const bool enable_phi_enforcement = fix_pkg->Param("enable_phi_enforcement"); + + if (enable_phi_enforcement) { + const Real enforced_phi_start_time = fix_pkg->Param("enforced_phi_start_time"); + const Real next_dphi_dt_update_time = + fix_pkg->Param("next_dphi_dt_update_time"); + if (t >= next_dphi_dt_update_time && t >= enforced_phi_start_time) { + TaskStatus status = net_field_totals->StartReduce(MPI_SUM); + return status; + } + } + + return TaskStatus::complete; +} + +TaskStatus NetFieldCheckReduce(MeshData *md, const Real t, const int stage, + AllReduce> *net_field_totals) { + if (stage != 1) { + return TaskStatus::complete; + } + + auto *pm = md->GetParentPointer(); + StateDescriptor *fix_pkg = pm->packages.Get("fixup").get(); + + const bool enable_phi_enforcement = fix_pkg->Param("enable_phi_enforcement"); + + if (enable_phi_enforcement) { + const Real enforced_phi_start_time = fix_pkg->Param("enforced_phi_start_time"); + const Real next_dphi_dt_update_time = + fix_pkg->Param("next_dphi_dt_update_time"); + if (t >= next_dphi_dt_update_time && t >= enforced_phi_start_time) { + TaskStatus status = net_field_totals->CheckReduce(); + return status; + } + } + + return TaskStatus::complete; +} + +TaskStatus UpdateNetFieldScaleControls(MeshData *md, const Real t, const Real dt, + const int stage, std::vector *vals0, + std::vector *vals1) { + if (stage != 1) { + return TaskStatus::complete; + } + auto *pm = md->GetParentPointer(); + StateDescriptor *fix_pkg = pm->packages.Get("fixup").get(); + + const bool enable_phi_enforcement = fix_pkg->Param("enable_phi_enforcement"); + + if (enable_phi_enforcement) { + const Real next_dphi_dt_update_time = + fix_pkg->Param("next_dphi_dt_update_time"); + if (t >= next_dphi_dt_update_time) { + // Real phi_factor = fix_pkg->Param("phi_factor"); + Real dphi_dt = fix_pkg->Param("dphi_dt"); + const Real enforced_phi = fix_pkg->Param("enforced_phi"); + const Real enforced_phi_timescale = fix_pkg->Param("enforced_phi_timescale"); + + Real phi_factor = (*vals1)[1] - (*vals0)[1]; + // Limit change of dphi_dt to some factor of the original value if is non-zero + const Real fiducial_dphi = + (enforced_phi - (*vals0)[1] / std::sqrt((*vals0)[0])) / enforced_phi_timescale; + const Real dphi_dt_change = fiducial_dphi - dphi_dt; + const Real reference_dphi_dt_change = + robust::sgn(dphi_dt_change) * std::fabs(dphi_dt); + if (dphi_dt != 0.) { + dphi_dt += std::clamp(dphi_dt_change, 3. / 4. * reference_dphi_dt_change, + 4. / 3. * reference_dphi_dt_change); + } else { + dphi_dt += dphi_dt_change; + } + + fix_pkg->UpdateParam("phi_factor", phi_factor); + fix_pkg->UpdateParam("dphi_dt", dphi_dt); + + // Update timing values + const Real enforced_phi_cadence = fix_pkg->Param("enforced_phi_cadence"); + fix_pkg->UpdateParam("next_dphi_dt_update_time", + next_dphi_dt_update_time + enforced_phi_cadence); + } + } + return TaskStatus::complete; +} + +TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, + const int stage, const bool fiducial, + const Real fiducial_factor) { + if (stage != 1 && fiducial) { + return TaskStatus::complete; + } + auto *pm = md->GetParentPointer(); + StateDescriptor *fix_pkg = pm->packages.Get("fixup").get(); + + const bool enable_phi_enforcement = fix_pkg->Param("enable_phi_enforcement"); + + if (enable_phi_enforcement) { + + const Real next_dphi_dt_update_time = + fix_pkg->Param("next_dphi_dt_update_time"); + + if (fiducial && !(t >= next_dphi_dt_update_time)) { + return TaskStatus::complete; + } + + namespace p = fluid_prim; + namespace c = fluid_cons; + const std::vector 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("derefine_poles"); + Real h = gpkg->Param("h"); + Real xt = gpkg->Param("xt"); + Real alpha = gpkg->Param("alpha"); + Real x0 = gpkg->Param("x0"); + Real smooth = gpkg->Param("smooth"); + const Real Rh = gpkg->Param("Rh"); + auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); + + Real phi_factor = + fix_pkg->Param("phi_factor") * fix_pkg->Param("dphi_dt") * dt; + if (fiducial) { + phi_factor = fiducial_factor; + } + + // Calculate hyperbola-based magnetic field configuration inside the event horizon + ParArrayND A("vector potential", pack.GetDim(5), jb.e + 2, ib.e + 2); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "Phoebus::Fixup::ModifyNetField::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 * std::sin(th); + const Real z = r * std::cos(th); + const Real a_hyp = Rh; + const Real b_hyp = 3. * a_hyp; + const Real x_hyp = a_hyp * std::sqrt(1. + std::pow(z / b_hyp, 2.)); + const Real q = (std::pow(x, 2) - std::pow(x_hyp, 2)) / std::pow(x_hyp, 2.); + if (x < x_hyp) { + A(b, j, i) = q; + } + }); + + // Modify B field + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "Phoebus::Fixup::ModifyNetField::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 * + (-(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 * + ((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)); + + SPACELOOP(ii) { + pack(b, cblo + ii, k, j, i) = pack(b, pblo + ii, k, j, i) * gammadet; + } + }); + } // enable_phi_enforcement + + return TaskStatus::complete; +} + +} // namespace fixup diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index ab772b7c..cb761d5b 100644 --- a/src/phoebus_driver.cpp +++ b/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,87 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) { } } + // Extra per-step user work + TaskRegion &sync_region_5 = tc.AddRegion(num_partitions); + int sync_region_5_dep_id; + net_field_totals.val.resize(2); // Mdot, Phi + net_field_totals_2.val.resize(2); + for (int i = 0; i < 2; i++) { + net_field_totals.val[i] = 0.; + net_field_totals_2.val[i] = 0.; + } + for (int i = 0; i < num_partitions; i++) { + sync_region_5_dep_id = 0; + auto &tl = sync_region_5[i]; + + auto &md = pmesh->mesh_data.GetOrAdd(stage_name[stage], i); + + // Begin tuning region that only evaluates occasionally and on first stage + // TODO(BRR) Setting scale factor of B field need only be done once per simulation + + // Evaluate current Mdot, Phi + auto sum_mdot_1 = tl.AddTask(none, fixup::SumMdotPhiForNetFieldScaling, md.get(), t, + stage, &(net_field_totals.val)); + + sync_region_5.AddRegionalDependencies(sync_region_5_dep_id, i, sum_mdot_1); + sync_region_5_dep_id++; + + TaskID start_reduce_1 = (i == 0 ? tl.AddTask(sum_mdot_1, fixup::NetFieldStartReduce, + md.get(), t, stage, &net_field_totals) + : none); + // Test the reduction until it completes + TaskID finish_reduce_1 = tl.AddTask(start_reduce_1, fixup::NetFieldCheckReduce, + md.get(), t, stage, &net_field_totals); + sync_region_5.AddRegionalDependencies(sync_region_5_dep_id, i, finish_reduce_1); + sync_region_5_dep_id++; + + auto mod_net = tl.AddTask(finish_reduce_1, fixup::ModifyNetField, md.get(), t, + beta * dt, stage, true, 1.); + + // Evaluate Mdot, Phi (only Phi changes) after modifying B field + auto sum_mdot_2 = tl.AddTask(mod_net, fixup::SumMdotPhiForNetFieldScaling, md.get(), + t, stage, &(net_field_totals_2.val)); + sync_region_5.AddRegionalDependencies(sync_region_5_dep_id, i, sum_mdot_2); + sync_region_5_dep_id++; + + TaskID start_reduce_2 = (i == 0 ? tl.AddTask(sum_mdot_2, fixup::NetFieldStartReduce, + md.get(), t, stage, &net_field_totals_2) + : none); + // Test the reduction until it completes + TaskID finish_reduce_2 = tl.AddTask(start_reduce_2, fixup::NetFieldCheckReduce, + md.get(), t, stage, &net_field_totals_2); + sync_region_5.AddRegionalDependencies(sync_region_5_dep_id, i, finish_reduce_2); + sync_region_5_dep_id++; + + // Remove artificial contribution to net field + auto mod_net_2 = tl.AddTask(finish_reduce_2, fixup::ModifyNetField, md.get(), t, + beta * dt, stage, true, -1.); + + auto set_scale = + tl.AddTask(mod_net_2, fixup::UpdateNetFieldScaleControls, md.get(), t, dt, stage, + &(net_field_totals.val), &(net_field_totals_2.val)); + + // End tuning region that only evaluates occasionally and on first stage + + // Update net field every stage of every timestep given current tuning parameters + auto update_netfield = tl.AddTask(set_scale, fixup::ModifyNetField, md.get(), t, + beta * dt, stage, false, 0.0); + } + // 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]); + + StateDescriptor *fix_pkg = pmb->packages.Get("fixup").get(); + const bool enable_phi_enforcement = fix_pkg->Param("enable_phi_enforcement"); + + if (enable_phi_enforcement) { + 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); diff --git a/src/phoebus_driver.hpp b/src/phoebus_driver.hpp index da0d20e9..9caf2486 100644 --- a/src/phoebus_driver.hpp +++ b/src/phoebus_driver.hpp @@ -45,6 +45,9 @@ class PhoebusDriver : public EvolutionDriver { AllReduce dNtot; AllReduce> particle_resolution; AllReduce particles_outstanding; + + AllReduce> net_field_totals; + AllReduce> net_field_totals_2; }; parthenon::Packages_t ProcessPackages(std::unique_ptr &pin); diff --git a/src/phoebus_utils/reduction.hpp b/src/phoebus_utils/reduction.hpp index cb7f555c..c548e773 100644 --- a/src/phoebus_utils/reduction.hpp +++ b/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 diff --git a/tst/regression/regression_test.py b/tst/regression/regression_test.py index 5dda3fa9..6f1450e3 100644 --- a/tst/regression/regression_test.py +++ b/tst/regression/regression_test.py @@ -37,6 +37,7 @@ # Utility functions # + # -- Compare two values up to some floating point tolerance def soft_equiv(val, ref, tol=1.0e-5): numerator = np.fabs(val - ref) @@ -48,6 +49,35 @@ def soft_equiv(val, ref, tol=1.0e-5): return True +# -- Read value of parameter in input file +def read_input_value(block, key, input_file): + with open(input_file, "r") as infile: + lines = infile.readlines() + for line in lines: + sline = line.strip() + + # Skip empty lines and comments + if len(sline) == 0 or sline[0] == "#": + continue + + # Check for block + elif sline[0] == "<": + current_block = sline.split("<")[1].split(">")[0] + continue + + # Ignore multi-value lines + elif len(sline.split("=")) != 2 or "," in sline or "&" in sline: + continue + + else: + current_key = sline.split("=")[0].strip() + + if block == current_block and key == current_key: + return sline.split("=")[1].strip() + + assert False, "block/key not found!" + + # -- Modify key in input file, add key (and block) if not present, write new file def modify_input(dict_key, value, input_file): key = dict_key.split("/")[-1] @@ -79,7 +109,6 @@ def modify_input(dict_key, value, input_file): continue else: - current_key = sline.split("=")[0].strip() current_value = sline.split("=")[1].strip() @@ -117,6 +146,7 @@ def modify_input(dict_key, value, input_file): # Common regression test tools # + # -- Configure and build phoebus with problem-specific options def build_code(geometry, use_gpu=False, build_type="Release"): if os.path.isdir(BUILD_DIR): @@ -216,6 +246,17 @@ def gold_comparison( compression_factor=1, tolerance=1.0e-5, ): + problem = read_input_value("phoebus", "problem", input_file) + print("\n=== GOLD COMPARISON TEST PROBLEM ===") + print(f"= problem: {problem}") + print(f"= executable: {executable}") + print(f"= geometry: {geometry}") + print(f"= use_gpu: {use_gpu}") + print(f"= use_mpiexec: {use_mpiexec}") + print(f"= build_type: {build_type}") + print(f"= compression: {compression_factor}") + print(f"= tolerance: {tolerance}") + print("====================================\n") if executable is None: executable = os.path.join(BUILD_DIR, "src", "phoebus")