mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-12-30 11:06:55 -06:00
602f3962f3
3x3 matrix block inversion in FieldMatrix is numerically unstable including version 2.5.0. Therefore the previous if clause was wrong as it activated the use of FieldMatrix already for 2.5.0 (the version in Debian stable). With this commit we use MatrixBlock for version 2.5.0.
620 lines
25 KiB
C++
620 lines
25 KiB
C++
/*
|
|
Copyright 2016 IRIS AS
|
|
|
|
This file is part of the Open Porous Media project (OPM).
|
|
|
|
OPM is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OPM is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#ifndef OPM_ISTLSOLVER_HEADER_INCLUDED
|
|
#define OPM_ISTLSOLVER_HEADER_INCLUDED
|
|
|
|
#include <opm/autodiff/AdditionalObjectDeleter.hpp>
|
|
#include <opm/autodiff/CPRPreconditioner.hpp>
|
|
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
|
|
#include <opm/autodiff/NewtonIterationUtilities.hpp>
|
|
#include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
|
|
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
|
|
#include <opm/autodiff/AutoDiffHelpers.hpp>
|
|
|
|
#include <opm/common/Exceptions.hpp>
|
|
#include <opm/core/linalg/ParallelIstlInformation.hpp>
|
|
#include <opm/common/utility/platform_dependent/disable_warnings.h>
|
|
|
|
#include <dune/istl/scalarproducts.hh>
|
|
#include <dune/istl/operators.hh>
|
|
#include <dune/istl/preconditioners.hh>
|
|
#include <dune/istl/solvers.hh>
|
|
#include <dune/istl/owneroverlapcopy.hh>
|
|
#include <dune/istl/paamg/amg.hh>
|
|
|
|
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
|
|
|
|
namespace Dune
|
|
{
|
|
namespace FMatrixHelp {
|
|
//! invert 4x4 Matrix without changing the original matrix
|
|
template <typename K>
|
|
static inline K invertMatrix (const FieldMatrix<K,4,4> &matrix, FieldMatrix<K,4,4> &inverse)
|
|
{
|
|
inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] -
|
|
matrix[1][1] * matrix[2][3] * matrix[3][2] -
|
|
matrix[2][1] * matrix[1][2] * matrix[3][3] +
|
|
matrix[2][1] * matrix[1][3] * matrix[3][2] +
|
|
matrix[3][1] * matrix[1][2] * matrix[2][3] -
|
|
matrix[3][1] * matrix[1][3] * matrix[2][2];
|
|
|
|
inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] +
|
|
matrix[1][0] * matrix[2][3] * matrix[3][2] +
|
|
matrix[2][0] * matrix[1][2] * matrix[3][3] -
|
|
matrix[2][0] * matrix[1][3] * matrix[3][2] -
|
|
matrix[3][0] * matrix[1][2] * matrix[2][3] +
|
|
matrix[3][0] * matrix[1][3] * matrix[2][2];
|
|
|
|
inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] -
|
|
matrix[1][0] * matrix[2][3] * matrix[3][1] -
|
|
matrix[2][0] * matrix[1][1] * matrix[3][3] +
|
|
matrix[2][0] * matrix[1][3] * matrix[3][1] +
|
|
matrix[3][0] * matrix[1][1] * matrix[2][3] -
|
|
matrix[3][0] * matrix[1][3] * matrix[2][1];
|
|
|
|
inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] +
|
|
matrix[1][0] * matrix[2][2] * matrix[3][1] +
|
|
matrix[2][0] * matrix[1][1] * matrix[3][2] -
|
|
matrix[2][0] * matrix[1][2] * matrix[3][1] -
|
|
matrix[3][0] * matrix[1][1] * matrix[2][2] +
|
|
matrix[3][0] * matrix[1][2] * matrix[2][1];
|
|
|
|
inverse[0][1]= -matrix[0][1] * matrix[2][2] * matrix[3][3] +
|
|
matrix[0][1] * matrix[2][3] * matrix[3][2] +
|
|
matrix[2][1] * matrix[0][2] * matrix[3][3] -
|
|
matrix[2][1] * matrix[0][3] * matrix[3][2] -
|
|
matrix[3][1] * matrix[0][2] * matrix[2][3] +
|
|
matrix[3][1] * matrix[0][3] * matrix[2][2];
|
|
|
|
inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] -
|
|
matrix[0][0] * matrix[2][3] * matrix[3][2] -
|
|
matrix[2][0] * matrix[0][2] * matrix[3][3] +
|
|
matrix[2][0] * matrix[0][3] * matrix[3][2] +
|
|
matrix[3][0] * matrix[0][2] * matrix[2][3] -
|
|
matrix[3][0] * matrix[0][3] * matrix[2][2];
|
|
|
|
inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] +
|
|
matrix[0][0] * matrix[2][3] * matrix[3][1] +
|
|
matrix[2][0] * matrix[0][1] * matrix[3][3] -
|
|
matrix[2][0] * matrix[0][3] * matrix[3][1] -
|
|
matrix[3][0] * matrix[0][1] * matrix[2][3] +
|
|
matrix[3][0] * matrix[0][3] * matrix[2][1];
|
|
|
|
inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] -
|
|
matrix[0][0] * matrix[2][2] * matrix[3][1] -
|
|
matrix[2][0] * matrix[0][1] * matrix[3][2] +
|
|
matrix[2][0] * matrix[0][2] * matrix[3][1] +
|
|
matrix[3][0] * matrix[0][1] * matrix[2][2] -
|
|
matrix[3][0] * matrix[0][2] * matrix[2][1];
|
|
|
|
inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] -
|
|
matrix[0][1] * matrix[1][3] * matrix[3][2] -
|
|
matrix[1][1] * matrix[0][2] * matrix[3][3] +
|
|
matrix[1][1] * matrix[0][3] * matrix[3][2] +
|
|
matrix[3][1] * matrix[0][2] * matrix[1][3] -
|
|
matrix[3][1] * matrix[0][3] * matrix[1][2];
|
|
|
|
inverse[1][2] = -matrix[0][0] * matrix[1][2] * matrix[3][3] +
|
|
matrix[0][0] * matrix[1][3] * matrix[3][2] +
|
|
matrix[1][0] * matrix[0][2] * matrix[3][3] -
|
|
matrix[1][0] * matrix[0][3] * matrix[3][2] -
|
|
matrix[3][0] * matrix[0][2] * matrix[1][3] +
|
|
matrix[3][0] * matrix[0][3] * matrix[1][2];
|
|
|
|
inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] -
|
|
matrix[0][0] * matrix[1][3] * matrix[3][1] -
|
|
matrix[1][0] * matrix[0][1] * matrix[3][3] +
|
|
matrix[1][0] * matrix[0][3] * matrix[3][1] +
|
|
matrix[3][0] * matrix[0][1] * matrix[1][3] -
|
|
matrix[3][0] * matrix[0][3] * matrix[1][1];
|
|
|
|
inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] +
|
|
matrix[0][0] * matrix[1][2] * matrix[3][1] +
|
|
matrix[1][0] * matrix[0][1] * matrix[3][2] -
|
|
matrix[1][0] * matrix[0][2] * matrix[3][1] -
|
|
matrix[3][0] * matrix[0][1] * matrix[1][2] +
|
|
matrix[3][0] * matrix[0][2] * matrix[1][1];
|
|
|
|
inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] +
|
|
matrix[0][1] * matrix[1][3] * matrix[2][2] +
|
|
matrix[1][1] * matrix[0][2] * matrix[2][3] -
|
|
matrix[1][1] * matrix[0][3] * matrix[2][2] -
|
|
matrix[2][1] * matrix[0][2] * matrix[1][3] +
|
|
matrix[2][1] * matrix[0][3] * matrix[1][2];
|
|
|
|
inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] -
|
|
matrix[0][0] * matrix[1][3] * matrix[2][2] -
|
|
matrix[1][0] * matrix[0][2] * matrix[2][3] +
|
|
matrix[1][0] * matrix[0][3] * matrix[2][2] +
|
|
matrix[2][0] * matrix[0][2] * matrix[1][3] -
|
|
matrix[2][0] * matrix[0][3] * matrix[1][2];
|
|
|
|
inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] +
|
|
matrix[0][0] * matrix[1][3] * matrix[2][1] +
|
|
matrix[1][0] * matrix[0][1] * matrix[2][3] -
|
|
matrix[1][0] * matrix[0][3] * matrix[2][1] -
|
|
matrix[2][0] * matrix[0][1] * matrix[1][3] +
|
|
matrix[2][0] * matrix[0][3] * matrix[1][1];
|
|
|
|
inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] -
|
|
matrix[0][0] * matrix[1][2] * matrix[2][1] -
|
|
matrix[1][0] * matrix[0][1] * matrix[2][2] +
|
|
matrix[1][0] * matrix[0][2] * matrix[2][1] +
|
|
matrix[2][0] * matrix[0][1] * matrix[1][2] -
|
|
matrix[2][0] * matrix[0][2] * matrix[1][1];
|
|
|
|
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];
|
|
|
|
// return identity for singular or nearly singular matrices.
|
|
if (std::abs(det) < 1e-40) {
|
|
for (int i = 0; i < 4; ++i){
|
|
inverse[i][i] = 1.0;
|
|
}
|
|
return 1.0;
|
|
}
|
|
K inv_det = 1.0 / det;
|
|
inverse *= inv_det;
|
|
|
|
return det;
|
|
}
|
|
} // end FMatrixHelp
|
|
|
|
namespace ISTLUtility {
|
|
|
|
//! invert matrix by calling FMatrixHelp::invert
|
|
template <typename K>
|
|
static inline void invertMatrix (FieldMatrix<K,1,1> &matrix)
|
|
{
|
|
FieldMatrix<K,1,1> A ( matrix );
|
|
FMatrixHelp::invertMatrix(A, matrix );
|
|
}
|
|
|
|
//! invert matrix by calling FMatrixHelp::invert
|
|
template <typename K>
|
|
static inline void invertMatrix (FieldMatrix<K,2,2> &matrix)
|
|
{
|
|
FieldMatrix<K,2,2> A ( matrix );
|
|
FMatrixHelp::invertMatrix(A, matrix );
|
|
}
|
|
|
|
//! invert matrix by calling FMatrixHelp::invert
|
|
template <typename K>
|
|
static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
|
|
{
|
|
FieldMatrix<K,3,3> A ( matrix );
|
|
FMatrixHelp::invertMatrix(A, matrix );
|
|
}
|
|
|
|
//! invert matrix by calling FMatrixHelp::invert
|
|
template <typename K>
|
|
static inline void invertMatrix (FieldMatrix<K,4,4> &matrix)
|
|
{
|
|
FieldMatrix<K,4,4> A ( matrix );
|
|
FMatrixHelp::invertMatrix(A, matrix );
|
|
}
|
|
|
|
//! invert matrix by calling matrix.invert
|
|
template <typename K, int n>
|
|
static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)
|
|
{
|
|
matrix.invert();
|
|
}
|
|
|
|
} // end ISTLUtility
|
|
|
|
template <class Scalar, int n, int m>
|
|
class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
|
|
{
|
|
public:
|
|
typedef Dune::FieldMatrix<Scalar, n, m> BaseType;
|
|
|
|
using BaseType :: operator= ;
|
|
using BaseType :: rows;
|
|
using BaseType :: cols;
|
|
explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {}
|
|
void invert()
|
|
{
|
|
ISTLUtility::invertMatrix( *this );
|
|
}
|
|
const BaseType& asBase() const { return static_cast< const BaseType& > (*this); }
|
|
BaseType& asBase() { return static_cast< BaseType& > (*this); }
|
|
};
|
|
|
|
template<class K, int n, int m>
|
|
void
|
|
print_row (std::ostream& s, const MatrixBlock<K,n,m>& A,
|
|
typename FieldMatrix<K,n,m>::size_type I,
|
|
typename FieldMatrix<K,n,m>::size_type J,
|
|
typename FieldMatrix<K,n,m>::size_type therow, int width,
|
|
int precision)
|
|
{
|
|
print_row(s, A.asBase(), I, J, therow, width, precision);
|
|
}
|
|
|
|
template<class K, int n, int m>
|
|
K& firstmatrixelement (MatrixBlock<K,n,m>& A)
|
|
{
|
|
return firstmatrixelement( A.asBase() );
|
|
}
|
|
|
|
|
|
|
|
template<typename Scalar, int n, int m>
|
|
struct MatrixDimension< MatrixBlock< Scalar, n, m > >
|
|
: public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType >
|
|
{
|
|
};
|
|
|
|
|
|
#if HAVE_UMFPACK
|
|
|
|
/// \brief UMFPack specialization for MatrixBlock to make AMG happy
|
|
///
|
|
/// Without this the empty default implementation would be used.
|
|
template<typename T, typename A, int n, int m>
|
|
class UMFPack<BCRSMatrix<MatrixBlock<T,n,m>, A> >
|
|
: public UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> >
|
|
{
|
|
typedef UMFPack<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
|
|
typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
|
|
|
|
public:
|
|
typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
|
|
|
|
UMFPack(const RealMatrix& matrix, int verbose, bool)
|
|
: Base(reinterpret_cast<const Matrix&>(matrix), verbose)
|
|
{}
|
|
};
|
|
#endif
|
|
|
|
#if HAVE_SUPERLU
|
|
|
|
/// \brief SuperLU specialization for MatrixBlock to make AMG happy
|
|
///
|
|
/// Without this the empty default implementation would be used.
|
|
template<typename T, typename A, int n, int m>
|
|
class SuperLU<BCRSMatrix<MatrixBlock<T,n,m>, A> >
|
|
: public SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> >
|
|
{
|
|
typedef SuperLU<BCRSMatrix<FieldMatrix<T,n,m>, A> > Base;
|
|
typedef BCRSMatrix<FieldMatrix<T,n,m>, A> Matrix;
|
|
|
|
public:
|
|
typedef BCRSMatrix<MatrixBlock<T,n,m>, A> RealMatrix;
|
|
|
|
SuperLU(const RealMatrix& matrix, int verbose, bool reuse=true)
|
|
: Base(reinterpret_cast<const Matrix&>(matrix), verbose, reuse)
|
|
{}
|
|
};
|
|
#endif
|
|
|
|
|
|
} // end namespace Dune
|
|
|
|
namespace Opm
|
|
{
|
|
/// This class solves the fully implicit black-oil system by
|
|
/// solving the reduced system (after eliminating well variables)
|
|
/// as a block-structured matrix (one block for all cell variables) for a fixed
|
|
/// number of cell variables np .
|
|
/// \tparam MatrixBlockType The type of the matrix block used.
|
|
/// \tparam VectorBlockType The type of the vector block used.
|
|
/// \tparam pressureIndex The index of the pressure component in the vector
|
|
/// vector block. It is used to guide the AMG coarsening.
|
|
/// Default is zero.
|
|
template < class MatrixBlockType, class VectorBlockType, int pressureIndex=0 >
|
|
class ISTLSolver : public NewtonIterationBlackoilInterface
|
|
{
|
|
typedef typename MatrixBlockType :: field_type Scalar;
|
|
|
|
typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
|
|
typedef Dune::BlockVector<VectorBlockType> Vector;
|
|
|
|
public:
|
|
typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType;
|
|
|
|
typedef NewtonIterationBlackoilInterface :: SolutionVector SolutionVector;
|
|
/// Construct a system solver.
|
|
/// \param[in] param parameters controlling the behaviour of the linear solvers
|
|
/// \param[in] parallelInformation In the case of a parallel run
|
|
/// with dune-istl the information about the parallelization.
|
|
ISTLSolver(const NewtonIterationBlackoilInterleavedParameters& param,
|
|
const boost::any& parallelInformation_arg=boost::any())
|
|
: iterations_( 0 ),
|
|
parallelInformation_(parallelInformation_arg),
|
|
isIORank_(isIORank(parallelInformation_arg)),
|
|
parameters_( param )
|
|
{
|
|
}
|
|
|
|
/// Construct a system solver.
|
|
/// \param[in] param ParameterGroup controlling the behaviour of the linear solvers
|
|
/// \param[in] parallelInformation In the case of a parallel run
|
|
/// with dune-istl the information about the parallelization.
|
|
ISTLSolver(const ParameterGroup& param,
|
|
const boost::any& parallelInformation_arg=boost::any())
|
|
: iterations_( 0 ),
|
|
parallelInformation_(parallelInformation_arg),
|
|
isIORank_(isIORank(parallelInformation_arg)),
|
|
parameters_( param )
|
|
{
|
|
}
|
|
|
|
// dummy method that is not implemented for this class
|
|
SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual&) const
|
|
{
|
|
OPM_THROW(std::logic_error,"This method is not implemented");
|
|
return SolutionVector();
|
|
}
|
|
|
|
/// Solve the system of linear equations Ax = b, with A being the
|
|
/// combined derivative matrix of the residual and b
|
|
/// being the residual itself.
|
|
/// \param[in] residual residual object containing A and b.
|
|
/// \return the solution x
|
|
|
|
/// \copydoc NewtonIterationBlackoilInterface::iterations
|
|
int iterations () const { return iterations_; }
|
|
|
|
/// \copydoc NewtonIterationBlackoilInterface::parallelInformation
|
|
const boost::any& parallelInformation() const { return parallelInformation_; }
|
|
|
|
public:
|
|
/// \brief construct the CPR preconditioner and the solver.
|
|
/// \tparam P The type of the parallel information.
|
|
/// \param parallelInformation the information about the parallelization.
|
|
template<int category=Dune::SolverCategory::sequential, class LinearOperator, class POrComm>
|
|
void constructPreconditionerAndSolve(LinearOperator& linearOperator,
|
|
Vector& x, Vector& istlb,
|
|
const POrComm& parallelInformation_arg,
|
|
Dune::InverseOperatorResult& result) const
|
|
{
|
|
// Construct scalar product.
|
|
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
|
|
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
|
|
SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
|
|
|
|
// Communicate if parallel.
|
|
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
|
|
|
|
#if FLOW_SUPPORT_AMG // activate AMG if either flow_ebos is used or UMFPack is not available
|
|
if( parameters_.linear_solver_use_amg_ )
|
|
{
|
|
typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType;
|
|
typedef typename CPRSelectorType::AMG AMG;
|
|
typedef typename CPRSelectorType::Operator MatrixOperator;
|
|
|
|
std::unique_ptr< AMG > amg;
|
|
std::unique_ptr< MatrixOperator > opA;
|
|
|
|
if( ! std::is_same< LinearOperator, MatrixOperator > :: value )
|
|
{
|
|
// create new operator in case linear operator and matrix operator differ
|
|
opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) );
|
|
}
|
|
|
|
const double relax = 1.0;
|
|
|
|
// Construct preconditioner.
|
|
constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
|
|
|
|
// Solve.
|
|
solve(linearOperator, x, istlb, *sp, *amg, result);
|
|
}
|
|
else
|
|
#endif
|
|
{
|
|
// Construct preconditioner.
|
|
auto precond = constructPrecond(linearOperator, parallelInformation_arg);
|
|
|
|
// Solve.
|
|
solve(linearOperator, x, istlb, *sp, *precond, result);
|
|
}
|
|
}
|
|
|
|
#if DUNE_VERSION_NEWER_REV(DUNE_ISTL, 2 , 5, 1)
|
|
// 3x3 matrix block inversion was unstable at least 2.3 until and including
|
|
// 2.5.0
|
|
typedef ParallelOverlappingILU0<Matrix,Vector,Vector> SeqPreconditioner;
|
|
#else
|
|
typedef ParallelOverlappingILU0<Dune::BCRSMatrix<Dune::MatrixBlock<typename Matrix::field_type,
|
|
Matrix::block_type::rows,
|
|
Matrix::block_type::cols> >,
|
|
Vector, Vector> SeqPreconditioner;
|
|
#endif
|
|
|
|
template <class Operator>
|
|
std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
|
|
{
|
|
const double relax = parameters_.ilu_relaxation_;
|
|
const int ilu_fillin = parameters_.ilu_fillin_level_;
|
|
std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax));
|
|
return precond;
|
|
}
|
|
|
|
#if HAVE_MPI
|
|
typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
|
|
#if DUNE_VERSION_NEWER_REV(DUNE_ISTL, 2 , 5, 1)
|
|
// 3x3 matrix block inversion was unstable from at least 2.3 until and
|
|
// including 2.5.0
|
|
typedef ParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
|
|
#else
|
|
typedef ParallelOverlappingILU0<Dune::BCRSMatrix<Dune::MatrixBlock<typename Matrix::field_type,
|
|
Matrix::block_type::rows,
|
|
Matrix::block_type::cols> >,
|
|
Vector, Vector, Comm> ParPreconditioner;
|
|
#endif
|
|
template <class Operator>
|
|
std::unique_ptr<ParPreconditioner>
|
|
constructPrecond(Operator& opA, const Comm& comm) const
|
|
{
|
|
typedef std::unique_ptr<ParPreconditioner> Pointer;
|
|
const double relax = parameters_.ilu_relaxation_;
|
|
return Pointer(new ParPreconditioner(opA.getmat(), comm, relax));
|
|
}
|
|
#endif
|
|
|
|
template <class LinearOperator, class MatrixOperator, class POrComm, class AMG >
|
|
void
|
|
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const
|
|
{
|
|
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, comm, amg );
|
|
}
|
|
|
|
|
|
template <class MatrixOperator, class POrComm, class AMG >
|
|
void
|
|
constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const
|
|
{
|
|
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax, comm, amg );
|
|
}
|
|
|
|
/// \brief Solve the system using the given preconditioner and scalar product.
|
|
template <class Operator, class ScalarProd, class Precond>
|
|
void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const
|
|
{
|
|
// TODO: Revise when linear solvers interface opm-core is done
|
|
// Construct linear solver.
|
|
// GMRes solver
|
|
int verbosity = ( isIORank_ ) ? parameters_.linear_solver_verbosity_ : 0;
|
|
|
|
if ( parameters_.newton_use_gmres_ ) {
|
|
Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond,
|
|
parameters_.linear_solver_reduction_,
|
|
parameters_.linear_solver_restart_,
|
|
parameters_.linear_solver_maxiter_,
|
|
verbosity);
|
|
// Solve system.
|
|
linsolve.apply(x, istlb, result);
|
|
}
|
|
else { // BiCGstab solver
|
|
Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond,
|
|
parameters_.linear_solver_reduction_,
|
|
parameters_.linear_solver_maxiter_,
|
|
verbosity);
|
|
// Solve system.
|
|
linsolve.apply(x, istlb, result);
|
|
}
|
|
}
|
|
|
|
|
|
/// Solve the linear system Ax = b, with A being the
|
|
/// combined derivative matrix of the residual and b
|
|
/// being the residual itself.
|
|
/// \param[in] A matrix A
|
|
/// \param[inout] x solution to be computed x
|
|
/// \param[in] b right hand side b
|
|
void solve(Matrix& A, Vector& x, Vector& b ) const
|
|
{
|
|
// Parallel version is deactivated until we figure out how to do it properly.
|
|
#if HAVE_MPI
|
|
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
|
|
{
|
|
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
|
|
const ParallelISTLInformation& info =
|
|
boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
|
|
Comm istlComm(info.communicator());
|
|
|
|
// Construct operator, scalar product and vectors needed.
|
|
typedef Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector,Comm> Operator;
|
|
Operator opA(A, istlComm);
|
|
solve( opA, x, b, istlComm );
|
|
}
|
|
else
|
|
#endif
|
|
{
|
|
// Construct operator, scalar product and vectors needed.
|
|
Dune::MatrixAdapter< Matrix, Vector, Vector> opA( A );
|
|
solve( opA, x, b );
|
|
}
|
|
}
|
|
|
|
/// Solve the linear system Ax = b, with A being the
|
|
/// combined derivative matrix of the residual and b
|
|
/// being the residual itself.
|
|
/// \param[in] A matrix A
|
|
/// \param[inout] x solution to be computed x
|
|
/// \param[in] b right hand side b
|
|
template <class Operator, class Comm >
|
|
void solve(Operator& opA, Vector& x, Vector& b, Comm& comm) const
|
|
{
|
|
Dune::InverseOperatorResult result;
|
|
// Parallel version is deactivated until we figure out how to do it properly.
|
|
#if HAVE_MPI
|
|
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
|
|
{
|
|
const size_t size = opA.getmat().N();
|
|
const ParallelISTLInformation& info =
|
|
boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
|
|
|
|
// As we use a dune-istl with block size np the number of components
|
|
// per parallel is only one.
|
|
info.copyValuesTo(comm.indexSet(), comm.remoteIndices(),
|
|
size, 1);
|
|
// Construct operator, scalar product and vectors needed.
|
|
constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, b, comm, result);
|
|
}
|
|
else
|
|
#endif
|
|
{
|
|
OPM_THROW(std::logic_error,"this method if for parallel solve only");
|
|
}
|
|
|
|
checkConvergence( result );
|
|
}
|
|
|
|
/// Solve the linear system Ax = b, with A being the
|
|
/// combined derivative matrix of the residual and b
|
|
/// being the residual itself.
|
|
/// \param[in] A matrix A
|
|
/// \param[inout] x solution to be computed x
|
|
/// \param[in] b right hand side b
|
|
template <class Operator>
|
|
void solve(Operator& opA, Vector& x, Vector& b ) const
|
|
{
|
|
Dune::InverseOperatorResult result;
|
|
// Construct operator, scalar product and vectors needed.
|
|
Dune::Amg::SequentialInformation info;
|
|
constructPreconditionerAndSolve(opA, x, b, info, result);
|
|
checkConvergence( result );
|
|
}
|
|
|
|
void checkConvergence( const Dune::InverseOperatorResult& result ) const
|
|
{
|
|
// store number of iterations
|
|
iterations_ = result.iterations;
|
|
|
|
// Check for failure of linear solver.
|
|
if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
|
|
const std::string msg("Convergence failure for linear solver.");
|
|
OPM_THROW_NOLOG(LinearSolverProblem, msg);
|
|
}
|
|
}
|
|
protected:
|
|
mutable int iterations_;
|
|
boost::any parallelInformation_;
|
|
bool isIORank_;
|
|
|
|
NewtonIterationBlackoilInterleavedParameters parameters_;
|
|
}; // end ISTLSolver
|
|
|
|
} // namespace Opm
|
|
#endif
|