2 #ifndef AMREX_MLCGSOLVER_H_
3 #define AMREX_MLCGSOLVER_H_
4 #include <AMReX_Config.H>
10 template <
typename MF>
37 int solve (MF& solnL,
const MF& rhsL,
RT eps_rel,
RT eps_abs);
58 [[nodiscard]]
RT dotxy (
const MF& r,
const MF& z,
bool local =
false);
59 [[nodiscard]]
RT norm_inf (
const MF& res,
bool local =
false);
61 int solve_cg (MF& solnL,
const MF& rhsL,
RT eps_rel,
RT eps_abs);
78 template <
typename MF>
80 : Lp(_lp), solver_type(_typ), mglev(_lp.NMGLevels(0)-1)
85 template <
typename MF>
89 if (solver_type == Type::BiCGStab) {
90 return solve_bicgstab(sol,rhs,eps_rel,eps_abs);
92 return solve_cg(sol,rhs,eps_rel,eps_abs);
96 template <
typename MF>
102 const int ncomp =
nComp(sol);
104 MF p = Lp.make(amrlev, mglev,
nGrowVect(sol));
105 MF
r = Lp.make(amrlev, mglev,
nGrowVect(sol));
109 MF rh = Lp.make(amrlev, mglev, nghost);
110 MF v = Lp.make(amrlev, mglev, nghost);
111 MF t = Lp.make(amrlev, mglev, nghost);
116 if ( initial_vec_zeroed ) {
119 sorig = Lp.make(amrlev, mglev, nghost);
128 Lp.normalize(amrlev, mglev,
r);
131 RT rnorm = norm_inf(
r);
132 const RT rnorm0 = rnorm;
136 amrex::Print() <<
"MLCGSolver_BiCGStab: Initial error (error0) = " << rnorm0 <<
'\n';
140 RT rho_1 = 0, alpha = 0, omega = 0;
142 if ( rnorm0 == 0 || rnorm0 < eps_abs )
147 <<
", rnorm = " << rnorm
148 <<
", eps_abs = " << eps_abs <<
'\n';
153 for (; iter <= maxiter; ++iter)
155 const RT rho = dotxy(rh,
r);
166 const RT beta = (rho/rho_1)*(alpha/omega);
167 Saxpy(p, -omega, v, 0, 0, ncomp, nghost);
168 Xpay(p, beta,
r, 0, 0, ncomp, nghost);
171 Lp.normalize(amrlev, mglev, v);
173 RT rhTv = dotxy(rh,v);
174 if ( rhTv !=
RT(0.0) )
182 Saxpy(sol, alpha, p, 0, 0, ncomp, nghost);
183 Saxpy(
r, -alpha, v, 0, 0, ncomp, nghost);
190 << std::setw(11) << iter
192 << rnorm/(rnorm0) <<
'\n';
195 if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) {
break; }
198 Lp.normalize(amrlev, mglev, t);
204 RT tvals[2] = { dotxy(t,t,
true), dotxy(t,
r,
true) };
210 if ( tvals[0] !=
RT(0.0) )
212 omega = tvals[1]/tvals[0];
218 Saxpy(sol, omega,
r, 0, 0, ncomp, nghost);
219 Saxpy(
r, -omega, t, 0, 0, ncomp, nghost);
226 << std::setw(11) << iter
228 << rnorm/(rnorm0) <<
'\n';
231 if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) {
break; }
242 amrex::Print() <<
"MLCGSolver_BiCGStab: Final: Iteration "
243 << std::setw(4) << iter
245 << rnorm/(rnorm0) <<
'\n';
248 if ( ret == 0 && rnorm > eps_rel*rnorm0 && rnorm > eps_abs)
256 if ( ( ret == 0 || ret == 8 ) && (rnorm < rnorm0) )
258 if ( !initial_vec_zeroed ) {
259 LocalAdd(sol, sorig, 0, 0, ncomp, nghost);
261 if (ret == 8) { ret = 9; }
266 if ( !initial_vec_zeroed ) {
267 LocalAdd(sol, sorig, 0, 0, ncomp, nghost);
274 template <
typename MF>
280 const int ncomp =
nComp(sol);
282 MF p = Lp.make(amrlev, mglev,
nGrowVect(sol));
285 MF
r = Lp.make(amrlev, mglev, nghost);
286 MF q = Lp.make(amrlev, mglev, nghost);
290 if ( initial_vec_zeroed ) {
293 sorig = Lp.make(amrlev, mglev, nghost);
301 RT rnorm = norm_inf(
r);
302 const RT rnorm0 = rnorm;
306 amrex::Print() <<
"MLCGSolver_CG: Initial error (error0) : " << rnorm0 <<
'\n';
313 if ( rnorm0 == 0 || rnorm0 < eps_abs )
317 <<
", rnorm = " << rnorm
318 <<
", eps_abs = " << eps_abs <<
'\n';
323 for (; iter <= maxiter; ++iter)
338 Xpay(p, beta,
r, 0, 0, ncomp, nghost);
358 <<
" alpha " << alpha <<
'\n';
360 Saxpy(sol, alpha, p, 0, 0, ncomp, nghost);
361 Saxpy(
r, -alpha, q, 0, 0, ncomp, nghost);
367 << std::setw(4) << iter
369 << rnorm/(rnorm0) <<
'\n';
372 if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) {
break; }
380 << std::setw(4) << iter
382 << rnorm/(rnorm0) <<
'\n';
385 if ( ret == 0 && rnorm > eps_rel*rnorm0 && rnorm > eps_abs )
393 if ( ( ret == 0 || ret == 8 ) && (rnorm < rnorm0) )
395 if ( !initial_vec_zeroed ) {
396 LocalAdd(sol, sorig, 0, 0, ncomp, nghost);
398 if (ret == 8) { ret = 9; }
403 if ( !initial_vec_zeroed ) {
404 LocalAdd(sol, sorig, 0, 0, ncomp, nghost);
411 template <
typename MF>
417 RT result = Lp.xdoty(amrlev, mglev,
r, z, local);
422 template <
typename MF>
426 int ncomp =
nComp(res);
#define BL_PROFILE_VAR_START(vname)
Definition: AMReX_BLProfiler.H:562
#define BL_PROFILE(a)
Definition: AMReX_BLProfiler.H:551
#define BL_PROFILE_VAR_STOP(vname)
Definition: AMReX_BLProfiler.H:563
#define BL_PROFILE_VAR(fname, vname)
Definition: AMReX_BLProfiler.H:560
#define BL_PROFILE_VAR_NS(fname, vname)
Definition: AMReX_BLProfiler.H:561
Definition: AMReX_MLCGSolver.H:12
RT dotxy(const MF &r, const MF &z, bool local=false)
Definition: AMReX_MLCGSolver.H:413
int verbose
Definition: AMReX_MLCGSolver.H:71
int iter
Definition: AMReX_MLCGSolver.H:74
MLLinOpT< MF > & Lp
Definition: AMReX_MLCGSolver.H:67
MLCGSolverT< MF > & operator=(const MLCGSolverT< MF > &rhs)=delete
void setSolver(Type _typ) noexcept
Definition: AMReX_MLCGSolver.H:28
bool getInitSolnZeroed() const
Definition: AMReX_MLCGSolver.H:53
void setVerbose(int _verbose)
Definition: AMReX_MLCGSolver.H:39
int getNGhost()
Definition: AMReX_MLCGSolver.H:56
const int mglev
Definition: AMReX_MLCGSolver.H:70
IntVect nghost
Definition: AMReX_MLCGSolver.H:73
int getNumIters() const noexcept
Definition: AMReX_MLCGSolver.H:63
int solve_cg(MF &solnL, const MF &rhsL, RT eps_rel, RT eps_abs)
Definition: AMReX_MLCGSolver.H:276
RT norm_inf(const MF &res, bool local=false)
Definition: AMReX_MLCGSolver.H:424
const int amrlev
Definition: AMReX_MLCGSolver.H:69
bool initial_vec_zeroed
Definition: AMReX_MLCGSolver.H:75
void setInitSolnZeroed(bool _sol_zeroed)
Definition: AMReX_MLCGSolver.H:52
int solve_bicgstab(MF &solnL, const MF &rhsL, RT eps_rel, RT eps_abs)
Definition: AMReX_MLCGSolver.H:98
MLCGSolverT(MLCGSolverT< MF > &&rhs)=delete
typename MLLinOpT< MF >::FAB FAB
Definition: AMReX_MLCGSolver.H:15
int getMaxIter() const
Definition: AMReX_MLCGSolver.H:43
int maxiter
Definition: AMReX_MLCGSolver.H:72
int solve(MF &solnL, const MF &rhsL, RT eps_rel, RT eps_abs)
Definition: AMReX_MLCGSolver.H:87
Type
Definition: AMReX_MLCGSolver.H:18
typename MLLinOpT< MF >::RT RT
Definition: AMReX_MLCGSolver.H:16
void setNGhost(int _nghost)
Definition: AMReX_MLCGSolver.H:55
MLCGSolverT(MLLinOpT< MF > &_lp, Type _typ=Type::BiCGStab)
Definition: AMReX_MLCGSolver.H:79
int getVerbose() const
Definition: AMReX_MLCGSolver.H:40
void setMaxIter(int _maxiter)
Definition: AMReX_MLCGSolver.H:42
Type solver_type
Definition: AMReX_MLCGSolver.H:68
MLCGSolverT(const MLCGSolverT< MF > &rhs)=delete
Definition: AMReX_MLLinOp.H:98
typename FabDataType< MF >::fab_type FAB
Definition: AMReX_MLLinOp.H:108
typename FabDataType< MF >::value_type RT
Definition: AMReX_MLLinOp.H:109
This class provides the user with a few print options.
Definition: AMReX_Print.H:35
void Sum(T &v, MPI_Comm comm)
Definition: AMReX_ParallelReduce.H:204
void Max(KeyValuePair< K, V > &vi, MPI_Comm comm)
Definition: AMReX_ParallelReduce.H:126
bool IOProcessor() noexcept
Is this CPU the I/O Processor? To get the rank number, call IOProcessorNumber()
Definition: AMReX_ParallelDescriptor.H:275
Definition: AMReX_Amr.cpp:49
MF::value_type norminf(MF const &mf, int scomp, int ncomp, IntVect const &nghost, bool local=false)
Definition: AMReX_FabArrayUtility.H:1682
int nComp(FabArrayBase const &fa)
void Saxpy(MF &dst, typename MF::value_type a, MF const &src, int scomp, int dcomp, int ncomp, IntVect const &nghost)
dst += a * src
Definition: AMReX_FabArrayUtility.H:1646
IntVect nGrowVect(FabArrayBase const &fa)
void LocalCopy(DMF &dst, SMF const &src, int scomp, int dcomp, int ncomp, IntVect const &nghost)
dst = src
Definition: AMReX_FabArrayUtility.H:1630
void LocalAdd(MF &dst, MF const &src, int scomp, int dcomp, int ncomp, IntVect const &nghost)
dst += src
Definition: AMReX_FabArrayUtility.H:1638
IntVectND< AMREX_SPACEDIM > IntVect
Definition: AMReX_BaseFwd.H:30
void Xpay(MF &dst, typename MF::value_type a, MF const &src, int scomp, int dcomp, int ncomp, IntVect const &nghost)
dst = src + a * dst
Definition: AMReX_FabArrayUtility.H:1654
void Warning(const std::string &msg)
Print out warning message to cerr.
Definition: AMReX.cpp:231
int verbose
Definition: AMReX_DistributionMapping.cpp:36
void setVal(MF &dst, typename MF::value_type val)
dst = val
Definition: AMReX_FabArrayUtility.H:1607