All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
ISTLSolver.hpp
1 /*
2  Copyright 2016 IRIS AS
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 3 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 
20 #ifndef OPM_ISTLSOLVER_HEADER_INCLUDED
21 #define OPM_ISTLSOLVER_HEADER_INCLUDED
22 
23 #include <opm/autodiff/AdditionalObjectDeleter.hpp>
24 #include <opm/autodiff/CPRPreconditioner.hpp>
25 #include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
26 #include <opm/autodiff/NewtonIterationUtilities.hpp>
27 #include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
28 #include <opm/autodiff/ParallelOverlappingILU0.hpp>
29 #include <opm/autodiff/AutoDiffHelpers.hpp>
30 
31 #include <opm/common/Exceptions.hpp>
32 #include <opm/core/linalg/ParallelIstlInformation.hpp>
33 #include <opm/common/utility/platform_dependent/disable_warnings.h>
34 
35 #include <dune/istl/scalarproducts.hh>
36 #include <dune/istl/operators.hh>
37 #include <dune/istl/preconditioners.hh>
38 #include <dune/istl/solvers.hh>
39 #include <dune/istl/owneroverlapcopy.hh>
40 #include <dune/istl/paamg/amg.hh>
41 
42 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
43 
44 namespace Dune
45 {
46 namespace FMatrixHelp {
48 template <typename K>
49 static inline K invertMatrix (const FieldMatrix<K,4,4> &matrix, FieldMatrix<K,4,4> &inverse)
50 {
51  inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] -
52  matrix[1][1] * matrix[2][3] * matrix[3][2] -
53  matrix[2][1] * matrix[1][2] * matrix[3][3] +
54  matrix[2][1] * matrix[1][3] * matrix[3][2] +
55  matrix[3][1] * matrix[1][2] * matrix[2][3] -
56  matrix[3][1] * matrix[1][3] * matrix[2][2];
57 
58  inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] +
59  matrix[1][0] * matrix[2][3] * matrix[3][2] +
60  matrix[2][0] * matrix[1][2] * matrix[3][3] -
61  matrix[2][0] * matrix[1][3] * matrix[3][2] -
62  matrix[3][0] * matrix[1][2] * matrix[2][3] +
63  matrix[3][0] * matrix[1][3] * matrix[2][2];
64 
65  inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] -
66  matrix[1][0] * matrix[2][3] * matrix[3][1] -
67  matrix[2][0] * matrix[1][1] * matrix[3][3] +
68  matrix[2][0] * matrix[1][3] * matrix[3][1] +
69  matrix[3][0] * matrix[1][1] * matrix[2][3] -
70  matrix[3][0] * matrix[1][3] * matrix[2][1];
71 
72  inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] +
73  matrix[1][0] * matrix[2][2] * matrix[3][1] +
74  matrix[2][0] * matrix[1][1] * matrix[3][2] -
75  matrix[2][0] * matrix[1][2] * matrix[3][1] -
76  matrix[3][0] * matrix[1][1] * matrix[2][2] +
77  matrix[3][0] * matrix[1][2] * matrix[2][1];
78 
79  inverse[0][1]= -matrix[0][1] * matrix[2][2] * matrix[3][3] +
80  matrix[0][1] * matrix[2][3] * matrix[3][2] +
81  matrix[2][1] * matrix[0][2] * matrix[3][3] -
82  matrix[2][1] * matrix[0][3] * matrix[3][2] -
83  matrix[3][1] * matrix[0][2] * matrix[2][3] +
84  matrix[3][1] * matrix[0][3] * matrix[2][2];
85 
86  inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] -
87  matrix[0][0] * matrix[2][3] * matrix[3][2] -
88  matrix[2][0] * matrix[0][2] * matrix[3][3] +
89  matrix[2][0] * matrix[0][3] * matrix[3][2] +
90  matrix[3][0] * matrix[0][2] * matrix[2][3] -
91  matrix[3][0] * matrix[0][3] * matrix[2][2];
92 
93  inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] +
94  matrix[0][0] * matrix[2][3] * matrix[3][1] +
95  matrix[2][0] * matrix[0][1] * matrix[3][3] -
96  matrix[2][0] * matrix[0][3] * matrix[3][1] -
97  matrix[3][0] * matrix[0][1] * matrix[2][3] +
98  matrix[3][0] * matrix[0][3] * matrix[2][1];
99 
100  inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] -
101  matrix[0][0] * matrix[2][2] * matrix[3][1] -
102  matrix[2][0] * matrix[0][1] * matrix[3][2] +
103  matrix[2][0] * matrix[0][2] * matrix[3][1] +
104  matrix[3][0] * matrix[0][1] * matrix[2][2] -
105  matrix[3][0] * matrix[0][2] * matrix[2][1];
106 
107  inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] -
108  matrix[0][1] * matrix[1][3] * matrix[3][2] -
109  matrix[1][1] * matrix[0][2] * matrix[3][3] +
110  matrix[1][1] * matrix[0][3] * matrix[3][2] +
111  matrix[3][1] * matrix[0][2] * matrix[1][3] -
112  matrix[3][1] * matrix[0][3] * matrix[1][2];
113 
114  inverse[1][2] = -matrix[0][0] * matrix[1][2] * matrix[3][3] +
115  matrix[0][0] * matrix[1][3] * matrix[3][2] +
116  matrix[1][0] * matrix[0][2] * matrix[3][3] -
117  matrix[1][0] * matrix[0][3] * matrix[3][2] -
118  matrix[3][0] * matrix[0][2] * matrix[1][3] +
119  matrix[3][0] * matrix[0][3] * matrix[1][2];
120 
121  inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] -
122  matrix[0][0] * matrix[1][3] * matrix[3][1] -
123  matrix[1][0] * matrix[0][1] * matrix[3][3] +
124  matrix[1][0] * matrix[0][3] * matrix[3][1] +
125  matrix[3][0] * matrix[0][1] * matrix[1][3] -
126  matrix[3][0] * matrix[0][3] * matrix[1][1];
127 
128  inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] +
129  matrix[0][0] * matrix[1][2] * matrix[3][1] +
130  matrix[1][0] * matrix[0][1] * matrix[3][2] -
131  matrix[1][0] * matrix[0][2] * matrix[3][1] -
132  matrix[3][0] * matrix[0][1] * matrix[1][2] +
133  matrix[3][0] * matrix[0][2] * matrix[1][1];
134 
135  inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] +
136  matrix[0][1] * matrix[1][3] * matrix[2][2] +
137  matrix[1][1] * matrix[0][2] * matrix[2][3] -
138  matrix[1][1] * matrix[0][3] * matrix[2][2] -
139  matrix[2][1] * matrix[0][2] * matrix[1][3] +
140  matrix[2][1] * matrix[0][3] * matrix[1][2];
141 
142  inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] -
143  matrix[0][0] * matrix[1][3] * matrix[2][2] -
144  matrix[1][0] * matrix[0][2] * matrix[2][3] +
145  matrix[1][0] * matrix[0][3] * matrix[2][2] +
146  matrix[2][0] * matrix[0][2] * matrix[1][3] -
147  matrix[2][0] * matrix[0][3] * matrix[1][2];
148 
149  inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] +
150  matrix[0][0] * matrix[1][3] * matrix[2][1] +
151  matrix[1][0] * matrix[0][1] * matrix[2][3] -
152  matrix[1][0] * matrix[0][3] * matrix[2][1] -
153  matrix[2][0] * matrix[0][1] * matrix[1][3] +
154  matrix[2][0] * matrix[0][3] * matrix[1][1];
155 
156  inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] -
157  matrix[0][0] * matrix[1][2] * matrix[2][1] -
158  matrix[1][0] * matrix[0][1] * matrix[2][2] +
159  matrix[1][0] * matrix[0][2] * matrix[2][1] +
160  matrix[2][0] * matrix[0][1] * matrix[1][2] -
161  matrix[2][0] * matrix[0][2] * matrix[1][1];
162 
163  K det = matrix[0][0] * inverse[0][0] + matrix[0][1] * inverse[1][0] + matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0];
164 
165  // return identity for singular or nearly singular matrices.
166  if (std::abs(det) < 1e-40) {
167  for (int i = 0; i < 4; ++i){
168  inverse[i][i] = 1.0;
169  }
170  return 1.0;
171  }
172  K inv_det = 1.0 / det;
173  inverse *= inv_det;
174 
175  return det;
176 }
177 } // end FMatrixHelp
178 
179 namespace ISTLUtility {
180 
182 template <typename K>
183 static inline void invertMatrix (FieldMatrix<K,1,1> &matrix)
184 {
185  FieldMatrix<K,1,1> A ( matrix );
186  FMatrixHelp::invertMatrix(A, matrix );
187 }
188 
190 template <typename K>
191 static inline void invertMatrix (FieldMatrix<K,2,2> &matrix)
192 {
193  FieldMatrix<K,2,2> A ( matrix );
194  FMatrixHelp::invertMatrix(A, matrix );
195 }
196 
198 template <typename K>
199 static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
200 {
201  FieldMatrix<K,3,3> A ( matrix );
202  FMatrixHelp::invertMatrix(A, matrix );
203 }
204 
206 template <typename K>
207 static inline void invertMatrix (FieldMatrix<K,4,4> &matrix)
208 {
209  FieldMatrix<K,4,4> A ( matrix );
210  FMatrixHelp::invertMatrix(A, matrix );
211 }
212 
214 template <typename K, int n>
215 static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)
216 {
217  matrix.invert();
218 }
219 
220 } // end ISTLUtility
221 
222 template <class Scalar, int n, int m>
223 class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
224 {
225 public:
226  typedef Dune::FieldMatrix<Scalar, n, m> BaseType;
227 
228  using BaseType :: operator= ;
229  using BaseType :: rows;
230  using BaseType :: cols;
231  explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {}
232  void invert()
233  {
234  ISTLUtility::invertMatrix( *this );
235  }
236  const BaseType& asBase() const { return static_cast< const BaseType& > (*this); }
237  BaseType& asBase() { return static_cast< BaseType& > (*this); }
238 };
239 
240 template<class K, int n, int m>
241 void
242 print_row (std::ostream& s, const MatrixBlock<K,n,m>& A,
243  typename FieldMatrix<K,n,m>::size_type I,
244  typename FieldMatrix<K,n,m>::size_type J,
245  typename FieldMatrix<K,n,m>::size_type therow, int width,
246  int precision)
247 {
248  print_row(s, A.asBase(), I, J, therow, width, precision);
249 }
250 
251 template<class K, int n, int m>
252 K& firstmatrixelement (MatrixBlock<K,n,m>& A)
253 {
254  return firstmatrixelement( A.asBase() );
255 }
256 
257 
258 
259 template<typename Scalar, int n, int m>
260 struct MatrixDimension< MatrixBlock< Scalar, n, m > >
261 : public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType >
262 {
263 };
264 
265 
266 #if HAVE_UMFPACK
267 
271 template<typename T, typename A, int n, int m>
272 class UMFPack<BCRSMatrix<MatrixBlock<T,n,m>, A> >
273  : public UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> >
274 {
275  typedef UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
276  typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
277 
278 public:
279  typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
280 
281  UMFPack(const RealMatrix& matrix, int verbose, bool)
282  : Base(reinterpret_cast<const Matrix&>(matrix), verbose)
283  {}
284 };
285 #endif
286 
287 #if HAVE_SUPERLU
288 
292 template<typename T, typename A, int n, int m>
293 class SuperLU<BCRSMatrix<MatrixBlock<T,n,m>, A> >
294  : public SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> >
295 {
296  typedef SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
297  typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
298 
299 public:
300  typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
301 
302  SuperLU(const RealMatrix& matrix, int verbose, bool reuse=true)
303  : Base(reinterpret_cast<const Matrix&>(matrix), verbose, reuse)
304  {}
305 };
306 #endif
307 
308 
309 } // end namespace Dune
310 
311 namespace Opm
312 {
322  template < class MatrixBlockType, class VectorBlockType, int pressureIndex=0 >
324  {
325  typedef typename MatrixBlockType :: field_type Scalar;
326 
327  typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
328  typedef Dune::BlockVector<VectorBlockType> Vector;
329 
330  public:
331  typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType;
332 
339  const boost::any& parallelInformation_arg=boost::any())
340  : iterations_( 0 ),
341  parallelInformation_(parallelInformation_arg),
342  isIORank_(isIORank(parallelInformation_arg)),
343  parameters_( param )
344  {
345  }
346 
351  ISTLSolver(const ParameterGroup& param,
352  const boost::any& parallelInformation_arg=boost::any())
353  : iterations_( 0 ),
354  parallelInformation_(parallelInformation_arg),
355  isIORank_(isIORank(parallelInformation_arg)),
356  parameters_( param )
357  {
358  }
359 
360  // dummy method that is not implemented for this class
362  {
363  OPM_THROW(std::logic_error,"This method is not implemented");
364  return SolutionVector();
365  }
366 
372 
374  int iterations () const { return iterations_; }
375 
377  const boost::any& parallelInformation() const { return parallelInformation_; }
378 
379  public:
383  template<int category=Dune::SolverCategory::sequential, class LinearOperator, class POrComm>
384  void constructPreconditionerAndSolve(LinearOperator& linearOperator,
385  Vector& x, Vector& istlb,
386  const POrComm& parallelInformation_arg,
387  Dune::InverseOperatorResult& result) const
388  {
389  // Construct scalar product.
390  typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
391  typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
392  SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
393 
394  // Communicate if parallel.
395  parallelInformation_arg.copyOwnerToAll(istlb, istlb);
396 
397 #if FLOW_SUPPORT_AMG // activate AMG if either flow_ebos is used or UMFPack is not available
398  if( parameters_.linear_solver_use_amg_ )
399  {
401  typedef typename CPRSelectorType::AMG AMG;
402  typedef typename CPRSelectorType::Operator MatrixOperator;
403 
404  std::unique_ptr< AMG > amg;
405  std::unique_ptr< MatrixOperator > opA;
406 
407  if( ! std::is_same< LinearOperator, MatrixOperator > :: value )
408  {
409  // create new operator in case linear operator and matrix operator differ
410  opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) );
411  }
412 
413  const double relax = 1.0;
414 
415  // Construct preconditioner.
416  constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
417 
418  // Solve.
419  solve(linearOperator, x, istlb, *sp, *amg, result);
420  }
421  else
422 #endif
423  {
424  // Construct preconditioner.
425  auto precond = constructPrecond(linearOperator, parallelInformation_arg);
426 
427  // Solve.
428  solve(linearOperator, x, istlb, *sp, *precond, result);
429  }
430  }
431 
432 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2 , 5)
433  typedef ParallelOverlappingILU0<Matrix,Vector,Vector> SeqPreconditioner;
434 #else
435  typedef ParallelOverlappingILU0<Dune::BCRSMatrix<Dune::MatrixBlock<typename Matrix::field_type,
436  Matrix::block_type::rows,
437  Matrix::block_type::cols> >,
438  Vector, Vector> SeqPreconditioner;
439 #endif
440 
441  template <class Operator>
442  std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
443  {
444  const double relax = parameters_.ilu_relaxation_;
445  const int ilu_fillin = parameters_.ilu_fillin_level_;
446  std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax));
447  return precond;
448  }
449 
450 #if HAVE_MPI
451  typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
452 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2 , 5)
453  typedef ParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
454 #else
455  typedef ParallelOverlappingILU0<Dune::BCRSMatrix<Dune::MatrixBlock<typename Matrix::field_type,
456  Matrix::block_type::rows,
457  Matrix::block_type::cols> >,
458  Vector, Vector, Comm> ParPreconditioner;
459 #endif
460  template <class Operator>
461  std::unique_ptr<ParPreconditioner>
462  constructPrecond(Operator& opA, const Comm& comm) const
463  {
464  typedef std::unique_ptr<ParPreconditioner> Pointer;
465  const double relax = parameters_.ilu_relaxation_;
466  return Pointer(new ParPreconditioner(opA.getmat(), comm, relax));
467  }
468 #endif
469 
470  template <class LinearOperator, class MatrixOperator, class POrComm, class AMG >
471  void
472  constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const
473  {
474  ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, comm, amg );
475  }
476 
477 
478  template <class MatrixOperator, class POrComm, class AMG >
479  void
480  constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const
481  {
482  ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax, comm, amg );
483  }
484 
486  template <class Operator, class ScalarProd, class Precond>
487  void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const
488  {
489  // TODO: Revise when linear solvers interface opm-core is done
490  // Construct linear solver.
491  // GMRes solver
492  int verbosity = ( isIORank_ ) ? parameters_.linear_solver_verbosity_ : 0;
493 
494  if ( parameters_.newton_use_gmres_ ) {
495  Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond,
496  parameters_.linear_solver_reduction_,
497  parameters_.linear_solver_restart_,
498  parameters_.linear_solver_maxiter_,
499  verbosity);
500  // Solve system.
501  linsolve.apply(x, istlb, result);
502  }
503  else { // BiCGstab solver
504  Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond,
505  parameters_.linear_solver_reduction_,
506  parameters_.linear_solver_maxiter_,
507  verbosity);
508  // Solve system.
509  linsolve.apply(x, istlb, result);
510  }
511  }
512 
513 
520  void solve(Matrix& A, Vector& x, Vector& b ) const
521  {
522  // Parallel version is deactivated until we figure out how to do it properly.
523 #if HAVE_MPI
524  if (parallelInformation_.type() == typeid(ParallelISTLInformation))
525  {
526  typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
527  const ParallelISTLInformation& info =
528  boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
529  Comm istlComm(info.communicator());
530 
531  // Construct operator, scalar product and vectors needed.
532  typedef Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector,Comm> Operator;
533  Operator opA(A, istlComm);
534  solve( opA, x, b, istlComm );
535  }
536  else
537 #endif
538  {
539  // Construct operator, scalar product and vectors needed.
540  Dune::MatrixAdapter< Matrix, Vector, Vector> opA( A );
541  solve( opA, x, b );
542  }
543  }
544 
551  template <class Operator, class Comm >
552  void solve(Operator& opA, Vector& x, Vector& b, Comm& comm) const
553  {
554  Dune::InverseOperatorResult result;
555  // Parallel version is deactivated until we figure out how to do it properly.
556 #if HAVE_MPI
557  if (parallelInformation_.type() == typeid(ParallelISTLInformation))
558  {
559  const size_t size = opA.getmat().N();
560  const ParallelISTLInformation& info =
561  boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
562 
563  // As we use a dune-istl with block size np the number of components
564  // per parallel is only one.
565  info.copyValuesTo(comm.indexSet(), comm.remoteIndices(),
566  size, 1);
567  // Construct operator, scalar product and vectors needed.
568  constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, b, comm, result);
569  }
570  else
571 #endif
572  {
573  OPM_THROW(std::logic_error,"this method if for parallel solve only");
574  }
575 
576  checkConvergence( result );
577  }
578 
585  template <class Operator>
586  void solve(Operator& opA, Vector& x, Vector& b ) const
587  {
588  Dune::InverseOperatorResult result;
589  // Construct operator, scalar product and vectors needed.
590  Dune::Amg::SequentialInformation info;
591  constructPreconditionerAndSolve(opA, x, b, info, result);
592  checkConvergence( result );
593  }
594 
595  void checkConvergence( const Dune::InverseOperatorResult& result ) const
596  {
597  // store number of iterations
598  iterations_ = result.iterations;
599 
600  // Check for failure of linear solver.
601  if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
602  const std::string msg("Convergence failure for linear solver.");
603  OPM_THROW_NOLOG(LinearSolverProblem, msg);
604  }
605  }
606  protected:
607  mutable int iterations_;
608  boost::any parallelInformation_;
609  bool isIORank_;
610 
611  NewtonIterationBlackoilInterleavedParameters parameters_;
612  }; // end ISTLSolver
613 
614 } // namespace Opm
615 #endif
ISTLSolver(const ParameterGroup &param, const boost::any &parallelInformation_arg=boost::any())
Construct a system solver.
Definition: ISTLSolver.hpp:351
LinearisedBlackoilResidual::ADB::V SolutionVector
Return type for linearSolve(). A simple, non-ad vector type.
Definition: NewtonIterationBlackoilInterface.hpp:35
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition: ParallelOverlappingILU0.hpp:35
Definition: ISTLSolver.hpp:223
ISTLSolver(const NewtonIterationBlackoilInterleavedParameters &param, const boost::any &parallelInformation_arg=boost::any())
Construct a system solver.
Definition: ISTLSolver.hpp:338
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition: ISTLSolver.hpp:323
void solve(Operator &opA, Vector &x, Vector &istlb, ScalarProd &sp, Precond &precond, Dune::InverseOperatorResult &result) const
Solve the system using the given preconditioner and scalar product.
Definition: ISTLSolver.hpp:487
void solve(Matrix &A, Vector &x, Vector &b) const
Solve the linear system Ax = b, with A being the combined derivative matrix of the residual and b bei...
Definition: ISTLSolver.hpp:520
void solve(Operator &opA, Vector &x, Vector &b) const
Solve the linear system Ax = b, with A being the combined derivative matrix of the residual and b bei...
Definition: ISTLSolver.hpp:586
SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual &) const
Solve the linear system Ax = b, with A being the combined derivative matrix of the residual and b bei...
Definition: ISTLSolver.hpp:361
bool isIORank(const boost::any &parallel_info)
Return true if this is a serial run, or rank zero on an MPI run.
Definition: NewtonIterationUtilities.cpp:294
void constructPreconditionerAndSolve(LinearOperator &linearOperator, Vector &x, Vector &istlb, const POrComm &parallelInformation_arg, Dune::InverseOperatorResult &result) const
construct the CPR preconditioner and the solver.
Definition: ISTLSolver.hpp:384
Residual structure of the fully implicit solver.
Definition: LinearisedBlackoilResidual.hpp:47
void solve(Operator &opA, Vector &x, Vector &b, Comm &comm) const
Solve the linear system Ax = b, with A being the combined derivative matrix of the residual and b bei...
Definition: ISTLSolver.hpp:552
const boost::any & parallelInformation() const
Get the information about the parallelization of the grid.
Definition: ISTLSolver.hpp:377
int iterations() const
Solve the system of linear equations Ax = b, with A being the combined derivative matrix of the resid...
Definition: ISTLSolver.hpp:374
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
This class carries all parameters for the NewtonIterationBlackoilInterleaved class.
Definition: NewtonIterationBlackoilInterleaved.hpp:37
A traits class for selecting the types of the preconditioner.
Definition: CPRPreconditioner.hpp:67