From e82a155e27cd18916ae6741c83a6c51a81715351 Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Thu, 30 Mar 2023 13:25:15 -0600 Subject: [PATCH 01/25] First attempt --- src/fixup/fixup.cpp | 118 ++++++++++++++++++++++++++++++++ src/fixup/fixup.hpp | 1 + src/phoebus_driver.cpp | 16 +++++ src/phoebus_utils/reduction.hpp | 9 ++- 4 files changed, 143 insertions(+), 1 deletion(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index b4380da4..b5a3ba34 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,18 @@ 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); + 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); + int enforced_phi_cadence = pin->GetOrAddInteger("fixup", "enforced_phi_cadence", 10); + params.Add("enforced_phi_cadence", enforced_phi_cadence); + 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 +659,110 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { return TaskStatus::complete; } +TaskStatus EndOfStepModify(MeshData *md, const Real dt) { + auto *pm = md->GetParentPointer(); + StateDescriptor *fix_pkg = pm->packages.Get("fixup").get(); + + const bool enable_phi_enforcement = fix_pkg->Param("enable_phi_enforcement"); + + 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"); + auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); + + if (enable_phi_enforcement) { + const Real enforced_phi = fix_pkg->Param("enforced_phi"); + const Real enforced_phi_timescale = fix_pkg->Param("enforced_phi_timescale"); + const int enforced_phi_cadence = fix_pkg->Param("enforced_phi_cadence"); + + // Measure phi + const Real phi_proc = History::ReduceMagneticFluxPhi(md); + const Real phi = reduction::Sum(phi_proc); + + // Calculate amount of phi to add this timestep + const Real dphi = (enforced_phi - phi)*dt/enforced_phi_timescale;//*enforced_phi_cadence; + + // Update B field + ParArrayND A("vector potential", pack.GetDim(5), jb.e + 1, ib.e + 1); + parthenon::par_for(parthenon::LoopPatternMDRange(), "Phoebus::Fixup::EndOfStepModify::EvaluateQ", + DevExecSpace(), + 0, pack.GetDim(5) - 1, jb.s + 1, jb.e, ib.s + 1, ib.e, + KOKKOS_LAMBDA(const int b, const int j, const int i) { + const auto &coords = pack.GetCoords(b); + const Real x1 = coords.Xc<1>(0, j, i); + const Real x2 = coords.Xc<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 = 2.; + const Real b_hyp = 60.; + 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; + } + } + ); + + 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); + printf("here! phi: %e\n", fiducial_phi); + exit(-1); + + // Rescale the phi contribution + +// pmb->par_for("Phoebus::Fixup::EndOfStepModify::EvaluateBField", +// kb.s, kb.e, jb.s, jb.e-1, ib.s, ib.e-1, +// KOKKOS_LAMBDA(const int k, const int j, const int i) { +// const real gammadet = geom.DetGamma(CellLocation::Cent, k, j, i); +// v(pblo, k, j, i) += (-(A(j, i) - A(j + 1, i) + A(j, i + 1) - A(j + 1, i + 1)) /(2.0 * coords.CellWidthFA(X2DIR, k, j, i) * gamdet)); +// v(pblo, k, j, i) += ((A(j, i) + A(j + 1, i) - A(j, i + 1) - A(j + 1, i + 1)) /(2.0 * coords.CellWidthFA(X1DIR, k, j, i) * gamdet)); +// } +// +// // Update conserved vars +// p2c(); +// ); + } + + + return TaskStatus::complete; +} + template TaskStatus ApplyFloors(T *rc) { auto *pm = rc->GetParentPointer().get(); diff --git a/src/fixup/fixup.hpp b/src/fixup/fixup.hpp index 6060454a..ff701b81 100644 --- a/src/fixup/fixup.hpp +++ b/src/fixup/fixup.hpp @@ -38,6 +38,7 @@ template TaskStatus ConservedToPrimitiveFixup(T *rc); template TaskStatus SourceFixup(T *rc); +TaskStatus EndOfStepModify(MeshData, const Real dt); static struct ConstantRhoSieFloor { } constant_rho_sie_floor_tag; diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index b05083eb..24ecfcd8 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -613,6 +613,22 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) { } } + // Communicate flux corrections and update independent data with fluxes and geometric + // sources + TaskRegion &sync_region_5 = tc.AddRegion(num_partitions); + PARTHENON_REQUIRE(num_partitions == 1, "Reductions don't work for multiple partitions?"); + if (stage == integrator->nstages) { + 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, dt); + + } + } + + return tc; } // namespace phoebus 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 From 32f4655365698445e1a712bb1b3a391281c48430 Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Tue, 18 Apr 2023 14:12:48 -0600 Subject: [PATCH 02/25] First attempt --- scripts/python/phoedf.py | 2 +- src/fixup/fixup.cpp | 192 ++++++++++++++++++++++++--------------- src/fixup/fixup.hpp | 2 +- src/phoebus_driver.cpp | 50 +++++++--- 4 files changed, 160 insertions(+), 86 deletions(-) diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index ae8ee698..8a0b5e9a 100644 --- a/scripts/python/phoedf.py +++ b/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): diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index b5a3ba34..8a0b1458 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -73,17 +73,19 @@ 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); + 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); + Real enforced_phi_timescale = + pin->GetOrAddReal("fixup", "enforced_phi_timescale", 1.e3); params.Add("enforced_phi_timescale", enforced_phi_timescale); - int enforced_phi_cadence = pin->GetOrAddInteger("fixup", "enforced_phi_cadence", 10); + Real enforced_phi_cadence = pin->GetOrAddInteger("fixup", "enforced_phi_cadence", 10); params.Add("enforced_phi_cadence", enforced_phi_cadence); if (enable_phi_enforcement) { PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::FMKS), - "Phi enforcement only supported for BH geometry!"); + "Phi enforcement only supported for BH geometry!"); } if (enable_mhd_floors && !enable_mhd) { @@ -659,7 +661,12 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) { return TaskStatus::complete; } -TaskStatus EndOfStepModify(MeshData *md, const Real dt) { +static Real dphi_dt = 0.; +static Real next_dphi_dt_update_time = 0.; +static Real phi_factor = 0.; + +TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, + bool last_stage) { auto *pm = md->GetParentPointer(); StateDescriptor *fix_pkg = pm->packages.Get("fixup").get(); @@ -681,84 +688,123 @@ TaskStatus EndOfStepModify(MeshData *md, const Real dt) { 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"); + 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"); auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); if (enable_phi_enforcement) { const Real enforced_phi = fix_pkg->Param("enforced_phi"); const Real enforced_phi_timescale = fix_pkg->Param("enforced_phi_timescale"); - const int enforced_phi_cadence = fix_pkg->Param("enforced_phi_cadence"); + const Real enforced_phi_cadence = fix_pkg->Param("enforced_phi_cadence"); + + // This only needs to be done at initialization + ParArrayND A("vector potential", pack.GetDim(5), jb.e + 1, ib.e + 1); + // Update B field + parthenon::par_for( + parthenon::LoopPatternMDRange(), "Phoebus::Fixup::EndOfStepModify::EvaluateQ", + DevExecSpace(), 0, pack.GetDim(5) - 1, jb.s + 1, jb.e, ib.s + 1, ib.e, + KOKKOS_LAMBDA(const int b, const int j, const int i) { + const auto &coords = pack.GetCoords(b); + const Real x1 = coords.Xc<1>(0, j, i); + const Real x2 = coords.Xc<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 = 2.; + const Real b_hyp = 60.; + 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; + } + }); - // Measure phi - const Real phi_proc = History::ReduceMagneticFluxPhi(md); - const Real phi = reduction::Sum(phi_proc); + // Recalculate dphi_dt occasionally + if (last_stage && t >= next_dphi_dt_update_time) { + next_dphi_dt_update_time += enforced_phi_cadence; - // Calculate amount of phi to add this timestep - const Real dphi = (enforced_phi - phi)*dt/enforced_phi_timescale;//*enforced_phi_cadence; + // Measure phi + const Real phi_proc = History::ReduceMagneticFluxPhi(md); + const Real phi = reduction::Sum(phi_proc); + printf("Original phi: %e\n", phi); - // Update B field - ParArrayND A("vector potential", pack.GetDim(5), jb.e + 1, ib.e + 1); - parthenon::par_for(parthenon::LoopPatternMDRange(), "Phoebus::Fixup::EndOfStepModify::EvaluateQ", - DevExecSpace(), - 0, pack.GetDim(5) - 1, jb.s + 1, jb.e, ib.s + 1, ib.e, - KOKKOS_LAMBDA(const int b, const int j, const int i) { - const auto &coords = pack.GetCoords(b); - const Real x1 = coords.Xc<1>(0, j, i); - const Real x2 = coords.Xc<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 = 2.; - const Real b_hyp = 60.; - 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; - } - } - ); + // Calculate amount of phi to add this timestep + const Real dphi = + (enforced_phi - phi) * dt / enforced_phi_timescale; //*enforced_phi_cadence; + dphi_dt = dphi / dt; + printf("dphi: %e\n", dphi); + printf("dphi_dt: %e\n", 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); - printf("here! phi: %e\n", fiducial_phi); - exit(-1); - - // Rescale the phi contribution - -// pmb->par_for("Phoebus::Fixup::EndOfStepModify::EvaluateBField", -// kb.s, kb.e, jb.s, jb.e-1, ib.s, ib.e-1, -// KOKKOS_LAMBDA(const int k, const int j, const int i) { -// const real gammadet = geom.DetGamma(CellLocation::Cent, k, j, i); -// v(pblo, k, j, i) += (-(A(j, i) - A(j + 1, i) + A(j, i + 1) - A(j + 1, i + 1)) /(2.0 * coords.CellWidthFA(X2DIR, k, j, i) * gamdet)); -// v(pblo, k, j, i) += ((A(j, i) + A(j + 1, i) - A(j, i + 1) - A(j + 1, i + 1)) /(2.0 * coords.CellWidthFA(X1DIR, k, j, i) * gamdet)); -// } -// -// // Update conserved vars -// p2c(); -// ); - } + 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); + printf("here! phi: %e\n", fiducial_phi); + phi_factor = dphi / (fiducial_phi - phi) / dt; + printf("phi_factor: %e\n", phi_factor); + // exit(-1); + + // 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 - 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) += + 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)); + }); + } return TaskStatus::complete; } diff --git a/src/fixup/fixup.hpp b/src/fixup/fixup.hpp index ff701b81..24767624 100644 --- a/src/fixup/fixup.hpp +++ b/src/fixup/fixup.hpp @@ -38,7 +38,7 @@ template TaskStatus ConservedToPrimitiveFixup(T *rc); template TaskStatus SourceFixup(T *rc); -TaskStatus EndOfStepModify(MeshData, const Real dt); +TaskStatus EndOfStepModify(MeshData*, const Real t, const Real dt, const bool last_stage); static struct ConstantRhoSieFloor { } constant_rho_sie_floor_tag; diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 0d3f45b9..12a3a736 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,27 @@ 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); @@ -614,20 +636,26 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) { } } - // Communicate flux corrections and update independent data with fluxes and geometric - // sources - TaskRegion &sync_region_5 = tc.AddRegion(num_partitions); - PARTHENON_REQUIRE(num_partitions == 1, "Reductions don't work for multiple partitions?"); - if (stage == integrator->nstages) { - for (int i = 0; i < num_partitions; i++) { - auto &tl = sync_region_5[i]; + //// 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 &md = pmesh->mesh_data.GetOrAdd(stage_name[stage], i); - auto end_mods = tl.AddTask(none, fixup::EndOfStepModify, md, dt); + // 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_3[i]; + // auto &sc = pmb->meshblock_data.Get(stage_name[stage]); + + // auto p2c = tl.AddTask(none, PrimitiveToConserved, sc); + //} return tc; From 884a4daca994fa462960a03b6fec796455f9e40f Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Tue, 18 Apr 2023 15:41:56 -0600 Subject: [PATCH 03/25] sing? --- external/singularity-eos | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/singularity-eos b/external/singularity-eos index e8c750dd..b4c9ce79 160000 --- a/external/singularity-eos +++ b/external/singularity-eos @@ -1 +1 @@ -Subproject commit e8c750dd8102b80b16de35a781e8f56a0b4a86b8 +Subproject commit b4c9ce7962cfeb7b46585e1ef323d9bbfde1101a From 191065882451a10cfae458052848ad34b9e60e2a Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Tue, 18 Apr 2023 15:59:19 -0600 Subject: [PATCH 04/25] This is better maybe it even works --- src/fixup/fixup.cpp | 38 +++++++++++++++++++++++++++++++------- 1 file changed, 31 insertions(+), 7 deletions(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index 8a0b1458..a71401d6 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -695,7 +695,7 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, Real smooth = gpkg->Param("smooth"); auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); - if (enable_phi_enforcement) { + if (enable_phi_enforcement) {// && t > 175) { const Real enforced_phi = fix_pkg->Param("enforced_phi"); const Real enforced_phi_timescale = fix_pkg->Param("enforced_phi_timescale"); const Real enforced_phi_cadence = fix_pkg->Param("enforced_phi_cadence"); @@ -708,15 +708,19 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, DevExecSpace(), 0, pack.GetDim(5) - 1, jb.s + 1, jb.e, ib.s + 1, ib.e, KOKKOS_LAMBDA(const int b, const int j, const int i) { const auto &coords = pack.GetCoords(b); - const Real x1 = coords.Xc<1>(0, j, i); - const Real x2 = coords.Xc<2>(0, j, i); + 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); + if (i == 10) { + printf("[%i %i %i] r: %e th: %e\n", 0,j,i,r,th); + } + const Real x = r * sin(th); const Real z = r * cos(th); - const Real a_hyp = 2.; - const Real b_hyp = 60.; + 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) { @@ -788,21 +792,41 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, }); } + // Temporary!! + phi_factor = 1.e5; + // Properly update B field + printf("fac: %e\n", phi_factor*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) { + 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); + if (i == 10) { + printf("[%i %i %i]\n", k, j, i); + printf(" before b: %e %e %e\n", pack(b, pblo, k, j, i), + pack(b, pblo+1, k, j, i), pack(b, pblo+2, 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)); + + if (i == 10) { + printf(" after b: %e %e %e\n", pack(b, pblo, k, j, i), + pack(b, pblo+1, k, j, i), pack(b, pblo+2, k, j, i)); + } + + if (i == 5 && j == 5) { + printf("b: %e %e %e\n", pack(b, pblo, k, j, i), pack(b, pblo+1, k, j, i), + pack(b, pblo+2, k, j, i)); + } }); } From 35ce8974705c30769135fe9549b2555ef6ce191d Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Wed, 19 Apr 2023 08:55:57 -0600 Subject: [PATCH 05/25] Fix bug in A --- src/fixup/fixup.cpp | 56 ++++++++++++++++++++++++--------------------- 1 file changed, 30 insertions(+), 26 deletions(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index a71401d6..098d1ff5 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -695,17 +695,18 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, Real smooth = gpkg->Param("smooth"); auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); - if (enable_phi_enforcement) {// && t > 175) { + if (enable_phi_enforcement){// && t > 175) { const Real enforced_phi = fix_pkg->Param("enforced_phi"); const Real enforced_phi_timescale = fix_pkg->Param("enforced_phi_timescale"); const Real enforced_phi_cadence = fix_pkg->Param("enforced_phi_cadence"); // This only needs to be done at initialization - ParArrayND A("vector potential", pack.GetDim(5), jb.e + 1, ib.e + 1); + ParArrayND 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 + 1, jb.e, ib.s + 1, ib.e, + //DevExecSpace(), 0, pack.GetDim(5) - 1, jb.s + 1, jb.e, ib.s + 1, ib.e, + 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); @@ -713,9 +714,9 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, const Real r = tr.bl_radius(x1); const Real th = tr.bl_theta(x1, x2); - if (i == 10) { - printf("[%i %i %i] r: %e th: %e\n", 0,j,i,r,th); - } + //if (i == 10) { + // printf("[%i %i %i] r: %e th: %e\n", 0,j,i,r,th); + //} const Real x = r * sin(th); const Real z = r * cos(th); @@ -726,6 +727,9 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, if (x < x_hyp) { A(b, j, i) = q; } + // if (i == 10) { + // printf("[%i %i %i] r: %e th: %e A = %e\n", 0,j,i,r,th, A(b,j,i)); + // } }); // Recalculate dphi_dt occasionally @@ -735,14 +739,14 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, // Measure phi const Real phi_proc = History::ReduceMagneticFluxPhi(md); const Real phi = reduction::Sum(phi_proc); - printf("Original phi: %e\n", phi); + //printf("Original phi: %e\n", phi); // Calculate amount of phi to add this timestep const Real dphi = (enforced_phi - phi) * dt / enforced_phi_timescale; //*enforced_phi_cadence; dphi_dt = dphi / dt; - printf("dphi: %e\n", dphi); - printf("dphi_dt: %e\n", dphi_dt); + //printf("dphi: %e\n", dphi); + //printf("dphi_dt: %e\n", dphi_dt); parthenon::par_for( parthenon::LoopPatternMDRange(), @@ -766,9 +770,9 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, // Measure change in phi const Real fiducial_phi_proc = History::ReduceMagneticFluxPhi(md); const Real fiducial_phi = reduction::Sum(fiducial_phi_proc); - printf("here! phi: %e\n", fiducial_phi); + // printf("here! phi: %e\n", fiducial_phi); phi_factor = dphi / (fiducial_phi - phi) / dt; - printf("phi_factor: %e\n", phi_factor); + // printf("phi_factor: %e\n", phi_factor); // exit(-1); // Remove the fiducial contribution @@ -793,10 +797,10 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, } // Temporary!! - phi_factor = 1.e5; + //phi_factor = 1.e5; // Properly update B field - printf("fac: %e\n", phi_factor*dt); + // printf("fac: %e\n", phi_factor*dt); parthenon::par_for( parthenon::LoopPatternMDRange(), "Phoebus::Fixup::EndOfStepModify::EvaluateBField", DevExecSpace(), 0, @@ -805,11 +809,11 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, const auto &coords = pack.GetCoords(b); const Real gammadet = geom.DetGamma(CellLocation::Cent, b, k, j, i); - if (i == 10) { - printf("[%i %i %i]\n", k, j, i); - printf(" before b: %e %e %e\n", pack(b, pblo, k, j, i), - pack(b, pblo+1, k, j, i), pack(b, pblo+2, k, j, i)); - } + // if (i == 10) { + // printf("[%i %i %i]\n", k, j, i); + //printf(" before b: %e %e %e\n", pack(b, pblo, k, j, i), + // pack(b, pblo+1, k, j, i), pack(b, pblo+2, 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)) / @@ -818,15 +822,15 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, 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)); - if (i == 10) { - printf(" after b: %e %e %e\n", pack(b, pblo, k, j, i), - pack(b, pblo+1, k, j, i), pack(b, pblo+2, k, j, i)); - } + // if (i == 10) { + // printf(" after b: %e %e %e\n", pack(b, pblo, k, j, i), + // pack(b, pblo+1, k, j, i), pack(b, pblo+2, k, j, i)); + // } - if (i == 5 && j == 5) { - printf("b: %e %e %e\n", pack(b, pblo, k, j, i), pack(b, pblo+1, k, j, i), - pack(b, pblo+2, k, j, i)); - } + // if (i == 5 && j == 5) { + // printf("b: %e %e %e\n", pack(b, pblo, k, j, i), pack(b, pblo+1, k, j, i), + /// pack(b, pblo+2, k, j, i)); + // } }); } From 9a5e2f49edda13a40ef9e974af7ad6564a4f19bc Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Wed, 19 Apr 2023 09:02:29 -0600 Subject: [PATCH 06/25] start time for phi adjustment --- src/fixup/fixup.cpp | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index 098d1ff5..959ffe28 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -81,8 +81,10 @@ std::shared_ptr Initialize(ParameterInput *pin) { 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->GetOrAddInteger("fixup", "enforced_phi_cadence", 10); + 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!"); @@ -695,10 +697,12 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, Real smooth = gpkg->Param("smooth"); auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); - if (enable_phi_enforcement){// && t > 175) { + if (enable_phi_enforcement) { const Real enforced_phi = fix_pkg->Param("enforced_phi"); const Real enforced_phi_timescale = fix_pkg->Param("enforced_phi_timescale"); const Real enforced_phi_cadence = fix_pkg->Param("enforced_phi_cadence"); + const Real enforced_phi_start_time = fix_pkg->Param("enforced_phi_start_time"); + if (t > enforced_phi_start_time) { // This only needs to be done at initialization ParArrayND A("vector potential", pack.GetDim(5), jb.e + 2, ib.e + 2); @@ -832,7 +836,8 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, /// pack(b, pblo+2, k, j, i)); // } }); - } + } // enforced_phi_start_time + } // enable_phi_enforcement return TaskStatus::complete; } From 96bce20217aeb7287b43af58a7baceb3d56f0441 Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Wed, 19 Apr 2023 14:31:28 -0600 Subject: [PATCH 07/25] phi not Phi --- src/fixup/fixup.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index 959ffe28..1ebea6af 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -740,9 +740,13 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, 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); + const Real phi = reduction::Sum(phi_proc) / std::sqrt(mdot); //printf("Original phi: %e\n", phi); // Calculate amount of phi to add this timestep @@ -773,7 +777,7 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, // Measure change in phi const Real fiducial_phi_proc = History::ReduceMagneticFluxPhi(md); - const Real fiducial_phi = reduction::Sum(fiducial_phi_proc); + const Real fiducial_phi = reduction::Sum(fiducial_phi_proc) / std::sqrt(mdot); // printf("here! phi: %e\n", fiducial_phi); phi_factor = dphi / (fiducial_phi - phi) / dt; // printf("phi_factor: %e\n", phi_factor); From e495489a02afece478e402ae3374449b4566acd5 Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Wed, 19 Apr 2023 14:31:45 -0600 Subject: [PATCH 08/25] Formatting --- src/fixup/fixup.cpp | 206 +++++++++++++++++++++-------------------- src/fixup/fixup.hpp | 3 +- src/phoebus_driver.cpp | 21 +++-- 3 files changed, 118 insertions(+), 112 deletions(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index 1ebea6af..1d38e28b 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -83,7 +83,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { 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.); + 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), @@ -704,13 +705,13 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, const Real enforced_phi_start_time = fix_pkg->Param("enforced_phi_start_time"); if (t > enforced_phi_start_time) { - // This only needs to be done at initialization - ParArrayND A("vector potential", pack.GetDim(5), jb.e + 2, ib.e + 2); + // This only needs to be done at initialization + ParArrayND 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 + 1, jb.e, ib.s + 1, ib.e, - DevExecSpace(), 0, pack.GetDim(5) - 1, jb.s, jb.e+1, ib.s, ib.e+1, + // DevExecSpace(), 0, pack.GetDim(5) - 1, jb.s + 1, jb.e, ib.s + 1, ib.e, + 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); @@ -718,130 +719,133 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, const Real r = tr.bl_radius(x1); const Real th = tr.bl_theta(x1, x2); - //if (i == 10) { + // if (i == 10) { // printf("[%i %i %i] r: %e th: %e\n", 0,j,i,r,th); //} 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 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; } - // if (i == 10) { - // printf("[%i %i %i] r: %e th: %e A = %e\n", 0,j,i,r,th, A(b,j,i)); - // } + // if (i == 10) { + // printf("[%i %i %i] r: %e th: %e A = %e\n", 0,j,i,r,th, A(b,j,i)); + // } }); - // 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); - //printf("Original phi: %e\n", phi); - - // Calculate amount of phi to add this timestep - const Real dphi = - (enforced_phi - phi) * dt / enforced_phi_timescale; //*enforced_phi_cadence; - dphi_dt = dphi / dt; - //printf("dphi: %e\n", dphi); - //printf("dphi_dt: %e\n", dphi_dt); - + // 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); + // printf("Original phi: %e\n", phi); + + // Calculate amount of phi to add this timestep + const Real dphi = + (enforced_phi - phi) * dt / enforced_phi_timescale; //*enforced_phi_cadence; + dphi_dt = dphi / dt; + // printf("dphi: %e\n", dphi); + // printf("dphi_dt: %e\n", 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); + // printf("here! phi: %e\n", fiducial_phi); + phi_factor = dphi / (fiducial_phi - phi) / dt; + // printf("phi_factor: %e\n", phi_factor); + // exit(-1); + + // 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; + }); + } + + // Temporary!! + // phi_factor = 1.e5; + + // Properly update B field + // printf("fac: %e\n", phi_factor*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, + 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); + + // if (i == 10) { + // printf("[%i %i %i]\n", k, j, i); + // printf(" before b: %e %e %e\n", pack(b, pblo, k, j, i), + // pack(b, pblo+1, k, j, i), pack(b, pblo+2, 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)); - // 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); - // printf("here! phi: %e\n", fiducial_phi); - phi_factor = dphi / (fiducial_phi - phi) / dt; - // printf("phi_factor: %e\n", phi_factor); - // exit(-1); + // if (i == 10) { + // printf(" after b: %e %e %e\n", pack(b, pblo, k, j, i), + // pack(b, pblo+1, k, j, i), pack(b, pblo+2, k, j, i)); + // } - // 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; + // if (i == 5 && j == 5) { + // printf("b: %e %e %e\n", pack(b, pblo, k, j, i), pack(b, pblo+1, k, j, + // i), + /// pack(b, pblo+2, k, j, i)); + // } }); - } - - // Temporary!! - //phi_factor = 1.e5; - - // Properly update B field - // printf("fac: %e\n", phi_factor*dt); - 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); - - // if (i == 10) { - // printf("[%i %i %i]\n", k, j, i); - //printf(" before b: %e %e %e\n", pack(b, pblo, k, j, i), - // pack(b, pblo+1, k, j, i), pack(b, pblo+2, 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)); - - // if (i == 10) { - // printf(" after b: %e %e %e\n", pack(b, pblo, k, j, i), - // pack(b, pblo+1, k, j, i), pack(b, pblo+2, k, j, i)); - // } - - // if (i == 5 && j == 5) { - // printf("b: %e %e %e\n", pack(b, pblo, k, j, i), pack(b, pblo+1, k, j, i), - /// pack(b, pblo+2, k, j, i)); - // } - }); - } // enforced_phi_start_time - } // enable_phi_enforcement + } // enforced_phi_start_time + } // enable_phi_enforcement return TaskStatus::complete; } diff --git a/src/fixup/fixup.hpp b/src/fixup/fixup.hpp index 24767624..be601563 100644 --- a/src/fixup/fixup.hpp +++ b/src/fixup/fixup.hpp @@ -38,7 +38,8 @@ template TaskStatus ConservedToPrimitiveFixup(T *rc); template TaskStatus SourceFixup(T *rc); -TaskStatus EndOfStepModify(MeshData*, const Real t, const Real dt, const bool last_stage); +TaskStatus EndOfStepModify(MeshData *, const Real t, const Real dt, + const bool last_stage); static struct ConstantRhoSieFloor { } constant_rho_sie_floor_tag; diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 12a3a736..45e063c0 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -354,14 +354,15 @@ 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?"); + 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); - + 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); @@ -637,19 +638,20 @@ 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++) { + // 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); + // 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++) { + // 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_3[i]; // auto &sc = pmb->meshblock_data.Get(stage_name[stage]); @@ -657,7 +659,6 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) { // auto p2c = tl.AddTask(none, PrimitiveToConserved, sc); //} - return tc; } // namespace phoebus From 2186e97b1d8977b88e551a2a8491df46aa126311 Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Wed, 19 Apr 2023 14:32:59 -0600 Subject: [PATCH 09/25] Update torus input file --- inputs/torus.pin | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) 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. From abfff6aba415dce3c90b5c8f855e0a8288118c67 Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Wed, 19 Apr 2023 14:38:58 -0600 Subject: [PATCH 10/25] Cleanup --- src/fixup/fixup.cpp | 38 +------------------------------------- src/phoebus_driver.cpp | 22 ---------------------- 2 files changed, 1 insertion(+), 59 deletions(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index 1d38e28b..b6e76f06 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -710,7 +710,6 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, // Update B field parthenon::par_for( parthenon::LoopPatternMDRange(), "Phoebus::Fixup::EndOfStepModify::EvaluateQ", - // DevExecSpace(), 0, pack.GetDim(5) - 1, jb.s + 1, jb.e, ib.s + 1, ib.e, 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); @@ -719,10 +718,6 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, const Real r = tr.bl_radius(x1); const Real th = tr.bl_theta(x1, x2); - // if (i == 10) { - // printf("[%i %i %i] r: %e th: %e\n", 0,j,i,r,th); - //} - const Real x = r * sin(th); const Real z = r * cos(th); const Real a_hyp = 1.2; @@ -732,9 +727,6 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, if (x < x_hyp) { A(b, j, i) = q; } - // if (i == 10) { - // printf("[%i %i %i] r: %e th: %e A = %e\n", 0,j,i,r,th, A(b,j,i)); - // } }); // Recalculate dphi_dt occasionally @@ -748,14 +740,10 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, // Measure phi const Real phi_proc = History::ReduceMagneticFluxPhi(md); const Real phi = reduction::Sum(phi_proc) / std::sqrt(mdot); - // printf("Original phi: %e\n", phi); // Calculate amount of phi to add this timestep - const Real dphi = - (enforced_phi - phi) * dt / enforced_phi_timescale; //*enforced_phi_cadence; + const Real dphi = (enforced_phi - phi) * dt / enforced_phi_timescale; dphi_dt = dphi / dt; - // printf("dphi: %e\n", dphi); - // printf("dphi_dt: %e\n", dphi_dt); parthenon::par_for( parthenon::LoopPatternMDRange(), @@ -779,10 +767,7 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, // Measure change in phi const Real fiducial_phi_proc = History::ReduceMagneticFluxPhi(md); const Real fiducial_phi = reduction::Sum(fiducial_phi_proc) / std::sqrt(mdot); - // printf("here! phi: %e\n", fiducial_phi); phi_factor = dphi / (fiducial_phi - phi) / dt; - // printf("phi_factor: %e\n", phi_factor); - // exit(-1); // Remove the fiducial contribution parthenon::par_for( @@ -805,11 +790,7 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, }); } - // Temporary!! - // phi_factor = 1.e5; - // Properly update B field - // printf("fac: %e\n", phi_factor*dt); parthenon::par_for( parthenon::LoopPatternMDRange(), "Phoebus::Fixup::EndOfStepModify::EvaluateBField", DevExecSpace(), 0, @@ -818,12 +799,6 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, const auto &coords = pack.GetCoords(b); const Real gammadet = geom.DetGamma(CellLocation::Cent, b, k, j, i); - // if (i == 10) { - // printf("[%i %i %i]\n", k, j, i); - // printf(" before b: %e %e %e\n", pack(b, pblo, k, j, i), - // pack(b, pblo+1, k, j, i), pack(b, pblo+2, 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)) / @@ -832,17 +807,6 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, 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)); - - // if (i == 10) { - // printf(" after b: %e %e %e\n", pack(b, pblo, k, j, i), - // pack(b, pblo+1, k, j, i), pack(b, pblo+2, k, j, i)); - // } - - // if (i == 5 && j == 5) { - // printf("b: %e %e %e\n", pack(b, pblo, k, j, i), pack(b, pblo+1, k, j, - // i), - /// pack(b, pblo+2, k, j, i)); - // } }); } // enforced_phi_start_time } // enable_phi_enforcement diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 45e063c0..69eb1c86 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -637,28 +637,6 @@ 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_3[i]; - // auto &sc = pmb->meshblock_data.Get(stage_name[stage]); - - // auto p2c = tl.AddTask(none, PrimitiveToConserved, sc); - //} - return tc; } // namespace phoebus From 6602af17aade043a92bac9e2af587321ec8f7fff Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Thu, 20 Apr 2023 13:50:24 -0600 Subject: [PATCH 11/25] Hopefully fix CI --- src/fixup/fixup.cpp | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index b6e76f06..7832340f 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -675,30 +675,30 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, const bool enable_phi_enforcement = fix_pkg->Param("enable_phi_enforcement"); - namespace p = fluid_prim; - namespace c = fluid_cons; - const std::vector vars({p::bfield, c::bfield}); - PackIndexMap imap; - auto pack = md->PackVariables(vars, imap); + if (enable_phi_enforcement) { + 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 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); + 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"); - auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); + 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"); + auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); - if (enable_phi_enforcement) { const Real enforced_phi = fix_pkg->Param("enforced_phi"); const Real enforced_phi_timescale = fix_pkg->Param("enforced_phi_timescale"); const Real enforced_phi_cadence = fix_pkg->Param("enforced_phi_cadence"); From 047a8ae6b88bc7696a9cb655cb9d618cdb0217d5 Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Thu, 27 Apr 2023 08:23:25 -0600 Subject: [PATCH 12/25] Have gold_comparison in regression_test provide more (any) information so CI failures are readable --- tst/regression/regression_test.py | 40 +++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/tst/regression/regression_test.py b/tst/regression/regression_test.py index 11e356af..ca13820d 100644 --- a/tst/regression/regression_test.py +++ b/tst/regression/regression_test.py @@ -46,6 +46,34 @@ def soft_equiv(val, ref, tol = 1.e-5): else: 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] @@ -192,6 +220,18 @@ def gold_comparison(variables, input_file, modified_inputs={}, executable=None, geometry='Minkowski', use_gpu=False, use_mpiexec=False, build_type='Release', upgold=False, compression_factor=1, tolerance=1.e-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') build_code(geometry, use_gpu, build_type) From 54529db5ff588ba918a5069316e16a22ae89cbc3 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Thu, 27 Apr 2023 10:06:34 -0600 Subject: [PATCH 13/25] Test of git with vscode --- src/fixup/fixup_c2p.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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. From 5f95e36f435959695a4fa1a0813808ed3bf3bb80 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Thu, 27 Apr 2023 15:59:08 -0600 Subject: [PATCH 14/25] Working on net field tasking --- src/CMakeLists.txt | 1 + src/fixup/fixup.cpp | 152 ++++++++++++++++++++++++++++++- src/fixup/fixup.hpp | 7 ++ src/fixup/fixup_netfield.cpp | 171 +++++++++++++++++++++++++++++++++++ src/phoebus_driver.cpp | 32 ++++++- src/phoebus_driver.hpp | 2 + 6 files changed, 359 insertions(+), 6 deletions(-) create mode 100644 src/fixup/fixup_netfield.cpp 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/fixup/fixup.cpp b/src/fixup/fixup.cpp index 7832340f..4308f917 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -85,11 +85,17 @@ std::shared_ptr Initialize(ParameterInput *pin) { 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); if (enable_phi_enforcement) { PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::FMKS), "Phi enforcement only supported for BH geometry!"); } + // 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; @@ -664,10 +670,6 @@ 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 *md, const Real t, const Real dt, bool last_stage) { auto *pm = md->GetParentPointer(); @@ -703,6 +705,11 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, const Real enforced_phi_timescale = fix_pkg->Param("enforced_phi_timescale"); const Real enforced_phi_cadence = fix_pkg->Param("enforced_phi_cadence"); const Real enforced_phi_start_time = fix_pkg->Param("enforced_phi_start_time"); + + Real phi_factor = fix_pkg->Param("phi_factor"); + Real dphi_dt = fix_pkg->Param("dphi_dt"); + Real next_dphi_dt_update_time = fix_pkg->Param("next_dphi_dt_update_time"); + if (t > enforced_phi_start_time) { // This only needs to be done at initialization @@ -788,6 +795,10 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, 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; }); + + fix_pkg->UpdateParam("phi_factor", phi_factor); + fix_pkg->UpdateParam("dphi_dt", dphi_dt); + fix_pkg->UpdateParam("next_dphi_dt_update_time", next_dphi_dt_update_time); } // Properly update B field @@ -814,6 +825,139 @@ TaskStatus EndOfStepModify(MeshData *md, const Real t, const Real dt, return TaskStatus::complete; } +// TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, +// std::vector *sums) { +// printf("%s:%i\n", __FILE__, __LINE__); +// 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) { +// const Real Mdot = History::ReduceMassAccretionRate(md); +// const Real Phi = History::ReduceMagneticFluxPhi(md); +// printf("Per-MD sum: Mdot: %e Phi: %e\n", Mdot, Phi); +// +// (*sums)[0] += Mdot; +// (*sums)[1] += Phi; +// } +// } +// printf("%s:%i\n", __FILE__, __LINE__); +// +// return TaskStatus::complete; +// } + +// TODO(BRR) separate task to update + +// TaskStatus UpdateNetFieldScaleControls(MeshData *md, const Real t, const Real dt, +// const Real Mdot, const Real Phi) { +// printf("%s:%i\n", __FILE__, __LINE__); +// auto *pm = md->GetParentPointer(); +// printf("%s:%i\n", __FILE__, __LINE__); +// return TaskStatus::complete; +// } +// +// TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, +// const Real Mdot, const Real Phi, const bool fiducial, +// const Real fiducial_factor) { +// printf("%s:%i\n", __FILE__, __LINE__); +// printf("Mdot: %e Phi: %e\n", Mdot, Phi); +// 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"); +// +// 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"); +// if (fiducial) { +// phi_factor = fiducial_factor; +// } +// +// if (t >= next_dphi_dt_update_time) { +// // next_dphi_dt_update_time += enforced_phi_cadence; +// +// // 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::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 * 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::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)); +// }); +// +// // TODO(BRR) optionally call p2c here? +// } // t >= next_dphi_dt_update_time +// } // enable_phie_enforcement +// +// exit(-1); +// return TaskStatus::complete; +// } + template TaskStatus ApplyFloors(T *rc) { auto *pm = rc->GetParentPointer().get(); diff --git a/src/fixup/fixup.hpp b/src/fixup/fixup.hpp index be601563..5b71427a 100644 --- a/src/fixup/fixup.hpp +++ b/src/fixup/fixup.hpp @@ -40,6 +40,13 @@ template TaskStatus SourceFixup(T *rc); TaskStatus EndOfStepModify(MeshData *, const Real t, const Real dt, const bool last_stage); +TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, + std::vector *sums); +TaskStatus ModifyNetField(MeshData *, const Real t, const Real dt, + std::vector *vals, const bool fiducial, + const Real fiducial_factor = 1.); +TaskStatus EvaluateNetFieldScaleControls(MeshData *, const Real t, const Real dt, + const Real Mdot, const Real Phi); static struct ConstantRhoSieFloor { } constant_rho_sie_floor_tag; diff --git a/src/fixup/fixup_netfield.cpp b/src/fixup/fixup_netfield.cpp new file mode 100644 index 00000000..ad92b9ab --- /dev/null +++ b/src/fixup/fixup_netfield.cpp @@ -0,0 +1,171 @@ +// © 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 "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, + std::vector *sums) { + printf("%s:%i\n", __FILE__, __LINE__); + 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) { + const Real Mdot = History::ReduceMassAccretionRate(md); + const Real Phi = History::ReduceMagneticFluxPhi(md); + printf("Per-MD sum: Mdot: %e Phi: %e\n", Mdot, Phi); + + (*sums)[0] += Mdot; + (*sums)[1] += Phi; + } + } + printf("%s:%i\n", __FILE__, __LINE__); + + return TaskStatus::complete; +} + +TaskStatus UpdateNetFieldScaleControls(MeshData *md, const Real t, const Real dt, + const Real Mdot, const Real Phi) { + printf("%s:%i\n", __FILE__, __LINE__); + auto *pm = md->GetParentPointer(); + printf("%s:%i\n", __FILE__, __LINE__); + return TaskStatus::complete; +} + +TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, + std::vector *vals, const bool fiducial, + const Real fiducial_factor) { + const Real Mdot = (*vals)[0]; + const Real Phi = (*vals)[1]; + const Real phi = Phi / std::sqrt(Mdot); + printf("%s:%i\n", __FILE__, __LINE__); + printf("Mdot: %e Phi: %e\n", Mdot, Phi); + 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"); + + 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"); + if (fiducial) { + phi_factor = fiducial_factor; + } + + if (t >= next_dphi_dt_update_time) { + // next_dphi_dt_update_time += enforced_phi_cadence; + + // 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::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 * 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::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)); + }); + + // TODO(BRR) optionally call p2c here? + } // t >= next_dphi_dt_update_time + } // enable_phie_enforcement + + exit(-1); + return TaskStatus::complete; +} + +} // namespace fixup \ No newline at end of file diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 69eb1c86..0e59d639 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -354,15 +354,43 @@ 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; PARTHENON_REQUIRE(num_partitions == 1, "Reductions don't work for multiple partitions?"); + net_field_totals.val.resize(2); // Mdot, Phi + for (int i = 0; i < 2; i++) { + net_field_totals.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); - auto end_mods = tl.AddTask(none, fixup::EndOfStepModify, md.get(), t, dt, - stage == integrator->nstages); + auto sum_mdot_1 = tl.AddTask(none, fixup::SumMdotPhiForNetFieldScaling, md.get(), t, + &(net_field_totals.val)); + + // TODO(BRR) StartVecReduce and FinishVecReduce + + sync_region_5.AddRegionalDependencies(sync_region_5_dep_id, i, sum_mdot_1); + sync_region_5_dep_id++; + + // TODO(BRR) only do these tasks if we are actually updating the net field controls! + TaskID start_reduce_1 = + (i == 0 ? tl.AddTask(sum_mdot_1, &AllReduce>::StartReduce, + &net_field_totals, MPI_SUM) + : none); + // Test the reduction until it completes + TaskID finish_reduce_1 = tl.AddTask( + start_reduce_1, &AllReduce>::CheckReduce, &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, dt, + &(net_field_totals.val), true, 1.); + + // 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); diff --git a/src/phoebus_driver.hpp b/src/phoebus_driver.hpp index e38105bf..81f9fd17 100644 --- a/src/phoebus_driver.hpp +++ b/src/phoebus_driver.hpp @@ -45,6 +45,8 @@ class PhoebusDriver : public EvolutionDriver { AllReduce dNtot; AllReduce> particle_resolution; AllReduce particles_outstanding; + + AllReduce> net_field_totals; }; parthenon::Packages_t ProcessPackages(std::unique_ptr &pin); From 29d7281c69b3e816ab9ff83e40af50b85659b986 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Fri, 28 Apr 2023 11:44:35 -0600 Subject: [PATCH 15/25] First pass seems to work with new machinery --- src/fixup/fixup.hpp | 21 +++- src/fixup/fixup_netfield.cpp | 219 +++++++++++++++++++++++++---------- src/phoebus_driver.cpp | 55 +++++++-- src/phoebus_driver.hpp | 1 + 4 files changed, 214 insertions(+), 82 deletions(-) diff --git a/src/fixup/fixup.hpp b/src/fixup/fixup.hpp index 5b71427a..095c5d86 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 { @@ -40,13 +42,20 @@ template TaskStatus SourceFixup(T *rc); TaskStatus EndOfStepModify(MeshData *, const Real t, const Real dt, const bool last_stage); -TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, +TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, const int stage, std::vector *sums); -TaskStatus ModifyNetField(MeshData *, const Real t, const Real dt, - std::vector *vals, const bool fiducial, - const Real fiducial_factor = 1.); +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 Real Mdot, const Real Phi); + 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_netfield.cpp b/src/fixup/fixup_netfield.cpp index ad92b9ab..94811ef2 100644 --- a/src/fixup/fixup_netfield.cpp +++ b/src/fixup/fixup_netfield.cpp @@ -35,7 +35,7 @@ using robust::ratio; namespace fixup { -TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, +TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, const int stage, std::vector *sums) { printf("%s:%i\n", __FILE__, __LINE__); auto *pm = md->GetParentPointer(); @@ -43,40 +43,124 @@ TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, 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); + printf("Per-MD sum: Mdot: %e Phi: %e\n", Mdot, Phi); + + (*sums)[0] += Mdot; + (*sums)[1] += Phi; + } + } + } + printf("%s:%i\n", __FILE__, __LINE__); + + 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) { - const Real Mdot = History::ReduceMassAccretionRate(md); - const Real Phi = History::ReduceMagneticFluxPhi(md); - printf("Per-MD sum: Mdot: %e Phi: %e\n", Mdot, Phi); + TaskStatus status = net_field_totals->StartReduce(MPI_SUM); + return status; + } + } - (*sums)[0] += Mdot; - (*sums)[1] += Phi; + 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; } } - printf("%s:%i\n", __FILE__, __LINE__); return TaskStatus::complete; } TaskStatus UpdateNetFieldScaleControls(MeshData *md, const Real t, const Real dt, - const Real Mdot, const Real Phi) { + const int stage, std::vector *vals0, + std::vector *vals1) { + if (stage != 1) { + return TaskStatus::complete; + } printf("%s:%i\n", __FILE__, __LINE__); 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 + if (dphi_dt > 0.) { + dphi_dt = std::clamp((enforced_phi - (*vals0)[1]) / enforced_phi_timescale, + 3. / 4. * dphi_dt, 4. / 3. * dphi_dt); + } else { + dphi_dt = (enforced_phi - (*vals0)[1]) / enforced_phi_timescale; + } + printf("dphi_dt: %e phi_factor: %e\n", dphi_dt, phi_factor); + + 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); + printf("next_dphi_dt_update_time = %e\n", next_dphi_dt_update_time); + } + } printf("%s:%i\n", __FILE__, __LINE__); + printf("vals0: %e %e\n", (*vals0)[0], (*vals0)[1]); + printf("vals1: %e %e\n", (*vals1)[0], (*vals1)[1]); return TaskStatus::complete; } TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, - std::vector *vals, const bool fiducial, + const int stage, const bool fiducial, const Real fiducial_factor) { - const Real Mdot = (*vals)[0]; - const Real Phi = (*vals)[1]; - const Real phi = Phi / std::sqrt(Mdot); - printf("%s:%i\n", __FILE__, __LINE__); - printf("Mdot: %e Phi: %e\n", Mdot, Phi); + if (stage != 1 && fiducial) { + return TaskStatus::complete; + } auto *pm = md->GetParentPointer(); StateDescriptor *fix_pkg = pm->packages.Get("fixup").get(); @@ -87,6 +171,10 @@ TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, 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}); @@ -111,60 +199,63 @@ TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, 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"); + Real phi_factor = + fix_pkg->Param("phi_factor") * fix_pkg->Param("dphi_dt") * dt; if (fiducial) { phi_factor = fiducial_factor; } + printf("phi_factor: %e\n", phi_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::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 * 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::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 * + (-(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; + } + }); + + if (!fiducial) { + // P2C + } + } // enable_phi_enforcement - if (t >= next_dphi_dt_update_time) { - // next_dphi_dt_update_time += enforced_phi_cadence; - - // 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::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 * 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::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)); - }); - - // TODO(BRR) optionally call p2c here? - } // t >= next_dphi_dt_update_time - } // enable_phie_enforcement - - exit(-1); return TaskStatus::complete; } diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 0e59d639..9909a6aa 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -358,8 +358,10 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) { PARTHENON_REQUIRE(num_partitions == 1, "Reductions don't work for multiple partitions?"); 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; @@ -367,27 +369,56 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) { auto &md = pmesh->mesh_data.GetOrAdd(stage_name[stage], i); - auto sum_mdot_1 = tl.AddTask(none, fixup::SumMdotPhiForNetFieldScaling, md.get(), t, - &(net_field_totals.val)); + // 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 - // TODO(BRR) StartVecReduce and FinishVecReduce + // 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++; - // TODO(BRR) only do these tasks if we are actually updating the net field controls! - TaskID start_reduce_1 = - (i == 0 ? tl.AddTask(sum_mdot_1, &AllReduce>::StartReduce, - &net_field_totals, MPI_SUM) - : none); + 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, &AllReduce>::CheckReduce, &net_field_totals); + 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, dt, - &(net_field_totals.val), true, 1.); + 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); // auto end_mods = tl.AddTask(none, fixup::EndOfStepModify, md.get(), t, dt, // stage == integrator->nstages); diff --git a/src/phoebus_driver.hpp b/src/phoebus_driver.hpp index 81f9fd17..ea05d5b3 100644 --- a/src/phoebus_driver.hpp +++ b/src/phoebus_driver.hpp @@ -47,6 +47,7 @@ class PhoebusDriver : public EvolutionDriver { AllReduce particles_outstanding; AllReduce> net_field_totals; + AllReduce> net_field_totals_2; }; parthenon::Packages_t ProcessPackages(std::unique_ptr &pin); From fce5dc12b52b91fb9bc7cc1918944744bcb4ca38 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Fri, 28 Apr 2023 12:11:40 -0600 Subject: [PATCH 16/25] Better dphi_dt damper --- src/fixup/fixup_netfield.cpp | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/src/fixup/fixup_netfield.cpp b/src/fixup/fixup_netfield.cpp index 94811ef2..77b4e31b 100644 --- a/src/fixup/fixup_netfield.cpp +++ b/src/fixup/fixup_netfield.cpp @@ -130,14 +130,24 @@ TaskStatus UpdateNetFieldScaleControls(MeshData *md, const Real t, const R const Real enforced_phi_timescale = fix_pkg->Param("enforced_phi_timescale"); Real phi_factor = (*vals1)[1] - (*vals0)[1]; - // Limit change of dphi_dt - if (dphi_dt > 0.) { - dphi_dt = std::clamp((enforced_phi - (*vals0)[1]) / enforced_phi_timescale, - 3. / 4. * dphi_dt, 4. / 3. * dphi_dt); + // Limit change of dphi_dt to some factor of the original value if is non-zero + const Real fiducial_dphi = (enforced_phi - (*vals0)[1]) / 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 = (enforced_phi - (*vals0)[1]) / enforced_phi_timescale; + dphi_dt += dphi_dt_change; } - printf("dphi_dt: %e phi_factor: %e\n", dphi_dt, phi_factor); + + // if (dphi_dt > 0. && dphi_dt * fiducial_dphi > 0.) { + // dphi_dt = std::clamp(fiducial_dphi, 3. / 4. * dphi_dt, 4. / 3. * dphi_dt); + // } else { + // dphi_dt = (enforced_phi - (*vals0)[1]) / enforced_phi_timescale; + // } + // printf("dphi_dt: %e phi_factor: %e\n", dphi_dt, phi_factor); fix_pkg->UpdateParam("phi_factor", phi_factor); fix_pkg->UpdateParam("dphi_dt", dphi_dt); From fd227af55c4f5c38f02287033c8b8c831266f96b Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Fri, 28 Apr 2023 14:08:24 -0600 Subject: [PATCH 17/25] history geom bugfix --- src/analysis/history.cpp | 4 ++++ src/analysis/history_utils.hpp | 28 ++++++++++++++-------------- src/fixup/fixup_netfield.cpp | 31 ++++++++++++++++++++++++++++++- 3 files changed, 48 insertions(+), 15 deletions(-) diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index 4f1c64e5..995f49fb 100644 --- a/src/analysis/history.cpp +++ b/src/analysis/history.cpp @@ -194,6 +194,7 @@ Real ReduceMagneticFluxPhi(MeshData *md) { auto pmb = md->GetParentPointer(); auto &pars = pmb->packages.Get("geometry")->AllParams(); const Real xh = pars.Get("xh"); + printf("xh: %e\n", xh); namespace c = fluid_cons; const std::vector vars({c::bfield}); @@ -224,6 +225,9 @@ Real ReduceMagneticFluxPhi(MeshData *md) { (xh - coords.Xc<1>(i + 1)) * m) * dx2 * dx3; + // printf("[%i %i %i %i] m: %e flux: %e b: %e lapse: %e\n", b, k, j, i, m, flux, + // pack(b, cb_lo, k, j, i), geom.Lapse(CellLocation::Cent, k, j, i)); + lresult += flux; } else { lresult += 0.0; 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_netfield.cpp b/src/fixup/fixup_netfield.cpp index 77b4e31b..23894e81 100644 --- a/src/fixup/fixup_netfield.cpp +++ b/src/fixup/fixup_netfield.cpp @@ -17,6 +17,7 @@ #include #include +#include #include "analysis/history.hpp" #include "fluid/con2prim_robust.hpp" @@ -131,7 +132,8 @@ TaskStatus UpdateNetFieldScaleControls(MeshData *md, const Real t, const R 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]) / enforced_phi_timescale; + 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); @@ -216,6 +218,9 @@ TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, } printf("phi_factor: %e\n", phi_factor); + // TODO(BRR) temporary + const int rank = parthenon::Globals::my_rank; + // 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( @@ -237,6 +242,10 @@ TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, if (x < x_hyp) { A(b, j, i) = q; } + // if (i == 10) { + // + // printf("[%i] A(%i %i %i) = %e\n", rank, b, j, i, A(b, j, i)); + // } }); // Modify B field @@ -256,6 +265,26 @@ TaskStatus ModifyNetField(MeshData *md, const Real t, const Real 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)); + pack(b, pblo, k, j, i) = 1.; + pack(b, pblo + 1, k, j, i) = 1.; + + if (i == 10 && j == 84 && rank == 0 && b == 0) { + printf("b: %e %e\n", pack(b, pblo, k, j, i), pack(b, pblo + 1, k, j, i)); + printf("A1: %e %e %e %e\n", A(b, j, i), A(b, j + 1, i), A(b, j, i + 1), + A(b, j + 1, i + 1)); + printf("A2: %e %e %e %e\n", A(b, j, i), A(b, j + 1, i), A(b, j, i + 1), + A(b, j + 1, i + 1)); + printf("fac: %e dx: %e %e gammadet: %e\n", phi_factor, + coords.CellWidthFA(X2DIR, k, j, i), coords.CellWidthFA(X1DIR, k, j, i), + gammadet); + // exit(-1); + } + + // if (i == 10) { + // printf("[%i] b(%i %i %i) = %e %e\n", rank, b, j, i, pack(b, pblo, k, j, i), + // pack(b, pblo + 1, k, j, i)); + // } + SPACELOOP(ii) { pack(b, cblo + ii, k, j, i) = pack(b, pblo + ii, k, j, i) * gammadet; } From a37c8fafc38605e3265e9b3330c58c1247687492 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Fri, 28 Apr 2023 15:09:35 -0600 Subject: [PATCH 18/25] working pre-formatting --- src/analysis/history.cpp | 1 - src/fixup/fixup_netfield.cpp | 44 ------------------------------------ src/phoebus_driver.cpp | 10 ++++---- 3 files changed, 6 insertions(+), 49 deletions(-) diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index 995f49fb..efad9606 100644 --- a/src/analysis/history.cpp +++ b/src/analysis/history.cpp @@ -194,7 +194,6 @@ Real ReduceMagneticFluxPhi(MeshData *md) { auto pmb = md->GetParentPointer(); auto &pars = pmb->packages.Get("geometry")->AllParams(); const Real xh = pars.Get("xh"); - printf("xh: %e\n", xh); namespace c = fluid_cons; const std::vector vars({c::bfield}); diff --git a/src/fixup/fixup_netfield.cpp b/src/fixup/fixup_netfield.cpp index 23894e81..630b3222 100644 --- a/src/fixup/fixup_netfield.cpp +++ b/src/fixup/fixup_netfield.cpp @@ -38,7 +38,6 @@ namespace fixup { TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, const int stage, std::vector *sums) { - printf("%s:%i\n", __FILE__, __LINE__); auto *pm = md->GetParentPointer(); StateDescriptor *fix_pkg = pm->packages.Get("fixup").get(); @@ -53,14 +52,12 @@ TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, const if (t >= next_dphi_dt_update_time && t >= enforced_phi_start_time) { const Real Mdot = History::ReduceMassAccretionRate(md); const Real Phi = History::ReduceMagneticFluxPhi(md); - printf("Per-MD sum: Mdot: %e Phi: %e\n", Mdot, Phi); (*sums)[0] += Mdot; (*sums)[1] += Phi; } } } - printf("%s:%i\n", __FILE__, __LINE__); return TaskStatus::complete; } @@ -115,7 +112,6 @@ TaskStatus UpdateNetFieldScaleControls(MeshData *md, const Real t, const R if (stage != 1) { return TaskStatus::complete; } - printf("%s:%i\n", __FILE__, __LINE__); auto *pm = md->GetParentPointer(); StateDescriptor *fix_pkg = pm->packages.Get("fixup").get(); @@ -144,13 +140,6 @@ TaskStatus UpdateNetFieldScaleControls(MeshData *md, const Real t, const R dphi_dt += dphi_dt_change; } - // if (dphi_dt > 0. && dphi_dt * fiducial_dphi > 0.) { - // dphi_dt = std::clamp(fiducial_dphi, 3. / 4. * dphi_dt, 4. / 3. * dphi_dt); - // } else { - // dphi_dt = (enforced_phi - (*vals0)[1]) / enforced_phi_timescale; - // } - // printf("dphi_dt: %e phi_factor: %e\n", dphi_dt, phi_factor); - fix_pkg->UpdateParam("phi_factor", phi_factor); fix_pkg->UpdateParam("dphi_dt", dphi_dt); @@ -158,12 +147,8 @@ TaskStatus UpdateNetFieldScaleControls(MeshData *md, const Real t, const R 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); - printf("next_dphi_dt_update_time = %e\n", next_dphi_dt_update_time); } } - printf("%s:%i\n", __FILE__, __LINE__); - printf("vals0: %e %e\n", (*vals0)[0], (*vals0)[1]); - printf("vals1: %e %e\n", (*vals1)[0], (*vals1)[1]); return TaskStatus::complete; } @@ -216,7 +201,6 @@ TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, if (fiducial) { phi_factor = fiducial_factor; } - printf("phi_factor: %e\n", phi_factor); // TODO(BRR) temporary const int rank = parthenon::Globals::my_rank; @@ -242,10 +226,6 @@ TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, if (x < x_hyp) { A(b, j, i) = q; } - // if (i == 10) { - // - // printf("[%i] A(%i %i %i) = %e\n", rank, b, j, i, A(b, j, i)); - // } }); // Modify B field @@ -265,34 +245,10 @@ TaskStatus ModifyNetField(MeshData *md, const Real t, const Real 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)); - pack(b, pblo, k, j, i) = 1.; - pack(b, pblo + 1, k, j, i) = 1.; - - if (i == 10 && j == 84 && rank == 0 && b == 0) { - printf("b: %e %e\n", pack(b, pblo, k, j, i), pack(b, pblo + 1, k, j, i)); - printf("A1: %e %e %e %e\n", A(b, j, i), A(b, j + 1, i), A(b, j, i + 1), - A(b, j + 1, i + 1)); - printf("A2: %e %e %e %e\n", A(b, j, i), A(b, j + 1, i), A(b, j, i + 1), - A(b, j + 1, i + 1)); - printf("fac: %e dx: %e %e gammadet: %e\n", phi_factor, - coords.CellWidthFA(X2DIR, k, j, i), coords.CellWidthFA(X1DIR, k, j, i), - gammadet); - // exit(-1); - } - - // if (i == 10) { - // printf("[%i] b(%i %i %i) = %e %e\n", rank, b, j, i, pack(b, pblo, k, j, i), - // pack(b, pblo + 1, k, j, i)); - // } - SPACELOOP(ii) { pack(b, cblo + ii, k, j, i) = pack(b, pblo + ii, k, j, i) * gammadet; } }); - - if (!fiducial) { - // P2C - } } // enable_phi_enforcement return TaskStatus::complete; diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index 9909a6aa..2325af79 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -419,9 +419,6 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int 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); - - // 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); @@ -430,7 +427,12 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) { auto &tl = async_region_4[i]; auto &sc = pmb->meshblock_data.Get(stage_name[stage]); - auto p2c = tl.AddTask(none, fluid::PrimitiveToConserved, sc.get()); + 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); From 36dd81124f54df4af6e72fba05467c220533eeb1 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Tue, 2 May 2023 14:26:14 -0600 Subject: [PATCH 19/25] Properly deal with phoedf geom issue --- scripts/python/phoedf.py | 32 +++++++++++++++++++++++++++++--- 1 file changed, 29 insertions(+), 3 deletions(-) diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index 239291af..9c736e36 100644 --- a/scripts/python/phoedf.py +++ b/scripts/python/phoedf.py @@ -63,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) @@ -72,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( @@ -81,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] > 1: + istart = 4 + iend = -4 + else: + istart = None + iend = None + if self.flatgcov.shape[-2] > 1: + jstart = 4 + jend = -4 + else: + jstart = None + jend = None + if self.flatgcov.shape[-3] > 1: + kstart = 4 + kend = -4 + 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]) From 9f7883b26e3819152f048f24d9859ae121e939ba Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Tue, 2 May 2023 14:29:09 -0600 Subject: [PATCH 20/25] Back to main versions of submodules --- external/parthenon | 2 +- external/singularity-eos | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/external/parthenon b/external/parthenon index eb16b502..04eceab8 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit eb16b502dc03d5288ef7ea06fb3399b480e17374 +Subproject commit 04eceab86606ae6dc38d1e2fcf1d2b0def501fc7 diff --git a/external/singularity-eos b/external/singularity-eos index b4c9ce79..e8c750dd 160000 --- a/external/singularity-eos +++ b/external/singularity-eos @@ -1 +1 @@ -Subproject commit b4c9ce7962cfeb7b46585e1ef323d9bbfde1101a +Subproject commit e8c750dd8102b80b16de35a781e8f56a0b4a86b8 From dec72ded2e2ab816c67b2f7aca0ecd27b0a17aa2 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Tue, 2 May 2023 14:34:37 -0600 Subject: [PATCH 21/25] More general ghost cell removal --- scripts/python/phoedf.py | 18 +++++++++--------- scripts/python/plot_torus.py | 2 -- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/scripts/python/phoedf.py b/scripts/python/phoedf.py index 9c736e36..2d627f5c 100644 --- a/scripts/python/phoedf.py +++ b/scripts/python/phoedf.py @@ -88,21 +88,21 @@ def flatten_indices(mu, nu): 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] > 1: - istart = 4 - iend = -4 + if self.flatgcov.shape[-1] != self.Nx1: + istart = self.NGhost + iend = -self.NGhost else: istart = None iend = None - if self.flatgcov.shape[-2] > 1: - jstart = 4 - jend = -4 + if self.flatgcov.shape[-2] != self.Nx2: + jstart = self.NGhost + jend = -self.NGhost else: jstart = None jend = None - if self.flatgcov.shape[-3] > 1: - kstart = 4 - kend = -4 + if self.flatgcov.shape[-3] != self.Nx3: + kstart = self.NGhost + kend = -self.NGhost else: kstart = None kend = None 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", From 4c490e198a99822d3881725a6586cb4fa651d6ef Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Tue, 2 May 2023 14:40:25 -0600 Subject: [PATCH 22/25] Remove unused variable --- src/fixup/fixup_netfield.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/fixup/fixup_netfield.cpp b/src/fixup/fixup_netfield.cpp index fa300980..c253145f 100644 --- a/src/fixup/fixup_netfield.cpp +++ b/src/fixup/fixup_netfield.cpp @@ -202,9 +202,6 @@ TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, phi_factor = fiducial_factor; } - // TODO(BRR) temporary - const int rank = parthenon::Globals::my_rank; - // 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( @@ -254,4 +251,4 @@ TaskStatus ModifyNetField(MeshData *md, const Real t, const Real dt, return TaskStatus::complete; } -} // namespace fixup \ No newline at end of file +} // namespace fixup From 780b9ea6c860ecf30cc8c859d3a4f768f8f3b3bc Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Wed, 3 May 2023 12:24:24 -0600 Subject: [PATCH 23/25] Remove unnecessary check --- src/phoebus_driver.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/phoebus_driver.cpp b/src/phoebus_driver.cpp index d9342457..cb761d5b 100644 --- a/src/phoebus_driver.cpp +++ b/src/phoebus_driver.cpp @@ -355,8 +355,6 @@ 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; - PARTHENON_REQUIRE(num_partitions == 1, - "Reductions don't work for multiple partitions?"); net_field_totals.val.resize(2); // Mdot, Phi net_field_totals_2.val.resize(2); for (int i = 0; i < 2; i++) { From 09ed3a508032152935332fde1df148de6a35b121 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Wed, 3 May 2023 12:27:13 -0600 Subject: [PATCH 24/25] No default value for enforced phi --- src/fixup/fixup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index bdc5b273..aa795f5f 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -76,7 +76,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { 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.); + 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); From c3bfe65c79480d7e09f270e5f04717ce766f3227 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Wed, 3 May 2023 13:14:44 -0600 Subject: [PATCH 25/25] Only set other phi enforcement values if active --- src/fixup/fixup.cpp | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/src/fixup/fixup.cpp b/src/fixup/fixup.cpp index aa795f5f..1ec75a90 100644 --- a/src/fixup/fixup.cpp +++ b/src/fixup/fixup.cpp @@ -76,26 +76,27 @@ std::shared_ptr Initialize(ParameterInput *pin) { bool enable_phi_enforcement = pin->GetOrAddBoolean("fixup", "enable_phi_enforcement", false); params.Add("enable_phi_enforcement", enable_phi_enforcement); - 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); 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); } - // 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;