Fillboundary

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<CopyMemory> m_localCopy;

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<int,Box> > isects;
 
     const std::vector<IntVect>& 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();

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;

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<RunOn::Host>(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();
 }

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 

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]);
 
     //