Chombo + EB  3.0
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: what the algorithm should consider "close to zero"
118  */
120
121  ///
122  /**
123  public member data: number of times the algorithm can restart
124  */
126
127  ///
128  /**
129  public member data: norm to be used when evaluation convergence.
130  0 is max norm, 1 is L(1), 2 is L(2) and so on.
131  */
133
134 };
135
136 // *******************************************************
137 // BiCGStabSolver Implementation
138 // *******************************************************
139
140 // For large elliptic problem, bottom smoother needed 67 iterations.
141 // Bumping imax to 80. (ndk)
142 template <class T>
144  :m_homogeneous(false),
145  m_op(NULL),
146  m_imax(80),
147  m_verbosity(3),
148  m_eps(1.0E-6),
149  m_reps(1.0E-12),
150  m_convergenceMetric(-1.0),
151  m_hang(1E-8),
152  m_exitStatus(-1),
153  m_small(1.0E-30),
154  m_numRestarts(5),
155  m_normType(2)
156 {
157 }
158
159 template <class T>
161 {
162  m_op = NULL;
163 }
164
165 template <class T>
166 void BiCGStabSolver<T>::define(LinearOp<T>* a_operator, bool a_homogeneous)
167 {
168  m_homogeneous = a_homogeneous;
169  m_op = a_operator;
170 }
171
172 template <class T>
173 void BiCGStabSolver<T>::solve(T& a_phi, const T& a_rhs)
174 {
175  CH_TIMERS("BiCGStabSolver::solve");
176
177  CH_TIMER("BiCGStabSolver::solve::Initialize",timeInitialize);
178  CH_TIMER("BiCGStabSolver::solve::MainLoop",timeMainLoop);
179  CH_TIMER("BiCGStabSolver::solve::Cleanup",timeCleanup);
180
181  CH_START(timeInitialize);
182
183  T r, r_tilde, e, p, p_tilde, s_tilde, t, v;
184
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);
193
194  int recount = 0;
195
196  CH_assert(m_op != NULL);
197
198  m_op->residual(r, a_phi, a_rhs, m_homogeneous);
199
200  m_op->assignLocal(r_tilde, r);
201  m_op->setToZero(e);
202  // (DFM 2/1/07) these next two need to be set to zero to prevent
203  // problems in the multilevel case
204  m_op->setToZero(p_tilde);
205  m_op->setToZero(s_tilde);
206
207  int i = 0;
208
209  // rho[0] = r_i , rho[1] = r_(i-1), etc.
210  Real rho[4]=
211  {
212  0,0,0,0
213  };
214  Real norm[2];
215  norm[0] = m_op->norm(r, m_normType);
216  Real initial_norm = norm[0];
217  Real initial_rnorm = norm[0];
218  norm[1] = norm[0];
219
220  Real alpha[2] =
221  {
222  0,0
223  };
224  Real beta[2] =
225  {
226  0,0
227  };
228  Real omega[2] =
229  {
230  0,0
231  };
232
233  bool init = true;
234  int restarts = 0;
235
236  if (m_verbosity >= 3)
237  {
238  pout() << " BiCGStab:: initial Residual norm = "
239  << initial_norm << "\n";
240  }
241
242  // if a convergence metric has been supplied, replace initial residual
243  // with the supplied convergence metric...
244  if (m_convergenceMetric > 0)
245  {
246  initial_norm = m_convergenceMetric;
247  }
248
249  CH_STOP(timeInitialize);
250
251  CH_START(timeMainLoop);
252  while ((i<m_imax && norm[0] > m_eps*norm[1]) && (norm[1] > 0))
253  {
254  i++;
255
256  norm[1] = norm[0];
257  alpha[1]= alpha[0];
258  beta[1] = beta[0];
259  omega[1]= omega[0];
260
261  if (m_verbosity >= 5)
262  {
263  pout() << " BiCGStab:: norm[0] = " << norm[0] << ", "
264  << "norm[1] = " << norm[1]
265  << "\n";
266  pout() << " BiCGStab:: alpha[0] = " << alpha[0] << ", "
267  << "alpha[1] = " << alpha[1]
268  << "\n";
269  pout() << " BiCGStab:: beta[0] = " << beta[0] << ", "
270  << "beta[1] = " << beta[1]
271  << "\n";
272  pout() << " BiCGStab:: omega[0] = " << omega[0] << ", "
273  << "omega[1] = " << omega[1]
274  << "\n";
275  }
276
277  rho[3] = rho[2];
278  rho[2] = rho[1];
279  rho[1] = m_op->dotProduct(r_tilde, r);
280
281  if (m_verbosity >= 5)
282  {
283  pout() << " BiCGStab:: rho[1] = " << rho[1] << ", "
284  << "rho[2] = " << rho[2] << ", "
285  << "rho[3] = " << rho[3]
286  << "\n";
287  }
288
289  if (rho[1] == 0.0)
290  {
291  // we are finished, we will not converge anymore
292  m_op->incr(a_phi, e, 1.0);
293
294  if (m_verbosity >= 3)
295  {
296  pout() << " BiCGStab:: rho = 0, returning"
297  << " -- Residual norm = "
298  << norm[0] << "\n";
299  }
300
301  CH_STOP(timeMainLoop);
302
303  CH_START(timeCleanup);
304
305  m_exitStatus = 2;
306  // need to call clear just in case
307  // this is not ideal -- maybe change the return
308  // to exit
309  m_op->clear(r);
310  m_op->clear(r_tilde);
311  m_op->clear(e);
312  m_op->clear(p);
313  m_op->clear(p_tilde);
314  m_op->clear(s_tilde);
315  m_op->clear(t);
316  m_op->clear(v);
317
318  CH_STOP(timeCleanup);
319
320  return;
321  }
322
323  if (init)
324  {
325  m_op->assignLocal(p, r);
326  init = false;
327  }
328  else
329  {
330  beta[1] = (rho[1]/rho[2])*(alpha[1]/omega[1]);
331  m_op->scale(p, beta[1]);
332  m_op->incr(p, v, -beta[1]*omega[1]);
333  m_op->incr(p, r, 1.0);
334  }
335
336  if (m_verbosity >= 5)
337  {
338  pout() << " BiCGStab:: beta[1] = " << beta[1]
339  << "\n";
340  }
341
342  m_op->preCond(p_tilde, p);
343  m_op->applyOp(v, p_tilde, true);
344  Real m = m_op->dotProduct(r_tilde, v);
345  alpha[0] = rho[1]/m;
346
347  if (m_verbosity >= 5)
348  {
349  pout() << " BiCGStab:: rho[1] = " << rho[1] << ", "
350  << "m = " << m << ", "
351  << "alpha[0] = " << alpha[0]
352  << "\n";
353  }
354
355  if (Abs(m) > m_small*Abs(rho[1]))
356  {
357  m_op->incr(r, v, -alpha[0]);
358  norm[0] = m_op->norm(r, m_normType);
359  m_op->incr(e, p_tilde, alpha[0]);
360  }
361  else
362  {
363  m_op->setToZero(r);
364  norm[0] = 0.0;
365  }
366
367  if (norm[0] > m_eps*initial_norm && norm[0] > m_reps*initial_rnorm)
368  {
369  m_op->preCond(s_tilde, r);
370  m_op->applyOp(t, s_tilde, true);
371  omega[0] = m_op->dotProduct(t, r)/m_op->dotProduct(t, t);
372  m_op->incr(e, s_tilde, omega[0]);
373  m_op->incr(r, t, -omega[0]);
374  norm[0] = m_op->norm(r, m_normType);
375  }
376
377  if (m_verbosity >= 4)
378  {
379  pout() << " BiCGStab:: iteration = " << i << ", error norm = " << norm[0] << ", rate = " << norm[1]/norm[0] << "\n";
380  }
381
382  if (norm[0] <= m_eps*initial_norm || norm[0] <= m_reps*initial_rnorm)
383  {
384  // converged to tolerance
385  m_exitStatus = 1;
386  break;
387  }
388
389  if (omega[0] == 0.0 || norm[0] > (1-m_hang)*norm[1])
390  {
391  if (recount == 0)
392  {
393  recount = 1;
394  }
395  else
396  {
397  recount = 0;
398  m_op->incr(a_phi, e, 1.0);
399
400  if (restarts == m_numRestarts)
401  {
402  if (m_verbosity >= 4)
403  {
404  pout() << " BiCGStab: max restarts reached" << endl;
405  pout() << " init norm = " << initial_norm << endl;
406  pout() << " final norm = " << norm[0] << endl;
407  }
408
409  CH_STOP(timeMainLoop);
410
411  CH_START(timeCleanup);
412
413  m_exitStatus = 3;
414
415  // need to call clear just in case
416  // this is not ideal -- maybe change the return
417  // to exit
418  m_op->clear(r);
419  m_op->clear(r_tilde);
420  m_op->clear(e);
421  m_op->clear(p);
422  m_op->clear(p_tilde);
423  m_op->clear(s_tilde);
424  m_op->clear(t);
425  m_op->clear(v);
426
427  CH_STOP(timeCleanup);
428
429  return;
430  }
431
432  {
433  CH_TIME("BiCGStabSolver::solve::Restart");
434
435  m_op->residual(r, a_phi, a_rhs, m_homogeneous);
436  norm[0] = m_op->norm(r, m_normType);
437  rho[1]=0.0; rho[1]=0.0; rho[2]=0.0; rho[3]=0.0;
438  alpha[0]=0; beta[0]=0; omega[0]=0;
439  m_op->assignLocal(r_tilde, r);
440  m_op->setToZero(e);
441
442  restarts++;
443  }
444
445  if (m_verbosity >= 4)
446  {
447  pout() << " BiCGStab:: restart = " << restarts << "\n";
448  }
449
450  init = true;
451  }
452  }
453  }
454  CH_STOP(timeMainLoop);
455
456  CH_START(timeCleanup);
457
458  if (m_verbosity >= 4)
459  {
460  pout() << " BiCGStab:: " << i << " iterations, final Residual norm = "
461  << norm[0] << "\n";
462  }
463
464  m_op->incr(a_phi, e, 1.0);
465
466  m_op->clear(r);
467  m_op->clear(r_tilde);
468  m_op->clear(e);
469  m_op->clear(p);
470  m_op->clear(p_tilde);
471  m_op->clear(s_tilde);
472  m_op->clear(t);
473  m_op->clear(v);
474
475  CH_STOP(timeCleanup);
476 }
477
478 template <class T>
480  Real a_tolerance)
481 {
482  m_convergenceMetric = a_metric;
483  m_eps = a_tolerance;
484 }
485
486 #include "NamespaceFooter.H"
487 #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:479
#define CH_TIMERS(name)
Definition: CH_Timer.H:70
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:78
Real m_reps
Definition: BiCGStabSolver.H:90
int m_exitStatus
Definition: BiCGStabSolver.H:113
#define CH_TIMER(name, tpointer)
Definition: CH_Timer.H:55
Real m_hang
Definition: BiCGStabSolver.H:103
int m_normType
Definition: BiCGStabSolver.H:132
#define CH_TIME(name)
Definition: CH_Timer.H:59
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:158
#define CH_STOP(tpointer)
Definition: CH_Timer.H:80
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