Skip to content

Commit 78268f1

Browse files
WeiqunZhangWeiqunZhang
authored
Reuse MLEBNodeFDLaplacian and MLCurlCurl (#4480)
Allow these two linear operators to be reused.
1 parent fc469a3 commit 78268f1

File tree

4 files changed

+105
-45
lines changed

4 files changed

+105
-45
lines changed

‎Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -82,6 +82,11 @@ public:
8282
const MF& b, BCMode bc_mode,
8383
const MF* crse_bcdata=nullptr) override;
8484

85+
[[nodiscard]] bool needsUpdate () const override {
86+
return (m_needs_update || MLLinOpT<Array<MultiFab,3>>::needsUpdate());
87+
}
88+
void update () override;
89+
8590
void prepareForSolve () override;
8691

8792
[[nodiscard]] bool isSingular (int /*amrlev*/) const override { return false; }
@@ -128,6 +133,8 @@ private:
128133
[[nodiscard]] CurlCurlDirichletInfo getDirichletInfo (int amrlev, int mglev) const;
129134
[[nodiscard]] CurlCurlSymmetryInfo getSymmetryInfo (int amrlev, int mglev) const;
130135

136+
void update_lusolver ();
137+
131138
RT m_alpha = std::numeric_limits<RT>::lowest();
132139
RT m_beta = std::numeric_limits<RT>::lowest();
133140

@@ -146,6 +153,7 @@ private:
146153
<LUSolver<AMREX_SPACEDIM*2,RT>>>>> m_lusolver;
147154
Vector<Vector<Array<std::unique_ptr<MultiFab>,3>>> m_bcoefs;
148155
bool m_use_pcg = false;
156+
bool m_needs_update = true;
149157
};
150158

151159
}

‎Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp‎

Lines changed: 21 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,13 +35,16 @@ void MLCurlCurl::define (const Vector<Geometry>& a_geom,
3535

3636
void MLCurlCurl::setScalars (RT a_alpha, RT a_beta) noexcept
3737
{
38+
m_needs_update = true;
3839
m_alpha = a_alpha;
3940
m_beta = a_beta;
4041
AMREX_ASSERT(m_beta > RT(0));
4142
}
4243

4344
void MLCurlCurl::setBeta (const Vector<Array<MultiFab const*,3>>& a_bcoefs)
4445
{
46+
m_needs_update = true;
47+
4548
Array<IntVect,3> ng;
4649
for (int idim = 0; idim < 3; ++idim) {
4750
ng[idim] = IntVect(1) - m_etype[idim]; // 1 ghost for cell direction, 0 for node
@@ -522,7 +525,7 @@ void MLCurlCurl::compresid (int amrlev, int mglev, MF& resid, MF const& b) const
522525
}
523526
}
524527

525-
void MLCurlCurl::prepareForSolve ()
528+
void MLCurlCurl::update_lusolver ()
526529
{
527530
#if (AMREX_SPACEDIM > 1)
528531
if (m_bcoefs[0][0][0] == nullptr) {
@@ -611,6 +614,11 @@ void MLCurlCurl::prepareForSolve ()
611614
#endif
612615
}
613616

617+
void MLCurlCurl::prepareForSolve ()
618+
{
619+
update_lusolver();
620+
}
621+
614622
Real MLCurlCurl::xdoty (int amrlev, int mglev, const MF& x, const MF& y,
615623
bool local) const
616624
{
@@ -942,4 +950,16 @@ CurlCurlSymmetryInfo MLCurlCurl::getSymmetryInfo (int amrlev, int mglev) const
942950
helper(2,1)))};
943951
}
944952

953+
void MLCurlCurl::update ()
954+
{
955+
if (MLLinOpT<Array<MultiFab,3>>::needsUpdate()) {
956+
MLLinOpT<Array<MultiFab,3>>::update();
957+
}
958+
959+
if (m_needs_update) {
960+
update_lusolver();
961+
m_needs_update = false;
962+
}
963+
}
964+
945965
}

‎Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.H‎

Lines changed: 13 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -90,6 +90,11 @@ public:
9090
void restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const final;
9191
void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const final;
9292

93+
[[nodiscard]] bool needsUpdate () const override {
94+
return (m_needs_update || MLNodeLinOp::needsUpdate());
95+
}
96+
void update () override;
97+
9398
void prepareForSolve () final;
9499
void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const final;
95100
void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs) const final;
@@ -123,10 +128,13 @@ private:
123128
GpuArray<Real,AMREX_SPACEDIM> m_sigma{{AMREX_D_DECL(1_rt,1_rt,1_rt)}};
124129
Vector<Vector<std::unique_ptr<MultiFab>>> m_sigma_mf;
125130
bool m_has_sigma_mf = false;
131+
bool m_needs_update = true;
126132
Real m_s_phi_eb = std::numeric_limits<Real>::lowest();
127133
Vector<MultiFab> m_phi_eb;
128134
int m_rz = false;
129135
Real m_rz_alpha = 0._rt;
136+
137+
void update_sigma ();
130138
};
131139

132140
#ifdef AMREX_USE_EB
@@ -142,9 +150,11 @@ MLEBNodeFDLaplacian::setEBDirichlet (F const& f)
142150
Geometry const& geom = m_geom[amrlev][0];
143151
auto const problo = geom.ProbLoArray();
144152
auto const cellsize = geom.CellSizeArray();
145-
m_phi_eb[amrlev].define(amrex::convert(m_grids[amrlev][0],IntVect(1)),
146-
m_dmap[amrlev][0], 1, 1);
147-
m_phi_eb[amrlev].setVal(0.0);
153+
if (m_phi_eb[amrlev].empty()) {
154+
m_phi_eb[amrlev].define(amrex::convert(m_grids[amrlev][0],IntVect(1)),
155+
m_dmap[amrlev][0], 1, 1);
156+
m_phi_eb[amrlev].setVal(0.0);
157+
}
148158
auto const& flags = factory->getMultiEBCellFlagFab();
149159
auto const& levset = factory->getLevelSet();
150160
#ifdef AMREX_USE_OMP

‎Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.cpp‎

Lines changed: 63 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,7 @@ MLEBNodeFDLaplacian::setSigma (Array<Real,AMREX_SPACEDIM> const& a_sigma) noexce
4343
void
4444
MLEBNodeFDLaplacian::setSigma (int amrlev, MultiFab const& a_sigma)
4545
{
46+
m_needs_update = true;
4647
m_has_sigma_mf = true;
4748
m_sigma_mf[amrlev][0] = std::make_unique<MultiFab>
4849
(this->m_grids[amrlev][0], this->m_dmap[amrlev][0], 1, 1, MFInfo{},
@@ -318,47 +319,7 @@ MLEBNodeFDLaplacian::prepareForSolve ()
318319
#endif
319320

320321
if (m_has_sigma_mf) {
321-
AMREX_D_TERM(m_sigma[0] = Real(1.0);,
322-
m_sigma[1] = Real(1.0);,
323-
m_sigma[2] = Real(1.0));
324-
AMREX_ALWAYS_ASSERT(this->m_num_amr_levels == 1);
325-
for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev) {
326-
for (int mglev = 1; mglev < this->m_num_mg_levels[amrlev]; ++mglev) {
327-
m_sigma_mf[amrlev][mglev] = std::make_unique<MultiFab>
328-
(this->m_grids[amrlev][mglev], this->m_dmap[amrlev][mglev], 1, 1,
329-
MFInfo{}, *(this->m_factory[amrlev][mglev]));
330-
IntVect const ratio = (amrlev > 0) ? IntVect (2)
331-
: this->mg_coarsen_ratio_vec[mglev-1];
332-
#ifdef AMREX_USE_EB
333-
amrex::EB_average_down
334-
#else
335-
amrex::average_down
336-
#endif
337-
(*m_sigma_mf[amrlev][mglev-1],
338-
*m_sigma_mf[amrlev][mglev], 0, 1, ratio);
339-
}
340-
341-
for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev) {
342-
auto const& geom = this->m_geom[amrlev][mglev];
343-
auto& sigma = *m_sigma_mf[amrlev][mglev];
344-
sigma.FillBoundary(geom.periodicity());
345-
346-
const Box& domain = geom.Domain();
347-
const auto lobc = LoBC();
348-
const auto hibc = HiBC();
349-
350-
MFItInfo mfi_info;
351-
if (Gpu::notInLaunchRegion()) { mfi_info.SetDynamic(true); }
352-
#ifdef AMREX_USE_OMP
353-
#pragma omp parallel if (Gpu::notInLaunchRegion())
354-
#endif
355-
for (MFIter mfi(sigma, mfi_info); mfi.isValid(); ++mfi)
356-
{
357-
Array4<Real> const& sfab = sigma.array(mfi);
358-
mlndlap_fillbc_cc<Real>(mfi.validbox(),sfab,domain,lobc,hibc);
359-
}
360-
}
361-
}
322+
update_sigma();
362323
}
363324
}
364325

@@ -792,4 +753,65 @@ MLEBNodeFDLaplacian::postSolve (Vector<MultiFab*> const& sol) const
792753
#endif
793754
}
794755

756+
void
757+
MLEBNodeFDLaplacian::update ()
758+
{
759+
if (MLNodeLinOp::needsUpdate()) {
760+
MLNodeLinOp::update();
761+
}
762+
763+
if (m_needs_update && m_has_sigma_mf) {
764+
update_sigma();
765+
}
766+
m_needs_update = false;
767+
}
768+
769+
void
770+
MLEBNodeFDLaplacian::update_sigma ()
771+
{
772+
AMREX_D_TERM(m_sigma[0] = Real(1.0);,
773+
m_sigma[1] = Real(1.0);,
774+
m_sigma[2] = Real(1.0));
775+
AMREX_ALWAYS_ASSERT(this->m_num_amr_levels == 1);
776+
for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev) {
777+
for (int mglev = 1; mglev < this->m_num_mg_levels[amrlev]; ++mglev) {
778+
if (m_sigma_mf[amrlev][mglev] == nullptr) {
779+
m_sigma_mf[amrlev][mglev] = std::make_unique<MultiFab>
780+
(this->m_grids[amrlev][mglev], this->m_dmap[amrlev][mglev], 1, 1,
781+
MFInfo{}, *(this->m_factory[amrlev][mglev]));
782+
}
783+
IntVect const ratio = (amrlev > 0) ? IntVect (2)
784+
: this->mg_coarsen_ratio_vec[mglev-1];
785+
#ifdef AMREX_USE_EB
786+
amrex::EB_average_down
787+
#else
788+
amrex::average_down
789+
#endif
790+
(*m_sigma_mf[amrlev][mglev-1],
791+
*m_sigma_mf[amrlev][mglev], 0, 1, ratio);
792+
}
793+
794+
for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev) {
795+
auto const& geom = this->m_geom[amrlev][mglev];
796+
auto& sigma = *m_sigma_mf[amrlev][mglev];
797+
sigma.FillBoundary(geom.periodicity());
798+
799+
const Box& domain = geom.Domain();
800+
const auto lobc = LoBC();
801+
const auto hibc = HiBC();
802+
803+
MFItInfo mfi_info;
804+
if (Gpu::notInLaunchRegion()) { mfi_info.SetDynamic(true); }
805+
#ifdef AMREX_USE_OMP
806+
#pragma omp parallel if (Gpu::notInLaunchRegion())
807+
#endif
808+
for (MFIter mfi(sigma, mfi_info); mfi.isValid(); ++mfi)
809+
{
810+
Array4<Real> const& sfab = sigma.array(mfi);
811+
mlndlap_fillbc_cc<Real>(mfi.validbox(),sfab,domain,lobc,hibc);
812+
}
813+
}
814+
}
815+
}
816+
795817
}

0 commit comments

Comments
 (0)