diff --git a/Src/Base/AMReX_FabArrayBase.H b/Src/Base/AMReX_FabArrayBase.H index 002d2a3586f..d45a416c5a9 100644 --- a/Src/Base/AMReX_FabArrayBase.H +++ b/Src/Base/AMReX_FabArrayBase.H @@ -143,6 +143,7 @@ public: */ bool is_cell_centered () const noexcept; + void setMultiGhost(bool a_multi_ghost) {m_multi_ghost = a_multi_ghost;} // These are provided for convenience to keep track of how many // ghost cells are up to date. The number of filled ghost cells @@ -466,6 +467,7 @@ public: int n_comp; mutable BDKey m_bdkey; IntVect n_filled; // Note that IntVect is zero by default. + bool m_multi_ghost = false; // // Tiling @@ -501,7 +503,7 @@ public: { FB (const FabArrayBase& fa, const IntVect& nghost, bool cross, const Periodicity& period, - bool enforce_periodicity_only); + bool enforce_periodicity_only, bool multi_ghost = false); ~FB (); IndexType m_typ; @@ -512,6 +514,7 @@ public: Periodicity m_period; // Long m_nuse; + bool m_multi_ghost = false; // #if ( defined(__CUDACC__) && (__CUDACC_VER_MAJOR__ >= 10) ) CudaGraph m_localCopy; diff --git a/Src/Base/AMReX_FabArrayBase.cpp b/Src/Base/AMReX_FabArrayBase.cpp index 4ede01f1c1b..4dcedc57559 100644 --- a/Src/Base/AMReX_FabArrayBase.cpp +++ b/Src/Base/AMReX_FabArrayBase.cpp @@ -626,11 +626,12 @@ FabArrayBase::getCPC (const IntVect& dstng, const FabArrayBase& src, const IntVe FabArrayBase::FB::FB (const FabArrayBase& fa, const IntVect& nghost, bool cross, const Periodicity& period, - bool enforce_periodicity_only) + bool enforce_periodicity_only, + bool multi_ghost) : m_typ(fa.boxArray().ixType()), m_crse_ratio(fa.boxArray().crseRatio()), m_ngrow(nghost), m_cross(cross), m_epo(enforce_periodicity_only), m_period(period), - m_nuse(0) + m_nuse(0), m_multi_ghost(multi_ghost) { BL_PROFILE("FabArrayBase::FB::FB()"); @@ -651,6 +652,8 @@ FabArrayBase::FB::FB (const FabArrayBase& fa, const IntVect& nghost, void FabArrayBase::FB::define_fb(const FabArrayBase& fa) { + AMREX_ASSERT(m_multi_ghost ? fa.nGrow() >= 2 : true); // must have >= 2 ghost nodes + AMREX_ASSERT(m_multi_ghost ? !m_period.isAnyPeriodic() : true); // this only works for non-periodic const int MyProc = ParallelDescriptor::MyProc(); const BoxArray& ba = fa.boxArray(); const DistributionMapping& dm = fa.DistributionMap(); @@ -661,6 +664,7 @@ FabArrayBase::FB::define_fb(const FabArrayBase& fa) const int nlocal = imap.size(); const IntVect& ng = m_ngrow; + const IntVect ng_ng = m_ngrow - 1; std::vector< std::pair > isects; const std::vector& pshifts = m_period.shiftIntVect(); @@ -671,7 +675,8 @@ FabArrayBase::FB::define_fb(const FabArrayBase& fa) { const int ksnd = imap[i]; const Box& vbx = ba[ksnd]; - + const Box& vbx_ng = amrex::grow(vbx,1); + for (auto pit=pshifts.cbegin(); pit!=pshifts.cend(); ++pit) { ba.intersections(vbx+(*pit), isects, false, ng); @@ -685,7 +690,20 @@ FabArrayBase::FB::define_fb(const FabArrayBase& fa) if (ParallelDescriptor::sameTeam(dst_owner)) { continue; // local copy will be dealt with later } else if (MyProc == dm[ksnd]) { - const BoxList& bl = amrex::boxDiff(bx, ba[krcv]); + BoxList bl = amrex::boxDiff(bx, ba[krcv]); + if (m_multi_ghost) + { + // In the case where ngrow>1, augment the send/rcv box list + // with boxes for overlapping ghost nodes. + const Box& ba_krcv = amrex::grow(ba[krcv],1); + const Box& dst_bx_ng = (amrex::grow(ba_krcv,ng_ng) & (vbx_ng + (*pit))); + const BoxList &bltmp = ba.complementIn(dst_bx_ng); + for (auto const& btmp : bltmp) + { + bl.join(amrex::boxDiff(btmp,ba_krcv)); + } + bl.simplify(); + } for (BoxList::const_iterator lit = bl.begin(); lit != bl.end(); ++lit) send_tags[dst_owner].push_back(CopyComTag(*lit, (*lit)-(*pit), krcv, ksnd)); } @@ -718,6 +736,7 @@ FabArrayBase::FB::define_fb(const FabArrayBase& fa) { const int krcv = imap[i]; const Box& vbx = ba[krcv]; + const Box& vbx_ng = amrex::grow(vbx,1); const Box& bxrcv = amrex::grow(vbx, ng); if (check_local) { @@ -740,7 +759,22 @@ FabArrayBase::FB::define_fb(const FabArrayBase& fa) const Box& dst_bx = isects[j].second - *pit; const int src_owner = dm[ksnd]; - const BoxList& bl = amrex::boxDiff(dst_bx, vbx); + BoxList bl = amrex::boxDiff(dst_bx, vbx); + + if (m_multi_ghost) + { + // In the case where ngrow>1, augment the send/rcv box list + // with boxes for overlapping ghost nodes. + Box ba_ksnd = ba[ksnd]; + ba_ksnd.grow(1); + const Box dst_bx_ng = (ba_ksnd & (bxrcv + (*pit))) - (*pit); + const BoxList &bltmp = ba.complementIn(dst_bx_ng); + for (auto const& btmp : bltmp) + { + bl.join(amrex::boxDiff(btmp,vbx_ng)); + } + bl.simplify(); + } for (BoxList::const_iterator lit = bl.begin(); lit != bl.end(); ++lit) { const Box& blbx = *lit; @@ -1057,6 +1091,7 @@ FabArrayBase::getFB (const IntVect& nghost, const Periodicity& period, it->second->m_crse_ratio == boxArray().crseRatio() && it->second->m_ngrow == nghost && it->second->m_cross == cross && + it->second->m_multi_ghost== m_multi_ghost && it->second->m_epo == enforce_periodicity_only && it->second->m_period == period ) { @@ -1067,7 +1102,7 @@ FabArrayBase::getFB (const IntVect& nghost, const Periodicity& period, } // Have to build a new one - FB* new_fb = new FB(*this, nghost, cross, period, enforce_periodicity_only); + FB* new_fb = new FB(*this, nghost, cross, period, enforce_periodicity_only,m_multi_ghost); #ifdef BL_PROFILE m_FBC_stats.bytes += new_fb->bytes(); diff --git a/Tutorials/LinearSolvers/MultiComponent/MCNodalLinOp.H b/Tutorials/LinearSolvers/MultiComponent/MCNodalLinOp.H index 17b0bad4421..ab19a8a1daf 100644 --- a/Tutorials/LinearSolvers/MultiComponent/MCNodalLinOp.H +++ b/Tutorials/LinearSolvers/MultiComponent/MCNodalLinOp.H @@ -97,9 +97,6 @@ public: private: // Usual buildMasks void buildMasks (); - // This function extends fillBoundary to account for the extra - // ghost nodes that are present. - void realFillBoundary(MultiFab &phi, const Geometry &geom) const; bool m_is_bottom_singular = false; bool m_masks_built = false; diff --git a/Tutorials/LinearSolvers/MultiComponent/MCNodalLinOp.cpp b/Tutorials/LinearSolvers/MultiComponent/MCNodalLinOp.cpp index 9899cfaa9b5..017b5dd4428 100644 --- a/Tutorials/LinearSolvers/MultiComponent/MCNodalLinOp.cpp +++ b/Tutorials/LinearSolvers/MultiComponent/MCNodalLinOp.cpp @@ -155,7 +155,8 @@ void MCNodalLinOp::Fsmooth (int amrlev, int mglev, amrex::MultiFab& a_x, const a } } amrex::Geometry geom = m_geom[amrlev][mglev]; - realFillBoundary(a_x,geom); + a_x.setMultiGhost(true); + a_x.FillBoundary(); nodalSync(amrlev, mglev, a_x); } @@ -403,6 +404,14 @@ void MCNodalLinOp::restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab // i,j,k == fine coordinates amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int I, int J, int K) { int i=2*I, j=2*J, k=2*K; +#if AMREX_SPACEDIM == 2 + cdata(I,J,K,n) = + (fdata(i-1,j-1,k,n) + fdata(i-1,j+1,k,n) + fdata(i+1,j-1,k,n) + fdata(i+1,j+1,k,n)) / 16.0 + + + (fdata(i-1,j,k,n) + fdata(i,j-1,k,n) + fdata(i+1,j,k,n) + fdata(i,j+1,k,n)) / 8.0 + + + fdata(i,j,k,n) / 4.0; +#elif AMREX_SPACEDIM == 3 cdata(I,J,K,n) = (fdata(i-1,j-1,k-1,n) + fdata(i-1,j-1,k+1,n) + fdata(i-1,j+1,k-1,n) + fdata(i-1,j+1,k+1,n) + fdata(i+1,j-1,k-1,n) + fdata(i+1,j-1,k+1,n) + fdata(i+1,j+1,k-1,n) + fdata(i+1,j+1,k+1,n)) / 64.0 @@ -415,6 +424,7 @@ void MCNodalLinOp::restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab fdata(i+1,j,k,n) + fdata(i,j+1,k,n) + fdata(i,j,k+1,n)) / 16.0 + fdata(i,j,k,n) / 8.0; +#endif }); } } @@ -424,7 +434,8 @@ void MCNodalLinOp::restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab } amrex::Geometry geom = m_geom[amrlev][cmglev]; - realFillBoundary(crse,geom); + crse.setMultiGhost(true); + crse.FillBoundary(); nodalSync(amrlev, cmglev, crse); } @@ -491,7 +502,8 @@ void MCNodalLinOp::interpolation (int amrlev, int fmglev, MultiFab& fine, const fine[mfi].plus(tmpfab,fine_bx,fine_bx,0,0,fine.nComp()); } amrex::Geometry geom = m_geom[amrlev][fmglev]; - realFillBoundary(fine,geom); + fine.setMultiGhost(true); + fine.FillBoundary(); nodalSync(amrlev, fmglev, fine); } @@ -504,27 +516,12 @@ void MCNodalLinOp::averageDownSolutionRHS (int camrlev, MultiFab& crse_sol, Mult if (isSingular(0)) amrex::Abort("Singular operators not supported!"); } -void MCNodalLinOp::realFillBoundary(MultiFab &phi, const Geometry &geom) const -{ - for (int i = 0; i < 2; i++) - { - MultiFab & mf = phi; - mf.FillBoundary(geom.periodicity()); - //const int ncomp = mf.nComp(); - const int ng1 = 1; - const int ng2 = 2; - MultiFab tmpmf(mf.boxArray(), mf.DistributionMap(), ncomp, ng1); - MultiFab::Copy(tmpmf, mf, 0, 0, ncomp, ng1); - mf.ParallelCopy (tmpmf, 0, 0, ncomp, ng1, ng2, geom.periodicity()); - } -} - void MCNodalLinOp::applyBC (int amrlev, int mglev, MultiFab& phi, BCMode, amrex::MLLinOp::StateMode , bool skip_fillboundary) const { BL_PROFILE("MCNodalLinOp::applyBC()"); const Geometry& geom = m_geom[amrlev][mglev]; - if (!skip_fillboundary) realFillBoundary(phi,geom); + if (!skip_fillboundary) {phi.setMultiGhost(true); phi.FillBoundary();} } void MCNodalLinOp::reflux (int crse_amrlev, @@ -629,9 +626,9 @@ void MCNodalLinOp::reflux (int crse_amrlev, // Sync up ghost nodes amrex::Geometry geom = m_geom[crse_amrlev][mglev]; - realFillBoundary(res,geom); + res.setMultiGhost(true); + res.FillBoundary(); nodalSync(crse_amrlev,mglev, res); - return; } void @@ -642,7 +639,8 @@ MCNodalLinOp::solutionResidual (int amrlev, MultiFab& resid, MultiFab& x, const apply(amrlev, mglev, resid, x, BCMode::Inhomogeneous, StateMode::Solution); MultiFab::Xpay(resid, -1.0, b, 0, 0, ncomp, 2); amrex::Geometry geom = m_geom[amrlev][mglev]; - realFillBoundary(resid,geom); + resid.setMultiGhost(true); + resid.FillBoundary(); } void @@ -653,5 +651,6 @@ MCNodalLinOp::correctionResidual (int amrlev, int mglev, MultiFab& resid, MultiF apply(amrlev, mglev, resid, x, BCMode::Homogeneous, StateMode::Correction); MultiFab::Xpay(resid, -1.0, b, 0, 0, ncomp, resid.nGrow()); amrex::Geometry geom = m_geom[amrlev][mglev]; - realFillBoundary(resid,geom); + resid.setMultiGhost(true); + resid.FillBoundary(); } diff --git a/Tutorials/LinearSolvers/MultiComponent/inputs b/Tutorials/LinearSolvers/MultiComponent/inputs index d93b29e9287..49945c89169 100644 --- a/Tutorials/LinearSolvers/MultiComponent/inputs +++ b/Tutorials/LinearSolvers/MultiComponent/inputs @@ -5,7 +5,6 @@ mesh.max_grid_size = 32 # Maximum grid size (default = very large --> no mlmg.fixed_iter = 1000 # Number of iterations before exiting gracefully mlmg.verbose = 2 # Verbosity of MLMG mlmg.cg_verbose = 0 # Verbosity of bottom solve -mlmg.max_coarsening_level = 100 # Max MG level mlmg.max_iter = 100 # Max number of iterations before error mlmg.max_fmg_iter = 100 # Number of F-cycles mlmg.agglomeration = 1 # Do agglomeration diff --git a/Tutorials/LinearSolvers/MultiComponent/main.cpp b/Tutorials/LinearSolvers/MultiComponent/main.cpp index eff9142faca..7477e2cc270 100644 --- a/Tutorials/LinearSolvers/MultiComponent/main.cpp +++ b/Tutorials/LinearSolvers/MultiComponent/main.cpp @@ -106,8 +106,8 @@ int main (int argc, char* argv[]) dmap.resize(mesh.nlevels); solution.resize(mesh.nlevels); rhs.resize(mesh.nlevels); - RealBox rb({AMREX_D_DECL(0.0,0.0,0.0)}, - {AMREX_D_DECL(1.0,1.0,1.0)}); + RealBox rb({AMREX_D_DECL(-0.5,-0.5,-0.5)}, + {AMREX_D_DECL(0.5,0.5,0.5)}); Geometry::Setup(&rb, 0); Box NDomain(IntVect{AMREX_D_DECL(0,0,0)}, IntVect{AMREX_D_DECL(mesh.nnodes,mesh.nnodes,mesh.nnodes)}, @@ -147,11 +147,18 @@ int main (int argc, char* argv[]) dmap [ilev].define(cgrids[ilev]); solution[ilev].define(ngrids[ilev], dmap[ilev], op.ncomp, nghost); solution[ilev].setVal(0.0); + solution[ilev].setMultiGhost(true); rhs [ilev].define(ngrids[ilev], dmap[ilev], op.ncomp, nghost); rhs [ilev].setVal(0.0); + rhs [ilev].setMultiGhost(true); Box domain(geom[ilev].Domain()); - const Real* DX = geom[ilev].CellSize(); + const Real AMREX_D_DECL( dx = geom[ilev].CellSize()[0], + dy = geom[ilev].CellSize()[1], + dz = geom[ilev].CellSize()[2]); + const Real AMREX_D_DECL( minx = geom[ilev].ProbLo()[0], + miny = geom[ilev].ProbLo()[1], + minz = geom[ilev].ProbLo()[2]); domain.convert(IntVect::TheNodeVector()); domain.grow(-1); // Shrink domain so we don't operate on any boundaries for (MFIter mfi(solution[ilev], TilingIfNotGPU()); mfi.isValid(); ++mfi) @@ -164,12 +171,17 @@ int main (int argc, char* argv[]) for (int n = 0; n < op.ncomp; n++) ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) { - Real x1 = i*DX[0], x2 = j*DX[1], x3 = k*DX[2]; + Real AMREX_D_DECL(x1 = i*dx + minx, + x2 = j*dy + miny, + x3 = k*dz + minz); - if (n==0) RHS(i,j,k,n) = x1*(1.0 - x1) * x2 * (1.0 - x2) * x3 * (1.0 - x3); + if (n==0) RHS(i,j,k,n) = AMREX_D_TERM( (x1-0.5)*(x1+0.5), + * (x2-0.5)*(x2+0.5), + * (x3-0.5)*(x3+0.5)); else RHS(i,j,k,n) = 0.0; }); } + rhs[ilev].FillBoundary(false,true); } // @@ -188,8 +200,8 @@ int main (int argc, char* argv[]) linop.setNComp(op.ncomp); linop.setCoeff(op.coeff); linop.define(geom,cgrids,dmap,info); - linop.setDomainBC({amrex::MLLinOp::BCType::Dirichlet,amrex::MLLinOp::BCType::Dirichlet,amrex::MLLinOp::BCType::Dirichlet}, - {amrex::MLLinOp::BCType::Dirichlet,amrex::MLLinOp::BCType::Dirichlet,amrex::MLLinOp::BCType::Dirichlet}); + linop.setDomainBC({AMREX_D_DECL(amrex::MLLinOp::BCType::Dirichlet,amrex::MLLinOp::BCType::Dirichlet,amrex::MLLinOp::BCType::Dirichlet)}, + {AMREX_D_DECL(amrex::MLLinOp::BCType::Dirichlet,amrex::MLLinOp::BCType::Dirichlet,amrex::MLLinOp::BCType::Dirichlet)}); for (int ilev = 0; ilev < mesh.nlevels; ++ilev) linop.setLevelBC(ilev,&solution[ilev]); //