27 #ifndef EWOMS_BICG_STAB_SOLVER_HH 28 #define EWOMS_BICG_STAB_SOLVER_HH 37 #include <opm/common/ErrorMacros.hpp> 38 #include <opm/common/Exceptions.hpp> 53 template <
class LinearOperator,
class Vector,
class Preconditioner>
57 typedef typename LinearOperator::field_type Scalar;
62 Dune::ScalarProduct<Vector>& scalarProduct)
63 : preconditioner_(preconditioner)
64 , convergenceCriterion_(convergenceCriterion)
65 , scalarProduct_(scalarProduct)
70 maxIterations_ = 1000;
78 { maxIterations_ = value; }
85 {
return maxIterations_; }
98 { verbosity_ = value; }
104 {
return verbosity_; }
124 const Scalar breakdownEps = std::numeric_limits<Scalar>::min() * Scalar(1e10);
132 report_.timer().
start();
146 preconditioner_.pre(x, r);
153 for (
unsigned i = 0; i < x.size(); ++i) {
154 const auto& u = x[i];
156 OPM_THROW(std::logic_error,
157 "The preconditioner is assumed not to modify the initial solution!");
163 report_.setConverged(
true);
164 return report_.converged();
167 if (verbosity_ > 0) {
168 std::cout <<
"-------- BiCGStabSolver --------" << std::endl;
176 const Vector& r0hat = *b_;
195 unsigned n = x.size();
197 for (; report_.iterations() < maxIterations_; report_.increment()) {
199 Scalar rho_i = scalarProduct_.dot(r0hat, r);
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);
214 for (
unsigned i = 0; i < n; ++i) {
228 preconditioner_.apply(y, p);
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)");
239 if (std::abs(alpha) <= breakdownEps)
240 OPM_THROW(Opm::NumericalProblem,
241 "Breakdown of the BiCGStab solver (stagnation detected)");
245 for (
unsigned i = 0; i < n; ++i) {
258 convergenceCriterion_.
update(h, y, s);
260 if (verbosity_ > 0) {
261 convergenceCriterion_.
print(report_.iterations() + 0.5);
262 std::cout <<
"-------- /BiCGStabSolver --------" << std::endl;
266 preconditioner_.post(x);
267 report_.setConverged(
true);
268 return report_.converged();
270 else if (convergenceCriterion_.
failed()) {
271 if (verbosity_ > 0) {
272 convergenceCriterion_.
print(report_.iterations() + 0.5);
273 std::cout <<
"-------- /BiCGStabSolver --------" << std::endl;
276 report_.setConverged(
false);
277 return report_.converged();
281 convergenceCriterion_.
print(report_.iterations() + 0.5);
285 preconditioner_.apply(z, s);
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)");
306 convergenceCriterion_.
update(x, z, r);
308 if (verbosity_ > 0) {
309 convergenceCriterion_.
print(1.0 + report_.iterations());
310 std::cout <<
"-------- /BiCGStabSolver --------" << std::endl;
313 preconditioner_.post(x);
314 report_.setConverged(
true);
315 return report_.converged();
317 else if (convergenceCriterion_.
failed()) {
318 if (verbosity_ > 0) {
319 convergenceCriterion_.
print(1.0 + report_.iterations());
320 std::cout <<
"-------- /BiCGStabSolver --------" << std::endl;
323 report_.setConverged(
false);
324 return report_.converged();
328 convergenceCriterion_.
print(1.0 + report_.iterations());
335 report_.setConverged(
false);
336 return report_.converged();
341 convergenceCriterion_ = &crit;
348 const LinearOperator* A_;
351 Preconditioner& preconditioner_;
352 ConvergenceCriterion& convergenceCriterion_;
353 Dune::ScalarProduct<Vector>& scalarProduct_;
356 unsigned maxIterations_;
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
Definition: baseauxiliarymodule.hh:37
virtual void printInitial(std::ostream &os OPM_UNUSED=std::cout) const
Prints the initial information about the convergence behaviour.
Definition: convergencecriterion.hh:137
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
void setRhs(const Vector *b)
Set the right hand side "b" 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 "x" vector.
Definition: bicgstabsolver.hh:121
unsigned maxIterations() const
Return the maximum number of iterations before we give up without achieving convergence.
Definition: bicgstabsolver.hh:84
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 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.
void setMaxIterations(unsigned value)
Set the maximum number of iterations before we give up without achieving convergence.
Definition: bicgstabsolver.hh:77
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 verbosity() const
Return the verbosity level of the linear solver.
Definition: bicgstabsolver.hh:103
A simple class which makes sure that a timer gets stopped if an exception is thrown.
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 "A" of the linear system.
Definition: bicgstabsolver.hh:109