Skip to content

Commit

Permalink
Fillboundary (AMReX-Codes#1171)
Browse files Browse the repository at this point in the history 
This pull request fixes the problem of edge nodes not getting filled properly when filling multiple layers of ghost nodes. This occurs when there are overlapping ghost nodes on a coarse/fine boundary. (See issue AMReX-Codes#568)

Summary of changes to AMreX:
- Added a `BoxArray::complementIn` function that takes a box _list_ and returns the complement box list. 
- Added an additional boolean argument `multi_ghost` to all `FillBoundary` functions that gets passed (eventually) to `define_fb`. In all cases, a default value of `false` is supplied.
- Modified `define_fb` to include the missing ghost nodes in the list of `send_tags`. (This uses the new `complementIn` function.)

Summery of changes to `LinearSolver/MultiComponent` tutorial:
- Removed the `realFillBoundary` function, replaced with `FillBoundary(false,true)
- Modified `main.cpp` so that it can be compiled and run (non-trivially) in 2D as well as 3D.

Diagram of the problem that this PR fixes:
![FillBoundaryIllustration](https://user-images.githubusercontent.com/6371567/88218141-d9982980-cc1c-11ea-84f7-194934dfed73.png)
  • Loading branch information
@bsrunnels @dwillcox
bsrunnels authored and dwillcox committed Oct 3, 2020
1 parent 88fb9cf commit 456a598
Show file tree
Hide file tree
Showing 6 changed files with 86 additions and 41 deletions.
5 changes: 4 additions & 1 deletion Src/Base/AMReX_FabArrayBase.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand Down
47 changes: 41 additions & 6 deletions Src/Base/AMReX_FabArrayBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()");

Expand All @@ -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();
Expand All @@ -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();
Expand All @@ -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);
Expand All @@ -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));
}
Expand Down Expand Up @@ -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) {
Expand All @@ -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;
Expand Down Expand Up @@ -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 )
{
Expand All @@ -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();
Expand Down
3 changes: 0 additions & 3 deletions Tutorials/LinearSolvers/MultiComponent/MCNodalLinOp.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
45 changes: 22 additions & 23 deletions Tutorials/LinearSolvers/MultiComponent/MCNodalLinOp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

Expand Down Expand Up @@ -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
Expand All @@ -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
});
}
}
Expand All @@ -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);
}

Expand Down Expand Up @@ -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);
}

Expand All @@ -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,
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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();
}
1 change: 0 additions & 1 deletion Tutorials/LinearSolvers/MultiComponent/inputs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
26 changes: 19 additions & 7 deletions Tutorials/LinearSolvers/MultiComponent/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)},
Expand Down Expand Up @@ -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)
Expand All @@ -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);
}

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

//
Expand Down

0 comments on commit 456a598

Please sign in to comment.