From ee1468df17b95452fecadb1f9e6605f61a4f7cb8 Mon Sep 17 00:00:00 2001 From: Brandon Barker Date: Mon, 18 Mar 2024 16:11:49 +0000 Subject: [PATCH 1/2] fix: no longer pull out MeshData pointer in torus history callbacks --- src/analysis/history.cpp | 30 +++++------------------------- src/analysis/history.hpp | 8 ++++---- src/fixup/fixup_netfield.cpp | 5 +++-- src/geometry/geometry.cpp | 18 ++++++++++-------- 4 files changed, 22 insertions(+), 39 deletions(-) diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp index d88e3d65..4ab05743 100644 --- a/src/analysis/history.cpp +++ b/src/analysis/history.cpp @@ -21,15 +21,11 @@ namespace History { -Real ReduceMassAccretionRate(MeshData *md) { +Real ReduceMassAccretionRate(MeshData *md, const Real xh) { const auto ib = md->GetBoundsI(IndexDomain::interior); const auto jb = md->GetBoundsJ(IndexDomain::interior); const auto kb = md->GetBoundsK(IndexDomain::interior); - Mesh *pmesh = md->GetMeshPointer(); - auto &pars = pmesh->packages.Get("geometry")->AllParams(); - const Real xh = pars.Get("xh"); - namespace p = fluid_prim; const std::vector vars({p::density::name(), p::velocity::name()}); @@ -70,15 +66,11 @@ Real ReduceMassAccretionRate(MeshData *md) { return result; } // mass accretion -Real ReduceJetEnergyFlux(MeshData *md) { +Real ReduceJetEnergyFlux(MeshData *md, const Real xh, const Real sigma_cutoff) { const auto ib = md->GetBoundsI(IndexDomain::interior); const auto jb = md->GetBoundsJ(IndexDomain::interior); const auto kb = md->GetBoundsK(IndexDomain::interior); - Mesh *pmesh = md->GetMeshPointer(); - auto &pars = pmesh->packages.Get("geometry")->AllParams(); - const Real xh = pars.Get("xh"); - namespace p = fluid_prim; const std::vector vars( {p::density::name(), p::bfield::name(), p::velocity::name()}); @@ -94,8 +86,6 @@ Real ReduceJetEnergyFlux(MeshData *md) { auto geom = Geometry::GetCoordinateSystem(md); - const Real sigma_cutoff = pmesh->packages.Get("fluid")->Param("sigma_cutoff"); - Real result = 0.0; parthenon::par_reduce( parthenon::LoopPatternMDRange(), "Phoebus History for Jet Energy Flux", @@ -130,15 +120,11 @@ Real ReduceJetEnergyFlux(MeshData *md) { } // JetEnergyFlux -Real ReduceJetMomentumFlux(MeshData *md) { +Real ReduceJetMomentumFlux(MeshData *md, const Real xh, const Real sigma_cutoff) { const auto ib = md->GetBoundsI(IndexDomain::interior); const auto jb = md->GetBoundsJ(IndexDomain::interior); const auto kb = md->GetBoundsK(IndexDomain::interior); - Mesh *pmesh = md->GetMeshPointer(); - auto &pars = pmesh->packages.Get("geometry")->AllParams(); - const Real xh = pars.Get("xh"); - namespace p = fluid_prim; const std::vector vars( {p::density::name(), p::bfield::name(), p::velocity::name()}); @@ -154,8 +140,6 @@ Real ReduceJetMomentumFlux(MeshData *md) { auto geom = Geometry::GetCoordinateSystem(md); - const Real sigma_cutoff = pmesh->packages.Get("fluid")->Param("sigma_cutoff"); - Real result = 0.0; parthenon::par_reduce( parthenon::LoopPatternMDRange(), "Phoebus History for Jet Momentum Flux", @@ -190,15 +174,11 @@ Real ReduceJetMomentumFlux(MeshData *md) { } // ReduceJetMomentumFlux -Real ReduceMagneticFluxPhi(MeshData *md) { +Real ReduceMagneticFluxPhi(MeshData *md, const Real xh) { const auto ib = md->GetBoundsI(IndexDomain::interior); const auto jb = md->GetBoundsJ(IndexDomain::interior); const auto kb = md->GetBoundsK(IndexDomain::interior); - Mesh *pmesh = md->GetMeshPointer(); - auto &pars = pmesh->packages.Get("geometry")->AllParams(); - const Real xh = pars.Get("xh"); - namespace c = fluid_cons; const std::vector vars({c::bfield::name()}); @@ -235,7 +215,7 @@ Real ReduceMagneticFluxPhi(MeshData *md) { }, result); return 0.5 * result; // 0.5 \int detg B^r dx2 dx3 -} // Phi +} // ReduceMagneticFluxPhi // SN analysis // ReduceLocalizationFunction is not used currently. However this function returns diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp index 8dbfb164..e74821bd 100644 --- a/src/analysis/history.hpp +++ b/src/analysis/history.hpp @@ -39,10 +39,10 @@ using namespace parthenon::package::prelude; namespace History { -Real ReduceMassAccretionRate(MeshData *md); -Real ReduceJetEnergyFlux(MeshData *md); -Real ReduceJetMomentumFlux(MeshData *md); -Real ReduceMagneticFluxPhi(MeshData *md); +Real ReduceMassAccretionRate(MeshData *md, const Real xh); +Real ReduceJetEnergyFlux(MeshData *md, const Real xh, const Real sigma_cutoff); +Real ReduceJetMomentumFlux(MeshData *md, const Real xh, const Real sigma_cutoff); +Real ReduceMagneticFluxPhi(MeshData *md, const Real xh); void ReduceLocalizationFunction(MeshData *md); template diff --git a/src/fixup/fixup_netfield.cpp b/src/fixup/fixup_netfield.cpp index 3195b0a7..2ea01e75 100644 --- a/src/fixup/fixup_netfield.cpp +++ b/src/fixup/fixup_netfield.cpp @@ -40,6 +40,7 @@ TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, const std::vector *sums) { Mesh *pm = md->GetMeshPointer(); StateDescriptor *fix_pkg = pm->packages.Get("fixup").get(); + const Real xh = pm->packages.Get("geometry")->Param("xh"); const bool enable_phi_enforcement = fix_pkg->Param("enable_phi_enforcement"); @@ -50,8 +51,8 @@ TaskStatus SumMdotPhiForNetFieldScaling(MeshData *md, const Real t, const 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); + const Real Mdot = History::ReduceMassAccretionRate(md, xh); + const Real Phi = History::ReduceMagneticFluxPhi(md, xh); (*sums)[0] += Mdot; (*sums)[1] += Phi; diff --git a/src/geometry/geometry.cpp b/src/geometry/geometry.cpp index 2ceff06d..48df5d43 100644 --- a/src/geometry/geometry.cpp +++ b/src/geometry/geometry.cpp @@ -62,6 +62,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { // Reductions const bool do_mhd = pin->GetOrAddBoolean("fluid", "mhd", false); const bool do_hydro = pin->GetBoolean("physics", "hydro"); + const Real xh = params.Get("xh"); + const Real sigma_cutoff = pin->GetOrAddReal("fluid", "sigma_cutoff", 1.0); if (params.hasKey("xh") && do_mhd && do_hydro) { auto HstSum = parthenon::UserHistoryOperation::sum; using History::ReduceJetEnergyFlux; @@ -71,20 +73,20 @@ std::shared_ptr Initialize(ParameterInput *pin) { using parthenon::HistoryOutputVar; parthenon::HstVar_list hst_vars = {}; - auto ReduceAccretionRate = [](MeshData *md) { - return ReduceMassAccretionRate(md); + auto ReduceAccretionRate = [xh](MeshData *md) { + return ReduceMassAccretionRate(md, xh); }; - auto ComputeJetEnergyFlux = [](MeshData *md) { - return ReduceJetEnergyFlux(md); + auto ComputeJetEnergyFlux = [xh, sigma_cutoff](MeshData *md) { + return ReduceJetEnergyFlux(md, xh, sigma_cutoff); }; - auto ComputeJetMomentumFlux = [](MeshData *md) { - return ReduceJetMomentumFlux(md); + auto ComputeJetMomentumFlux = [xh, sigma_cutoff](MeshData *md) { + return ReduceJetMomentumFlux(md, xh, sigma_cutoff); }; - auto ComputeMagneticFluxPhi = [](MeshData *md) { - return ReduceMagneticFluxPhi(md); + auto ComputeMagneticFluxPhi = [xh](MeshData *md) { + return ReduceMagneticFluxPhi(md, xh); }; hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceAccretionRate, "mdot")); From 4a74790dfa7563d357e7ee1b1727019098821d30 Mon Sep 17 00:00:00 2001 From: Brandon Barker Date: Mon, 18 Mar 2024 17:04:55 +0000 Subject: [PATCH 2/2] fix: pull xh in the correct place. now runs without FMKS.. --- src/geometry/geometry.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/geometry/geometry.cpp b/src/geometry/geometry.cpp index 48df5d43..746c4b45 100644 --- a/src/geometry/geometry.cpp +++ b/src/geometry/geometry.cpp @@ -62,9 +62,9 @@ std::shared_ptr Initialize(ParameterInput *pin) { // Reductions const bool do_mhd = pin->GetOrAddBoolean("fluid", "mhd", false); const bool do_hydro = pin->GetBoolean("physics", "hydro"); - const Real xh = params.Get("xh"); - const Real sigma_cutoff = pin->GetOrAddReal("fluid", "sigma_cutoff", 1.0); if (params.hasKey("xh") && do_mhd && do_hydro) { + const Real xh = params.Get("xh"); + const Real sigma_cutoff = pin->GetOrAddReal("fluid", "sigma_cutoff", 1.0); auto HstSum = parthenon::UserHistoryOperation::sum; using History::ReduceJetEnergyFlux; using History::ReduceJetMomentumFlux;