Skip to content

Commit b78921a

Browse files
WeiqunZhangWeiqunZhang
authored
Fix solvability issue in the nodal solver RAP approach (#2783)
In the RAP approach of the nodal solver, the RHS at Neumann boundaries only includes the integral inside the domain. That is it's only half of what its "physical" value is. The usual way of subtracting a constant from the RHS does not work for RAP. We should only subtract half of the constant offset at Neumann boundaries. The computation of the constant offset needed for the solvability fix is also affected by the way how RHS is computed at Neumann boundaries. For EB, the computation of the constant offset is an integral of the RHS multiplied by the volume fraction, and the subtraction also needs to be weighted by the volume fraction.
1 parent 806108c commit b78921a

File tree

7 files changed

+255
-18
lines changed

7 files changed

+255
-18
lines changed

‎Src/LinearSolvers/MLMG/AMReX_MLLinOp.H‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -231,6 +231,8 @@ public:
231231
virtual void applyInhomogNeumannTerm (int /*amrlev*/, MultiFab& /*rhs*/) const {}
232232
virtual void applyOverset (int /*amlev*/, MultiFab& /*rhs*/) const {}
233233
virtual void scaleRHS (int /*amrlev*/, MultiFab& /*rhs*/) const {}
234+
virtual Real getSolvabilityOffset (int /*amrlev*/, int /*mglev*/, MultiFab const& /*rhs*/) const { return 0._rt; } // Only nodal solvers need it
235+
virtual void fixSolvabilityByOffset (int /*amrlev*/, int /*mglev*/, MultiFab& /*rhs*/, Real /*offset*/) const {} // Only nodal solvers need it
234236

235237
virtual void prepareForSolve () = 0;
236238
virtual bool isSingular (int amrlev) const = 0;

‎Src/LinearSolvers/MLMG/AMReX_MLMG.H‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -156,7 +156,6 @@ public:
156156
void computeVolInv ();
157157
void makeSolvable ();
158158
void makeSolvable (int amrlev, int mglev, MultiFab& mf);
159-
Real getNodalSum (int amrlev, int mglev, MultiFab& mf) const;
160159

161160
#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
162161
void bottomSolveWithHypre (MultiFab& x, const MultiFab& b);

‎Src/LinearSolvers/MLMG/AMReX_MLMG.cpp‎

Lines changed: 4 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -1776,12 +1776,12 @@ MLMG::makeSolvable ()
17761776
else
17771777
{
17781778
AMREX_ASSERT_WITH_MESSAGE(ncomp==1, "ncomp > 1 not supported for singular nodal problem");
1779-
Real offset = getNodalSum(0, 0, rhs[0]);
1779+
Real offset = linop.getSolvabilityOffset(0, 0, rhs[0]);
17801780
if (verbose >= 4) {
17811781
amrex::Print() << "MLMG: Subtracting " << offset << " from rhs\n";
17821782
}
17831783
for (int alev = 0; alev < namrlevs; ++alev) {
1784-
rhs[alev].plus(-offset, 0, 1);
1784+
linop.fixSolvabilityByOffset(alev, 0, rhs[alev], offset);
17851785
}
17861786
}
17871787
}
@@ -1833,27 +1833,15 @@ MLMG::makeSolvable (int amrlev, int mglev, MultiFab& mf)
18331833
else
18341834
{
18351835
AMREX_ASSERT_WITH_MESSAGE(ncomp==1, "ncomp > 1 not supported for singular nodal problem");
1836-
Real offset = getNodalSum(amrlev, mglev, mf);
1836+
Real offset = linop.getSolvabilityOffset(amrlev, mglev, mf);
18371837
if (verbose >= 4) {
18381838
amrex::Print() << "MLMG: Subtracting " << offset << " on level (" << amrlev << ", "
18391839
<< mglev << ")\n";
18401840
}
1841-
mf.plus(-offset, 0, 1);
1841+
linop.fixSolvabilityByOffset(amrlev, mglev, mf, offset);
18421842
}
18431843
}
18441844

1845-
Real
1846-
MLMG::getNodalSum (int amrlev, int mglev, MultiFab& mf) const
1847-
{
1848-
MultiFab one(mf.boxArray(), mf.DistributionMap(), 1, 0, MFInfo(), mf.Factory());
1849-
one.setVal(Real(1.0));
1850-
const bool local = true;
1851-
Real s1 = linop.xdoty(amrlev, mglev, mf, one, local);
1852-
Real s2 = linop.xdoty(amrlev, mglev, one, one, local);
1853-
ParallelAllReduce::Sum<Real>({s1,s2}, ParallelContext::CommunicatorSub());
1854-
return s1/s2;
1855-
}
1856-
18571845
#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
18581846
void
18591847
MLMG::bottomSolveWithHypre (MultiFab& x, const MultiFab& b)

‎Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.H‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -117,6 +117,8 @@ public :
117117
virtual void getFluxes (const Vector<MultiFab*>& a_flux,
118118
const Vector<MultiFab*>& a_sol) const final override;
119119
virtual void unimposeNeumannBC (int amrlev, MultiFab& rhs) const final override;
120+
virtual Real getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rhs) const override;
121+
virtual void fixSolvabilityByOffset (int amrlev, int mglev, MultiFab& rhs, Real offset) const override;
120122

121123
virtual void compGrad (int /*amrlev*/, const Array<MultiFab*,AMREX_SPACEDIM>& /*grad*/,
122124
MultiFab& /*sol*/, Location /*loc*/) const final override {

‎Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp‎

Lines changed: 213 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -171,6 +171,219 @@ MLNodeLaplacian::unimposeNeumannBC (int amrlev, MultiFab& rhs) const
171171
}
172172
}
173173

174+
Real
175+
MLNodeLaplacian::getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rhs) const
176+
{
177+
amrex::ignore_unused(amrlev);
178+
AMREX_ASSERT(amrlev==0);
179+
AMREX_ASSERT(mglev+1==m_num_mg_levels[0] || mglev==0);
180+
181+
if (m_coarsening_strategy == CoarseningStrategy::RAP) {
182+
#ifdef AMREX_USE_EB
183+
auto factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][0].get());
184+
if (factory && !factory->isAllRegular()) {
185+
if (mglev > 0) {
186+
return 0._rt;
187+
} else {
188+
const MultiFab& vfrac = factory->getVolFrac();
189+
const auto& vfrac_ma = vfrac.const_arrays();
190+
191+
Box dom = Geom(amrlev,mglev).Domain();
192+
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
193+
if (m_lobc[0][idim] != LinOpBCType::Neumann &&
194+
m_lobc[0][idim] != LinOpBCType::inflow)
195+
{
196+
dom.growLo(idim, 10);
197+
}
198+
if (m_hibc[0][idim] != LinOpBCType::Neumann &&
199+
m_hibc[0][idim] != LinOpBCType::inflow)
200+
{
201+
dom.growHi(idim, 10);
202+
}
203+
}
204+
205+
const auto& mask = (mglev+1 == m_num_mg_levels[0]) ? m_bottom_dot_mask : m_coarse_dot_mask;
206+
const auto& mask_ma = mask.const_arrays();
207+
const auto& rhs_ma = rhs.const_arrays();
208+
auto r = ParReduce(TypeList<ReduceOpSum,ReduceOpSum>{}, TypeList<Real,Real>{},
209+
rhs, IntVect(0),
210+
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
211+
-> GpuTuple<Real,Real>
212+
{
213+
Real scale = 0.0_rt;
214+
#if (AMREX_SPACEDIM == 3)
215+
int const koff = 1;
216+
Real const fac = 0.125_rt;
217+
#else
218+
int const koff = 0;
219+
Real const fac = 0.25_rt;
220+
#endif
221+
for (int kc = k-koff; kc <= k; ++kc) {
222+
for (int jc = j-1 ; jc <= j; ++jc) {
223+
for (int ic = i-1 ; ic <= i; ++ic) {
224+
if (dom.contains(ic,jc,kc)) {
225+
scale += vfrac_ma[box_no](ic,jc,kc) * fac;
226+
}
227+
}}}
228+
return { mask_ma[box_no](i,j,k) * rhs_ma[box_no](i,j,k),
229+
mask_ma[box_no](i,j,k) * scale };
230+
});
231+
232+
Real s1 = amrex::get<0>(r);
233+
Real s2 = amrex::get<1>(r);
234+
ParallelAllReduce::Sum<Real>({s1,s2}, ParallelContext::CommunicatorSub());
235+
return s1/s2;
236+
}
237+
} else
238+
#endif
239+
{
240+
Box nddom = amrex::surroundingNodes(Geom(amrlev,mglev).Domain());
241+
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
242+
if (m_lobc[0][idim] != LinOpBCType::Neumann &&
243+
m_lobc[0][idim] != LinOpBCType::inflow)
244+
{
245+
nddom.growLo(idim, 10); // so that the test in ParReduce will faill
246+
}
247+
if (m_hibc[0][idim] != LinOpBCType::Neumann &&
248+
m_hibc[0][idim] != LinOpBCType::inflow)
249+
{
250+
nddom.growHi(idim, 10);
251+
}
252+
}
253+
254+
const auto& mask = (mglev+1 == m_num_mg_levels[0]) ? m_bottom_dot_mask : m_coarse_dot_mask;
255+
const auto& mask_ma = mask.const_arrays();
256+
const auto& rhs_ma = rhs.const_arrays();
257+
auto r = ParReduce(TypeList<ReduceOpSum,ReduceOpSum>{}, TypeList<Real,Real>{},
258+
rhs, IntVect(0),
259+
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
260+
-> GpuTuple<Real,Real>
261+
{
262+
Real scale = 1.0_rt;
263+
if (i == nddom.smallEnd(0) ||
264+
i == nddom.bigEnd(0)) {
265+
scale *= 0.5_rt;
266+
}
267+
#if (AMREX_SPACEDIM >= 2)
268+
if (j == nddom.smallEnd(1) ||
269+
j == nddom.bigEnd(1)) {
270+
scale *= 0.5_rt;
271+
}
272+
#endif
273+
#if (AMREX_SPACEDIM == 3)
274+
if (k == nddom.smallEnd(2) ||
275+
k == nddom.bigEnd(2)) {
276+
scale *= 0.5_rt;
277+
}
278+
#endif
279+
return { mask_ma[box_no](i,j,k) * rhs_ma[box_no](i,j,k),
280+
mask_ma[box_no](i,j,k) * scale };
281+
});
282+
283+
Real s1 = amrex::get<0>(r);
284+
Real s2 = amrex::get<1>(r);
285+
ParallelAllReduce::Sum<Real>({s1,s2}, ParallelContext::CommunicatorSub());
286+
return s1/s2;
287+
}
288+
} else {
289+
return MLNodeLinOp::getSolvabilityOffset(amrlev, mglev, rhs);
290+
}
291+
}
292+
293+
void
294+
MLNodeLaplacian::fixSolvabilityByOffset (int amrlev, int mglev, MultiFab& rhs, Real offset) const
295+
{
296+
if (m_coarsening_strategy == CoarseningStrategy::RAP) {
297+
#ifdef AMREX_USE_EB
298+
auto factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][0].get());
299+
if (factory && !factory->isAllRegular()) {
300+
if (mglev == 0) {
301+
const MultiFab& vfrac = factory->getVolFrac();
302+
const auto& vfrac_ma = vfrac.const_arrays();
303+
304+
Box dom = Geom(amrlev,mglev).Domain();
305+
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
306+
if (m_lobc[0][idim] != LinOpBCType::Neumann &&
307+
m_lobc[0][idim] != LinOpBCType::inflow)
308+
{
309+
dom.growLo(idim, 10);
310+
}
311+
if (m_hibc[0][idim] != LinOpBCType::Neumann &&
312+
m_hibc[0][idim] != LinOpBCType::inflow)
313+
{
314+
dom.growHi(idim, 10);
315+
}
316+
}
317+
318+
auto const& rhs_ma = rhs.arrays();
319+
ParallelFor(rhs, IntVect(0),
320+
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
321+
{
322+
Real scale = 0.0_rt;
323+
#if (AMREX_SPACEDIM == 3)
324+
int const koff = 1;
325+
Real const fac = 0.125_rt;
326+
#else
327+
int const koff = 0;
328+
Real const fac = 0.25_rt;
329+
#endif
330+
for (int kc = k-koff; kc <= k; ++kc) {
331+
for (int jc = j-1 ; jc <= j; ++jc) {
332+
for (int ic = i-1 ; ic <= i; ++ic) {
333+
if (dom.contains(ic,jc,kc)) {
334+
scale += vfrac_ma[box_no](ic,jc,kc) * fac;
335+
}
336+
}}}
337+
rhs_ma[box_no](i,j,k) -= offset * scale;
338+
});
339+
}
340+
} else
341+
#endif
342+
{
343+
Box nddom = amrex::surroundingNodes(Geom(amrlev,mglev).Domain());
344+
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
345+
if (m_lobc[0][idim] != LinOpBCType::Neumann &&
346+
m_lobc[0][idim] != LinOpBCType::inflow)
347+
{
348+
nddom.growLo(idim, 10); // so that the test in ParReduce will faill
349+
}
350+
if (m_hibc[0][idim] != LinOpBCType::Neumann &&
351+
m_hibc[0][idim] != LinOpBCType::inflow)
352+
{
353+
nddom.growHi(idim, 10);
354+
}
355+
}
356+
357+
auto const& rhs_ma = rhs.arrays();
358+
ParallelFor(rhs, IntVect(0),
359+
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
360+
{
361+
Real scale = 1.0_rt;
362+
if (i == nddom.smallEnd(0) ||
363+
i == nddom.bigEnd(0)) {
364+
scale *= 0.5_rt;
365+
}
366+
#if (AMREX_SPACEDIM >= 2)
367+
if (j == nddom.smallEnd(1) ||
368+
j == nddom.bigEnd(1)) {
369+
scale *= 0.5_rt;
370+
}
371+
#endif
372+
#if (AMREX_SPACEDIM == 3)
373+
if (k == nddom.smallEnd(2) ||
374+
k == nddom.bigEnd(2)) {
375+
scale *= 0.5_rt;
376+
}
377+
#endif
378+
rhs_ma[box_no](i,j,k) -= offset * scale;
379+
});
380+
}
381+
Gpu::streamSynchronize();
382+
} else {
383+
rhs.plus(-offset, 0, 1);
384+
}
385+
}
386+
174387
void
175388
MLNodeLaplacian::setSigma (int amrlev, const MultiFab& a_sigma)
176389
{

‎Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -67,7 +67,10 @@ public:
6767
amrex::Abort("AMReX_MLNodeLinOp::fillSolutionBC::How did we get here?");
6868
}
6969

70-
virtual void applyInhomogNeumannTerm (int armlev, MultiFab& rhs) const override;
70+
virtual void applyInhomogNeumannTerm (int amrlev, MultiFab& rhs) const override;
71+
72+
virtual Real getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rhs) const override;
73+
virtual void fixSolvabilityByOffset (int amrlev, int mglev, MultiFab& rhs, Real offset) const override;
7174

7275
virtual void prepareForSolve () override {}
7376

‎Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.cpp‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -182,6 +182,36 @@ MLNodeLinOp::applyInhomogNeumannTerm (int /*amrlev*/, MultiFab& /*rhs*/) const
182182
{
183183
}
184184

185+
Real
186+
MLNodeLinOp::getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rhs) const
187+
{
188+
amrex::ignore_unused(amrlev);
189+
AMREX_ASSERT(amrlev==0);
190+
AMREX_ASSERT(mglev+1==m_num_mg_levels[0] || mglev==0);
191+
const auto& mask = (mglev+1 == m_num_mg_levels[0]) ? m_bottom_dot_mask : m_coarse_dot_mask;
192+
const auto& mask_ma = mask.const_arrays();
193+
const auto& rhs_ma = rhs.const_arrays();
194+
auto r = ParReduce(TypeList<ReduceOpSum,ReduceOpSum>{}, TypeList<Real,Real>{},
195+
rhs, IntVect(0),
196+
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
197+
-> GpuTuple<Real,Real>
198+
{
199+
return { mask_ma[box_no](i,j,k) * rhs_ma[box_no](i,j,k),
200+
mask_ma[box_no](i,j,k) };
201+
});
202+
203+
Real s1 = amrex::get<0>(r);
204+
Real s2 = amrex::get<1>(r);
205+
ParallelAllReduce::Sum<Real>({s1,s2}, ParallelContext::CommunicatorSub());
206+
return s1/s2;
207+
}
208+
209+
void
210+
MLNodeLinOp::fixSolvabilityByOffset (int /*amrlev*/, int /*mglev*/, MultiFab& rhs, Real offset) const
211+
{
212+
rhs.plus(-offset, 0, 1);
213+
}
214+
185215
namespace {
186216

187217
void MLNodeLinOp_set_dot_mask (MultiFab& dot_mask, iMultiFab const& omask, Geometry const& geom,

0 commit comments

Comments
 (0)