Skip to content

Commit

Permalink
Solvability of bottom solver: Follow-up on #2783 (#2801)
Browse files Browse the repository at this point in the history 
In #2783, the solvability fix at the bottom level of EB nodal projection
solver was disabled because we don't have sufficient geometry information.
This broke a number of tests.  It appears that we should still do the best
we can (i.e., use the old approach at the bottom level).
  • Loading branch information
@WeiqunZhang
WeiqunZhang authored May 31, 2022
1 parent 1dff173 commit f24582d
Showing 1 changed file with 77 additions and 83 deletions.
160 changes: 77 additions & 83 deletions Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,59 +181,55 @@ MLNodeLaplacian::getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rh
if (m_coarsening_strategy == CoarseningStrategy::RAP) {
#ifdef AMREX_USE_EB
auto factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][0].get());
if (factory && !factory->isAllRegular()) {
if (mglev > 0) {
return 0._rt;
} else {
const MultiFab& vfrac = factory->getVolFrac();
const auto& vfrac_ma = vfrac.const_arrays();
if (mglev == 0 && factory && !factory->isAllRegular()) {
const MultiFab& vfrac = factory->getVolFrac();
const auto& vfrac_ma = vfrac.const_arrays();

Box dom = Geom(amrlev,mglev).Domain();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (m_lobc[0][idim] != LinOpBCType::Neumann &&
m_lobc[0][idim] != LinOpBCType::inflow)
{
dom.growLo(idim, 10);
}
if (m_hibc[0][idim] != LinOpBCType::Neumann &&
m_hibc[0][idim] != LinOpBCType::inflow)
{
dom.growHi(idim, 10);
}
Box dom = Geom(amrlev,mglev).Domain();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (m_lobc[0][idim] != LinOpBCType::Neumann &&
m_lobc[0][idim] != LinOpBCType::inflow)
{
dom.growLo(idim, 10);
}
if (m_hibc[0][idim] != LinOpBCType::Neumann &&
m_hibc[0][idim] != LinOpBCType::inflow)
{
dom.growHi(idim, 10);
}
}

const auto& mask = (mglev+1 == m_num_mg_levels[0]) ? m_bottom_dot_mask : m_coarse_dot_mask;
const auto& mask_ma = mask.const_arrays();
const auto& rhs_ma = rhs.const_arrays();
auto r = ParReduce(TypeList<ReduceOpSum,ReduceOpSum>{}, TypeList<Real,Real>{},
rhs, IntVect(0),
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
-> GpuTuple<Real,Real>
{
Real scale = 0.0_rt;
const auto& mask = (mglev+1 == m_num_mg_levels[0]) ? m_bottom_dot_mask : m_coarse_dot_mask;
const auto& mask_ma = mask.const_arrays();
const auto& rhs_ma = rhs.const_arrays();
auto r = ParReduce(TypeList<ReduceOpSum,ReduceOpSum>{}, TypeList<Real,Real>{},
rhs, IntVect(0),
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
-> GpuTuple<Real,Real>
{
Real scale = 0.0_rt;
#if (AMREX_SPACEDIM == 3)
int const koff = 1;
Real const fac = 0.125_rt;
int const koff = 1;
Real const fac = 0.125_rt;
#else
int const koff = 0;
Real const fac = 0.25_rt;
int const koff = 0;
Real const fac = 0.25_rt;
#endif
for (int kc = k-koff; kc <= k; ++kc) {
for (int jc = j-1 ; jc <= j; ++jc) {
for (int ic = i-1 ; ic <= i; ++ic) {
if (dom.contains(ic,jc,kc)) {
scale += vfrac_ma[box_no](ic,jc,kc) * fac;
}
}}}
return { mask_ma[box_no](i,j,k) * rhs_ma[box_no](i,j,k),
mask_ma[box_no](i,j,k) * scale };
});

Real s1 = amrex::get<0>(r);
Real s2 = amrex::get<1>(r);
ParallelAllReduce::Sum<Real>({s1,s2}, ParallelContext::CommunicatorSub());
return s1/s2;
}
for (int kc = k-koff; kc <= k; ++kc) {
for (int jc = j-1 ; jc <= j; ++jc) {
for (int ic = i-1 ; ic <= i; ++ic) {
if (dom.contains(ic,jc,kc)) {
scale += vfrac_ma[box_no](ic,jc,kc) * fac;
}
}}}
return { mask_ma[box_no](i,j,k) * rhs_ma[box_no](i,j,k),
mask_ma[box_no](i,j,k) * scale };
});

Real s1 = amrex::get<0>(r);
Real s2 = amrex::get<1>(r);
ParallelAllReduce::Sum<Real>({s1,s2}, ParallelContext::CommunicatorSub());
return s1/s2;
} else
#endif
{
Expand Down Expand Up @@ -296,47 +292,45 @@ MLNodeLaplacian::fixSolvabilityByOffset (int amrlev, int mglev, MultiFab& rhs, R
if (m_coarsening_strategy == CoarseningStrategy::RAP) {
#ifdef AMREX_USE_EB
auto factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][0].get());
if (factory && !factory->isAllRegular()) {
if (mglev == 0) {
const MultiFab& vfrac = factory->getVolFrac();
const auto& vfrac_ma = vfrac.const_arrays();

Box dom = Geom(amrlev,mglev).Domain();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (m_lobc[0][idim] != LinOpBCType::Neumann &&
m_lobc[0][idim] != LinOpBCType::inflow)
{
dom.growLo(idim, 10);
}
if (m_hibc[0][idim] != LinOpBCType::Neumann &&
m_hibc[0][idim] != LinOpBCType::inflow)
{
dom.growHi(idim, 10);
}
if (mglev == 0 && factory && !factory->isAllRegular()) {
const MultiFab& vfrac = factory->getVolFrac();
const auto& vfrac_ma = vfrac.const_arrays();

Box dom = Geom(amrlev,mglev).Domain();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (m_lobc[0][idim] != LinOpBCType::Neumann &&
m_lobc[0][idim] != LinOpBCType::inflow)
{
dom.growLo(idim, 10);
}
if (m_hibc[0][idim] != LinOpBCType::Neumann &&
m_hibc[0][idim] != LinOpBCType::inflow)
{
dom.growHi(idim, 10);
}
}

auto const& rhs_ma = rhs.arrays();
ParallelFor(rhs, IntVect(0),
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
{
Real scale = 0.0_rt;
auto const& rhs_ma = rhs.arrays();
ParallelFor(rhs, IntVect(0),
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
{
Real scale = 0.0_rt;
#if (AMREX_SPACEDIM == 3)
int const koff = 1;
Real const fac = 0.125_rt;
int const koff = 1;
Real const fac = 0.125_rt;
#else
int const koff = 0;
Real const fac = 0.25_rt;
int const koff = 0;
Real const fac = 0.25_rt;
#endif
for (int kc = k-koff; kc <= k; ++kc) {
for (int jc = j-1 ; jc <= j; ++jc) {
for (int ic = i-1 ; ic <= i; ++ic) {
if (dom.contains(ic,jc,kc)) {
scale += vfrac_ma[box_no](ic,jc,kc) * fac;
}
}}}
rhs_ma[box_no](i,j,k) -= offset * scale;
});
}
for (int kc = k-koff; kc <= k; ++kc) {
for (int jc = j-1 ; jc <= j; ++jc) {
for (int ic = i-1 ; ic <= i; ++ic) {
if (dom.contains(ic,jc,kc)) {
scale += vfrac_ma[box_no](ic,jc,kc) * fac;
}
}}}
rhs_ma[box_no](i,j,k) -= offset * scale;
});
} else
#endif
{
Expand Down

0 comments on commit f24582d

Please sign in to comment.