Chombo + EB  3.2
BiCGStabSolver.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
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"
19 #include "NamespaceHeader.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  m_op->setToZero(r); // added by petermc, 26 Nov 2013, to zero out ghosts
198 
199  m_op->residual(r, a_phi, a_rhs, m_homogeneous);
200 
201  m_op->assignLocal(r_tilde, r);
202  m_op->setToZero(e);
203  // (DFM 2/1/07) these next two need to be set to zero to prevent
204  // problems in the multilevel case
205  m_op->setToZero(p_tilde);
206  m_op->setToZero(s_tilde);
207 
208  int i = 0;
209 
210  // rho[0] = r_i , rho[1] = r_(i-1), etc.
211  Real rho[4]=
212  {
213  0,0,0,0
214  };
215  Real norm[2];
216  norm[0] = m_op->norm(r, m_normType);
217  Real initial_norm = norm[0];
218  Real initial_rnorm = norm[0];
219  norm[1] = norm[0];
220 
221  Real alpha[2] =
222  {
223  0,0
224  };
225  Real beta[2] =
226  {
227  0,0
228  };
229  Real omega[2] =
230  {
231  0,0
232  };
233 
234  bool init = true;
235  int restarts = 0;
236 
237  if (m_verbosity >= 5)
238  {
239  pout() << " BiCGStab:: initial Residual norm = "
240  << initial_norm << "\n";
241  }
242 
243  // if a convergence metric has been supplied, replace initial residual
244  // with the supplied convergence metric...
245  if (m_convergenceMetric > 0)
246  {
247  initial_norm = m_convergenceMetric;
248  }
249 
250  CH_STOP(timeInitialize);
251 
252  CH_START(timeMainLoop);
253  while ((i<m_imax && norm[0] > m_eps*norm[1]) && (norm[1] > 0))
254  {
255  i++;
256 
257  norm[1] = norm[0];
258  alpha[1]= alpha[0];
259  beta[1] = beta[0];
260  omega[1]= omega[0];
261 
262  if (m_verbosity >= 5)
263  {
264  pout() << " BiCGStab:: norm[0] = " << norm[0] << ", "
265  << "norm[1] = " << norm[1]
266  << "\n";
267  pout() << " BiCGStab:: alpha[0] = " << alpha[0] << ", "
268  << "alpha[1] = " << alpha[1]
269  << "\n";
270  pout() << " BiCGStab:: beta[0] = " << beta[0] << ", "
271  << "beta[1] = " << beta[1]
272  << "\n";
273  pout() << " BiCGStab:: omega[0] = " << omega[0] << ", "
274  << "omega[1] = " << omega[1]
275  << "\n";
276  }
277 
278  rho[3] = rho[2];
279  rho[2] = rho[1];
280  rho[1] = m_op->dotProduct(r_tilde, r);
281 
282  if (m_verbosity >= 5)
283  {
284  pout() << " BiCGStab:: rho[1] = " << rho[1] << ", "
285  << "rho[2] = " << rho[2] << ", "
286  << "rho[3] = " << rho[3]
287  << "\n";
288  }
289 
290  if (rho[1] == 0.0)
291  {
292  // we are finished, we will not converge anymore
293  m_op->incr(a_phi, e, 1.0);
294 
295  if (m_verbosity >= 5)
296  {
297  pout() << " BiCGStab:: rho = 0, returning"
298  << " -- Residual norm = "
299  << norm[0] << "\n";
300  }
301 
302  CH_STOP(timeMainLoop);
303 
304  CH_START(timeCleanup);
305 
306  m_exitStatus = 2;
307  // need to call clear just in case
308  // this is not ideal -- maybe change the return
309  // to exit
310  m_op->clear(r);
311  m_op->clear(r_tilde);
312  m_op->clear(e);
313  m_op->clear(p);
314  m_op->clear(p_tilde);
315  m_op->clear(s_tilde);
316  m_op->clear(t);
317  m_op->clear(v);
318 
319  CH_STOP(timeCleanup);
320 
321  return;
322  }
323 
324  if (init)
325  {
326  m_op->assignLocal(p, r);
327  init = false;
328  }
329  else
330  {
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);
335  }
336 
337  if (m_verbosity >= 5)
338  {
339  pout() << " BiCGStab:: beta[1] = " << beta[1]
340  << "\n";
341  }
342 
343  m_op->preCond(p_tilde, p);
344  m_op->setToZero(v); // added by petermc, 27 Nov 2013, to zero out ghosts
345  m_op->applyOp(v, p_tilde, true);
346  Real m = m_op->dotProduct(r_tilde, v);
347  alpha[0] = rho[1]/m;
348 
349  if (m_verbosity >= 5)
350  {
351  pout() << " BiCGStab:: rho[1] = " << rho[1] << ", "
352  << "m = " << m << ", "
353  << "alpha[0] = " << alpha[0]
354  << "\n";
355  }
356 
357  if (Abs(m) > m_small*Abs(rho[1]))
358  {
359  m_op->incr(r, v, -alpha[0]);
360  norm[0] = m_op->norm(r, m_normType);
361  m_op->incr(e, p_tilde, alpha[0]);
362  }
363  else
364  {
365  m_op->setToZero(r);
366  norm[0] = 0.0;
367  }
368 
369  if (m_verbosity >= 4)
370  {
371  pout() << "norm[0] = " << norm[0]
372  << " initial_norm = " << initial_norm
373  << " initial_rnorm = " << initial_rnorm << endl;
374  }
375 
376  if (norm[0] > m_eps*initial_norm && norm[0] > m_reps*initial_rnorm)
377  {
378  m_op->preCond(s_tilde, r);
379  m_op->setToZero(t); // added by petermc, 27 Nov 2013, to zero out ghosts
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]);
384  norm[0] = m_op->norm(r, m_normType);
385  }
386 
387  if (m_verbosity >= 4)
388  {
389  pout() << " BiCGStab:: iteration = " << i << ", error norm = " << norm[0] << ", rate = " << norm[1]/norm[0] << "\n";
390  }
391 
392  if (norm[0] <= m_eps*initial_norm || norm[0] <= m_reps*initial_rnorm)
393  {
394  // converged to tolerance
395  m_exitStatus = 1;
396  break;
397  }
398 
399  if (omega[0] == 0.0 || norm[0] > (1-m_hang)*norm[1])
400  {
401  if (recount == 0)
402  {
403  recount = 1;
404  }
405  else
406  {
407  recount = 0;
408  m_op->incr(a_phi, e, 1.0);
409 
410  if (restarts == m_numRestarts)
411  {
412  if (m_verbosity >= 4)
413  {
414  pout() << " BiCGStab: max restarts reached" << endl;
415  pout() << " init norm = " << initial_norm << endl;
416  pout() << " final norm = " << norm[0] << endl;
417  }
418 
419  CH_STOP(timeMainLoop);
420 
421  CH_START(timeCleanup);
422 
423  m_exitStatus = 3;
424 
425  // need to call clear just in case
426  // this is not ideal -- maybe change the return
427  // to exit
428  m_op->clear(r);
429  m_op->clear(r_tilde);
430  m_op->clear(e);
431  m_op->clear(p);
432  m_op->clear(p_tilde);
433  m_op->clear(s_tilde);
434  m_op->clear(t);
435  m_op->clear(v);
436 
437  CH_STOP(timeCleanup);
438 
439  return;
440  }
441 
442  {
443  CH_TIME("BiCGStabSolver::solve::Restart");
444 
445  m_op->residual(r, a_phi, a_rhs, m_homogeneous);
446  norm[0] = m_op->norm(r, m_normType);
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);
450  m_op->setToZero(e);
451 
452  restarts++;
453  }
454 
455  if (m_verbosity >= 4)
456  {
457  pout() << " BiCGStab:: restart = " << restarts << "\n";
458  }
459 
460  init = true;
461  }
462  }
463  }
464  CH_STOP(timeMainLoop);
465 
466  CH_START(timeCleanup);
467 
468  if (m_verbosity >= 4)
469  {
470  pout() << " BiCGStab:: " << i << " iterations, final Residual norm = "
471  << norm[0] << "\n";
472  }
473 
474  m_op->incr(a_phi, e, 1.0);
475 
476  m_op->clear(r);
477  m_op->clear(r_tilde);
478  m_op->clear(e);
479  m_op->clear(p);
480  m_op->clear(p_tilde);
481  m_op->clear(s_tilde);
482  m_op->clear(t);
483  m_op->clear(v);
484 
485  CH_STOP(timeCleanup);
486 }
487 
488 template <class T>
490  Real a_tolerance)
491 {
492  m_convergenceMetric = a_metric;
493  m_eps = a_tolerance;
494 }
495 
496 #include "NamespaceFooter.H"
497 #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: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