All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
bicgstabsolver.hh
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef EWOMS_BICG_STAB_SOLVER_HH
28 #define EWOMS_BICG_STAB_SOLVER_HH
29 
30 #include "convergencecriterion.hh"
32 #include "linearsolverreport.hh"
33 
34 #include <ewoms/common/timer.hh>
36 
37 #include <opm/common/ErrorMacros.hpp>
38 #include <opm/common/Exceptions.hpp>
39 
40 #include <memory>
41 
42 namespace Ewoms {
43 namespace Linear {
53 template <class LinearOperator, class Vector, class Preconditioner>
55 {
57  typedef typename LinearOperator::field_type Scalar;
58 
59 public:
60  BiCGStabSolver(Preconditioner& preconditioner,
61  ConvergenceCriterion& convergenceCriterion,
62  Dune::ScalarProduct<Vector>& scalarProduct)
63  : preconditioner_(preconditioner)
64  , convergenceCriterion_(convergenceCriterion)
65  , scalarProduct_(scalarProduct)
66  {
67  A_ = nullptr;
68  b_ = nullptr;
69 
70  maxIterations_ = 1000;
71  }
72 
77  void setMaxIterations(unsigned value)
78  { maxIterations_ = value; }
79 
84  unsigned maxIterations() const
85  { return maxIterations_; }
86 
97  void setVerbosity(unsigned value)
98  { verbosity_ = value; }
99 
103  unsigned verbosity() const
104  { return verbosity_; }
105 
109  void setLinearOperator(const LinearOperator* A)
110  { A_ = A; }
111 
115  void setRhs(const Vector* b)
116  { b_ = b; }
117 
121  bool apply(Vector& x)
122  {
123  // epsilon used for detecting breakdowns
124  const Scalar breakdownEps = std::numeric_limits<Scalar>::min() * Scalar(1e10);
125 
126  // start the stop watch for the solution proceedure, but make sure that it is
127  // turned off regardless of how we leave the stadium. (i.e., that the timer gets
128  // stopped in case exceptions are thrown as well as if the method returns
129  // regularly.)
130  report_.reset();
131  Ewoms::TimerGuard reportTimerGuard(report_.timer());
132  report_.timer().start();
133 
134  // preconditioned stabilized biconjugate gradient method
135  //
136  // See https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method,
137  // (article date: December 19, 2016)
138 
139  // set the initial solution to the zero vector
140  x = 0.0;
141 
142  // prepare the preconditioner. to allow some optimizations, we assume that the
143  // preconditioner does not change the initial solution x if the initial solution
144  // is a zero vector.
145  Vector r = *b_;
146  preconditioner_.pre(x, r);
147 
148 #ifndef NDEBUG
149  // ensure that the preconditioner does not change the initial solution. since
150  // this is a debugging check, we don't care if it does not work properly in
151  // parallel. (because this goes wrong, it should be considered to be a bug in the
152  // code anyway.)
153  for (unsigned i = 0; i < x.size(); ++i) {
154  const auto& u = x[i];
155  if (u*u != 0.0)
156  OPM_THROW(std::logic_error,
157  "The preconditioner is assumed not to modify the initial solution!");
158  }
159 #endif // NDEBUG
160 
161  convergenceCriterion_.setInitial(x, r);
162  if (convergenceCriterion_.converged()) {
163  report_.setConverged(true);
164  return report_.converged();
165  }
166 
167  if (verbosity_ > 0) {
168  std::cout << "-------- BiCGStabSolver --------" << std::endl;
169  convergenceCriterion_.printInitial();
170  }
171 
172  // r0 = b - Ax (i.e., r -= A*x_0 = b, because x_0 == 0)
173  //A_->applyscaleadd(/*alpha=*/-1.0, x, r);
174 
175  // r0hat = r0
176  const Vector& r0hat = *b_;
177 
178  // rho0 = alpha = omega0 = 1
179  Scalar rho = 1.0;
180  Scalar alpha = 1.0;
181  Scalar omega = 1.0;
182 
183  // v_0 = p_0 = 0;
184  Vector v(r);
185  v = 0.0;
186  Vector p(v);
187 
188  // create all the temporary vectors which we need. Be aware that some of them
189  // actually point to the same object because they are not needed at the same time!
190  Vector y(x);
191  Vector& h(x);
192  Vector& s(r);
193  Vector z(x);
194  Vector& t(y);
195  unsigned n = x.size();
196 
197  for (; report_.iterations() < maxIterations_; report_.increment()) {
198  // rho_i = (r0hat,r_(i-1))
199  Scalar rho_i = scalarProduct_.dot(r0hat, r);
200 
201  // beta = (rho_i/rho_(i-1))*(alpha/omega_(i-1))
202  if (std::abs(rho) <= breakdownEps || std::abs(omega) <= breakdownEps)
203  OPM_THROW(Opm::NumericalProblem,
204  "Breakdown of the BiCGStab solver (division by zero)");
205  Scalar beta = (rho_i/rho)*(alpha/omega);
206 
207  // make rho correspond to the current iteration (i.e., forget rho_(i-1))
208  rho = rho_i;
209 
210  // this loop conflates the following operations:
211  //
212  // p_i = r_(i-1) + beta*(p_(i-1) - omega_(i-1)*v_(i-1))
213  // y = p
214  for (unsigned i = 0; i < n; ++i) {
215  // p_i = r_(i-1) + beta*(p_(i-1) - omega_(i-1)*v_(i-1))
216  auto tmp = v[i];
217  tmp *= omega;
218  tmp -= p[i];
219  tmp *= -beta;
220  p[i] = r[i];
221  p[i] += tmp;
222 
223  // y = p; not required because the precontioner overwrites y anyway...
224  // y[i] = p[i];
225  }
226 
227  // y = K^-1 * p_i
228  preconditioner_.apply(y, p);
229 
230  // v_i = A*y
231  A_->apply(y, v);
232 
233  // alpha = rho_i/(r0hat,v_i)
234  Scalar denom = scalarProduct_.dot(r0hat, v);
235  if (std::abs(denom) <= breakdownEps)
236  OPM_THROW(Opm::NumericalProblem,
237  "Breakdown of the BiCGStab solver (division by zero)");
238  alpha = rho_i/denom;
239  if (std::abs(alpha) <= breakdownEps)
240  OPM_THROW(Opm::NumericalProblem,
241  "Breakdown of the BiCGStab solver (stagnation detected)");
242 
243  // h = x_(i-1) + alpha*y
244  // s = r_(i-1) - alpha*v_i
245  for (unsigned i = 0; i < n; ++i) {
246  auto tmp = y[i];
247  tmp *= alpha;
248  tmp += x[i];
249  h[i] = tmp;
250 
251  //s[i] = r[i]; // not necessary because r and s are the same object
252  tmp = v[i];
253  tmp *= alpha;
254  s[i] -= tmp;
255  }
256 
257  // do convergence check and print terminal output
258  convergenceCriterion_.update(/*curSol=*/h, /*delta=*/y, s);
259  if (convergenceCriterion_.converged()) {
260  if (verbosity_ > 0) {
261  convergenceCriterion_.print(report_.iterations() + 0.5);
262  std::cout << "-------- /BiCGStabSolver --------" << std::endl;
263  }
264 
265  // x = h; // not necessary because x and h are the same object
266  preconditioner_.post(x);
267  report_.setConverged(true);
268  return report_.converged();
269  }
270  else if (convergenceCriterion_.failed()) {
271  if (verbosity_ > 0) {
272  convergenceCriterion_.print(report_.iterations() + 0.5);
273  std::cout << "-------- /BiCGStabSolver --------" << std::endl;
274  }
275 
276  report_.setConverged(false);
277  return report_.converged();
278  }
279 
280  if (verbosity_ > 1)
281  convergenceCriterion_.print(report_.iterations() + 0.5);
282 
283  // z = K^-1*s
284  z = s;
285  preconditioner_.apply(z, s);
286 
287  // t = Az
288  t = z;
289  A_->apply(z, t);
290 
291  // omega_i = (t*s)/(t*t)
292  denom = scalarProduct_.dot(t, t);
293  if (std::abs(denom) <= breakdownEps)
294  OPM_THROW(Opm::NumericalProblem,
295  "Breakdown of the BiCGStab solver (division by zero)");
296  omega = scalarProduct_.dot(t, s)/denom;
297  if (std::abs(omega) <= breakdownEps)
298  OPM_THROW(Opm::NumericalProblem,
299  "Breakdown of the BiCGStab solver (stagnation detected)");
300 
301  // x_i = h + omega_i*z
302  // x = h; // not necessary because x and h are the same object
303  x.axpy(/*a=*/omega, /*y=*/z);
304 
305  // do convergence check and print terminal output
306  convergenceCriterion_.update(/*curSol=*/x, /*delta=*/z, r);
307  if (convergenceCriterion_.converged()) {
308  if (verbosity_ > 0) {
309  convergenceCriterion_.print(1.0 + report_.iterations());
310  std::cout << "-------- /BiCGStabSolver --------" << std::endl;
311  }
312 
313  preconditioner_.post(x);
314  report_.setConverged(true);
315  return report_.converged();
316  }
317  else if (convergenceCriterion_.failed()) {
318  if (verbosity_ > 0) {
319  convergenceCriterion_.print(1.0 + report_.iterations());
320  std::cout << "-------- /BiCGStabSolver --------" << std::endl;
321  }
322 
323  report_.setConverged(false);
324  return report_.converged();
325  }
326 
327  if (verbosity_ > 1)
328  convergenceCriterion_.print(1.0 + report_.iterations());
329 
330  // r_i = s - omega*t
331  // r = s; // not necessary because r and s are the same object
332  r.axpy(/*a=*/-omega, /*y=*/t);
333  }
334 
335  report_.setConverged(false);
336  return report_.converged();
337  }
338 
339  void setConvergenceCriterion(ConvergenceCriterion& crit)
340  {
341  convergenceCriterion_ = &crit;
342  }
343 
344  const Ewoms::Linear::SolverReport& report() const
345  { return report_; }
346 
347 private:
348  const LinearOperator* A_;
349  const Vector* b_;
350 
351  Preconditioner& preconditioner_;
352  ConvergenceCriterion& convergenceCriterion_;
353  Dune::ScalarProduct<Vector>& scalarProduct_;
355 
356  unsigned maxIterations_;
357  unsigned verbosity_;
358 };
359 
360 } // namespace Linear
361 } // namespace Ewoms
362 
363 #endif
Base class for all convergence criteria which only defines an virtual API.
Definition: convergencecriterion.hh:57
A simple class which makes sure that a timer gets stopped if an exception is thrown.
Definition: timerguard.hh:40
void start()
Start counting the time resources used by the simulation.
Definition: timer.hh:62
virtual void printInitial(std::ostream &os OPM_UNUSED=std::cout) const
Prints the initial information about the convergence behaviour.
Definition: convergencecriterion.hh:137
void setRhs(const Vector *b)
Set the right hand side &quot;b&quot; of the linear system.
Definition: bicgstabsolver.hh:115
virtual void setInitial(const Vector &curSol, const Vector &curResid)=0
Set the initial solution of the linear system of equations.
Collects summary information about the execution of the linear solver.
Definition: linearsolverreport.hh:41
bool apply(Vector &x)
Run the stabilized BiCG solver and store the result into the &quot;x&quot; vector.
Definition: bicgstabsolver.hh:121
void setVerbosity(unsigned value)
Set the verbosity level of the linear solver.
Definition: bicgstabsolver.hh:97
Collects summary information about the execution of the linear solver.
virtual bool failed() const
Returns true if the convergence criterion cannot be met anymore because the solver has broken down...
Definition: convergencecriterion.hh:117
unsigned maxIterations() const
Return the maximum number of iterations before we give up without achieving convergence.
Definition: bicgstabsolver.hh:84
virtual void update(const Vector &curSol, const Vector &changeIndicator, const Vector &curResid)=0
Update the internal members of the convergence criterion with the current solution.
Provides an encapsulation to measure the system time.
unsigned verbosity() const
Return the verbosity level of the linear solver.
Definition: bicgstabsolver.hh:103
void setMaxIterations(unsigned value)
Set the maximum number of iterations before we give up without achieving convergence.
Definition: bicgstabsolver.hh:77
A simple class which makes sure that a timer gets stopped if an exception is thrown.
virtual void print(Scalar iter OPM_UNUSED, std::ostream &os OPM_UNUSED=std::cout) const
Prints the information about the convergence behaviour for the current iteration. ...
Definition: convergencecriterion.hh:148
Implements a preconditioned stabilized BiCG linear solver.
Definition: bicgstabsolver.hh:54
virtual bool converged() const =0
Returns true if and only if the convergence criterion is met.
void setLinearOperator(const LinearOperator *A)
Set the matrix &quot;A&quot; of the linear system.
Definition: bicgstabsolver.hh:109