13 #ifndef _BICGSTABSOLVER_H_ 14 #define _BICGSTABSOLVER_H_ 19 #include "NamespaceHeader.H" 48 virtual void solve(T& a_phi,
const T& a_rhs);
177 CH_TIMER(
"BiCGStabSolver::solve::Initialize",timeInitialize);
178 CH_TIMER(
"BiCGStabSolver::solve::MainLoop",timeMainLoop);
179 CH_TIMER(
"BiCGStabSolver::solve::Cleanup",timeCleanup);
183 T r, r_tilde, e, p, p_tilde, s_tilde, t, v;
185 m_op->create(r, a_rhs);
186 m_op->create(r_tilde, a_rhs);
187 m_op->create(e, a_phi);
188 m_op->create(p, a_rhs);
189 m_op->create(p_tilde, a_phi);
190 m_op->create(s_tilde, a_phi);
191 m_op->create(t, a_rhs);
192 m_op->create(v, a_rhs);
201 m_op->assignLocal(r_tilde, r);
205 m_op->setToZero(p_tilde);
206 m_op->setToZero(s_tilde);
217 Real initial_norm = norm[0];
218 Real initial_rnorm = norm[0];
239 pout() <<
" BiCGStab:: initial Residual norm = " 240 << initial_norm <<
"\n";
253 while ((i<
m_imax && norm[0] >
m_eps*norm[1]) && (norm[1] > 0))
264 pout() <<
" BiCGStab:: norm[0] = " << norm[0] <<
", " 265 <<
"norm[1] = " << norm[1]
267 pout() <<
" BiCGStab:: alpha[0] = " << alpha[0] <<
", " 268 <<
"alpha[1] = " << alpha[1]
270 pout() <<
" BiCGStab:: beta[0] = " << beta[0] <<
", " 271 <<
"beta[1] = " << beta[1]
273 pout() <<
" BiCGStab:: omega[0] = " << omega[0] <<
", " 274 <<
"omega[1] = " << omega[1]
280 rho[1] =
m_op->dotProduct(r_tilde, r);
284 pout() <<
" BiCGStab:: rho[1] = " << rho[1] <<
", " 285 <<
"rho[2] = " << rho[2] <<
", " 286 <<
"rho[3] = " << rho[3]
293 m_op->incr(a_phi, e, 1.0);
297 pout() <<
" BiCGStab:: rho = 0, returning" 298 <<
" -- Residual norm = " 311 m_op->clear(r_tilde);
314 m_op->clear(p_tilde);
315 m_op->clear(s_tilde);
326 m_op->assignLocal(p, r);
331 beta[1] = (rho[1]/rho[2])*(alpha[1]/omega[1]);
332 m_op->scale(p, beta[1]);
333 m_op->incr(p, v, -beta[1]*omega[1]);
334 m_op->incr(p, r, 1.0);
339 pout() <<
" BiCGStab:: beta[1] = " << beta[1]
343 m_op->preCond(p_tilde, p);
345 m_op->applyOp(v, p_tilde,
true);
346 Real m =
m_op->dotProduct(r_tilde, v);
351 pout() <<
" BiCGStab:: rho[1] = " << rho[1] <<
", " 352 <<
"m = " << m <<
", " 353 <<
"alpha[0] = " << alpha[0]
359 m_op->incr(r, v, -alpha[0]);
361 m_op->incr(e, p_tilde, alpha[0]);
371 pout() <<
"norm[0] = " << norm[0]
372 <<
" initial_norm = " << initial_norm
373 <<
" initial_rnorm = " << initial_rnorm << endl;
376 if (norm[0] >
m_eps*initial_norm && norm[0] >
m_reps*initial_rnorm)
378 m_op->preCond(s_tilde, r);
380 m_op->applyOp(t, s_tilde,
true);
381 omega[0] =
m_op->dotProduct(t, r)/
m_op->dotProduct(t, t);
382 m_op->incr(e, s_tilde, omega[0]);
383 m_op->incr(r, t, -omega[0]);
389 pout() <<
" BiCGStab:: iteration = " << i <<
", error norm = " << norm[0] <<
", rate = " << norm[1]/norm[0] <<
"\n";
392 if (norm[0] <=
m_eps*initial_norm || norm[0] <=
m_reps*initial_rnorm)
399 if (omega[0] == 0.0 || norm[0] > (1-
m_hang)*norm[1])
408 m_op->incr(a_phi, e, 1.0);
414 pout() <<
" BiCGStab: max restarts reached" << endl;
415 pout() <<
" init norm = " << initial_norm << endl;
416 pout() <<
" final norm = " << norm[0] << endl;
429 m_op->clear(r_tilde);
432 m_op->clear(p_tilde);
433 m_op->clear(s_tilde);
443 CH_TIME(
"BiCGStabSolver::solve::Restart");
447 rho[1]=0.0; rho[1]=0.0; rho[2]=0.0; rho[3]=0.0;
448 alpha[0]=0; beta[0]=0; omega[0]=0;
449 m_op->assignLocal(r_tilde, r);
457 pout() <<
" BiCGStab:: restart = " << restarts <<
"\n";
470 pout() <<
" BiCGStab:: " << i <<
" iterations, final Residual norm = " 474 m_op->incr(a_phi, e, 1.0);
477 m_op->clear(r_tilde);
480 m_op->clear(p_tilde);
481 m_op->clear(s_tilde);
496 #include "NamespaceFooter.H" std::ostream & pout()
Use this in place of std::cout for program output.
virtual void setConvergenceMetrics(Real a_metric, Real a_tolerance)
Set a convergence metric, along with solver tolerance, if desired.
Definition: BiCGStabSolver.H:489
#define CH_TIMERS(name)
Definition: CH_Timer.H:133
Real m_convergenceMetric
Definition: BiCGStabSolver.H:97
Real m_eps
Definition: BiCGStabSolver.H:84
Definition: LinearSolver.H:28
int m_imax
Definition: BiCGStabSolver.H:71
#define CH_assert(cond)
Definition: CHArray.H:37
LinearOp< T > * m_op
Definition: BiCGStabSolver.H:65
int m_verbosity
Definition: BiCGStabSolver.H:78
Real norm(const BoxLayoutData< FArrayBox > &A, const Interval &interval, const int &p=2)
#define CH_START(tpointer)
Definition: CH_Timer.H:145
Real m_reps
Definition: BiCGStabSolver.H:90
int m_exitStatus
Definition: BiCGStabSolver.H:113
#define CH_TIMER(name, tpointer)
Definition: CH_Timer.H:63
Real m_hang
Definition: BiCGStabSolver.H:103
int m_normType
Definition: BiCGStabSolver.H:132
#define CH_TIME(name)
Definition: CH_Timer.H:82
virtual ~BiCGStabSolver()
Definition: BiCGStabSolver.H:160
double Real
Definition: REAL.H:33
int m_numRestarts
Definition: BiCGStabSolver.H:125
T Abs(const T &a_a)
Definition: Misc.H:53
virtual void solve(T &a_phi, const T &a_rhs)
solve the equation.
Definition: BiCGStabSolver.H:173
Definition: LinearSolver.H:156
#define CH_STOP(tpointer)
Definition: CH_Timer.H:150
virtual void define(LinearOp< T > *a_op, bool a_homogeneous)
Definition: BiCGStabSolver.H:166
BiCGStabSolver()
Definition: BiCGStabSolver.H:143
virtual void setHomogeneous(bool a_homogeneous)
Definition: BiCGStabSolver.H:34
Real m_small
Definition: BiCGStabSolver.H:119
bool m_homogeneous
Definition: BiCGStabSolver.H:59
Definition: BiCGStabSolver.H:26