20 #ifndef OPM_ISTLSOLVER_HEADER_INCLUDED 21 #define OPM_ISTLSOLVER_HEADER_INCLUDED 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> 31 #include <opm/common/Exceptions.hpp> 32 #include <opm/core/linalg/ParallelIstlInformation.hpp> 33 #include <opm/common/utility/platform_dependent/disable_warnings.h> 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> 42 #include <opm/common/utility/platform_dependent/reenable_warnings.h> 46 namespace FMatrixHelp {
49 static inline K invertMatrix (
const FieldMatrix<K,4,4> &matrix, FieldMatrix<K,4,4> &inverse)
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
166 if (std::abs(det) < 1e-40) {
167 for (
int i = 0; i < 4; ++i){
172 K inv_det = 1.0 / det;
179 namespace ISTLUtility {
182 template <
typename K>
183 static inline void invertMatrix (FieldMatrix<K,1,1> &matrix)
185 FieldMatrix<K,1,1> A ( matrix );
186 FMatrixHelp::invertMatrix(A, matrix );
190 template <
typename K>
191 static inline void invertMatrix (FieldMatrix<K,2,2> &matrix)
193 FieldMatrix<K,2,2> A ( matrix );
194 FMatrixHelp::invertMatrix(A, matrix );
198 template <
typename K>
199 static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
201 FieldMatrix<K,3,3> A ( matrix );
202 FMatrixHelp::invertMatrix(A, matrix );
206 template <
typename K>
207 static inline void invertMatrix (FieldMatrix<K,4,4> &matrix)
209 FieldMatrix<K,4,4> A ( matrix );
210 FMatrixHelp::invertMatrix(A, matrix );
214 template <
typename K,
int n>
215 static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)
222 template <
class Scalar,
int n,
int m>
226 typedef Dune::FieldMatrix<Scalar, n, m> BaseType;
228 using BaseType :: operator= ;
229 using BaseType :: rows;
230 using BaseType :: cols;
231 explicit MatrixBlock(
const Scalar scalar = 0 ) : BaseType( scalar ) {}
234 ISTLUtility::invertMatrix( *
this );
236 const BaseType& asBase()
const {
return static_cast< const BaseType&
> (*this); }
237 BaseType& asBase() {
return static_cast< BaseType&
> (*this); }
240 template<
class K,
int n,
int m>
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,
248 print_row(s, A.asBase(), I, J, therow, width, precision);
251 template<
class K,
int n,
int m>
252 K& firstmatrixelement (MatrixBlock<K,n,m>& A)
254 return firstmatrixelement( A.asBase() );
259 template<
typename Scalar,
int n,
int m>
261 :
public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType >
271 template<
typename T,
typename A,
int n,
int m>
273 :
public UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> >
275 typedef UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
276 typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
279 typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
281 UMFPack(
const RealMatrix& matrix,
int verbose,
bool)
282 : Base(reinterpret_cast<const Matrix&>(matrix), verbose)
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> >
296 typedef SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
297 typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
300 typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
302 SuperLU(
const RealMatrix& matrix,
int verbose,
bool reuse=
true)
303 : Base(reinterpret_cast<const Matrix&>(matrix), verbose, reuse)
322 template <
class MatrixBlockType,
class VectorBlockType,
int pressureIndex=0 >
325 typedef typename MatrixBlockType :: field_type Scalar;
327 typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
328 typedef Dune::BlockVector<VectorBlockType> Vector;
331 typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType;
339 const boost::any& parallelInformation_arg=boost::any())
341 parallelInformation_(parallelInformation_arg),
342 isIORank_(
isIORank(parallelInformation_arg)),
352 const boost::any& parallelInformation_arg=boost::any())
354 parallelInformation_(parallelInformation_arg),
355 isIORank_(
isIORank(parallelInformation_arg)),
363 OPM_THROW(std::logic_error,
"This method is not implemented");
383 template<
int category=Dune::SolverCategory::sequential,
class LinearOperator,
class POrComm>
385 Vector& x, Vector& istlb,
386 const POrComm& parallelInformation_arg,
387 Dune::InverseOperatorResult& result)
const 390 typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
391 typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
392 SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
395 parallelInformation_arg.copyOwnerToAll(istlb, istlb);
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_ )
401 typedef typename CPRSelectorType::AMG AMG;
402 typedef typename CPRSelectorType::Operator MatrixOperator;
404 std::unique_ptr< AMG > amg;
405 std::unique_ptr< MatrixOperator > opA;
407 if( ! std::is_same< LinearOperator, MatrixOperator > :: value )
410 opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) );
413 const double relax = 1.0;
416 constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
419 solve(linearOperator, x, istlb, *sp, *amg, result);
425 auto precond = constructPrecond(linearOperator, parallelInformation_arg);
428 solve(linearOperator, x, istlb, *sp, *precond, result);
432 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2 , 5) 436 Matrix::block_type::rows,
437 Matrix::block_type::cols> >,
438 Vector, Vector> SeqPreconditioner;
441 template <
class Operator>
442 std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA,
const Dune::Amg::SequentialInformation&)
const 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));
451 typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
452 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2 , 5) 453 typedef ParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
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;
460 template <
class Operator>
461 std::unique_ptr<ParPreconditioner>
462 constructPrecond(Operator& opA,
const Comm& comm)
const 464 typedef std::unique_ptr<ParPreconditioner> Pointer;
465 const double relax = parameters_.ilu_relaxation_;
466 return Pointer(
new ParPreconditioner(opA.getmat(), comm, relax));
470 template <
class LinearOperator,
class MatrixOperator,
class POrComm,
class AMG >
472 constructAMGPrecond(LinearOperator& ,
const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA,
const double relax )
const 474 ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, comm, amg );
478 template <
class MatrixOperator,
class POrComm,
class AMG >
480 constructAMGPrecond(MatrixOperator& opA,
const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&,
const double relax )
const 482 ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax, comm, amg );
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 492 int verbosity = ( isIORank_ ) ? parameters_.linear_solver_verbosity_ : 0;
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_,
501 linsolve.apply(x, istlb, result);
504 Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond,
505 parameters_.linear_solver_reduction_,
506 parameters_.linear_solver_maxiter_,
509 linsolve.apply(x, istlb, result);
520 void solve(Matrix& A, Vector& x, Vector& b )
const 524 if (parallelInformation_.type() ==
typeid(ParallelISTLInformation))
526 typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
527 const ParallelISTLInformation& info =
528 boost::any_cast<
const ParallelISTLInformation&>( parallelInformation_);
529 Comm istlComm(info.communicator());
532 typedef Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector,Comm> Operator;
533 Operator opA(A, istlComm);
534 solve( opA, x, b, istlComm );
540 Dune::MatrixAdapter< Matrix, Vector, Vector> opA( A );
551 template <
class Operator,
class Comm >
552 void solve(Operator& opA, Vector& x, Vector& b, Comm& comm)
const 554 Dune::InverseOperatorResult result;
557 if (parallelInformation_.type() ==
typeid(ParallelISTLInformation))
559 const size_t size = opA.getmat().N();
560 const ParallelISTLInformation& info =
561 boost::any_cast<
const ParallelISTLInformation&>( parallelInformation_);
565 info.copyValuesTo(comm.indexSet(), comm.remoteIndices(),
568 constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, b, comm, result);
573 OPM_THROW(std::logic_error,
"this method if for parallel solve only");
576 checkConvergence( result );
585 template <
class Operator>
586 void solve(Operator& opA, Vector& x, Vector& b )
const 588 Dune::InverseOperatorResult result;
590 Dune::Amg::SequentialInformation info;
592 checkConvergence( result );
595 void checkConvergence(
const Dune::InverseOperatorResult& result )
const 598 iterations_ = result.iterations;
601 if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
602 const std::string msg(
"Convergence failure for linear solver.");
603 OPM_THROW_NOLOG(LinearSolverProblem, msg);
607 mutable int iterations_;
608 boost::any parallelInformation_;
611 NewtonIterationBlackoilInterleavedParameters parameters_;
ISTLSolver(const ParameterGroup ¶m, const boost::any ¶llelInformation_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
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
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
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
Definition: ISTLSolver.hpp:223
Definition: ISTLSolver.hpp:44
ISTLSolver(const NewtonIterationBlackoilInterleavedParameters ¶m, const boost::any ¶llelInformation_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 &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
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
void constructPreconditionerAndSolve(LinearOperator &linearOperator, Vector &x, Vector &istlb, const POrComm ¶llelInformation_arg, Dune::InverseOperatorResult &result) const
construct the CPR preconditioner and the solver.
Definition: ISTLSolver.hpp:384
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
bool isIORank(const boost::any ¶llel_info)
Return true if this is a serial run, or rank zero on an MPI run.
Definition: NewtonIterationUtilities.cpp:294
Residual structure of the fully implicit solver.
Definition: LinearisedBlackoilResidual.hpp:47
const boost::any & parallelInformation() const
Get the information about the parallelization of the grid.
Definition: ISTLSolver.hpp:377
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
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