OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Allow user to choose modified ILU0 decomposition.

Browse Source 
Using the parameter ilu_milu=true|false (default=false) the user can now choose
to use the modified ILU0 decomposition. If selected values will
not be dropped for nonzero entries but added to the diagonal of U.
This approach will result in A*e = L*U*e for vector e with all entries
beging 1.
This commit is contained in:
Markus Blatt
2018-05-17 12:01:02 +02:00
parent 5333a2430f
commit 0bae240a42
4 changed files with 146 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files

59
opm/autodiff/CPRPreconditioner.hpp
View File

@@ -64,6 +64,7 @@ struct CPRParameter
double cpr_relax_;
double cpr_solver_tol_;
int cpr_ilu_n_;
bool cpr_ilu_milu_;
int cpr_max_ell_iter_;
bool cpr_use_amg_;
bool cpr_use_bicgstab_;
@@ -80,6 +81,7 @@ struct CPRParameter
cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_);
cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_);
cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_);
cpr_ilu_milu_ = param.getDefault("ilu_milu", cpr_ilu_milu_);
cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_);
cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_);
cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_);
@@ -92,6 +94,7 @@ struct CPRParameter
cpr_relax_ = 1.0;
cpr_solver_tol_ = 1e-2;
cpr_ilu_n_ = 0;
cpr_ilu_milu_ = false;
cpr_max_ell_iter_ = 25;
cpr_use_amg_ = true;
cpr_use_bicgstab_ = true;
@@ -103,6 +106,31 @@ struct CPRParameter
namespace ISTLUtility
{
template<class T>
void setILUParameters(Opm::ParallelOverlappingILU0Args<T>& args,
const CPRParameter& params)
{
args.setN(params.cpr_ilu_n_);
args.setMilu(params.cpr_ilu_milu_);
}
template<class T>
void setILUParameters(Opm::ParallelOverlappingILU0Args<T>& args,
bool milu, int n=0)
{
args.setN(n);
args.setMilu(milu);
}
template<class S, class P>
void setILUParameters(S&, const P&)
{}
template<class S, class P>
void setILUParameters(S&, bool, int)
{}
///
/// \brief A traits class for selecting the types of the preconditioner.
///
@@ -167,9 +195,9 @@ struct CPRSelector<M,X,Y,Dune::Amg::SequentialInformation>
//! \param relax The relaxation factor to use.
template<class M, class X, class C>
std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
createILU0Ptr(const M& A, const C& comm, double relax)
createILU0Ptr(const M& A, const C& comm, double relax, bool milu)
{
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >(A, comm, relax);
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >(A, comm, relax, milu);
}
//! \brief Creates and initializes a shared pointer to an ILUn preconditioner.
//! \param A The matrix of the linear system to solve.
@@ -177,25 +205,28 @@ createILU0Ptr(const M& A, const C& comm, double relax)
//! \param relax The relaxation factor to use.
template<class M, class X, class C>
std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
createILUnPtr(const M& A, const C& comm, int ilu_n, double relax)
createILUnPtr(const M& A, const C& comm, int ilu_n, double relax, bool milu)
{
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >( A, comm, ilu_n, relax);
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >( A, comm, ilu_n, relax, milu );
}
/// \brief Creates the elliptic preconditioner (ILU0)
/// \param Ae The matrix of the elliptic system.
/// \param relax The relaxation parameter for ILU0.
/// \param milu If true, the modified ilu approach is used. Dropped elements
/// will get added to the diagonal of U to preserve the row sum
/// for constant vectors (Ae = LUe).
/// \param comm The object describing the parallelization information and communication.
template<class M, class X=typename M::range_type, class P>
typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
createEllipticPreconditionerPointer(const M& Ae, double relax,
const P& comm)
bool milu, const P& comm)
{
typedef typename CPRSelector<M,X,X,P >
::EllipticPreconditioner ParallelPreconditioner;
typedef typename CPRSelector<M,X,X,P>
::EllipticPreconditionerPointer EllipticPreconditionerPointer;
return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax));
return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax, milu));
}
template < class C, class Op, class P, class S, std::size_t index >
@@ -220,13 +251,14 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm,
SmootherArgs smootherArgs;
smootherArgs.iterations = 1;
smootherArgs.relaxationFactor = relax;
setILUParameters(smootherArgs, params);
amgPtr.reset( new AMG( params, opA, criterion, smootherArgs, comm ) );
}
template < class C, class Op, class P, class AMG >
inline void
createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::unique_ptr< AMG >& amgPtr)
createAMGPreconditionerPointer(Op& opA, const double relax, const bool milu, const P& comm, std::unique_ptr< AMG >& amgPtr)
{
// TODO: revise choice of parameters
int coarsenTarget=1200;
@@ -243,6 +275,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::
SmootherArgs smootherArgs;
smootherArgs.iterations = 1;
smootherArgs.relaxationFactor = relax;
setILUParameters(smootherArgs, milu);
amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) );
}
@@ -254,7 +287,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::
// \param amgPtr The unique_ptr to be filled (return)
template < int pressureIndex=0, class Op, class P, class AMG >
inline void
createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std::unique_ptr< AMG >& amgPtr )
createAMGPreconditionerPointer( Op& opA, const double relax, const bool milu, const P& comm, std::unique_ptr< AMG >& amgPtr )
{
// type of matrix
typedef typename Op::matrix_type M;
@@ -268,7 +301,7 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
// The coarsening criterion used in the AMG
typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
createAMGPreconditionerPointer<Criterion>(opA, relax, comm, amgPtr);
createAMGPreconditionerPointer<Criterion>(opA, relax, milu, comm, amgPtr);
}
} // end namespace ISTLUtility
@@ -370,10 +403,10 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
// create the preconditioner for the whole system.
if( param_.cpr_ilu_n_ == 0 ) {
pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, comm, param_.cpr_relax_ );
pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, comm, param_.cpr_relax_, param_.cpr_ilu_milu_ );
}
else {
pre_ = ISTLUtility::createILUnPtr<M,X>( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_);
pre_ = ISTLUtility::createILUnPtr<M,X>( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_, param_.cpr_ilu_milu_);
}
}
@@ -531,11 +564,11 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
{
if( amg )
{
ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, comm, amg_ );
ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, param_.cpr_ilu_milu_, comm, amg_ );
}
else
{
precond_ = ISTLUtility::createEllipticPreconditionerPointer<M,X>( Ae_, param_.cpr_relax_, comm);
precond_ = ISTLUtility::createEllipticPreconditionerPointer<M,X>( Ae_, param_.cpr_relax_, param_.cpr_ilu_milu_, comm);
}
}
};

32
opm/autodiff/ISTLSolver.hpp
View File

@@ -440,6 +440,7 @@ namespace Detail
}
const double relax = parameters_.ilu_relaxation_;
const bool ilu_milu = parameters_.ilu_milu_;
if ( parameters_.use_cpr_ )
{
using Matrix = typename MatrixOperator::matrix_type;
@@ -452,7 +453,7 @@ namespace Detail
std::unique_ptr< AMG > amg;
// Construct preconditioner.
Criterion crit(15, 2000);
constructAMGPrecond<Criterion>( linearOperator, parallelInformation_arg, amg, opA, relax );
constructAMGPrecond<Criterion>( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu );
// Solve.
solve(linearOperator, x, istlb, *sp, *amg, result);
@@ -463,7 +464,7 @@ namespace Detail
std::unique_ptr< AMG > amg;
// Construct preconditioner.
constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu );
// Solve.
solve(linearOperator, x, istlb, *sp, *amg, result);
@@ -495,7 +496,8 @@ namespace Detail
{
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));
const bool ilu_milu = parameters_.ilu_milu_;
std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax, ilu_milu));
return precond;
}
@@ -517,38 +519,44 @@ namespace Detail
{
typedef std::unique_ptr<ParPreconditioner> Pointer;
const double relax = parameters_.ilu_relaxation_;
return Pointer(new ParPreconditioner(opA.getmat(), comm, relax));
const bool ilu_milu = parameters_.ilu_milu_;
return Pointer(new ParPreconditioner(opA.getmat(), comm, relax, ilu_milu));
}
#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
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, const bool milu) const
{
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, comm, amg );
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, milu, 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
constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax,
const bool milu) const
{
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax, comm, amg );
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax,
milu, comm, amg );
}
template <class C, 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
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax,
const bool milu ) const
{
ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax, comm, amg, parameters_ );
ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax,
comm, amg, parameters_ );
}
template <class C, 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
constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax, const bool milu ) const
{
ISTLUtility::template createAMGPreconditionerPointer<C>( opA, relax, comm, amg, parameters_ );
ISTLUtility::template createAMGPreconditionerPointer<C>( opA, relax, milu,
comm, amg, parameters_ );
}
/// \brief Solve the system using the given preconditioner and scalar product.
template <class Operator, class ScalarProd, class Precond>

3
opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
View File

@@ -44,6 +44,7 @@ namespace Opm
int linear_solver_restart_;
int linear_solver_verbosity_;
int ilu_fillin_level_;
bool ilu_milu_;
bool newton_use_gmres_;
bool require_full_sparsity_pattern_;
bool ignoreConvergenceFailure_;
@@ -69,6 +70,7 @@ namespace Opm
linear_solver_use_amg_ = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ );
ilu_relaxation_ = param.getDefault("ilu_relaxation", ilu_relaxation_ );
ilu_fillin_level_ = param.getDefault("ilu_fillin_level", ilu_fillin_level_ );
ilu_milu_ = param.getDefault("ilu_milu", ilu_milu_);
// Check whether to use cpr approach
const std::string cprSolver = "cpr";
@@ -89,6 +91,7 @@ namespace Opm
linear_solver_use_amg_ = false;
ilu_fillin_level_ = 0;
ilu_relaxation_ = 0.9;
ilu_milu_ = false;
}
};

93
opm/autodiff/ParallelOverlappingILU0.hpp
View File

@@ -37,6 +37,34 @@ namespace Opm
template<class Matrix, class Domain, class Range, class ParallelInfo = Dune::Amg::SequentialInformation>
class ParallelOverlappingILU0;
template<class F>
class ParallelOverlappingILU0Args
: public Dune::Amg::DefaultSmootherArgs<F>
{
public:
ParallelOverlappingILU0Args(bool milu = false)
: milu_(milu)
{}
void setMilu(bool milu)
{
milu_ = milu;
}
bool getMilu() const
{
return milu_;
}
void setN(int n)
{
n_ = n;
}
int getN() const
{
return n_;
}
private:
bool milu_;
int n_;
};
} // end namespace Opm
namespace Dune
@@ -45,6 +73,13 @@ namespace Dune
namespace Amg
{
template<class M, class X, class Y, class C>
struct SmootherTraits<Opm::ParallelOverlappingILU0<M,X,Y,C> >
{
using Arguments = Opm::ParallelOverlappingILU0Args<typename M::field_type>;
};
/// \brief Tells AMG how to construct the Opm::ParallelOverlappingILU0 smoother
/// \tparam Matrix The type of the Matrix.
/// \tparam Domain The type of the Vector representing the domain.
@@ -54,17 +89,18 @@ namespace Amg
template<class Matrix, class Domain, class Range, class ParallelInfo>
struct ConstructionTraits<Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> >
{
typedef Dune::SeqILU0<Matrix,Domain,Range> T;
typedef Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> T;
typedef DefaultParallelConstructionArgs<T,ParallelInfo> Arguments;
typedef ConstructionTraits<T> SeqConstructionTraits;
static inline Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* construct(Arguments& args)
{
return new Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>(args.getMatrix(),
args.getComm(),
args.getArgs().relaxationFactor);
return new T(args.getMatrix(),
args.getComm(),
args.getArgs().getN(),
args.getArgs().relaxationFactor,
args.getArgs().getMilu());
}
static inline void deconstruct(Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* bp)
static inline void deconstruct(T* bp)
{
delete bp;
}
@@ -330,10 +366,14 @@ public:
\param A The matrix to operate on.
\param n ILU fill in level (for testing). This does not work in parallel.
\param w The relaxation factor.
\param milu If true, the modified ilu approach is used. Dropped elements
will get added to the diagonal of U to preserve the row sum
for constant vectors (Ae = LUe).
*/
template<class BlockType, class Alloc>
ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const int n, const field_type w )
const int n, const field_type w,
bool milu)
: lower_(),
upper_(),
inv_(),
@@ -342,7 +382,7 @@ public:
{
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), n );
init( reinterpret_cast<const Matrix&>(A), n, milu );
}
/*! \brief Constructor gets all parameters to operate the prec.
@@ -350,10 +390,14 @@ public:
\param comm communication object, e.g. Dune::OwnerOverlapCopyCommunication
\param n ILU fill in level (for testing). This does not work in parallel.
\param w The relaxation factor.
\param milu If true, the modified ilu approach is used. Dropped elements
will get added to the diagonal of U to preserve the row sum
for constant vectors (Ae = LUe).
*/
template<class BlockType, class Alloc>
ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const ParallelInfo& comm, const int n, const field_type w )
const ParallelInfo& comm, const int n, const field_type w,
bool milu)
: lower_(),
upper_(),
inv_(),
@@ -362,7 +406,7 @@ public:
{
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), n );
init( reinterpret_cast<const Matrix&>(A), n, milu );
}
/*! \brief Constructor.
@@ -370,11 +414,14 @@ public:
Constructor gets all parameters to operate the prec.
\param A The matrix to operate on.
\param w The relaxation factor.
\param milu If true, the modified ilu approach is used. Dropped elements
will get added to the diagonal of U to preserve the row sum
for constant vectors (Ae = LUe).
*/
template<class BlockType, class Alloc>
ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const field_type w)
: ParallelOverlappingILU0( A, 0, w )
const field_type w, bool milu)
: ParallelOverlappingILU0( A, 0, w, milu )
{
}
@@ -384,10 +431,14 @@ public:
\param A The matrix to operate on.
\param comm communication object, e.g. Dune::OwnerOverlapCopyCommunication
\param w The relaxation factor.
\param milu If true, the modified ilu approach is used. Dropped elements
will get added to the diagonal of U to preserve the row sum
for constant vectors (Ae = LUe).
*/
template<class BlockType, class Alloc>
ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const ParallelInfo& comm, const field_type w)
const ParallelInfo& comm, const field_type w,
bool milu)
: lower_(),
upper_(),
inv_(),
@@ -396,7 +447,7 @@ public:
{
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), 0 );
init( reinterpret_cast<const Matrix&>(A), 0, milu );
}
/*!
@@ -490,7 +541,7 @@ public:
}
protected:
void init( const Matrix& A, const int iluIteration )
void init( const Matrix& A, const int iluIteration, bool milu )
{
// (For older DUNE versions the communicator might be
// invalid if redistribution in AMG happened on the coarset level.
@@ -519,11 +570,21 @@ protected:
if( iluIteration == 0 ) {
// create ILU-0 decomposition
ILU.reset( new Matrix( A ) );
bilu0_decomposition( *ILU );
if ( milu )
{
detail::milu0_decomposition ( *ILU );
}else
{
bilu0_decomposition( *ILU );
}
}
else {
// create ILU-n decomposition
ILU.reset( new Matrix( A.N(), A.M(), Matrix::row_wise) );
if ( milu )
{
OpmLog::warning("MILU_N_NOT_SUPPORTED", "MILU variant only supported for zero fill-in");
}
bilu_decomposition( A, iluIteration, *ILU );
}
}