From 865ca53cdaf5cabd58c3eb8a84d624e006450490 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Fri, 14 Oct 2022 17:59:50 +0000 Subject: [PATCH 1/5] use FillPatcher class, update InterpHooks --- extern/amrex | 2 +- src/RadhydroSimulation.hpp | 48 ++++++++++++------------- src/simulation.hpp | 72 +++++++++++++++++++++++--------------- 3 files changed, 68 insertions(+), 54 deletions(-) diff --git a/extern/amrex b/extern/amrex index 1bc4e4eb5..1ad414466 160000 --- a/extern/amrex +++ b/extern/amrex @@ -1 +1 @@ -Subproject commit 1bc4e4eb5a25f4bdf9933695ead86f17dfdee9ed +Subproject commit 1ad4144668b0656d42950be92936073c64c56db7 diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp index e4c40cef4..072cbfc79 100644 --- a/src/RadhydroSimulation.hpp +++ b/src/RadhydroSimulation.hpp @@ -163,11 +163,11 @@ template class RadhydroSimulation : public AMRSimulation @@ -555,7 +555,7 @@ void RadhydroSimulation::FixupState(int lev) template void RadhydroSimulation::FillPatch(int lev, amrex::Real time, amrex::MultiFab &mf, int icomp, - int ncomp) { + int ncomp, FillPatchType fptype) { BL_PROFILE("AMRSimulation::FillPatch()"); amrex::Vector cmf; @@ -573,50 +573,48 @@ void RadhydroSimulation::FillPatch(int lev, amrex::Real time, GetData(lev - 1, time, cmf, ctime); } - FillPatchWithData(lev, time, mf, cmf, ctime, fmf, ftime, icomp, ncomp, + FillPatchWithData(lev, time, mf, cmf, ctime, fmf, ftime, icomp, ncomp, fptype, PreInterpState, PostInterpState); } template -void RadhydroSimulation::PreInterpState(amrex::FArrayBox &fab, amrex::Box const &box, - int scomp, int ncomp) +void RadhydroSimulation::PreInterpState(amrex::MultiFab &mf, int scomp, int ncomp) { BL_PROFILE("RadhydroSimulation::PreInterpState()"); - auto const &cons = fab.array(); - amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - const auto rho = cons(i, j, k, HydroSystem::density_index); - const auto px = cons(i, j, k, HydroSystem::x1Momentum_index); - const auto py = cons(i, j, k, HydroSystem::x2Momentum_index); - const auto pz = cons(i, j, k, HydroSystem::x3Momentum_index); - const auto Etot = cons(i, j, k, HydroSystem::energy_index); + auto const &cons = mf.arrays(); + amrex::ParallelFor(mf, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) { + const auto rho = cons[bx](i, j, k, HydroSystem::density_index); + const auto px = cons[bx](i, j, k, HydroSystem::x1Momentum_index); + const auto py = cons[bx](i, j, k, HydroSystem::x2Momentum_index); + const auto pz = cons[bx](i, j, k, HydroSystem::x3Momentum_index); + const auto Etot = cons[bx](i, j, k, HydroSystem::energy_index); const auto kinetic_energy = (px*px + py*py + pz*pz) / (2.0*rho); // replace hydro total energy with specific internal energy (SIE) const auto e = (Etot - kinetic_energy) / rho; - cons(i, j, k, HydroSystem::energy_index) = e; + cons[bx](i, j, k, HydroSystem::energy_index) = e; }); } template -void RadhydroSimulation::PostInterpState(amrex::FArrayBox &fab, amrex::Box const &box, - int scomp, int ncomp) +void RadhydroSimulation::PostInterpState(amrex::MultiFab &mf, int scomp, int ncomp) { BL_PROFILE("RadhydroSimulation::PostInterpState()"); - auto const &cons = fab.array(); - amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - const auto rho = cons(i, j, k, HydroSystem::density_index); - const auto px = cons(i, j, k, HydroSystem::x1Momentum_index); - const auto py = cons(i, j, k, HydroSystem::x2Momentum_index); - const auto pz = cons(i, j, k, HydroSystem::x3Momentum_index); - const auto e = cons(i, j, k, HydroSystem::energy_index); + auto const &cons = mf.arrays(); + amrex::ParallelFor(mf, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) { + const auto rho = cons[bx](i, j, k, HydroSystem::density_index); + const auto px = cons[bx](i, j, k, HydroSystem::x1Momentum_index); + const auto py = cons[bx](i, j, k, HydroSystem::x2Momentum_index); + const auto pz = cons[bx](i, j, k, HydroSystem::x3Momentum_index); + const auto e = cons[bx](i, j, k, HydroSystem::energy_index); const auto Eint = rho * e; const auto kinetic_energy = (px*px + py*py + pz*pz) / (2.0*rho); // recompute hydro total energy from Eint + KE const auto Etot = Eint + kinetic_energy; - cons(i, j, k, HydroSystem::energy_index) = Etot; + cons[bx](i, j, k, HydroSystem::energy_index) = Etot; }); } diff --git a/src/simulation.hpp b/src/simulation.hpp index 5eeb87e80..dd30eaa64 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -34,6 +34,7 @@ #include "AMReX_Extension.H" #include "AMReX_FArrayBox.H" #include "AMReX_FillPatchUtil.H" +#include "AMReX_FillPatcher.H" #include "AMReX_FluxRegister.H" #include "AMReX_GpuQualifiers.H" #include "AMReX_INT.H" @@ -75,6 +76,8 @@ using namespace conduit; using namespace ascent; #endif +enum class FillPatchType { fillpatch_class, fillpatch_function }; + // Main simulation class; solvers should inherit from this template class AMRSimulation : public amrex::AmrCore { public: @@ -165,12 +168,13 @@ template class AMRSimulation : public amrex::AmrCore { amrex::Vector &coarseTime, amrex::Vector &fineData, amrex::Vector &fineTime, int icomp, - int ncomp, PreInterpHook const &pre_interp, + int ncomp, FillPatchType fptype, + PreInterpHook const &pre_interp, PostInterpHook const &post_interp); - static void InterpHookNone(amrex::FArrayBox &fab, amrex::Box const &box, int scomp, int ncomp); + static void InterpHookNone(amrex::MultiFab &mf, int scomp, int ncomp); virtual void FillPatch(int lev, amrex::Real time, amrex::MultiFab &mf, int icomp, - int ncomp); + int ncomp, FillPatchType fptype); void FillCoarsePatch(int lev, amrex::Real time, amrex::MultiFab &mf, int icomp, int ncomp); void GetData(int lev, amrex::Real time, @@ -224,6 +228,9 @@ template class AMRSimulation : public amrex::AmrCore { // the reflux operation amrex::Vector> flux_reg_; + // This is for fillpatch during timestepping, but not for regridding. + amrex::Vector>> fillpatcher_; + // Nghost = number of ghost cells for each array int nghost_ = 4; // PPM needs nghost >= 3, PPM+flattening needs nghost >= 4 int ncomp_cc_ = 0; // = number of components (conserved variables) for each array @@ -287,6 +294,7 @@ void AMRSimulation::initialize( state_old_cc_.resize(nlevs_max); max_signal_speed_.resize(nlevs_max); flux_reg_.resize(nlevs_max + 1); + fillpatcher_.resize(nlevs_max + 1); cellUpdatesEachLevel_.resize(nlevs_max, 0); BCs_cc_ = boundaryConditions; @@ -846,6 +854,8 @@ auto AMRSimulation::timeStepWithSubcycling(int lev, amrex::Real time, AverageDownTo(lev); // average lev+1 down to lev FixupState(lev); // fix any unphysical states created by reflux or averaging + + fillpatcher_[lev+1].reset(); // because the data on lev have changed. } return stepsLeft; @@ -920,7 +930,7 @@ void AMRSimulation::RemakeLevel( amrex::MultiFab old_state(ba, dm, ncomp, nghost); amrex::MultiFab max_signal_speed(ba, dm, 1, nghost); - FillPatch(level, time, new_state, 0, ncomp); + FillPatch(level, time, new_state, 0, ncomp, FillPatchType::fillpatch_function); std::swap(new_state, state_new_cc_[level]); std::swap(old_state, state_old_cc_[level]); @@ -945,11 +955,11 @@ void AMRSimulation::ClearLevel(int level) { state_old_cc_[level].clear(); max_signal_speed_[level].clear(); flux_reg_[level].reset(nullptr); + fillpatcher_[level].reset(nullptr); } template -void AMRSimulation::InterpHookNone( - amrex::FArrayBox &fab, amrex::Box const &box, int scomp, int ncomp) +void AMRSimulation::InterpHookNone(amrex::MultiFab &mf, int scomp, int ncomp) { // do nothing } @@ -961,7 +971,7 @@ void AMRSimulation::InterpHookNone( template void AMRSimulation::FillPatch(int lev, amrex::Real time, amrex::MultiFab &mf, int icomp, - int ncomp) { + int ncomp, FillPatchType fptype) { BL_PROFILE("AMRSimulation::FillPatch()"); amrex::Vector cmf; @@ -979,7 +989,7 @@ void AMRSimulation::FillPatch(int lev, amrex::Real time, GetData(lev - 1, time, cmf, ctime); } - FillPatchWithData(lev, time, mf, cmf, ctime, fmf, ftime, icomp, ncomp, + FillPatchWithData(lev, time, mf, cmf, ctime, fmf, ftime, icomp, ncomp, fptype, InterpHookNone, InterpHookNone); } @@ -1083,7 +1093,8 @@ void AMRSimulation::fillBoundaryConditions(amrex::MultiFab &S_filled, } FillPatchWithData(lev, time, S_filled, coarseData, coarseTime, fineData, - fineTime, 0, S_filled.nComp(), pre_interp, post_interp); + fineTime, 0, S_filled.nComp(), FillPatchType::fillpatch_class, + pre_interp, post_interp); } else { // level 0 // fill internal and periodic boundaries, ignoring corners (cross=true) // (there is no performance benefit for this in practice) @@ -1122,9 +1133,21 @@ void AMRSimulation::FillPatchWithData( amrex::Vector &coarseTime, amrex::Vector &fineData, amrex::Vector &fineTime, int icomp, int ncomp, + FillPatchType fptype, PreInterpHook const &pre_interp, PostInterpHook const &post_interp) { BL_PROFILE("AMRSimulation::FillPatchWithData()"); + // use CellConservativeLinear interpolation onto fine grid + amrex::Interpolater *mapper = &amrex::cell_cons_interp; + + if (fptype == FillPatchType::fillpatch_class) { + if (fillpatcher_[lev] == nullptr) { + fillpatcher_[lev] = std::make_unique>( + boxArray(lev), DistributionMap(lev), Geom(lev), boxArray(lev - 1), DistributionMap(lev - 1), + Geom(lev - 1), mf.nGrowVect(), mf.nComp(), mapper); + } + } + // create functor to fill ghost zones at domain boundaries // (note that domain boundaries may be present at any refinement level) amrex::GpuBndryFuncFab> boundaryFunctor( @@ -1142,19 +1165,18 @@ void AMRSimulation::FillPatchWithData( coarsePhysicalBoundaryFunctor(geom[lev - 1], BCs_cc_, boundaryFunctor); - // use CellConservativeLinear interpolation onto fine grid - amrex::Interpolater *mapper = &amrex::cell_cons_interp; - // amrex::MFInterpolater *mapper = &amrex::mf_cell_cons_interp; - // amrex::MFInterpolater *mapper = &amrex::mf_pc_interp; - // copies interior zones, fills ghost zones with space-time interpolated // data - amrex::FillPatchTwoLevels(mf, time, coarseData, coarseTime, fineData, - fineTime, 0, icomp, ncomp, geom[lev - 1], - geom[lev], coarsePhysicalBoundaryFunctor, 0, - finePhysicalBoundaryFunctor, 0, refRatio(lev - 1), - mapper, BCs_cc_, 0, - pre_interp, post_interp); + if (fptype == FillPatchType::fillpatch_class) { + fillpatcher_[lev]->fill(mf, mf.nGrowVect(), time, coarseData, coarseTime, fineData, fineTime, 0, icomp, + ncomp, coarsePhysicalBoundaryFunctor, 0, finePhysicalBoundaryFunctor, 0, BCs_cc_, 0, + pre_interp, post_interp); + } else { + amrex::FillPatchTwoLevels(mf, time, coarseData, coarseTime, fineData, fineTime, 0, icomp, ncomp, + geom[lev - 1], geom[lev], coarsePhysicalBoundaryFunctor, 0, + finePhysicalBoundaryFunctor, 0, refRatio(lev - 1), mapper, BCs_cc_, 0, pre_interp, + post_interp); + } } } @@ -1206,16 +1228,10 @@ void AMRSimulation::GetData(int lev, amrex::Real time, data.clear(); datatime.clear(); - const amrex::Real teps = - (tNew_[lev] - tOld_[lev]) * 1.e-3; // generous roundoff error threshold - - if (time > tNew_[lev] - teps && - time < tNew_[lev] + teps) { // if time == tNew_[lev] within roundoff + if (amrex::almostEqual(time, tNew_[lev], 5)) { // if time == tNew_[lev] within roundoff data.push_back(&state_new_cc_[lev]); datatime.push_back(tNew_[lev]); - } else if (time > tOld_[lev] - teps && - time < - tOld_[lev] + teps) { // if time == tOld_[lev] within roundoff + } else if (amrex::almostEqual(time, tOld_[lev], 5)) { // if time == tOld_[lev] within roundoff data.push_back(&state_old_cc_[lev]); datatime.push_back(tOld_[lev]); } else { // otherwise return both old and new states for interpolation From caf046670f19d0e9d1de3cbbdd8e127b261c5f2e Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Mon, 17 Oct 2022 14:48:43 +0000 Subject: [PATCH 2/5] test FillPatcher PR --- extern/amrex | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/amrex b/extern/amrex index 1ad414466..fd083c1c0 160000 --- a/extern/amrex +++ b/extern/amrex @@ -1 +1 @@ -Subproject commit 1ad4144668b0656d42950be92936073c64c56db7 +Subproject commit fd083c1c0de5658e7243e09d14c3f3961bd5d977 From 9dc8867eb53ec6c736e83d55c13405a6e20697d9 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Mon, 17 Oct 2022 15:27:14 +0000 Subject: [PATCH 3/5] use FillPatcher class --- src/simulation.hpp | 75 ++++++++++++++++++++++++---------------------- 1 file changed, 39 insertions(+), 36 deletions(-) diff --git a/src/simulation.hpp b/src/simulation.hpp index f8f5d5e47..b45b9657e 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -161,7 +161,8 @@ template class AMRSimulation : public amrex::AmrCore { void fillBoundaryConditions(amrex::MultiFab &S_filled, amrex::MultiFab &state, int lev, amrex::Real time, PreInterpHook const &pre_interp, - PostInterpHook const&post_interp); + PostInterpHook const&post_interp, + FillPatchType fptype = FillPatchType::fillpatch_class); template void FillPatchWithData(int lev, amrex::Real time, amrex::MultiFab &mf, @@ -971,6 +972,31 @@ void AMRSimulation::InterpHookNone(amrex::MultiFab &mf, int scomp, in // do nothing } + +template struct setBoundaryFunctor { + AMREX_GPU_DEVICE void + operator()(const amrex::IntVect &iv, amrex::Array4 const &dest, + const int &dcomp, const int &numcomp, + amrex::GeometryData const &geom, const amrex::Real &time, + const amrex::BCRec *bcr, int bcomp, const int &orig_comp) const { + AMRSimulation::setCustomBoundaryConditions( + iv, dest, dcomp, numcomp, geom, time, bcr, bcomp, orig_comp); + } +}; + +template +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +AMRSimulation::setCustomBoundaryConditions( + const amrex::IntVect &iv, amrex::Array4 const &dest, int dcomp, + int numcomp, amrex::GeometryData const &geom, const amrex::Real time, + const amrex::BCRec *bcr, int bcomp, int orig_comp) { + // user should implement if needed using template specialization + // (This is only called when amrex::BCType::ext_dir is set for a given + // boundary.) + + // set boundary condition for cell 'iv' +} + // Compute a new multifab 'mf' by copying in state from valid region and filling // ghost cells // NOTE: This implementation is only used by AdvectionSimulation. @@ -1031,37 +1057,12 @@ void AMRSimulation::MakeNewLevelFromScratch( // check that state_new_cc_[lev] is properly filled AMREX_ALWAYS_ASSERT(!state_new_cc_[level].contains_nan(0, ncomp)); - // fill ghost zones + // fill ghost zones (needed for some refinement criteria) fillBoundaryConditions(state_new_cc_[level], state_new_cc_[level], level, time, - InterpHookNone, InterpHookNone); + InterpHookNone, InterpHookNone, FillPatchType::fillpatch_function); // copy to state_old_cc_ (including ghost zones) - state_old_cc_[level].ParallelCopy(state_new_cc_[level], 0, 0, ncomp, nghost, - nghost); -} - -template struct setBoundaryFunctor { - AMREX_GPU_DEVICE void - operator()(const amrex::IntVect &iv, amrex::Array4 const &dest, - const int &dcomp, const int &numcomp, - amrex::GeometryData const &geom, const amrex::Real &time, - const amrex::BCRec *bcr, int bcomp, const int &orig_comp) const { - AMRSimulation::setCustomBoundaryConditions( - iv, dest, dcomp, numcomp, geom, time, bcr, bcomp, orig_comp); - } -}; - -template -AMREX_GPU_DEVICE AMREX_FORCE_INLINE void -AMRSimulation::setCustomBoundaryConditions( - const amrex::IntVect &iv, amrex::Array4 const &dest, int dcomp, - int numcomp, amrex::GeometryData const &geom, const amrex::Real time, - const amrex::BCRec *bcr, int bcomp, int orig_comp) { - // user should implement if needed using template specialization - // (This is only called when amrex::BCType::ext_dir is set for a given - // boundary.) - - // set boundary condition for cell 'iv' + state_old_cc_[level].ParallelCopy(state_new_cc_[level], 0, 0, ncomp, nghost, nghost); } template @@ -1071,7 +1072,8 @@ void AMRSimulation::fillBoundaryConditions(amrex::MultiFab &S_filled, int const lev, amrex::Real const time, PreInterpHook const &pre_interp, - PostInterpHook const &post_interp) { + PostInterpHook const &post_interp, + FillPatchType fptype) { BL_PROFILE("AMRSimulation::fillBoundaryConditions()"); // On a single level, any periodic boundaries are filled first @@ -1100,7 +1102,7 @@ void AMRSimulation::fillBoundaryConditions(amrex::MultiFab &S_filled, } FillPatchWithData(lev, time, S_filled, coarseData, coarseTime, fineData, - fineTime, 0, S_filled.nComp(), FillPatchType::fillpatch_class, + fineTime, 0, S_filled.nComp(), fptype, pre_interp, post_interp); } else { // level 0 // fill internal and periodic boundaries, ignoring corners (cross=true) @@ -1150,8 +1152,8 @@ void AMRSimulation::FillPatchWithData( if (fptype == FillPatchType::fillpatch_class) { if (fillpatcher_[lev] == nullptr) { fillpatcher_[lev] = std::make_unique>( - boxArray(lev), DistributionMap(lev), Geom(lev), boxArray(lev - 1), DistributionMap(lev - 1), - Geom(lev - 1), mf.nGrowVect(), mf.nComp(), mapper); + grids[lev], dmap[lev], geom[lev], grids[lev - 1], dmap[lev - 1], geom[lev - 1], + mf.nGrowVect(), mf.nComp(), mapper); } } @@ -1175,9 +1177,10 @@ void AMRSimulation::FillPatchWithData( // copies interior zones, fills ghost zones with space-time interpolated // data if (fptype == FillPatchType::fillpatch_class) { - fillpatcher_[lev]->fill(mf, mf.nGrowVect(), time, coarseData, coarseTime, fineData, fineTime, 0, icomp, - ncomp, coarsePhysicalBoundaryFunctor, 0, finePhysicalBoundaryFunctor, 0, BCs_cc_, 0, - pre_interp, post_interp); + fillpatcher_[lev]->fill(mf, mf.nGrowVect(), time, + coarseData, coarseTime, fineData, fineTime, 0, icomp, ncomp, + coarsePhysicalBoundaryFunctor, 0, finePhysicalBoundaryFunctor, 0, + BCs_cc_, 0, pre_interp, post_interp); } else { amrex::FillPatchTwoLevels(mf, time, coarseData, coarseTime, fineData, fineTime, 0, icomp, ncomp, geom[lev - 1], geom[lev], coarsePhysicalBoundaryFunctor, 0, From a201da7a5efadc045090c1d4bbcc2acd96d9091b Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Wed, 26 Oct 2022 19:59:02 +0000 Subject: [PATCH 4/5] update amrex submodule --- extern/amrex | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/amrex b/extern/amrex index fd083c1c0..91f8a63f2 160000 --- a/extern/amrex +++ b/extern/amrex @@ -1 +1 @@ -Subproject commit fd083c1c0de5658e7243e09d14c3f3961bd5d977 +Subproject commit 91f8a63f2187a9565c76e96ff0d92b7e16cb154c From a43cb425a68af6237746d99570d20121c8852b9d Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Wed, 26 Oct 2022 22:05:24 +0000 Subject: [PATCH 5/5] update amrex submodule --- extern/amrex | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/amrex b/extern/amrex index 91f8a63f2..b3e0a62ba 160000 --- a/extern/amrex +++ b/extern/amrex @@ -1 +1 @@ -Subproject commit 91f8a63f2187a9565c76e96ff0d92b7e16cb154c +Subproject commit b3e0a62ba4d8c66b7cc40ab439b94835a5f4247c