Skip to content

Commit

Permalink
Fixes for the obvious reduction bugs, & a new name collision with pac…
Browse files Browse the repository at this point in the history 
…kage names
  • Loading branch information
@bprather
bprather committed Aug 18, 2023
1 parent 2f9fbdd commit 5ba8f0d
Show file tree
Hide file tree
Showing 4 changed files with 14 additions and 35 deletions.
11 changes: 6 additions & 5 deletions kharma/emhd/emhd.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,15 +122,16 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
// Only enable limits internally if we're actually doing EMHD
params.Add("enable_emhd_limits", enable_emhd_limits);

Metadata::AddUserFlag("EMHD");

Metadata::AddUserFlag("EMHDVar");

// General options for primitive and conserved scalar variables in ImEx driver
// EMHD is supported only with imex driver and implicit evolution,
// synchronizing primitive variables
Metadata m_con = Metadata({Metadata::Real, Metadata::Cell, Metadata::Independent, Metadata::GetUserFlag("Implicit"),
Metadata::WithFluxes, Metadata::Conserved, Metadata::GetUserFlag("EMHD")});
Metadata::WithFluxes, Metadata::Conserved, Metadata::GetUserFlag("EMHDVar")});
Metadata m_prim = Metadata({Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::GetUserFlag("Implicit"),
Metadata::Restart, Metadata::FillGhost, Metadata::GetUserFlag("Primitive"), Metadata::GetUserFlag("EMHD")});
Metadata::Restart, Metadata::FillGhost, Metadata::GetUserFlag("Primitive"), Metadata::GetUserFlag("EMHDVar")});

// Heat conduction
if (conduction) {
Expand Down Expand Up @@ -184,7 +185,7 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
// auto pmb = rc->GetBlockPointer();

// PackIndexMap prims_map, cons_map;
// auto U_E = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("EMHD"), Metadata::Conserved}, cons_map);
// auto U_E = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("EMHDVar"), Metadata::Conserved}, cons_map);
// auto P = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("Primitive")}, prims_map);
// const VarMap m_p(prims_map, false), m_u(cons_map, true);

Expand Down Expand Up @@ -216,7 +217,7 @@ void BlockPtoU(MeshBlockData<Real> *rc, IndexDomain domain, bool coarse)
auto pmb = rc->GetBlockPointer();

PackIndexMap prims_map, cons_map;
auto U_E = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("EMHD"), Metadata::Conserved}, cons_map);
auto U_E = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("EMHDVar"), Metadata::Conserved}, cons_map);
auto P = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("Primitive")}, prims_map);
const VarMap m_p(prims_map, false), m_u(cons_map, true);

Expand Down
2 changes: 1 addition & 1 deletion kharma/grmhd/grmhd.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ void FillOutput(MeshBlock *pmb, ParameterInput *pin);

/**
* Diagnostics performed after each step.
* Currently finds any negative flags or 0/NaN values in ctop
* Currently just looks for negative density/internal energy
*/
TaskStatus PostStepDiagnostics(const SimTime& tm, MeshData<Real> *rc);
}
22 changes: 0 additions & 22 deletions kharma/prob/post_initialize.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,8 +92,6 @@ void KHARMA::SeedAndNormalizeB(ParameterInput *pin, std::shared_ptr<MeshData<Rea
const bool use_b_cd = pmesh->packages.AllPackages().count("B_CD");
const int verbose = pmesh->packages.Get("Globals")->Param<int>("verbose");

fprintf(stderr, "0.5");

Flag("SeedBField");
// Seed the magnetic field on each block
for (auto &pmb : pmesh->block_list) {
Expand All @@ -108,8 +106,6 @@ void KHARMA::SeedAndNormalizeB(ParameterInput *pin, std::shared_ptr<MeshData<Rea
}
EndFlag();

fprintf(stderr, "0.9");

// Then, if we're in a torus problem or we explicitly ask for it,
// normalize the magnetic field according to the density
auto prob = pin->GetString("parthenon/job", "problem_id");
Expand All @@ -127,11 +123,8 @@ void KHARMA::SeedAndNormalizeB(ParameterInput *pin, std::shared_ptr<MeshData<Rea
// Calculate current beta_min value
Real bsq_max, p_max, beta_min;
if (beta_calc_legacy) {
fprintf(stderr, "1");
bsq_max = MPIReduce_once(MaxBsq(md.get()), MPI_MAX);
fprintf(stderr, "2");
p_max = MPIReduce_once(MaxPressure(md.get()), MPI_MAX);
fprintf(stderr, "3");
beta_min = p_max / (0.5 * bsq_max);
} else {
beta_min = MPIReduce_once(MinBeta(md.get()), MPI_MIN);
Expand Down Expand Up @@ -194,36 +187,29 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
// If your problem requires custom boundary conditions, these should be implemented
// with the problem and assigned to the relevant functions in the "Boundaries" package.

fprintf(stderr, "0.0");
auto &md = pmesh->mesh_data.Get();

auto& pkgs = pmesh->packages.AllPackages();

fprintf(stderr, "0.1");
// Magnetic field operations
if (pin->GetString("b_field", "solver") != "none") {
// If we need to seed a field based on the problem's fluid initialization...
if (pin->GetOrAddString("b_field", "type", "none") != "none" && !is_restart) {
// B field init is not stencil-1, needs boundaries sync'd.
// FreezeDirichlet ensures any Dirichlet conditions aren't overwritten by zeros
fprintf(stderr, "0.2");
KBoundaries::FreezeDirichlet(md);
KHARMADriver::SyncAllBounds(md);

fprintf(stderr, "0.3");
// Then init B field on each block...
KHARMA::SeedAndNormalizeB(pin, md);
}
fprintf(stderr, "4");

// Regardless, if evolving a field we should print max(divB)
// divB is not stencil-1 and we may not have run the above.
// If we did, we still need another sync, so it works out
KBoundaries::FreezeDirichlet(md);
KHARMADriver::SyncAllBounds(md);

fprintf(stderr, "5");

if (pkgs.count("B_FluxCT")) {
B_FluxCT::PrintGlobalMaxDivB(md.get());
} else if (pkgs.count("B_CT")) {
Expand All @@ -233,8 +219,6 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
}
}

fprintf(stderr, "6");

// Add any hotspots.
// Note any other modifications made when restarting should be made around here
if (pin->GetOrAddBoolean("blob", "add_blob", false)) {
Expand All @@ -245,17 +229,13 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
}
}

fprintf(stderr, "7");

// Any extra cleanup & init especially when restarting
if (is_restart) {
// Parthenon restores all parameters (global vars) when restarting,
// but KHARMA needs a few (currently one) reset instead
KHARMA::ResetGlobals(pin, pmesh);
}

fprintf(stderr, "8");

// Clean the B field if we've introduced a divergence somewhere
// We call this function any time the package is loaded:
// if we decided to load it in kharma.cpp, we need to clean.
Expand All @@ -270,8 +250,6 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
B_Cleanup::CleanupDivergence(md);
}

fprintf(stderr, "9");

// Finally, synchronize boundary values.
// Freeze any Dirichlet physical boundaries as they are now, after cleanup/sync/etc.
KBoundaries::FreezeDirichlet(md);
Expand Down
14 changes: 7 additions & 7 deletions kharma/reductions/reductions_impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,10 +73,10 @@ void Reductions::Start(MeshData<Real> *md, int channel, T val, MPI_Op op)
auto& pars = md->GetMeshPointer()->packages.Get("Reductions")->AllParams();
auto *reduce_pool = pars.GetMutable<std::vector<Reduce<T>>>(pool_name);
while (reduce_pool->size() <= channel) reduce_pool->push_back(Reduce<T>());
auto& vector_int_reduce = (*reduce_pool)[channel];
auto& reduce = (*reduce_pool)[channel];
// Fill with flags
vector_int_reduce.val = val;
vector_int_reduce.StartReduce(0, op);
reduce.val = val;
reduce.StartReduce(0, op);
}
template<typename T>
void Reductions::StartToAll(MeshData<Real> *md, int channel, T val, MPI_Op op)
Expand All @@ -86,10 +86,10 @@ void Reductions::StartToAll(MeshData<Real> *md, int channel, T val, MPI_Op op)
auto& pars = md->GetMeshPointer()->packages.Get("Reductions")->AllParams();
auto *allreduce_pool = pars.GetMutable<std::vector<AllReduce<T>>>(pool_name);
while (allreduce_pool->size() <= channel) allreduce_pool->push_back(AllReduce<T>());
auto& vector_int_reduce = (*allreduce_pool)[channel];
auto& reduce = (*allreduce_pool)[channel];
// Fill with flags
vector_int_reduce.val = val;
vector_int_reduce.StartReduce(op);
reduce.val = val;
reduce.StartReduce(op);
}

// MPI reduction checks
Expand All @@ -111,7 +111,7 @@ T Reductions::CheckOnAll(MeshData<Real> *md, int channel)
// Get the relevant reducer and result
const std::string pool_name = GetPoolName<T, true>();
auto& pars = md->GetMeshPointer()->packages.Get("Reductions")->AllParams();
auto *reduce_pool = pars.GetMutable<std::vector<Reduce<T>>>(pool_name);
auto *reduce_pool = pars.GetMutable<std::vector<AllReduce<T>>>(pool_name);
auto& reducer = (*reduce_pool)[channel];

while (reducer.CheckReduce() == TaskStatus::incomplete);
Expand Down

0 comments on commit 5ba8f0d

Please sign in to comment.