From 2a3827f2826f5e1f31c47b68638caa4106467ace Mon Sep 17 00:00:00 2001 From: Ben Prather Date: Wed, 22 Nov 2023 16:27:17 -0700 Subject: [PATCH] Update imex driver and boundaries to keep divB under SMR --- kharma/boundaries/boundaries.cpp | 19 ++++++++++++------- kharma/driver/imex_step.cpp | 9 +++++---- kharma/driver/kharma_step.cpp | 4 ++-- 3 files changed, 19 insertions(+), 13 deletions(-) diff --git a/kharma/boundaries/boundaries.cpp b/kharma/boundaries/boundaries.cpp index cd75ee9f..0edd8fbb 100644 --- a/kharma/boundaries/boundaries.cpp +++ b/kharma/boundaries/boundaries.cpp @@ -265,9 +265,12 @@ void KBoundaries::ApplyBoundary(std::shared_ptr> &rc, IndexD EndFlag(); // If we're syncing EMFs and in spherical, explicitly zero polar faces - // TODO allow any other face? + // Since we manipulate the j coord, we'd overstep coarse bufs auto& emfpack = rc->PackVariables(std::vector{"B_CT.emf"}); - if (bdir == X2DIR && pmb->coords.coords.is_spherical() && emfpack.GetDim(4) > 0) { + if (bdir == X2DIR && + pmb->coords.coords.is_spherical() && + emfpack.GetDim(4) > 0 && + !coarse) { Flag("BoundaryEdge_"+bname); for (TE el : {TE::E1, TE::E3}) { int off = (binner) ? 1 : -1; @@ -283,7 +286,9 @@ void KBoundaries::ApplyBoundary(std::shared_ptr> &rc, IndexD // Zero/invert X2 faces at polar X2 boundary auto fpack = rc->PackVariables({Metadata::Face, Metadata::FillGhost}); - if (bdir == X2DIR && pmb->coords.coords.is_spherical() && fpack.GetDim(4) > 0) { + if (bdir == X2DIR && + pmb->coords.coords.is_spherical() && + fpack.GetDim(4) > 0) { Flag("BoundaryFace_"+bname); // Zero face fluxes auto b = KDomain::GetRange(rc, domain, coarse); @@ -298,7 +303,7 @@ void KBoundaries::ApplyBoundary(std::shared_ptr> &rc, IndexD pmb->par_for_bndry( "invert_F2_" + bname, IndexRange{0, fpack.GetDim(4)-1}, domain, F2, coarse, KOKKOS_LAMBDA (const int &v, const int &k, const int &j, const int &i) { - fpack(F2, 0, k, j, i) *= -1; + fpack(F2, v, k, j, i) *= -1; } ); EndFlag(); @@ -333,7 +338,7 @@ void KBoundaries::ApplyBoundary(std::shared_ptr> &rc, IndexD const auto &range = (bdir == 1) ? bounds.GetBoundsI(IndexDomain::interior) : (bdir == 2 ? bounds.GetBoundsJ(IndexDomain::interior) : bounds.GetBoundsK(IndexDomain::interior)); - const int ref = BoundaryIsInner(domain) ? range.s : range.e; + const int ref = binner ? range.s : range.e; pmb->par_for_bndry( "outflow_EMHD", IndexRange{0,EMHDg.GetDim(4)-1}, domain, CC, coarse, KOKKOS_LAMBDA (const int &v, const int &k, const int &j, const int &i) { @@ -384,9 +389,9 @@ void KBoundaries::ApplyBoundary(std::shared_ptr> &rc, IndexD // TODO there should be a set of B field wrappers that dispatch this auto pkgs = pmb->packages.AllPackages(); if (pkgs.count("B_FluxCT")) { - B_FluxCT::BlockUtoP(rc.get(), IndexDomain::entire); + B_FluxCT::BlockUtoP(rc.get(), domain, coarse); } else if (pkgs.count("B_CT")) { - B_CT::BlockUtoP(rc.get(), IndexDomain::entire); + B_CT::BlockUtoP(rc.get(), domain, coarse); } Flux::BlockPtoU(rc.get(), domain, coarse); } else { diff --git a/kharma/driver/imex_step.cpp b/kharma/driver/imex_step.cpp index 0b179327..0ec5815c 100644 --- a/kharma/driver/imex_step.cpp +++ b/kharma/driver/imex_step.cpp @@ -153,15 +153,16 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta // If we're in AMR, correct fluxes from neighbors auto t_flux_bounds = t_fix_flux; if (pmesh->multilevel || use_b_ct) { - auto t_load_send_flux = tl.AddTask(t_fix_flux, parthenon::LoadAndSendFluxCorrections, md_sub_step_init); - auto t_recv_flux = tl.AddTask(t_load_send_flux, parthenon::ReceiveFluxCorrections, md_sub_step_init); - t_flux_bounds = tl.AddTask(t_recv_flux, parthenon::SetFluxCorrections, md_sub_step_init); + auto t_emf = t_flux_bounds; if (use_b_ct) { // Pull out a container of only EMF to synchronize auto &md_emf_only = pmesh->mesh_data.AddShallow("EMF", std::vector{"B_CT.emf"}); // TODO this gets weird if we partition auto t_emf_local = tl.AddTask(t_flux_bounds, B_CT::CalculateEMF, md_sub_step_init.get()); - t_flux_bounds = KHARMADriver::AddBoundarySync(t_emf_local, tl, md_sub_step_init); + t_emf = KHARMADriver::AddBoundarySync(t_emf_local, tl, md_emf_only); } + auto t_load_send_flux = tl.AddTask(t_emf, parthenon::LoadAndSendFluxCorrections, md_sub_step_init); + auto t_recv_flux = tl.AddTask(t_load_send_flux, parthenon::ReceiveFluxCorrections, md_sub_step_init); + t_flux_bounds = tl.AddTask(t_recv_flux, parthenon::SetFluxCorrections, md_sub_step_init); } // Apply the fluxes to calculate a change in cell-centered values "md_flux_src" diff --git a/kharma/driver/kharma_step.cpp b/kharma/driver/kharma_step.cpp index 58c1dea2..f59e7976 100644 --- a/kharma/driver/kharma_step.cpp +++ b/kharma/driver/kharma_step.cpp @@ -165,9 +165,9 @@ TaskCollection KHARMADriver::MakeDefaultTaskCollection(BlockList_t &blocks, int auto t_emf = t_flux_bounds; if (use_b_ct) { // Pull out a container of only EMF to synchronize - auto &md_b_ct = pmesh->mesh_data.AddShallow("B_CT", std::vector{"B_CT.emf"}); // TODO this gets weird if we partition + auto &md_emf_only = pmesh->mesh_data.AddShallow("EMF", std::vector{"B_CT.emf"}); // TODO this gets weird if we partition auto t_emf_local = tl.AddTask(t_flux_bounds, B_CT::CalculateEMF, md_sub_step_init.get()); - t_emf = KHARMADriver::AddBoundarySync(t_emf_local, tl, md_b_ct); + t_emf = KHARMADriver::AddBoundarySync(t_emf_local, tl, md_emf_only); } auto t_load_send_flux = tl.AddTask(t_emf, parthenon::LoadAndSendFluxCorrections, md_sub_step_init); auto t_recv_flux = tl.AddTask(t_load_send_flux, parthenon::ReceiveFluxCorrections, md_sub_step_init);