2#ifndef AMREX_MLCGSOLVER_H_
3#define AMREX_MLCGSOLVER_H_
4#include <AMReX_Config.H>
37 int solve (MF& solnL,
const MF& rhsL,
RT eps_rel,
RT eps_abs);
59 [[nodiscard]]
RT dotxy (
const MF&
r,
const MF& z,
bool local =
false);
60 [[nodiscard]]
RT norm_inf (
const MF& res,
bool local =
false);
62 int solve_cg (MF& solnL,
const MF& rhsL,
RT eps_rel,
RT eps_abs);
82 : Lp(_lp), solver_type(_typ), mglev(_lp.NMGLevels(0)-1)
91 if (solver_type == Type::BiCGStab) {
92 return solve_bicgstab(sol,rhs,eps_rel,eps_abs);
94 return solve_cg(sol,rhs,eps_rel,eps_abs);
104 const int ncomp =
nComp(sol);
106 MF p = Lp.make(amrlev, mglev,
nGrowVect(sol));
107 MF
r = Lp.make(amrlev, mglev,
nGrowVect(sol));
111 MF rh = Lp.make(amrlev, mglev, nghost);
112 MF v = Lp.make(amrlev, mglev, nghost);
113 MF t = Lp.make(amrlev, mglev, nghost);
118 if ( initial_vec_zeroed ) {
121 sorig = Lp.make(amrlev, mglev, nghost);
130 Lp.normalize(amrlev, mglev,
r);
133 RT rnorm = norm_inf(
r);
134 const RT rnorm0 = rnorm;
138 amrex::Print() << print_ident <<
"MLCGSolver_BiCGStab: Initial error (error0) = " << rnorm0 <<
'\n';
142 RT rho_1 = 0, alpha = 0, omega = 0;
144 if ( rnorm0 == 0 || rnorm0 < eps_abs )
148 amrex::Print() << print_ident <<
"MLCGSolver_BiCGStab: niter = 0,"
149 <<
", rnorm = " << rnorm
150 <<
", eps_abs = " << eps_abs <<
'\n';
155 for (; iter <= maxiter; ++iter)
157 const RT rho = dotxy(rh,
r);
168 const RT beta = (rho/rho_1)*(alpha/omega);
169 Saxpy(p, -omega, v, 0, 0, ncomp, nghost);
170 Xpay(p, beta,
r, 0, 0, ncomp, nghost);
173 Lp.normalize(amrlev, mglev, v);
175 RT rhTv = dotxy(rh,v);
176 if ( rhTv !=
RT(0.0) )
184 Saxpy(sol, alpha, p, 0, 0, ncomp, nghost);
185 Saxpy(
r, -alpha, v, 0, 0, ncomp, nghost);
191 amrex::Print() << print_ident <<
"MLCGSolver_BiCGStab: Half Iter "
192 << std::setw(11) << iter
194 << rnorm/(rnorm0) <<
'\n';
197 if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) {
break; }
200 Lp.normalize(amrlev, mglev, t);
206 RT tvals[2] = { dotxy(t,t,
true), dotxy(t,
r,
true) };
212 if ( tvals[0] !=
RT(0.0) )
214 omega = tvals[1]/tvals[0];
220 Saxpy(sol, omega,
r, 0, 0, ncomp, nghost);
221 Saxpy(
r, -omega, t, 0, 0, ncomp, nghost);
227 amrex::Print() << print_ident <<
"MLCGSolver_BiCGStab: Iteration "
228 << std::setw(11) << iter
230 << rnorm/(rnorm0) <<
'\n';
233 if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) {
break; }
244 amrex::Print() << print_ident <<
"MLCGSolver_BiCGStab: Final: Iteration "
245 << std::setw(4) << iter
247 << rnorm/(rnorm0) <<
'\n';
250 if ( ret == 0 && rnorm > eps_rel*rnorm0 && rnorm > eps_abs)
258 if ( ( ret == 0 || ret == 8 ) && (rnorm < rnorm0) )
260 if ( !initial_vec_zeroed ) {
261 LocalAdd(sol, sorig, 0, 0, ncomp, nghost);
263 if (ret == 8) { ret = 9; }
268 if ( !initial_vec_zeroed ) {
269 LocalAdd(sol, sorig, 0, 0, ncomp, nghost);
276template <
typename MF>
282 const int ncomp =
nComp(sol);
284 MF p = Lp.make(amrlev, mglev,
nGrowVect(sol));
287 MF
r = Lp.make(amrlev, mglev, nghost);
288 MF q = Lp.make(amrlev, mglev, nghost);
292 if ( initial_vec_zeroed ) {
295 sorig = Lp.make(amrlev, mglev, nghost);
303 RT rnorm = norm_inf(
r);
304 const RT rnorm0 = rnorm;
308 amrex::Print() << print_ident <<
"MLCGSolver_CG: Initial error (error0) : " << rnorm0 <<
'\n';
315 if ( rnorm0 == 0 || rnorm0 < eps_abs )
318 amrex::Print() << print_ident <<
"MLCGSolver_CG: niter = 0,"
319 <<
", rnorm = " << rnorm
320 <<
", eps_abs = " << eps_abs <<
'\n';
325 for (; iter <= maxiter; ++iter)
340 Xpay(p, beta,
r, 0, 0, ncomp, nghost);
360 <<
" alpha " << alpha <<
'\n';
362 Saxpy(sol, alpha, p, 0, 0, ncomp, nghost);
363 Saxpy(
r, -alpha, q, 0, 0, ncomp, nghost);
368 amrex::Print() << print_ident <<
"MLCGSolver_cg: Iteration"
369 << std::setw(4) << iter
371 << rnorm/(rnorm0) <<
'\n';
374 if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) {
break; }
381 amrex::Print() << print_ident <<
"MLCGSolver_cg: Final Iteration"
382 << std::setw(4) << iter
384 << rnorm/(rnorm0) <<
'\n';
387 if ( ret == 0 && rnorm > eps_rel*rnorm0 && rnorm > eps_abs )
395 if ( ( ret == 0 || ret == 8 ) && (rnorm < rnorm0) )
397 if ( !initial_vec_zeroed ) {
398 LocalAdd(sol, sorig, 0, 0, ncomp, nghost);
400 if (ret == 8) { ret = 9; }
405 if ( !initial_vec_zeroed ) {
406 LocalAdd(sol, sorig, 0, 0, ncomp, nghost);
413template <
typename MF>
419 RT result = Lp.xdoty(amrlev, mglev,
r, z, local);
424template <
typename MF>
428 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:415
int verbose
Definition AMReX_MLCGSolver.H:72
int iter
Definition AMReX_MLCGSolver.H:75
MLLinOpT< MF > & Lp
Definition AMReX_MLCGSolver.H:68
void setSolver(Type _typ) noexcept
Definition AMReX_MLCGSolver.H:28
bool getInitSolnZeroed() const
Definition AMReX_MLCGSolver.H:54
void setVerbose(int _verbose)
Definition AMReX_MLCGSolver.H:39
int getNGhost()
Definition AMReX_MLCGSolver.H:57
const int mglev
Definition AMReX_MLCGSolver.H:71
IntVect nghost
Definition AMReX_MLCGSolver.H:74
MLCGSolverT< MF > & operator=(const MLCGSolverT< MF > &rhs)=delete
int getNumIters() const noexcept
Definition AMReX_MLCGSolver.H:64
int solve_cg(MF &solnL, const MF &rhsL, RT eps_rel, RT eps_abs)
Definition AMReX_MLCGSolver.H:278
RT norm_inf(const MF &res, bool local=false)
Definition AMReX_MLCGSolver.H:426
const int amrlev
Definition AMReX_MLCGSolver.H:70
bool initial_vec_zeroed
Definition AMReX_MLCGSolver.H:76
void setInitSolnZeroed(bool _sol_zeroed)
Definition AMReX_MLCGSolver.H:53
int solve_bicgstab(MF &solnL, const MF &rhsL, RT eps_rel, RT eps_abs)
Definition AMReX_MLCGSolver.H:100
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:73
void setPrintIdentation(std::string s)
Definition AMReX_MLCGSolver.H:45
int solve(MF &solnL, const MF &rhsL, RT eps_rel, RT eps_abs)
Definition AMReX_MLCGSolver.H:89
std::string print_ident
Definition AMReX_MLCGSolver.H:77
Type
Definition AMReX_MLCGSolver.H:18
typename MLLinOpT< MF >::RT RT
Definition AMReX_MLCGSolver.H:16
void setNGhost(int _nghost)
Definition AMReX_MLCGSolver.H:56
MLCGSolverT(MLLinOpT< MF > &_lp, Type _typ=Type::BiCGStab)
Definition AMReX_MLCGSolver.H:81
int getVerbose() const
Definition AMReX_MLCGSolver.H:40
void setMaxIter(int _maxiter)
Definition AMReX_MLCGSolver.H:42
Type solver_type
Definition AMReX_MLCGSolver.H:69
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:1883
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:1847
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:1831
void LocalAdd(MF &dst, MF const &src, int scomp, int dcomp, int ncomp, IntVect const &nghost)
dst += src
Definition AMReX_FabArrayUtility.H:1839
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:1855
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:1808