From 5ba8f0dc6fe4c9b59b8b5dbf67a168b7ed261e9e Mon Sep 17 00:00:00 2001 From: Ben Prather Date: Thu, 17 Aug 2023 22:49:21 -0600 Subject: [PATCH] Fixes for the obvious reduction bugs, & a new name collision with package names --- kharma/emhd/emhd.cpp | 11 ++++++----- kharma/grmhd/grmhd.hpp | 2 +- kharma/prob/post_initialize.cpp | 22 ---------------------- kharma/reductions/reductions_impl.hpp | 14 +++++++------- 4 files changed, 14 insertions(+), 35 deletions(-) diff --git a/kharma/emhd/emhd.cpp b/kharma/emhd/emhd.cpp index 34153ce0..3150a198 100644 --- a/kharma/emhd/emhd.cpp +++ b/kharma/emhd/emhd.cpp @@ -122,15 +122,16 @@ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr

Initialize(ParameterInput *pin, std::shared_ptr

GetBlockPointer(); // PackIndexMap prims_map, cons_map; -// auto U_E = rc->PackVariables(std::vector{Metadata::GetUserFlag("EMHD"), Metadata::Conserved}, cons_map); +// auto U_E = rc->PackVariables(std::vector{Metadata::GetUserFlag("EMHDVar"), Metadata::Conserved}, cons_map); // auto P = rc->PackVariables(std::vector{Metadata::GetUserFlag("Primitive")}, prims_map); // const VarMap m_p(prims_map, false), m_u(cons_map, true); @@ -216,7 +217,7 @@ void BlockPtoU(MeshBlockData *rc, IndexDomain domain, bool coarse) auto pmb = rc->GetBlockPointer(); PackIndexMap prims_map, cons_map; - auto U_E = rc->PackVariables(std::vector{Metadata::GetUserFlag("EMHD"), Metadata::Conserved}, cons_map); + auto U_E = rc->PackVariables(std::vector{Metadata::GetUserFlag("EMHDVar"), Metadata::Conserved}, cons_map); auto P = rc->PackVariables(std::vector{Metadata::GetUserFlag("Primitive")}, prims_map); const VarMap m_p(prims_map, false), m_u(cons_map, true); diff --git a/kharma/grmhd/grmhd.hpp b/kharma/grmhd/grmhd.hpp index 51736a2f..53878118 100644 --- a/kharma/grmhd/grmhd.hpp +++ b/kharma/grmhd/grmhd.hpp @@ -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 *rc); } diff --git a/kharma/prob/post_initialize.cpp b/kharma/prob/post_initialize.cpp index 6b338019..5944bbe1 100644 --- a/kharma/prob/post_initialize.cpp +++ b/kharma/prob/post_initialize.cpp @@ -92,8 +92,6 @@ void KHARMA::SeedAndNormalizeB(ParameterInput *pin, std::shared_ptrpackages.AllPackages().count("B_CD"); const int verbose = pmesh->packages.Get("Globals")->Param("verbose"); - fprintf(stderr, "0.5"); - Flag("SeedBField"); // Seed the magnetic field on each block for (auto &pmb : pmesh->block_list) { @@ -108,8 +106,6 @@ void KHARMA::SeedAndNormalizeB(ParameterInput *pin, std::shared_ptrGetString("parthenon/job", "problem_id"); @@ -127,11 +123,8 @@ void KHARMA::SeedAndNormalizeB(ParameterInput *pin, std::shared_ptrmesh_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. @@ -222,8 +210,6 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart) 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")) { @@ -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)) { @@ -245,8 +229,6 @@ 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, @@ -254,8 +236,6 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart) 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. @@ -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); diff --git a/kharma/reductions/reductions_impl.hpp b/kharma/reductions/reductions_impl.hpp index 0575cd88..df42b1b8 100644 --- a/kharma/reductions/reductions_impl.hpp +++ b/kharma/reductions/reductions_impl.hpp @@ -73,10 +73,10 @@ void Reductions::Start(MeshData *md, int channel, T val, MPI_Op op) auto& pars = md->GetMeshPointer()->packages.Get("Reductions")->AllParams(); auto *reduce_pool = pars.GetMutable>>(pool_name); while (reduce_pool->size() <= channel) reduce_pool->push_back(Reduce()); - 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 void Reductions::StartToAll(MeshData *md, int channel, T val, MPI_Op op) @@ -86,10 +86,10 @@ void Reductions::StartToAll(MeshData *md, int channel, T val, MPI_Op op) auto& pars = md->GetMeshPointer()->packages.Get("Reductions")->AllParams(); auto *allreduce_pool = pars.GetMutable>>(pool_name); while (allreduce_pool->size() <= channel) allreduce_pool->push_back(AllReduce()); - 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 @@ -111,7 +111,7 @@ T Reductions::CheckOnAll(MeshData *md, int channel) // Get the relevant reducer and result const std::string pool_name = GetPoolName(); auto& pars = md->GetMeshPointer()->packages.Get("Reductions")->AllParams(); - auto *reduce_pool = pars.GetMutable>>(pool_name); + auto *reduce_pool = pars.GetMutable>>(pool_name); auto& reducer = (*reduce_pool)[channel]; while (reducer.CheckReduce() == TaskStatus::incomplete);