Skip to content

Commit

Permalink
fix: no longer pull out MeshData pointer in torus history callbacks
Browse files Browse the repository at this point in the history
  • Loading branch information
@AstroBarker
AstroBarker committed Mar 18, 2024
1 parent f877cdf commit ee1468d
Show file tree
Hide file tree
Showing 4 changed files with 22 additions and 39 deletions.
30 changes: 5 additions & 25 deletions src/analysis/history.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,15 +21,11 @@

namespace History {

Real ReduceMassAccretionRate(MeshData<Real> *md) {
Real ReduceMassAccretionRate(MeshData<Real> *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<Real>("xh");

namespace p = fluid_prim;
const std::vector<std::string> vars({p::density::name(), p::velocity::name()});

Expand Down Expand Up @@ -70,15 +66,11 @@ Real ReduceMassAccretionRate(MeshData<Real> *md) {
return result;
} // mass accretion

Real ReduceJetEnergyFlux(MeshData<Real> *md) {
Real ReduceJetEnergyFlux(MeshData<Real> *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<Real>("xh");

namespace p = fluid_prim;
const std::vector<std::string> vars(
{p::density::name(), p::bfield::name(), p::velocity::name()});
Expand All @@ -94,8 +86,6 @@ Real ReduceJetEnergyFlux(MeshData<Real> *md) {

auto geom = Geometry::GetCoordinateSystem(md);

const Real sigma_cutoff = pmesh->packages.Get("fluid")->Param<Real>("sigma_cutoff");

Real result = 0.0;
parthenon::par_reduce(
parthenon::LoopPatternMDRange(), "Phoebus History for Jet Energy Flux",
Expand Down Expand Up @@ -130,15 +120,11 @@ Real ReduceJetEnergyFlux(MeshData<Real> *md) {

} // JetEnergyFlux

Real ReduceJetMomentumFlux(MeshData<Real> *md) {
Real ReduceJetMomentumFlux(MeshData<Real> *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<Real>("xh");

namespace p = fluid_prim;
const std::vector<std::string> vars(
{p::density::name(), p::bfield::name(), p::velocity::name()});
Expand All @@ -154,8 +140,6 @@ Real ReduceJetMomentumFlux(MeshData<Real> *md) {

auto geom = Geometry::GetCoordinateSystem(md);

const Real sigma_cutoff = pmesh->packages.Get("fluid")->Param<Real>("sigma_cutoff");

Real result = 0.0;
parthenon::par_reduce(
parthenon::LoopPatternMDRange(), "Phoebus History for Jet Momentum Flux",
Expand Down Expand Up @@ -190,15 +174,11 @@ Real ReduceJetMomentumFlux(MeshData<Real> *md) {

} // ReduceJetMomentumFlux

Real ReduceMagneticFluxPhi(MeshData<Real> *md) {
Real ReduceMagneticFluxPhi(MeshData<Real> *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<Real>("xh");

namespace c = fluid_cons;
const std::vector<std::string> vars({c::bfield::name()});

Expand Down Expand Up @@ -235,7 +215,7 @@ Real ReduceMagneticFluxPhi(MeshData<Real> *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
Expand Down
8 changes: 4 additions & 4 deletions src/analysis/history.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,10 +39,10 @@ using namespace parthenon::package::prelude;

namespace History {

Real ReduceMassAccretionRate(MeshData<Real> *md);
Real ReduceJetEnergyFlux(MeshData<Real> *md);
Real ReduceJetMomentumFlux(MeshData<Real> *md);
Real ReduceMagneticFluxPhi(MeshData<Real> *md);
Real ReduceMassAccretionRate(MeshData<Real> *md, const Real xh);
Real ReduceJetEnergyFlux(MeshData<Real> *md, const Real xh, const Real sigma_cutoff);
Real ReduceJetMomentumFlux(MeshData<Real> *md, const Real xh, const Real sigma_cutoff);
Real ReduceMagneticFluxPhi(MeshData<Real> *md, const Real xh);
void ReduceLocalizationFunction(MeshData<Real> *md);

template <typename Reducer_t>
Expand Down
5 changes: 3 additions & 2 deletions src/fixup/fixup_netfield.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ TaskStatus SumMdotPhiForNetFieldScaling(MeshData<Real> *md, const Real t, const
std::vector<Real> *sums) {
Mesh *pm = md->GetMeshPointer();
StateDescriptor *fix_pkg = pm->packages.Get("fixup").get();
const Real xh = pm->packages.Get("geometry")->Param<Real>("xh");

const bool enable_phi_enforcement = fix_pkg->Param<bool>("enable_phi_enforcement");

Expand All @@ -50,8 +51,8 @@ TaskStatus SumMdotPhiForNetFieldScaling(MeshData<Real> *md, const Real t, const
const Real next_dphi_dt_update_time =
fix_pkg->Param<Real>("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;
Expand Down
18 changes: 10 additions & 8 deletions src/geometry/geometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,8 @@ std::shared_ptr<StateDescriptor> 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<Real>("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;
Expand All @@ -71,20 +73,20 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
using parthenon::HistoryOutputVar;
parthenon::HstVar_list hst_vars = {};

auto ReduceAccretionRate = [](MeshData<Real> *md) {
return ReduceMassAccretionRate(md);
auto ReduceAccretionRate = [xh](MeshData<Real> *md) {
return ReduceMassAccretionRate(md, xh);
};

auto ComputeJetEnergyFlux = [](MeshData<Real> *md) {
return ReduceJetEnergyFlux(md);
auto ComputeJetEnergyFlux = [xh, sigma_cutoff](MeshData<Real> *md) {
return ReduceJetEnergyFlux(md, xh, sigma_cutoff);
};

auto ComputeJetMomentumFlux = [](MeshData<Real> *md) {
return ReduceJetMomentumFlux(md);
auto ComputeJetMomentumFlux = [xh, sigma_cutoff](MeshData<Real> *md) {
return ReduceJetMomentumFlux(md, xh, sigma_cutoff);
};

auto ComputeMagneticFluxPhi = [](MeshData<Real> *md) {
return ReduceMagneticFluxPhi(md);
auto ComputeMagneticFluxPhi = [xh](MeshData<Real> *md) {
return ReduceMagneticFluxPhi(md, xh);
};

hst_vars.emplace_back(HistoryOutputVar(HstSum, ReduceAccretionRate, "mdot"));
Expand Down

0 comments on commit ee1468d

Please sign in to comment.