Skip to content

Commit e4d1756

Browse files
eebassoeebasso
committed
Initial commit
1 parent fa3743f commit e4d1756

File tree

10 files changed

+5187
-2
lines changed

10 files changed

+5187
-2
lines changed

‎Src/Base/AMReX_FabArray.H‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -229,6 +229,7 @@ template <class DFAB, class SFAB,
229229
void
230230
Copy (FabArray<DFAB>& dst, FabArray<SFAB> const& src, int srccomp, int dstcomp, int numcomp, const IntVect& nghost)
231231
{
232+
amrex::Print() << "amrex::Copy start" << "\n";
232233
BL_PROFILE("amrex::Copy()");
233234

234235
using DT = typename DFAB::value_type;
@@ -263,6 +264,7 @@ Copy (FabArray<DFAB>& dst, FabArray<SFAB> const& src, int srccomp, int dstcomp,
263264
}
264265
}
265266
}
267+
amrex::Print() << "amrex::Copy end" << "\n";
266268
}
267269

268270
template <class FAB,
@@ -1765,7 +1767,9 @@ void
17651767
FabArray<FAB>::LocalCopy (FabArray<SFAB> const& src, int scomp, int dcomp, int ncomp,
17661768
IntVect const& nghost)
17671769
{
1770+
amrex::Print() << "FabArray<FAB>::LocalCopy start" << "\n";
17681771
amrex::Copy(*this, src, scomp, dcomp, ncomp, nghost);
1772+
amrex::Print() << "FabArray<FAB>::LocalCopy end" << "\n";
17691773
}
17701774

17711775
template <class FAB>
@@ -2454,6 +2458,7 @@ FabArray<FAB>::setVal (value_type val,
24542458
int ncomp,
24552459
const IntVect& nghost)
24562460
{
2461+
amrex::Print() << "FabArray<FAB>::setVal(val,comp,ncomp,nghost) start" << "\n";
24572462
BL_ASSERT(nghost.allGE(IntVect::TheZeroVector()) && nghost.allLE(n_grow));
24582463
BL_ASSERT(comp+ncomp <= n_comp);
24592464

@@ -2486,6 +2491,7 @@ FabArray<FAB>::setVal (value_type val,
24862491
});
24872492
}
24882493
}
2494+
amrex::Print() << "FabArray<FAB>::setVal(val,comp,ncomp,nghost) end" << "\n";
24892495
}
24902496

24912497
template <class FAB>

‎Src/Base/AMReX_Loop.H‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include <AMReX_Config.H>
44

55
#include <AMReX_Box.H>
6+
#include <AMReX_Print.H>
67

78
namespace amrex {
89

@@ -187,6 +188,7 @@ void LoopOnCpu (Box const& bx, int ncomp, F&& f) noexcept
187188
template <class F>
188189
void LoopConcurrentOnCpu (Box const& bx, F&& f) noexcept
189190
{
191+
// amrex::Print() << "LoopConcurrentOnCpu start" << "\n";
190192
const auto lo = amrex::lbound(bx);
191193
const auto hi = amrex::ubound(bx);
192194
for (int k = lo.z; k <= hi.z; ++k) {
@@ -195,6 +197,7 @@ void LoopConcurrentOnCpu (Box const& bx, F&& f) noexcept
195197
for (int i = lo.x; i <= hi.x; ++i) {
196198
f(i,j,k);
197199
}}}
200+
// amrex::Print() << "LoopConcurrentOnCpu end" << "\n";
198201
}
199202

200203
template <class F>

‎Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H‎

Lines changed: 38 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,12 @@ public:
4545
void setNGhost(int _nghost) {nghost = IntVect(_nghost);}
4646
[[nodiscard]] int getNGhost() {return nghost[0];}
4747

48+
void setAMRLevel (int _amrlev) { amrlev = _amrlev; }
49+
[[nodiscard]] int getAMRLevel () const { return amrlev; }
50+
51+
void setMGLevel (int _mglev) { mglev = _mglev; }
52+
[[nodiscard]] int getMGLevel () const { return mglev; }
53+
4854
[[nodiscard]] RT dotxy (const MF& r, const MF& z, bool local = false);
4955
[[nodiscard]] RT norm_inf (const MF& res, bool local = false);
5056
int solve_bicgstab (MF& solnL, const MF& rhsL, RT eps_rel, RT eps_abs);
@@ -56,8 +62,8 @@ private:
5662

5763
MLLinOpT<MF>& Lp;
5864
Type solver_type;
59-
const int amrlev = 0;
60-
const int mglev;
65+
int amrlev = 0;
66+
int mglev;
6167
int verbose = 0;
6268
int maxiter = 100;
6369
IntVect nghost = IntVect(0);
@@ -87,14 +93,18 @@ int
8793
MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
8894
{
8995
BL_PROFILE("MLCGSolver::bicgstab");
96+
amrex::Print() << "MLCGSolver::solve_bicgstab start" << "\n";
9097

9198
const int ncomp = sol.nComp();
9299

100+
amrex::Print() << "MLCGSolver_BiCGStab: make ph and sh" << "\n";
93101
MF ph = Lp.make(amrlev, mglev, sol.nGrowVect());
94102
MF sh = Lp.make(amrlev, mglev, sol.nGrowVect());
103+
amrex::Print() << "MLCGSolver_BiCGStab: ph.setVal(0) and sh.setVal(0)" << "\n";
95104
ph.setVal(RT(0.0));
96105
sh.setVal(RT(0.0));
97106

107+
amrex::Print() << "MLCGSolver_BiCGStab: make sorig,p,r,s,rh,v,t" << "\n";
98108
MF sorig = Lp.make(amrlev, mglev, nghost);
99109
MF p = Lp.make(amrlev, mglev, nghost);
100110
MF r = Lp.make(amrlev, mglev, nghost);
@@ -103,16 +113,22 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
103113
MF v = Lp.make(amrlev, mglev, nghost);
104114
MF t = Lp.make(amrlev, mglev, nghost);
105115

116+
amrex::Print() << "MLCGSolver_BiCGStab: Lp.correctionResidual(r,sol,rhs)" << "\n";
106117
Lp.correctionResidual(amrlev, mglev, r, sol, rhs, MLLinOpT<MF>::BCMode::Homogeneous);
107118

119+
amrex::Print() << "MLCGSolver_BiCGStab: Lp.normalize(r)" << "\n";
108120
// Then normalize
109121
Lp.normalize(amrlev, mglev, r);
110122

123+
amrex::Print() << "MLCGSolver_BiCGStab: LocalCopy sorig from sol" << "\n";
111124
sorig.LocalCopy(sol,0,0,ncomp,nghost);
125+
amrex::Print() << "MLCGSolver_BiCGStab: LocalCopy rh from r" << "\n";
112126
rh.LocalCopy (r ,0,0,ncomp,nghost);
113127

128+
amrex::Print() << "MLCGSolver_BiCGStab: sol.setVal(0)" << "\n";
114129
sol.setVal(RT(0.0));
115130

131+
amrex::Print() << "MLCGSolver_BiCGStab: rnorm = norm_inf(r)" << "\n";
116132
RT rnorm = norm_inf(r);
117133
const RT rnorm0 = rnorm;
118134

@@ -137,25 +153,33 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
137153

138154
for (; iter <= maxiter; ++iter)
139155
{
156+
amrex::Print() << "solve_bicgstab: start iter = " << iter << "\n";
157+
amrex::Print() << "Calculate rho" << "\n";
140158
const RT rho = dotxy(rh,r);
141159
if ( rho == 0 )
142160
{
143161
ret = 1; break;
144162
}
145163
if ( iter == 1 )
146164
{
165+
amrex::Print() << "Local copy p from r" << "\n";
147166
p.LocalCopy(r,0,0,ncomp,nghost);
148167
}
149168
else
150169
{
170+
amrex::Print() << "Calculate p" << "\n";
151171
const RT beta = (rho/rho_1)*(alpha/omega);
152172
MF::Saxpy(p, -omega, v, 0, 0, ncomp, nghost); // p += -omega*v
153173
MF::Xpay(p, beta, r, 0, 0, ncomp, nghost); // p = r + beta*p
154174
}
175+
amrex::Print() << "Local copy ph from p" << "\n";
155176
ph.LocalCopy(p,0,0,ncomp,nghost);
177+
amrex::Print() << "Lp.apply(v,ph)" << "\n";
156178
Lp.apply(amrlev, mglev, v, ph, MLLinOpT<MF>::BCMode::Homogeneous, MLLinOpT<MF>::StateMode::Correction);
179+
amrex::Print() << "Lp.normalize(v)" << "\n";
157180
Lp.normalize(amrlev, mglev, v);
158181

182+
amrex::Print() << "Calculate rhTv" << "\n";
159183
RT rhTv = dotxy(rh,v);
160184
if ( rhTv != RT(0.0) )
161185
{
@@ -165,9 +189,12 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
165189
{
166190
ret = 2; break;
167191
}
192+
amrex::Print() << "Saxpy sol += alpha*ph" << "\n";
168193
MF::Saxpy(sol, alpha, ph, 0, 0, ncomp, nghost); // sol += alpha * ph
194+
amrex::Print() << "Lincomb s = r - alpha*v" << "\n";
169195
MF::LinComb(s, RT(1.0), r, 0, -alpha, v, 0, 0, ncomp, nghost); // s = r - alpha * v
170196

197+
amrex::Print() << "rnorm = norm_inf(s)" << "\n";
171198
rnorm = norm_inf(s);
172199

173200
if ( verbose > 2 && ParallelDescriptor::IOProcessor() )
@@ -180,14 +207,18 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
180207

181208
if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) { break; }
182209

210+
amrex::Print() << "LocalCopy sh from s" << "\n";
183211
sh.LocalCopy(s,0,0,ncomp,nghost);
212+
amrex::Print() << "Lp.apply(t,sh)" << "\n";
184213
Lp.apply(amrlev, mglev, t, sh, MLLinOpT<MF>::BCMode::Homogeneous, MLLinOpT<MF>::StateMode::Correction);
214+
amrex::Print() << "Lp.normalize(t)" << "\n";
185215
Lp.normalize(amrlev, mglev, t);
186216
//
187217
// This is a little funky. I want to elide one of the reductions
188218
// in the following two dotxy()s. We do that by calculating the "local"
189219
// values and then reducing the two local values at the same time.
190220
//
221+
amrex::Print() << "tvals[2]" << "\n";
191222
RT tvals[2] = { dotxy(t,t,true), dotxy(t,s,true) };
192223

193224
BL_PROFILE_VAR("MLCGSolver::ParallelAllReduce", blp_par);
@@ -202,9 +233,12 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
202233
{
203234
ret = 3; break;
204235
}
236+
amrex::Print() << "Saxpy sol += omega * sh" << "\n";
205237
MF::Saxpy(sol, omega, sh, 0, 0, ncomp, nghost); // sol += omega * sh
238+
amrex::Print() << "Lincomb r = s - omega * t" << "\n";
206239
MF::LinComb(r, RT(1.0), s, 0, -omega, t, 0, 0, ncomp, nghost); // r = s - omega * t
207240

241+
amrex::Print() << "rnorm = norm_inf(r)" << "\n";
208242
rnorm = norm_inf(r);
209243

210244
if ( verbose > 2 )
@@ -222,6 +256,7 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
222256
ret = 4; break;
223257
}
224258
rho_1 = rho;
259+
amrex::Print() << "solve_bicgstab: end iter = " << iter << "\n";
225260
}
226261

227262
if ( verbose > 0 )
@@ -250,6 +285,7 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
250285
sol.LocalAdd(sorig, 0, 0, ncomp, nghost);
251286
}
252287

288+
amrex::Print() << "MLCGSolver::solve_bicgstab end" << "\n";
253289
return ret;
254290
}
255291

‎Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -850,6 +850,7 @@ MLEBABecLap::compGrad (int amrlev, const Array<MultiFab*,AMREX_SPACEDIM>& grad,
850850
void
851851
MLEBABecLap::normalize (int amrlev, int mglev, MultiFab& mf) const
852852
{
853+
amrex::Print() << "MLEBABecLap::normalize start" << "\n";
853854
const MultiFab& acoef = m_a_coeffs[amrlev][mglev];
854855
AMREX_D_TERM(const MultiFab& bxcoef = m_b_coeffs[amrlev][mglev][0];,
855856
const MultiFab& bycoef = m_b_coeffs[amrlev][mglev][1];,
@@ -938,6 +939,7 @@ MLEBABecLap::normalize (int amrlev, int mglev, MultiFab& mf) const
938939
});
939940
}
940941
}
942+
amrex::Print() << "MLEBABecLap::normalize end" << "\n";
941943
}
942944

943945
void

‎Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_F.cpp‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ void
1212
MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const
1313
{
1414
BL_PROFILE("MLEBABecLap::Fapply()");
15+
amrex::Print() << "MLEBABecLap::Fapply: start" << "\n";
1516

1617
const MultiFab& acoef = m_a_coeffs[amrlev][mglev];
1718
AMREX_D_TERM(const MultiFab& bxcoef = m_b_coeffs[amrlev][mglev][0];,
@@ -32,6 +33,8 @@ MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) c
3233
const MultiCutFab* bcent = (factory) ? &(factory->getBndryCent()) : nullptr;
3334
const auto *const ccent = (factory) ? &(factory->getCentroid()) : nullptr;
3435

36+
amrex::Print() << "MLEBABecLap::Fapply: step 1" << "\n";
37+
3538
const bool is_eb_dirichlet = isEBDirichlet();
3639
const bool is_eb_inhomog = m_is_eb_inhomog;
3740

@@ -44,6 +47,8 @@ MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) c
4447

4548
const Box& domain_box = m_geom[amrlev][mglev].Domain();
4649

50+
amrex::Print() << "MLEBABecLap::Fapply: step 2" << "\n";
51+
4752
AMREX_D_TERM(
4853
const int domlo_x = domain_box.smallEnd(0);
4954
const int domhi_x = domain_box.bigEnd(0);,
@@ -58,12 +63,15 @@ MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) c
5863
const bool extdir_z = !(m_geom[amrlev][mglev].isPeriodic(2)););
5964

6065
MFItInfo mfi_info;
66+
amrex::Print() << "MLEBABecLap::Fapply: step 3" << "\n";
6167
if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
6268
#ifdef AMREX_USE_OMP
6369
#pragma omp parallel if (Gpu::notInLaunchRegion())
6470
#endif
71+
amrex::Print() << "MLEBABecLap::Fapply: step 4" << "\n";
6572
for (MFIter mfi(out, mfi_info); mfi.isValid(); ++mfi)
6673
{
74+
amrex::Print() << "MLEBABecLap::Fapply: MFIter loop start" << "\n";
6775
const Box& bx = mfi.tilebox();
6876
Array4<Real const> const& xfab = in.const_array(mfi);
6977
Array4<Real> const& yfab = out.array(mfi);
@@ -74,19 +82,23 @@ MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) c
7482

7583
auto fabtyp = (flags) ? (*flags)[mfi].getType(bx) : FabType::regular;
7684

85+
amrex::Print() << "MLEBABecLap::Fapply: Check FabType" << "\n";
7786
if (fabtyp == FabType::covered) {
87+
amrex::Print() << "MLEBABecLap::Fapply: fabtyp == FabType::covered" << "\n";
7888
AMREX_HOST_DEVICE_PARALLEL_FOR_4D( bx, ncomp, i, j, k, n,
7989
{
8090
yfab(i,j,k,n) = 0.0;
8191
});
8292
} else if (fabtyp == FabType::regular) {
93+
amrex::Print() << "MLEBABecLap::Fapply: fabtyp == FabType::regular" << "\n";
8394
AMREX_HOST_DEVICE_PARALLEL_FOR_4D( bx, ncomp, i, j, k, n,
8495
{
8596
mlabeclap_adotx(i,j,k,n, yfab, xfab, afab,
8697
AMREX_D_DECL(bxfab,byfab,bzfab),
8798
dxinvarr, ascalar, bscalar);
8899
});
89100
} else {
101+
amrex::Print() << "MLEBABecLap::Fapply: fabtyp == FabType::singlevalued" << "\n";
90102
Array4<int const> const& ccmfab = ccmask.const_array(mfi);
91103
Array4<EBCellFlag const> const& flagfab = flags->const_array(mfi);
92104
Array4<Real const> const& vfracfab = vfrac->const_array(mfi);
@@ -121,6 +133,7 @@ MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) c
121133
AMREX_D_DECL(extdir_x, extdir_y, extdir_z));
122134
amrex::ignore_unused(ccfab);
123135
#else
136+
amrex::Print() << "MLEBABecLap::Fapply mlebabeclap_adotx_centroid" << "\n";
124137
AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bx, tbx,
125138
{
126139
mlebabeclap_adotx_centroid(tbx, yfab, xfab, afab, AMREX_D_DECL(bxfab,byfab,bzfab),
@@ -136,6 +149,7 @@ MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) c
136149
});
137150
#endif
138151
} else {
152+
amrex::Print() << "MLEBABecLap::Fapply mlebabeclap_adotx" << "\n";
139153
AMREX_LAUNCH_HOST_DEVICE_LAMBDA ( bx, tbx,
140154
{
141155
mlebabeclap_adotx(tbx, yfab, xfab, afab, AMREX_D_DECL(bxfab,byfab,bzfab),
@@ -151,6 +165,7 @@ MLEBABecLap::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) c
151165
}
152166
}
153167
}
168+
amrex::Print() << "MLEBABecLap::Fapply end" << "\n";
154169
}
155170

156171
void

‎Src/LinearSolvers/MLMG/AMReX_MLMG.H‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@ public:
2626

2727
using BCMode = typename MLLinOpT<MF>::BCMode;
2828
using Location = typename MLLinOpT<MF>::Location;
29+
using StateMode = typename MLLinOpT<MF>::StateMode;
2930

3031
using BottomSolver = amrex::BottomSolver;
3132
enum class CFStrategy : int {none,ghostnodes};
@@ -174,6 +175,8 @@ public:
174175
void mgVcycle (int amrlev, int mglev);
175176
void mgFcycle ();
176177

178+
// void smooth (int amrlev, int mglev, MF& sol, const MF& rhs, bool skip_fillboundary=false);
179+
177180
void bottomSolve ();
178181
void NSolve (MLMGT<MF>& a_solver, MF& a_sol, MF& a_rhs);
179182
void actualBottomSolve ();
@@ -346,6 +349,8 @@ MLMGT<MF>::solve (const Vector<AMF*>& a_sol, const Vector<AMF const*>& a_rhs,
346349
{
347350
BL_PROFILE("MLMG::solve()");
348351

352+
verbose = 5;
353+
349354
if constexpr (std::is_same<AMF,MultiFab>()) {
350355
if (checkpoint_file != nullptr) {
351356
checkPoint(a_sol, a_rhs, a_tol_rel, a_tol_abs, checkpoint_file);
@@ -1370,6 +1375,20 @@ MLMGT<MF>::mgFcycle ()
13701375
}
13711376
}
13721377

1378+
// template <typename MF>
1379+
// void
1380+
// MLMGT<MF>::smooth (int amrlev, int mglev, MF& sol, const MF& rhs, bool skip_fillboundary)
1381+
// {
1382+
// linop.smooth(amrlev, mglev, sol, rhs, skip_fillboundary);
1383+
// MLCGSolverT<MF> cg_solver(linop);
1384+
// cg_solver.setAMRLevel(amrlev);
1385+
// cg_solver.setMGLevel(mglev);
1386+
// cg_solver.setSolver(MLCGSolverT<MF>::Type::BiCGStab);
1387+
// cg_solver.setVerbose(5);
1388+
// cg_solver.setMaxIter(2);
1389+
// int ret = cg_solver.solve(sol, rhs, RT(-1.0), RT(-1.0));
1390+
// }
1391+
13731392
// At the true bottom of the coarset AMR level.
13741393
// in : Residual (res) as b
13751394
// out : Correction (cor) as x

‎Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -839,6 +839,7 @@ void
839839
MLNodeLaplacian::normalize (int amrlev, int mglev, MultiFab& mf) const
840840
{
841841
BL_PROFILE("MLNodeLaplacian::normalize()");
842+
amrex::Print() << "MLNodeLaplacian::normalize start" << "\n";
842843

843844
if (m_sigma[0][0][0] == nullptr) { return; }
844845

@@ -927,6 +928,7 @@ MLNodeLaplacian::normalize (int amrlev, int mglev, MultiFab& mf) const
927928
}
928929
}
929930
}
931+
amrex::Print() << "MLNodeLaplacian::normalize end" << "\n";
930932
}
931933

932934
void

0 commit comments

Comments
 (0)