Chombo + EB + MF  3.2
BiCGStabSolver.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
8  */
9 #endif
10
11 // BVS, June 18, 2003
12
13 #ifndef _BICGSTABSOLVER_H_
14 #define _BICGSTABSOLVER_H_
15
16 #include "LinearSolver.H"
17 #include "parstream.H"
18 #include "CH_Timer.H"
20
21 ///
22 /**
23  Elliptic solver using the BiCGStab algorithm.
24  */
25 template <class T>
26 class BiCGStabSolver : public LinearSolver<T>
27 {
28 public:
29
31
32  virtual ~BiCGStabSolver();
33
34  virtual void setHomogeneous(bool a_homogeneous)
35  {
36  m_homogeneous = a_homogeneous;
37  }
38
39  ///
40  /**
41  define the solver. a_op is the linear operator.
42  a_homogeneous is whether the solver uses homogeneous boundary
43  conditions.
44  */
45  virtual void define(LinearOp<T>* a_op, bool a_homogeneous);
46
47  ///solve the equation.
48  virtual void solve(T& a_phi, const T& a_rhs);
49
50  ///
51  virtual void setConvergenceMetrics(Real a_metric,
52  Real a_tolerance);
53
54  ///
55  /**
56  public member data: whether the solver is restricted to
57  homogeneous boundary conditions
58  */
60
61  ///
62  /**
63  public member data: operator to solve.
64  */
66
67  ///
68  /**
69  public member data: maximum number of iterations
70  */
71  int m_imax;
72
73  ///
74  /**
75  public member data: how much screen out put the user wants.
76  set = 0 for no output.
77  */
79
80  ///
81  /**
82  public member data: solver tolerance
83  */
85
86  ///
87  /**
88  public member data: relative solver tolerance
89  */
91
92  ///
93  /**
94  public member data: solver convergence metric -- if negative, use
95  initial residual; if positive, then use m_convergenceMetric
96  */
98
99  ///
100  /**
101  public member data: minium norm of solution should change per iterations
102  */
104
105  ///
106  /**
107  public member data:
108  set = -1 if solver exited for an unknown reason
109  set = 1 if solver converged to tolerance
110  set = 2 if rho = 0
111  set = 3 if max number of restarts was reached
112  */
114
115  ///
116  /**
117  public member data: iteration number. This is a member so that it can be
118  queried after a solve, for example.
119  */
120  int m_iter;
121
122  ///
123  /**
124  public member data: current residual norm. This is a member so that
125  it can be queried after a solve, for example
126  */
128
129  ///
130  /**
131  public member data: initial residual norm. This is a member so that
132  it can be queried after a solve, for example
133  */
135
136  ///
137  /**
138  public member data: what the algorithm should consider "close to zero"
139  */
141
142  ///
143  /**
144  public member data: number of times the algorithm can restart
145  */
147
148  ///
149  /**
150  public member data: norm to be used when evaluation convergence.
151  0 is max norm, 1 is L(1), 2 is L(2) and so on.
152  */
154
155 };
156
157
158
159 // *******************************************************
160 // BiCGStabSolver Implementation
161 // *******************************************************
162
163 // For large elliptic problem, bottom smoother needed 67 iterations.
164 // Bumping imax to 80. (ndk)
165 template <class T>
167  :m_homogeneous(false),
168  m_op(NULL),
169  m_imax(80),
170  m_verbosity(3),
171  m_eps(1.0E-6),
172  m_reps(1.0E-12),
173  m_convergenceMetric(-1.0),
174  m_hang(1E-8),
175  m_exitStatus(-1),
176  m_iter(-1),
177  m_residual(-1e300),
178  m_initial_residual(-1e300),
179  m_small(1.0E-30),
180  m_numRestarts(5),
181  m_normType(2)
182 {
183 }
184
185 template <class T>
187 {
188  m_op = NULL;
189 }
190
191 template <class T>
192 void BiCGStabSolver<T>::define(LinearOp<T>* a_operator, bool a_homogeneous)
193 {
194  m_homogeneous = a_homogeneous;
195  m_op = a_operator;
196 }
197
198 template <class T>
199 void BiCGStabSolver<T>::solve(T& a_phi, const T& a_rhs)
200 {
201  CH_TIMERS("BiCGStabSolver::solve");
202
203  CH_TIMER("BiCGStabSolver::solve::Initialize",timeInitialize);
204  CH_TIMER("BiCGStabSolver::solve::MainLoop",timeMainLoop);
205  CH_TIMER("BiCGStabSolver::solve::Cleanup",timeCleanup);
206
207  CH_START(timeInitialize);
208
209  T r, r_tilde, e, p, p_tilde, s_tilde, t, v;
210
211  m_op->create(r, a_rhs);
212  m_op->create(r_tilde, a_rhs);
213  m_op->create(e, a_phi);
214  m_op->create(p, a_rhs);
215  m_op->create(p_tilde, a_phi);
216  m_op->create(s_tilde, a_phi);
217  m_op->create(t, a_rhs);
218  m_op->create(v, a_rhs);
219
220  int recount = 0;
221
222  CH_assert(m_op != NULL);
223  m_op->setToZero(r); // added by petermc, 26 Nov 2013, to zero out ghosts
224
225  m_op->residual(r, a_phi, a_rhs, m_homogeneous);
226
227  m_op->assignLocal(r_tilde, r);
228  m_op->setToZero(e);
229  // (DFM 2/1/07) these next two need to be set to zero to prevent
230  // problems in the multilevel case
231  m_op->setToZero(p_tilde);
232  m_op->setToZero(s_tilde);
233
234  m_iter = 0;
235
236  // rho[0] = r_i , rho[1] = r_(i-1), etc.
237  Real rho[4]=
238  {
239  0,0,0,0
240  };
241  Real norm[2];
242  norm[0] = m_op->norm(r, m_normType);
243  Real initial_norm = norm[0];
244  m_initial_residual = norm[0];
245  m_residual = norm[0];
246  norm[1] = norm[0];
247
248  Real alpha[2] =
249  {
250  0,0
251  };
252  Real beta[2] =
253  {
254  0,0
255  };
256  Real omega[2] =
257  {
258  0,0
259  };
260
261  bool init = true;
262  int restarts = 0;
263
264  if (m_verbosity >= 5)
265  {
266  pout() << " BiCGStab:: initial Residual norm = "
267  << initial_norm << "\n";
268  }
269
270  // if a convergence metric has been supplied, replace initial residual
271  // with the supplied convergence metric...
272  if (m_convergenceMetric > 0)
273  {
274  initial_norm = m_convergenceMetric;
275  }
276
277  CH_STOP(timeInitialize);
278
279  CH_START(timeMainLoop);
280  while ((m_iter<m_imax && norm[0] > m_eps*norm[1]) && (norm[1] > 0))
281  {
282  m_iter++;
283
284  norm[1] = norm[0];
285  alpha[1]= alpha[0];
286  beta[1] = beta[0];
287  omega[1]= omega[0];
288
289  if (m_verbosity >= 5)
290  {
291  pout() << " BiCGStab:: norm[0] = " << norm[0] << ", "
292  << "norm[1] = " << norm[1]
293  << "\n";
294  pout() << " BiCGStab:: alpha[0] = " << alpha[0] << ", "
295  << "alpha[1] = " << alpha[1]
296  << "\n";
297  pout() << " BiCGStab:: beta[0] = " << beta[0] << ", "
298  << "beta[1] = " << beta[1]
299  << "\n";
300  pout() << " BiCGStab:: omega[0] = " << omega[0] << ", "
301  << "omega[1] = " << omega[1]
302  << "\n";
303  }
304
305  rho[3] = rho[2];
306  rho[2] = rho[1];
307  rho[1] = m_op->dotProduct(r_tilde, r);
308
309  if (m_verbosity >= 5)
310  {
311  pout() << " BiCGStab:: rho[1] = " << rho[1] << ", "
312  << "rho[2] = " << rho[2] << ", "
313  << "rho[3] = " << rho[3]
314  << "\n";
315  }
316
317  if (rho[1] == 0.0)
318  {
319  // we are finished, we will not converge anymore
320  m_op->incr(a_phi, e, 1.0);
321
322  if (m_verbosity >= 5)
323  {
324  pout() << " BiCGStab:: rho = 0, returning"
325  << " -- Residual norm = "
326  << norm[0] << "\n";
327  }
328
329  m_residual = norm[0];
330
331  CH_STOP(timeMainLoop);
332
333  CH_START(timeCleanup);
334
335  m_exitStatus = 2;
336  // need to call clear just in case
337  // this is not ideal -- maybe change the return
338  // to exit
339  m_op->clear(r);
340  m_op->clear(r_tilde);
341  m_op->clear(e);
342  m_op->clear(p);
343  m_op->clear(p_tilde);
344  m_op->clear(s_tilde);
345  m_op->clear(t);
346  m_op->clear(v);
347
348  CH_STOP(timeCleanup);
349
350  return;
351  }
352
353  if (init)
354  {
355  m_op->assignLocal(p, r);
356  init = false;
357  }
358  else
359  {
360  beta[1] = (rho[1]/rho[2])*(alpha[1]/omega[1]);
361  m_op->scale(p, beta[1]);
362  m_op->incr(p, v, -beta[1]*omega[1]);
363  m_op->incr(p, r, 1.0);
364  }
365
366  if (m_verbosity >= 5)
367  {
368  pout() << " BiCGStab:: beta[1] = " << beta[1]
369  << "\n";
370  }
371
372  m_op->preCond(p_tilde, p);
373  m_op->setToZero(v); // added by petermc, 27 Nov 2013, to zero out ghosts
374  m_op->applyOp(v, p_tilde, true);
375  Real m = m_op->dotProduct(r_tilde, v);
376  alpha[0] = rho[1]/m;
377
378  if (m_verbosity >= 5)
379  {
380  pout() << " BiCGStab:: rho[1] = " << rho[1] << ", "
381  << "m = " << m << ", "
382  << "alpha[0] = " << alpha[0]
383  << "\n";
384  }
385
386  if (Abs(m) > m_small*Abs(rho[1]))
387  {
388  m_op->incr(r, v, -alpha[0]);
389  norm[0] = m_op->norm(r, m_normType);
390  m_op->incr(e, p_tilde, alpha[0]);
391  }
392  else
393  {
394  m_op->setToZero(r);
395  norm[0] = 0.0;
396  }
397
398  if (m_verbosity >= 4)
399  {
400  pout() << "norm[0] = " << norm[0]
401  << " initial_norm = " << initial_norm
402  << " initial_rnorm = " << m_initial_residual << endl;
403  }
404  m_residual = norm[0];
405
406  if (norm[0] > m_eps*initial_norm && norm[0] > m_reps*m_initial_residual)
407  {
408  m_op->preCond(s_tilde, r);
409  m_op->setToZero(t); // added by petermc, 27 Nov 2013, to zero out ghosts
410  m_op->applyOp(t, s_tilde, true);
411  omega[0] = m_op->dotProduct(t, r)/m_op->dotProduct(t, t);
412  m_op->incr(e, s_tilde, omega[0]);
413  m_op->incr(r, t, -omega[0]);
414  norm[0] = m_op->norm(r, m_normType);
415  }
416
417  if (m_verbosity >= 4)
418  {
419  pout() << " BiCGStab:: iteration = " << m_iter << ", error norm = " << norm[0] << ", rate = " << norm[1]/norm[0] << "\n";
420  }
421  m_residual = norm[0];
422
423  if (norm[0] <= m_eps*initial_norm || norm[0] <= m_reps*m_initial_residual)
424  {
425  // converged to tolerance
426  m_exitStatus = 1;
427  break;
428  }
429
430  if (omega[0] == 0.0 || norm[0] > (1-m_hang)*norm[1])
431  {
432  if (recount == 0)
433  {
434  recount = 1;
435  }
436  else
437  {
438  recount = 0;
439  m_op->incr(a_phi, e, 1.0);
440
441  if (restarts == m_numRestarts)
442  {
443  if (m_verbosity >= 4)
444  {
445  pout() << " BiCGStab: max restarts reached" << endl;
446  pout() << " init norm = " << initial_norm << endl;
447  pout() << " final norm = " << norm[0] << endl;
448  }
449
450  CH_STOP(timeMainLoop);
451
452  CH_START(timeCleanup);
453
454  m_exitStatus = 3;
455
456  // need to call clear just in case
457  // this is not ideal -- maybe change the return
458  // to exit
459  m_op->clear(r);
460  m_op->clear(r_tilde);
461  m_op->clear(e);
462  m_op->clear(p);
463  m_op->clear(p_tilde);
464  m_op->clear(s_tilde);
465  m_op->clear(t);
466  m_op->clear(v);
467
468  CH_STOP(timeCleanup);
469
470  return;
471  }
472
473  {
474  CH_TIME("BiCGStabSolver::solve::Restart");
475
476  m_op->residual(r, a_phi, a_rhs, m_homogeneous);
477  norm[0] = m_op->norm(r, m_normType);
478  rho[1]=0.0; rho[1]=0.0; rho[2]=0.0; rho[3]=0.0;
479  alpha[0]=0; beta[0]=0; omega[0]=0;
480  m_op->assignLocal(r_tilde, r);
481  m_op->setToZero(e);
482
483  restarts++;
484  }
485
486  if (m_verbosity >= 4)
487  {
488  pout() << " BiCGStab:: restart = " << restarts << "\n";
489  }
490
491  init = true;
492  }
493  }
494  }
495  CH_STOP(timeMainLoop);
496
497  CH_START(timeCleanup);
498
499  if (m_verbosity >= 4)
500  {
501  pout() << " BiCGStab:: " << m_iter << " iterations, final Residual norm = "
502  << norm[0] << "\n";
503  }
504
505  m_residual = norm[0];
506
507  m_op->incr(a_phi, e, 1.0);
508
509  m_op->clear(r);
510  m_op->clear(r_tilde);
511  m_op->clear(e);
512  m_op->clear(p);
513  m_op->clear(p_tilde);
514  m_op->clear(s_tilde);
515  m_op->clear(t);
516  m_op->clear(v);
517
518  CH_STOP(timeCleanup);
519 }
520
521 template <class T>
523  Real a_tolerance)
524 {
525  m_convergenceMetric = a_metric;
526  m_eps = a_tolerance;
527 }
528
529 #include "NamespaceFooter.H"
530 #endif /*_BICGSTABSOLVER_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:522
#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 m_residual
Definition: BiCGStabSolver.H:127
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:153
#define CH_TIME(name)
Definition: CH_Timer.H:82
virtual ~BiCGStabSolver()
Definition: BiCGStabSolver.H:186
double Real
Definition: REAL.H:33
int m_numRestarts
Definition: BiCGStabSolver.H:146
int m_iter
Definition: BiCGStabSolver.H:120
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:199
Definition: LinearSolver.H:156
#define CH_STOP(tpointer)
Definition: CH_Timer.H:150
Real m_initial_residual
Definition: BiCGStabSolver.H:134
virtual void define(LinearOp< T > *a_op, bool a_homogeneous)
Definition: BiCGStabSolver.H:192
BiCGStabSolver()
Definition: BiCGStabSolver.H:166
virtual void setHomogeneous(bool a_homogeneous)
Definition: BiCGStabSolver.H:34
Real m_small
Definition: BiCGStabSolver.H:140
bool m_homogeneous
Definition: BiCGStabSolver.H:59
Definition: BiCGStabSolver.H:26