Skip to content

Commit

Permalink
Neumann BC at coarse/fine interface for cell-centered solvers (#3926)
Browse files Browse the repository at this point in the history 
For cell-centered linear solvers, if the lowest AMR level in the solver
is not AMR level 0, the coarse/fine interface was previously assumed to
be Dirichlet. In this commit, this has been relaxed to allow for
homogeneous Neumann BC at the coarse/fine interface.
  • Loading branch information
@WeiqunZhang
WeiqunZhang committed May 9, 2024
1 parent b082d3e commit a89b465
Show file tree
Hide file tree
Showing 6 changed files with 175 additions and 53 deletions.
38 changes: 35 additions & 3 deletions Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.H
View file
Expand Up @@ -726,7 +726,7 @@ MLABecLaplacianT<MF>::update_singular_flags ()
// For now this assumes that overset regions are treated as Dirichlet bc's
if (this->m_domain_covered[alev] && !this->m_overset_mask[alev][0])
{
if (m_a_scalar == 0.0)
if (m_a_scalar == Real(0.0))
{
m_is_singular[alev] = true;
}
Expand All @@ -739,6 +739,37 @@ MLABecLaplacianT<MF>::update_singular_flags ()
}
}
}

if (!m_is_singular[0] && this->m_needs_coarse_data_for_bc &&
this->m_coarse_fine_bc_type == LinOpBCType::Neumann)
{
AMREX_ASSERT(this->m_overset_mask[0][0] == nullptr);

bool lev0_a_is_zero = false;
if (m_a_scalar == Real(0.0)) {
lev0_a_is_zero = true;
} else {
RT asum = m_a_coeffs[0].back().sum(0,IntVect(0));
RT amax = m_a_coeffs[0].back().norminf(0,1,IntVect(0));
bool a_is_almost_zero = std::abs(asum) <= amax * RT(1.e-12);
if (a_is_almost_zero) { lev0_a_is_zero = true; }
}

if (lev0_a_is_zero) {
auto bbox = this->m_grids[0][0].minimalBox();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (this->m_lobc[0][idim] == LinOpBCType::Dirichlet) {
bbox.growLo(idim,1);
}
if (this->m_hibc[0][idim] == LinOpBCType::Dirichlet) {
bbox.growHi(idim,1);
}
}
if (this->m_geom[0][0].Domain().contains(bbox)) {
m_is_singular[0] = true;
}
}
}
}

template <typename MF>
Expand Down Expand Up @@ -1217,8 +1248,9 @@ MLABecLaplacianT<MF>::makeNLinOp (int /*grid_size*/) const
if (this->needsCoarseDataForBC())
{
const Real* dx0 = this->m_geom[0][0].CellSize();
const Real fac = Real(0.5)*this->m_coarse_data_crse_ratio;
RealVect cbloc {AMREX_D_DECL(dx0[0]*fac, dx0[1]*fac, dx0[2]*fac)};
RealVect fac(this->m_coarse_data_crse_ratio);
fac *= Real(0.5);
RealVect cbloc {AMREX_D_DECL(dx0[0]*fac[0], dx0[1]*fac[1], dx0[2]*fac[2])};
nop->setCoarseFineBCLocation(cbloc);
}

Expand Down
47 changes: 30 additions & 17 deletions Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H
View file
Expand Up @@ -156,9 +156,10 @@ protected:
void setLOBndryConds (const Geometry& geom, const Real* dx,
const Vector<Array<BCType,AMREX_SPACEDIM> >& lobc,
const Vector<Array<BCType,AMREX_SPACEDIM> >& hibc,
int ratio, const RealVect& interior_bloc,
IntVect const& ratio, const RealVect& interior_bloc,
const Array<Real,AMREX_SPACEDIM>& domain_bloc_lo,
const Array<Real,AMREX_SPACEDIM>& domain_bloc_hi);
const Array<Real,AMREX_SPACEDIM>& domain_bloc_hi,
LinOpBCType crse_fine_bc_type);

const Vector<BCTuple>& bndryConds (const MFIter& mfi) const noexcept {
return bcond[mfi];
Expand Down Expand Up @@ -289,9 +290,10 @@ MLCellLinOpT<MF>::BndryCondLoc::
setLOBndryConds (const Geometry& geom, const Real* dx,
const Vector<Array<BCType,AMREX_SPACEDIM> >& lobc,
const Vector<Array<BCType,AMREX_SPACEDIM> >& hibc,
int ratio, const RealVect& interior_bloc,
IntVect const& ratio, const RealVect& interior_bloc,
const Array<Real,AMREX_SPACEDIM>& domain_bloc_lo,
const Array<Real,AMREX_SPACEDIM>& domain_bloc_hi)
const Array<Real,AMREX_SPACEDIM>& domain_bloc_hi,
LinOpBCType crse_fine_bc_type)
{
const Box& domain = geom.Domain();

Expand All @@ -308,7 +310,8 @@ setLOBndryConds (const Geometry& geom, const Real* dx,
lobc[icomp], hibc[icomp],
dx, ratio, interior_bloc,
domain_bloc_lo, domain_bloc_hi,
geom.isPeriodicArray());
geom.isPeriodicArray(),
crse_fine_bc_type);
}
}

Expand Down Expand Up @@ -513,58 +516,61 @@ MLCellLinOpT<MF>::setLevelBC (int amrlev, const MF* a_levelbcdata, const MF* rob
}
const MF& bcdata = (a_levelbcdata == nullptr) ? zero : *a_levelbcdata;

int br_ref_ratio = -1;
IntVect br_ref_ratio(-1);

if (amrlev == 0)
{
if (this->needsCoarseDataForBC())
{
AMREX_ALWAYS_ASSERT(!this->hasHiddenDimension());
br_ref_ratio = this->m_coarse_data_crse_ratio > 0 ? this->m_coarse_data_crse_ratio : 2;
if (m_crse_sol_br[amrlev] == nullptr && br_ref_ratio > 0)
br_ref_ratio = this->m_coarse_data_crse_ratio.allGT(0) ? this->m_coarse_data_crse_ratio : IntVect(2);
if (m_crse_sol_br[amrlev] == nullptr && br_ref_ratio.allGT(0))
{
const int in_rad = 0;
const int out_rad = 1;
const int extent_rad = 2;
const int crse_ratio = br_ref_ratio;
const IntVect crse_ratio = br_ref_ratio;
BoxArray cba = this->m_grids[amrlev][0];
cba.coarsen(crse_ratio);
m_crse_sol_br[amrlev] = std::make_unique<BndryRegisterT<MF>>
(cba, this->m_dmap[amrlev][0], in_rad, out_rad, extent_rad, ncomp);
}
if (this->m_coarse_data_for_bc != nullptr) {
AMREX_ALWAYS_ASSERT(this->m_coarse_data_crse_ratio > 0);
AMREX_ALWAYS_ASSERT(this->m_coarse_data_crse_ratio.allGT(0));
const Box& cbx = amrex::coarsen(this->m_geom[0][0].Domain(), this->m_coarse_data_crse_ratio);
m_crse_sol_br[amrlev]->copyFrom(*this->m_coarse_data_for_bc, 0, 0, 0, ncomp,
this->m_geom[0][0].periodicity(cbx));
} else {
m_crse_sol_br[amrlev]->setVal(RT(0.0));
}
m_bndry_sol[amrlev]->setBndryValues(*m_crse_sol_br[amrlev], 0,
bcdata, 0, 0, ncomp, IntVect(br_ref_ratio));
bcdata, 0, 0, ncomp, br_ref_ratio);
br_ref_ratio = this->m_coarse_data_crse_ratio;
}
else
{
m_bndry_sol[amrlev]->setPhysBndryValues(bcdata,0,0,ncomp);
br_ref_ratio = 1;
br_ref_ratio = IntVect(1);
}
}
else
{
m_bndry_sol[amrlev]->setPhysBndryValues(bcdata,0,0,ncomp);
br_ref_ratio = this->m_amr_ref_ratio[amrlev-1];
br_ref_ratio = IntVect(this->m_amr_ref_ratio[amrlev-1]);
}

m_bndry_sol[amrlev]->setLOBndryConds(this->m_lobc, this->m_hibc, br_ref_ratio, this->m_coarse_bc_loc);
auto crse_fine_bc_type = (amrlev == 0) ? this->m_coarse_fine_bc_type : LinOpBCType::Dirichlet;
m_bndry_sol[amrlev]->setLOBndryConds(this->m_lobc, this->m_hibc, br_ref_ratio,
this->m_coarse_bc_loc, crse_fine_bc_type);

const Real* dx = this->m_geom[amrlev][0].CellSize();
for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
{
m_bcondloc[amrlev][mglev]->setLOBndryConds(this->m_geom[amrlev][mglev], dx,
this->m_lobc, this->m_hibc,
br_ref_ratio, this->m_coarse_bc_loc,
this->m_domain_bloc_lo, this->m_domain_bloc_hi);
this->m_domain_bloc_lo, this->m_domain_bloc_hi,
crse_fine_bc_type);
}

if (this->hasRobinBC()) {
Expand Down Expand Up @@ -1911,6 +1917,8 @@ MLCellLinOpT<MF>::computeVolInv () const
m_volinv[amrlev].resize(this->NMGLevels(amrlev));
}

AMREX_ASSERT(this->m_coarse_fine_bc_type == LinOpBCType::Dirichlet || ! this->hasHiddenDimension());

// We don't need to compute for every level

auto f = [&] (int amrlev, int mglev) {
Expand All @@ -1928,8 +1936,13 @@ MLCellLinOpT<MF>::computeVolInv () const
else
#endif
{
m_volinv[amrlev][mglev]
= RT(1.0 / this->compactify(this->Geom(amrlev,mglev).Domain()).d_numPts());
if (this->m_coarse_fine_bc_type == LinOpBCType::Dirichlet) {
m_volinv[amrlev][mglev]
= RT(1.0 / this->compactify(this->Geom(amrlev,mglev).Domain()).d_numPts());
} else {
m_volinv[amrlev][mglev]
= RT(1.0 / this->m_grids[amrlev][mglev].d_numPts());
}
}
};

Expand Down
50 changes: 42 additions & 8 deletions Src/LinearSolvers/MLMG/AMReX_MLLinOp.H
View file
Expand Up @@ -184,20 +184,32 @@ public:
* coarsest AMR level of the solve come from a coarser level (e.g. the
* base AMR level of the solve is > 0 and does not cover the entire
* domain), we must explicitly provide the coarser data. Boundary
* conditions from a coarser level are always Dirichlet. The MultiFab
* conditions from a coarser level are Dirichlet by default. The MultiFab
* crse does not need to have ghost cells and is at a coarser resolution
* than the coarsest AMR level of the solve; it is used to supply
* (interpolated) boundary conditions for the solve. NOTE: If this is
* called, it must be called before `setLevelBC`. If crse is nullptr,
* then the bc values are assumed to be zero.
* then the bc values are assumed to be zero. The coarse/fine BC type
* can be changed to homogeneous Neumann by the bc_type argument. In that
* case, use nullptr for the crse argument.
*
* \param crse the coarse AMR level data
* \param crse_ratio the coarsening ratio between fine and coarse AMR levels.
* \param bc_type optional. It's Dirichlet by default, and can be Neumann.
*/
void setCoarseFineBC (const MF* crse, int crse_ratio) noexcept;
void setCoarseFineBC (const MF* crse, int crse_ratio,
LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;

void setCoarseFineBC (const MF* crse, IntVect const& crse_ratio,
LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;

template <typename AMF, std::enable_if_t<!std::is_same_v<MF,AMF>,int> = 0>
void setCoarseFineBC (const AMF* crse, int crse_ratio,
LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;

template <typename AMF, std::enable_if_t<!std::is_same_v<MF,AMF>,int> = 0>
void setCoarseFineBC (const AMF* crse, int crse_ratio) noexcept;
void setCoarseFineBC (const AMF* crse, IntVect const& crse_ratio,
LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;

/**
* \brief Set boundary conditions for given level. For cell-centered
Expand Down Expand Up @@ -612,7 +624,8 @@ protected:
Array<Real, AMREX_SPACEDIM> m_domain_bloc_hi {{AMREX_D_DECL(0._rt,0._rt,0._rt)}};

bool m_needs_coarse_data_for_bc = false;
int m_coarse_data_crse_ratio = -1;
LinOpBCType m_coarse_fine_bc_type = LinOpBCType::Dirichlet;
IntVect m_coarse_data_crse_ratio = IntVect(-1);
RealVect m_coarse_bc_loc;
const MF* m_coarse_data_for_bc = nullptr;
MF m_coarse_data_for_bc_raii;
Expand Down Expand Up @@ -889,7 +902,7 @@ MLLinOpT<MF>::defineGrids (const Vector<Geometry>& a_geom,
Box aggbox;
bool aggable = false;

if (info.do_agglomeration)
if (m_grids[0][0].size() > 1 && info.do_agglomeration)
{
if (m_domain_covered[0])
{
Expand Down Expand Up @@ -1404,23 +1417,44 @@ MLLinOpT<MF>::setDomainBCLoc (const Array<Real,AMREX_SPACEDIM>& lo_bcloc,

template <typename MF>
void
MLLinOpT<MF>::setCoarseFineBC (const MF* crse, int crse_ratio) noexcept
MLLinOpT<MF>::setCoarseFineBC (const MF* crse, int crse_ratio,
LinOpBCType bc_type) noexcept
{
setCoarseFineBC(crse, IntVect(crse_ratio), bc_type);
}

template <typename MF>
void
MLLinOpT<MF>::setCoarseFineBC (const MF* crse, IntVect const& crse_ratio,
LinOpBCType bc_type) noexcept
{
m_coarse_data_for_bc = crse;
m_coarse_data_crse_ratio = crse_ratio;
m_coarse_fine_bc_type = bc_type;
}

template <typename MF>
template <typename AMF, std::enable_if_t<!std::is_same_v<MF,AMF>,int>>
void
MLLinOpT<MF>::setCoarseFineBC (const AMF* crse, int crse_ratio,
LinOpBCType bc_type) noexcept
{
setCoarseFineBC(crse, IntVect(crse_ratio), bc_type);
}

template <typename MF>
template <typename AMF, std::enable_if_t<!std::is_same_v<MF,AMF>,int>>
void
MLLinOpT<MF>::setCoarseFineBC (const AMF* crse, int crse_ratio) noexcept
MLLinOpT<MF>::setCoarseFineBC (const AMF* crse, IntVect const& crse_ratio,
LinOpBCType bc_type) noexcept
{
m_coarse_data_for_bc_raii = MF(crse->boxArray(), crse->DistributionMap(),
crse->nComp(), crse->nGrowVect());
m_coarse_data_for_bc_raii.LocalCopy(*crse, 0, 0, crse->nComp(),
crse->nGrowVect());
m_coarse_data_for_bc = &m_coarse_data_for_bc_raii;
m_coarse_data_crse_ratio = crse_ratio;
m_coarse_fine_bc_type = bc_type;
}

template <typename MF>
Expand Down
22 changes: 12 additions & 10 deletions Src/LinearSolvers/MLMG/AMReX_MLMG.H
View file
Expand Up @@ -1893,11 +1893,12 @@ MLMGT<MF>::bottomSolveWithHypre (MF& x, const MF& b)
hypre_bndry = std::make_unique<MLMGBndryT<MF>>(ba, dm, ncomp, geom);
hypre_bndry->setHomogValues();
const Real* dx = linop.m_geom[0][0].CellSize();
int crse_ratio = linop.m_coarse_data_crse_ratio > 0 ? linop.m_coarse_data_crse_ratio : 1;
RealVect bclocation(AMREX_D_DECL(0.5*dx[0]*crse_ratio,
0.5*dx[1]*crse_ratio,
0.5*dx[2]*crse_ratio));
hypre_bndry->setLOBndryConds(linop.m_lobc, linop.m_hibc, -1, bclocation);
IntVect crse_ratio = linop.m_coarse_data_crse_ratio.allGT(0) ? linop.m_coarse_data_crse_ratio : IntVect(1);
RealVect bclocation(AMREX_D_DECL(0.5*dx[0]*crse_ratio[0],
0.5*dx[1]*crse_ratio[1],
0.5*dx[2]*crse_ratio[2]));
hypre_bndry->setLOBndryConds(linop.m_lobc, linop.m_hibc, IntVect(-1), bclocation,
linop.m_coarse_fine_bc_type);
}

// IJ interface understands absolute tolerance API of hypre
Expand Down Expand Up @@ -1947,11 +1948,12 @@ MLMGT<MF>::bottomSolveWithPETSc (MF& x, const MF& b)
petsc_bndry = std::make_unique<MLMGBndryT<MF>>(ba, dm, ncomp, geom);
petsc_bndry->setHomogValues();
const Real* dx = linop.m_geom[0][0].CellSize();
int crse_ratio = linop.m_coarse_data_crse_ratio > 0 ? linop.m_coarse_data_crse_ratio : 1;
RealVect bclocation(AMREX_D_DECL(0.5*dx[0]*crse_ratio,
0.5*dx[1]*crse_ratio,
0.5*dx[2]*crse_ratio));
petsc_bndry->setLOBndryConds(linop.m_lobc, linop.m_hibc, -1, bclocation);
auto crse_ratio = linop.m_coarse_data_crse_ratio.allGT(0) ? linop.m_coarse_data_crse_ratio : IntVect(1);
RealVect bclocation(AMREX_D_DECL(0.5*dx[0]*crse_ratio[0],
0.5*dx[1]*crse_ratio[1],
0.5*dx[2]*crse_ratio[2]));
petsc_bndry->setLOBndryConds(linop.m_lobc, linop.m_hibc, IntVect(-1), bclocation,
linop.m_coarse_fine_bc_type);
}
petsc_solver->solve(x, b, bottom_reltol, Real(-1.), bottom_maxiter, *petsc_bndry,
linop.getMaxOrder());
Expand Down

0 comments on commit a89b465

Please sign in to comment.