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

Merge pull request #1757 from atgeirr/improve_cpr_mod

Browse Source 
Improve CPR solver
This commit is contained in:
Atgeirr Flø Rasmussen
2019-03-21 11:19:10 +01:00
committed by GitHub
parent a7362abb04 19b78b78c0
commit f18477ce71
7 changed files with 470 additions and 146 deletions
Download Patch File Download Diff File Expand all files Collapse all files

180
opm/autodiff/BlackoilAmg.hpp
View File

@@ -21,6 +21,7 @@
#include <ewoms/linear/matrixblock.hh>
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
#include <opm/autodiff/FlowLinearSolverParameters.hpp>
#include <opm/autodiff/CPRPreconditioner.hpp>
#include <dune/istl/paamg/twolevelmethod.hh>
#include <dune/istl/paamg/aggregates.hh>
@@ -41,6 +42,27 @@ class UnSymmetricCriterion;
}
}
namespace Opm
{
namespace Amg
{
template<int Row, int Column>
class Element
{
public:
enum { /* @brief We preserve the sign.*/
is_sign_preserving = true
};
template<class M>
typename M::field_type operator()(const M& m) const
{
return m[Row][Column];
}
};
} // namespace Amg
} // namespace Opm
namespace Dune
{
@@ -67,6 +89,17 @@ Dune::MatrixAdapter<M,X,Y> createOperator(const Dune::MatrixAdapter<M,X,Y>&, con
return Dune::MatrixAdapter<M,X,Y>(matrix);
}
/**
* \brief Creates a MatrixAdapter as an operator, storing it in a unique_ptr.
*
* The first argument is used to specify the return type using function overloading.
* \param matrix The matrix to wrap.
*/
template<class M, class X, class Y, class T>
std::unique_ptr< Dune::MatrixAdapter<M,X,Y> > createOperatorPtr(const Dune::MatrixAdapter<M,X,Y>&, const M& matrix, const T&)
{
return std::make_unique< Dune::MatrixAdapter<M,X,Y> >(matrix);
}
/**
* \brief Creates an OverlappingSchwarzOperator as an operator.
*
@@ -87,30 +120,27 @@ Dune::OverlappingSchwarzOperator<M,X,Y,T> createOperator(const Dune::Overlapping
//! Sedimentary Basin Simulations, 2003.
//! \param op The operator that stems from the discretization.
//! \param comm The communication objecte describing the data distribution.
//! \param pressureIndex The index of the pressure in the matrix block
//! \param pressureEqnIndex The index of the pressure in the matrix block
//! \retun A pair of the scaled matrix and the associated operator-
template<class Operator, class Communication>
template<class Operator, class Communication, class Vector>
std::tuple<std::unique_ptr<typename Operator::matrix_type>, Operator>
scaleMatrixQuasiImpes(const Operator& op, const Communication& comm,
std::size_t pressureIndex)
scaleMatrixDRS(const Operator& op, const Communication& comm,
std::size_t pressureEqnIndex, const Vector& weights, const Opm::CPRParameter& param)
{
using Matrix = typename Operator::matrix_type;
using Block = typename Matrix::block_type;
using BlockVector = typename Vector::block_type;
std::unique_ptr<Matrix> matrix(new Matrix(op.getmat()));
for ( auto& row : *matrix )
{
for ( auto& block : row )
{
for ( std::size_t i = 0; i < Block::rows; i++ )
{
if ( i != pressureIndex )
{
for(std::size_t j=0; j < Block::cols; j++)
{
block[pressureIndex][j] += block[i][j];
}
}
if (param.cpr_use_drs_) {
const auto endi = matrix->end();
for (auto i = matrix->begin(); i != endi; ++i) {
const BlockVector& bw = weights[i.index()];
const auto endj = (*i).end();
for (auto j = (*i).begin(); j != endj; ++j) {
Block& block = *j;
BlockVector& bvec = block[pressureEqnIndex];
// should introduce limits which also change the weights
block.mtv(bw, bvec);
}
}
}
@@ -122,20 +152,16 @@ scaleMatrixQuasiImpes(const Operator& op, const Communication& comm,
//! See section 3.2.3 of Scheichl, Masson: Decoupling and Block Preconditioning for
//! Sedimentary Basin Simulations, 2003.
//! \param vector The vector to scale
//! \param pressureIndex The index of the pressure in the matrix block
//! \param pressureEqnIndex The index of the pressure in the matrix block
template<class Vector>
void scaleVectorQuasiImpes(Vector& vector, std::size_t pressureIndex)
void scaleVectorDRS(Vector& vector, std::size_t pressureEqnIndex, const Opm::CPRParameter& param, const Vector& weights)
{
using Block = typename Vector::block_type;
for ( auto& block: vector)
{
for ( std::size_t i = 0; i < Block::dimension; i++ )
{
if ( i != pressureIndex )
{
block[pressureIndex] += block[i];
}
if (param.cpr_use_drs_) {
for (std::size_t j = 0; j < vector.size(); ++j) {
Block& block = vector[j];
const Block& bw = weights[j];
block[pressureEqnIndex] = bw.dot(block);
}
}
}
@@ -284,23 +310,23 @@ struct ScalarType<Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Du
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<typename ScalarType<B>::value>, Dune::Amg::FirstDiagonal> >;
};
template<class C, std::size_t COMPONENT_INDEX>
template<class C, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
struct OneComponentCriterionType
{};
template<class B, class N, std::size_t COMPONENT_INDEX>
struct OneComponentCriterionType<Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<B>,N> >,COMPONENT_INDEX>
template<class B, class N, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
struct OneComponentCriterionType<Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<B>,N> >, COMPONENT_INDEX, VARIABLE_INDEX>
{
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<B>, Dune::Amg::Diagonal<COMPONENT_INDEX> > >;
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<B>, Opm::Amg::Element<COMPONENT_INDEX, VARIABLE_INDEX> > >;
};
template<class B, class N, std::size_t COMPONENT_INDEX>
struct OneComponentCriterionType<Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<B>,N> >,COMPONENT_INDEX>
template<class B, class N, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
struct OneComponentCriterionType<Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<B>,N> >, COMPONENT_INDEX, VARIABLE_INDEX>
{
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<B>, Dune::Amg::Diagonal<COMPONENT_INDEX> > >;
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<B>, Opm::Amg::Element<COMPONENT_INDEX, VARIABLE_INDEX> > >;
};
template<class Operator, class Criterion, class Communication, std::size_t COMPONENT_INDEX>
template<class Operator, class Criterion, class Communication, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
class OneComponentAggregationLevelTransferPolicy;
@@ -401,9 +427,12 @@ private:
// Linear solver parameters
const double tolerance = param_->cpr_solver_tol_;
const int maxit = param_->cpr_max_ell_iter_;
const int verbosity = ( param_->cpr_solver_verbose_ &&
comm_.communicator().rank()==0 ) ? 1 : 0;
if ( param_->cpr_use_bicgstab_ )
int verbosity = 0;
if (comm_.communicator().rank() == 0) {
verbosity = param_->cpr_solver_verbose_;
}
if ( param_->cpr_ell_solvetype_ == 0)
{
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
Dune::BiCGSTABSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp, *prec,
@@ -428,7 +457,7 @@ private:
}
#endif
}
else
else if (param_->cpr_ell_solvetype_ == 1)
{
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
Dune::CGSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp, *prec,
@@ -453,6 +482,36 @@ private:
}
#endif
}
else
{
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
Dune::LoopSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp, *prec,
tolerance, maxit, verbosity);
solver.apply(x,b,res);
#else
if ( !amg_ )
{
Dune::LoopSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp,
reinterpret_cast<Smoother&>(*prec),
tolerance, maxit, verbosity);
solver.apply(x,b,res);
}
else
{
Dune::LoopSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp,
reinterpret_cast<AMGType&>(*prec),
tolerance, maxit, verbosity);
solver.apply(x,b,res);
}
#endif
}
// Warn if unknown options.
if (param_->cpr_ell_solvetype_ > 2 && comm_.communicator().rank() == 0) {
OpmLog::warning("cpr_ell_solver_type_unknown", "Unknown CPR elliptic solver type specification, using LoopSolver.");
}
#if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
delete sp;
@@ -645,7 +704,7 @@ void buildCoarseSparseMatrix(M& coarseMatrix, G& fineGraph, const V& visitedMap,
* @tparam Criterion The criterion that describes the aggregation procedure.
* @tparam Communication The class that describes the communication pattern.
*/
template<class Operator, class Criterion, class Communication, std::size_t COMPONENT_INDEX>
template<class Operator, class Criterion, class Communication, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
class OneComponentAggregationLevelTransferPolicy
: public Dune::Amg::LevelTransferPolicy<Operator, typename Detail::ScalarType<Operator>::value>
{
@@ -756,7 +815,7 @@ public:
for ( auto col = row.begin(), cend = row.end(); col != cend; ++col, ++coarseCol )
{
assert( col.index() == coarseCol.index() );
*coarseCol = (*col)[COMPONENT_INDEX][COMPONENT_INDEX];
*coarseCol = (*col)[COMPONENT_INDEX][VARIABLE_INDEX];
}
++coarseRow;
}
@@ -786,7 +845,7 @@ public:
const auto& j = (*aggregatesMap_)[entry.index()];
if ( j != AggregatesMap::ISOLATED )
{
(*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][COMPONENT_INDEX];
(*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][VARIABLE_INDEX];
}
}
}
@@ -873,18 +932,20 @@ private:
* \brief An algebraic twolevel or multigrid approach for solving blackoil (supports CPR with and without AMG)
*
* This preconditioner first decouples the component used for coarsening using a simple scaling
* approach (e.g. Scheichl, Masson 2013,\see scaleMatrixQuasiImpes). Then it constructs the first
* coarse level system, either by simply extracting the coupling between the components at COMPONENT_INDEX
* in the matrix blocks or by extracting them and applying aggregation to them directly. This coarse level
* approach (e.g. Scheichl, Masson 2013,\see scaleMatrixDRS). Then it constructs the
* coarse level system. The coupling is defined by the weights corresponding to the element located at
* (COMPONENT_INDEX, VARIABLE_INDEX) in the block matrix. Then the coarse level system is constructed
* either by extracting these elements, or by applying aggregation to them directly. This coarse level
* can be solved either by AMG or by ILU. The preconditioner is configured using CPRParameter.
* \tparam O The type of the operator (encapsulating a BCRSMatrix).
* \tparam S The type of the smoother.
* \tparam C The type of coarsening criterion to use.
* \tparam P The type of the class describing the parallelization.
* \tparam COMPONENT_INDEX The index of the component to use for coarsening (usually the pressure).
* \tparam COMPONENT_INDEX The index of the component to use for coarsening (usually water).
* \tparam VARIABLE_INDEX The index of the variable to use for coarsening (usually pressure).
*/
template<typename O, typename S, typename C,
typename P, std::size_t COMPONENT_INDEX>
typename P, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
class BlackoilAmg
: public Dune::Preconditioner<typename O::domain_type, typename O::range_type>
{
@@ -905,13 +966,14 @@ protected:
using CoarseOperator = typename Detail::ScalarType<Operator>::value;
using CoarseSmoother = typename Detail::ScalarType<Smoother>::value;
using FineCriterion =
typename Detail::OneComponentCriterionType<Criterion,COMPONENT_INDEX>::value;
typename Detail::OneComponentCriterionType<Criterion, COMPONENT_INDEX, VARIABLE_INDEX>::value;
using CoarseCriterion = typename Detail::ScalarType<Criterion>::value;
using LevelTransferPolicy =
OneComponentAggregationLevelTransferPolicy<Operator,
FineCriterion,
Communication,
COMPONENT_INDEX>;
COMPONENT_INDEX,
VARIABLE_INDEX>;
using CoarseSolverPolicy =
Detail::OneStepAMGCoarseSolverPolicy<CoarseOperator,
CoarseSmoother,
@@ -944,18 +1006,19 @@ public:
* \param comm The information about the parallelization.
*/
BlackoilAmg(const CPRParameter& param,
const typename TwoLevelMethod::FineDomainType& weights,
const Operator& fineOperator, const Criterion& criterion,
const SmootherArgs& smargs, const Communication& comm)
: param_(param),
scaledMatrixOperator_(Detail::scaleMatrixQuasiImpes(fineOperator, comm,
COMPONENT_INDEX)),
weights_(weights),
scaledMatrixOperator_(Detail::scaleMatrixDRS(fineOperator, comm,
COMPONENT_INDEX, weights, param)),
smoother_(Detail::constructSmoother<Smoother>(std::get<1>(scaledMatrixOperator_),
smargs, comm)),
levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_),
coarseSolverPolicy_(&param, smargs, criterion),
twoLevelMethod_(std::get<1>(scaledMatrixOperator_), smoother_,
levelTransferPolicy_,
coarseSolverPolicy_, 0, 1)
levelTransferPolicy_, coarseSolverPolicy_, 0, 1)
{}
void pre(typename TwoLevelMethod::FineDomainType& x,
@@ -973,11 +1036,12 @@ public:
const typename TwoLevelMethod::FineRangeType& d)
{
auto scaledD = d;
Detail::scaleVectorQuasiImpes(scaledD, COMPONENT_INDEX);
Detail::scaleVectorDRS(scaledD, COMPONENT_INDEX, param_, weights_);
twoLevelMethod_.apply(v, scaledD);
}
private:
const CPRParameter& param_;
const typename TwoLevelMethod::FineDomainType& weights_;
std::tuple<std::unique_ptr<Matrix>, Operator> scaledMatrixOperator_;
std::shared_ptr<Smoother> smoother_;
LevelTransferPolicy levelTransferPolicy_;
@@ -997,7 +1061,7 @@ namespace ISTLUtility
/// \tparam C The type of the coarsening criterion to use.
/// \tparam index The pressure index.
////
template<class M, class X, class Y, class P, class C, std::size_t index>
template<class M, class X, class Y, class P, class C, std::size_t pressureEqnIndex, std::size_t pressureVarIndex>
struct BlackoilAmgSelector
{
using Criterion = C;
@@ -1005,7 +1069,7 @@ struct BlackoilAmgSelector
using ParallelInformation = typename Selector::ParallelInformation;
using Operator = typename Selector::Operator;
using Smoother = typename Selector::EllipticPreconditioner;
using AMG = BlackoilAmg<Operator,Smoother,Criterion,ParallelInformation,index>;
using AMG = BlackoilAmg<Operator, Smoother, Criterion, ParallelInformation, pressureEqnIndex, pressureVarIndex>;
};
} // end namespace ISTLUtility
} // end namespace Opm

29
opm/autodiff/CPRPreconditioner.hpp
View File

@@ -55,9 +55,15 @@ namespace Opm
{
template<typename O, typename S, typename C,
typename P, std::size_t COMPONENT_INDEX>
typename P, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
class BlackoilAmg;
namespace Amg
{
template<int Row, int Column>
class Element;
}
namespace ISTLUtility
{
@@ -183,18 +189,23 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax, milu));
}
template < class C, class Op, class P, class S, std::size_t index >
template < class C, class Op, class P, class S, std::size_t PressureEqnIndex, std::size_t PressureVarIndex, class Vector>
inline void
createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm,
std::unique_ptr< BlackoilAmg<Op,S,C,P,index> >& amgPtr,
const CPRParameter& params)
std::unique_ptr< BlackoilAmg<Op,S,C,P,PressureEqnIndex,PressureVarIndex> >& amgPtr,
const CPRParameter& params,
const Vector& weights)
{
using AMG = BlackoilAmg<Op,S,C,P,index>;
using AMG = BlackoilAmg<Op,S,C,P,PressureEqnIndex,PressureVarIndex>;
int verbosity = 0;
if (comm.communicator().rank() == 0) {
verbosity = params.cpr_solver_verbose_;
}
// TODO: revise choice of parameters
int coarsenTarget=1200;
using Criterion = C;
Criterion criterion(15, coarsenTarget);
criterion.setDebugLevel( 0 ); // no debug information, 1 for printing hierarchy information
criterion.setDebugLevel( verbosity ); // no debug information, 1 for printing hierarchy information
criterion.setDefaultValuesIsotropic(2);
criterion.setNoPostSmoothSteps( 1 );
criterion.setNoPreSmoothSteps( 1 );
@@ -207,7 +218,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm,
smootherArgs.relaxationFactor = relax;
setILUParameters(smootherArgs, params);
amgPtr.reset( new AMG( params, opA, criterion, smootherArgs, comm ) );
amgPtr.reset( new AMG( params, weights, opA, criterion, smootherArgs, comm ) );
}
template < class C, class Op, class P, class AMG >
@@ -239,7 +250,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const MILU_VARIANT m
/// \param relax The relaxation parameter for ILU0.
/// \param comm The object describing the parallelization information and communication.
// \param amgPtr The unique_ptr to be filled (return)
template < int pressureIndex=0, class Op, class P, class AMG >
template < int PressureEqnIndex, int PressureVarIndex, class Op, class P, class AMG >
inline void
createAMGPreconditionerPointer( Op& opA, const double relax, const MILU_VARIANT milu, const P& comm, std::unique_ptr< AMG >& amgPtr )
{
@@ -247,7 +258,7 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const MILU_VARIANT
typedef typename Op::matrix_type M;
// The coupling metric used in the AMG
typedef Dune::Amg::Diagonal<pressureIndex> CouplingMetric;
typedef Opm::Amg::Element<PressureEqnIndex, PressureVarIndex> CouplingMetric;
// The coupling criterion used in the AMG
typedef Dune::Amg::SymmetricCriterion<M, CouplingMetric> CritBase;

90
opm/autodiff/FlowLinearSolverParameters.hpp
View File

@@ -59,6 +59,13 @@ NEW_PROP_TAG(UseAmg);
NEW_PROP_TAG(UseCpr);
NEW_PROP_TAG(LinearSolverBackend);
NEW_PROP_TAG(PreconditionerAddWellContributions);
NEW_PROP_TAG(SystemStrategy);
NEW_PROP_TAG(ScaleLinearSystem);
NEW_PROP_TAG(CprSolverVerbose);
NEW_PROP_TAG(CprUseDrs);
NEW_PROP_TAG(CprMaxEllIter);
NEW_PROP_TAG(CprEllSolvetype);
NEW_PROP_TAG(CprReuseSetup);
SET_SCALAR_PROP(FlowIstlSolverParams, LinearSolverReduction, 1e-2);
SET_SCALAR_PROP(FlowIstlSolverParams, IluRelaxation, 0.9);
@@ -76,6 +83,13 @@ SET_BOOL_PROP(FlowIstlSolverParams, UseAmg, false);
SET_BOOL_PROP(FlowIstlSolverParams, UseCpr, false);
SET_TYPE_PROP(FlowIstlSolverParams, LinearSolverBackend, Opm::ISTLSolverEbos<TypeTag>);
SET_BOOL_PROP(FlowIstlSolverParams, PreconditionerAddWellContributions, false);
SET_STRING_PROP(FlowIstlSolverParams, SystemStrategy, "none");
SET_BOOL_PROP(FlowIstlSolverParams, ScaleLinearSystem, false);
SET_INT_PROP(FlowIstlSolverParams, CprSolverVerbose, 0);
SET_BOOL_PROP(FlowIstlSolverParams, CprUseDrs, false);
SET_INT_PROP(FlowIstlSolverParams, CprMaxEllIter, 20);
SET_INT_PROP(FlowIstlSolverParams, CprEllSolvetype, 0);
SET_INT_PROP(FlowIstlSolverParams, CprReuseSetup, 0);
@@ -95,47 +109,31 @@ namespace Opm
MILU_VARIANT cpr_ilu_milu_;
bool cpr_ilu_redblack_;
bool cpr_ilu_reorder_sphere_;
bool cpr_use_drs_;
int cpr_max_ell_iter_;
int cpr_ell_solvetype_;
bool cpr_use_amg_;
bool cpr_use_bicgstab_;
bool cpr_solver_verbose_;
int cpr_solver_verbose_;
bool cpr_pressure_aggregation_;
int cpr_reuse_setup_;
CPRParameter() { reset(); }
CPRParameter( const ParameterGroup& param)
{
// reset values to default
reset();
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_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_);
cpr_ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_);
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_);
cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_);
cpr_pressure_aggregation_ = param.getDefault("cpr_pressure_aggregation", cpr_pressure_aggregation_);
std::string milu("ILU");
cpr_ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu));
}
void reset()
{
cpr_relax_ = 1.0;
cpr_solver_tol_ = 1e-2;
cpr_ilu_n_ = 0;
cpr_ilu_milu_ = MILU_VARIANT::ILU;
cpr_ilu_redblack_ = false;
cpr_ilu_reorder_sphere_ = true;
cpr_max_ell_iter_ = 25;
cpr_ell_solvetype_ = 0;
cpr_use_drs_ = false;
cpr_use_amg_ = true;
cpr_use_bicgstab_ = true;
cpr_solver_verbose_ = false;
cpr_solver_verbose_ = 0;
cpr_pressure_aggregation_ = false;
cpr_reuse_setup_ = 0;
}
};
@@ -160,6 +158,8 @@ namespace Opm
bool ignoreConvergenceFailure_;
bool linear_solver_use_amg_;
bool use_cpr_;
std::string system_strategy_;
bool scale_linear_system_;
template <class TypeTag>
void init()
@@ -179,6 +179,13 @@ namespace Opm
ignoreConvergenceFailure_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure);
linear_solver_use_amg_ = EWOMS_GET_PARAM(TypeTag, bool, UseAmg);
use_cpr_ = EWOMS_GET_PARAM(TypeTag, bool, UseCpr);
system_strategy_ = EWOMS_GET_PARAM(TypeTag, std::string, SystemStrategy);
scale_linear_system_ = EWOMS_GET_PARAM(TypeTag, bool, ScaleLinearSystem);
cpr_solver_verbose_ = EWOMS_GET_PARAM(TypeTag, int, CprSolverVerbose);
cpr_use_drs_ = EWOMS_GET_PARAM(TypeTag, bool, CprUseDrs);
cpr_max_ell_iter_ = EWOMS_GET_PARAM(TypeTag, int, CprMaxEllIter);
cpr_ell_solvetype_ = EWOMS_GET_PARAM(TypeTag, int, CprEllSolvetype);
cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup);
}
template <class TypeTag>
@@ -198,37 +205,16 @@ namespace Opm
EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseAmg, "Use AMG as the linear solver's preconditioner");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseCpr, "Use CPR as the linear solver's preconditioner");
EWOMS_REGISTER_PARAM(TypeTag, std::string, SystemStrategy, "Strategy for reformulating and scaling linear system (none: no scaling -- should not be used with CPR, original: use weights that are equivalent to no scaling -- should not be used with CPR, simple: form pressure equation as simple sum of conservation equations, quasiimpes: form pressure equation based on diagonal block, trueimpes: form pressure equation based on linearization of accumulation term)");
EWOMS_REGISTER_PARAM(TypeTag, bool, ScaleLinearSystem, "Scale linear system according to equation scale and primary variable types");
EWOMS_REGISTER_PARAM(TypeTag, int, CprSolverVerbose, "Verbosity of cpr solver (0: silent, 1: print summary of inner linear solver, 2: print extensive information about inner linear solve, including setup information)");
EWOMS_REGISTER_PARAM(TypeTag, bool, CprUseDrs, "Use dynamic row sum using weights");
EWOMS_REGISTER_PARAM(TypeTag, int, CprMaxEllIter, "MaxIterations of the elliptic pressure part of the cpr solver");
EWOMS_REGISTER_PARAM(TypeTag, int, CprEllSolvetype, "Solver type of elliptic pressure solve (0: bicgstab, 1: cg, 2: only amg preconditioner)");
EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse Amg Setup");
}
FlowLinearSolverParameters() { reset(); }
// read values from parameter class
FlowLinearSolverParameters( const ParameterGroup& param )
: CPRParameter(param)
{
// set default parameters
reset();
// read parameters (using previsouly set default values)
newton_use_gmres_ = param.getDefault("newton_use_gmres", newton_use_gmres_ );
linear_solver_reduction_ = param.getDefault("linear_solver_reduction", linear_solver_reduction_ );
linear_solver_maxiter_ = param.getDefault("linear_solver_maxiter", linear_solver_maxiter_);
linear_solver_restart_ = param.getDefault("linear_solver_restart", linear_solver_restart_);
linear_solver_verbosity_ = param.getDefault("linear_solver_verbosity", linear_solver_verbosity_);
require_full_sparsity_pattern_ = param.getDefault("require_full_sparsity_pattern", require_full_sparsity_pattern_);
ignoreConvergenceFailure_ = param.getDefault("linear_solver_ignoreconvergencefailure", ignoreConvergenceFailure_);
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_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_);
ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_);
std::string milu("ILU");
ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu));
// Check whether to use cpr approach
const std::string cprSolver = "cpr";
use_cpr_ = ( param.getDefault("solver_approach", std::string()) == cprSolver );
}
// set default values
void reset()

309
opm/autodiff/ISTLSolverEbos.hpp
View File

@@ -1,5 +1,6 @@
/*
Copyright 2016 IRIS AS
Copyright 2019 Equinor ASA
This file is part of the Open Porous Media project (OPM).
@@ -31,6 +32,7 @@
#include <opm/common/Exceptions.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
#include <ewoms/common/parametersystem.hh>
#include <ewoms/common/propertysystem.hh>
@@ -59,6 +61,17 @@ END_PROPERTIES
namespace Opm
{
template <class DenseMatrix>
DenseMatrix transposeDenseMatrix(const DenseMatrix& M)
{
DenseMatrix tmp;
for (int i = 0; i < M.rows; ++i)
for (int j = 0; j < M.cols; ++j)
tmp[j][i] = M[i][j];
return tmp;
}
//=====================================================================
// Implementation for ISTL-matrix based operator
//=====================================================================
@@ -162,14 +175,10 @@ protected:
/// 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 TypeTag>
class ISTLSolverEbos
{
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
@@ -178,10 +187,17 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename SparseMatrixAdapter::IstlMatrix Matrix;
enum { pressureIndex = Indices::pressureSwitchIdx };
typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
typedef typename Vector::block_type BlockVector;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
// Due to miscibility oil <-> gas the water eqn is the one we can replace with a pressure equation.
enum { pressureEqnIndex = BlackOilDefaultIndexTraits::waterCompIdx };
enum { pressureVarIndex = Indices::pressureSwitchIdx };
static const int numEq = Indices::numEq;
public:
typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType;
@@ -208,19 +224,69 @@ protected:
matrix_for_preconditioner_.reset();
}
void prepare(const SparseMatrixAdapter& M, const Vector& b) {
void prepare(const SparseMatrixAdapter& M, Vector& b)
{
matrix_.reset(new Matrix(M.istlMatrix()));
rhs_ = &b;
this->scaleSystem();
}
void setResidual(Vector& b) {
rhs_ = &b;
void scaleSystem()
{
const bool matrix_cont_added = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
if (matrix_cont_added) {
bool form_cpr = true;
if (parameters_.system_strategy_ == "quasiimpes") {
weights_ = getQuasiImpesWeights();
} else if (parameters_.system_strategy_ == "trueimpes") {
weights_ = getStorageWeights();
} else if (parameters_.system_strategy_ == "simple") {
BlockVector bvec(1.0);
weights_ = getSimpleWeights(bvec);
} else if (parameters_.system_strategy_ == "original") {
BlockVector bvec(0.0);
bvec[pressureEqnIndex] = 1;
weights_ = getSimpleWeights(bvec);
} else {
if (parameters_.system_strategy_ != "none") {
OpmLog::warning("unknown_system_strategy", "Unknown linear solver system strategy: '" + parameters_.system_strategy_ + "', applying 'none' strategy.");
}
form_cpr = false;
}
if (parameters_.scale_linear_system_) {
// also scale weights
this->scaleEquationsAndVariables(weights_);
}
if (form_cpr && !(parameters_.cpr_use_drs_)) {
scaleMatrixAndRhs(weights_);
}
if (weights_.size() == 0) {
// if weights are not set cpr_use_drs_=false;
parameters_.cpr_use_drs_ = false;
}
} else {
if (parameters_.use_cpr_ && parameters_.cpr_use_drs_) {
OpmLog::warning("DRS_DISABLE", "Disabling DRS as matrix does not contain well contributions");
}
parameters_.cpr_use_drs_ = false;
if (parameters_.scale_linear_system_) {
// also scale weights
this->scaleEquationsAndVariables(weights_);
}
}
}
void setResidual(Vector& /* b */) {
// rhs_ = &b; // Must be handled in prepare() instead.
}
void getResidual(Vector& b) const {
b = *rhs_;
}
void setMatrix(const SparseMatrixAdapter& M) {
matrix_ = &M.istlMatrix();
void setMatrix(const SparseMatrixAdapter& /* M */) {
// matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
}
bool solve(Vector& x) {
@@ -245,13 +311,17 @@ protected:
}
else
{
const WellModel& wellModel = simulator_.problem().wellModel();
typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false > Operator;
Operator opA(*matrix_, *matrix_, wellModel);
solve( opA, x, *rhs_ );
}
return converged_;
if (parameters_.scale_linear_system_) {
scaleSolution(x);
}
return converged_;
}
@@ -313,11 +383,11 @@ protected:
if ( parameters_.use_cpr_ )
{
using Matrix = typename MatrixOperator::matrix_type;
using CouplingMetric = Dune::Amg::Diagonal<pressureIndex>;
using CouplingMetric = Opm::Amg::Element<pressureEqnIndex, pressureVarIndex>;
using CritBase = Dune::Amg::SymmetricCriterion<Matrix, CouplingMetric>;
using Criterion = Dune::Amg::CoarsenCriterion<CritBase>;
using AMG = typename ISTLUtility
::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureIndex >::AMG;
::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureEqnIndex, pressureVarIndex >::AMG;
std::unique_ptr< AMG > amg;
// Construct preconditioner.
@@ -351,13 +421,13 @@ protected:
}
// 3x3 matrix block inversion was unstable at least 2.3 until and including
// 2.5.0. There may still be some issue with the 4x4 matrix block inversion
// we therefore still use the block inversion in OPM
// 3x3 matrix block inversion was unstable at least 2.3 until and including
// 2.5.0. There may still be some issue with the 4x4 matrix block inversion
// we therefore still use the block inversion in OPM
typedef ParallelOverlappingILU0<Dune::BCRSMatrix<Dune::MatrixBlock<typename Matrix::field_type,
Matrix::block_type::rows,
Matrix::block_type::cols> >,
Vector, Vector> SeqPreconditioner;
Matrix::block_type::rows,
Matrix::block_type::cols> >,
Vector, Vector> SeqPreconditioner;
template <class Operator>
@@ -401,7 +471,7 @@ protected:
void
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, const MILU_VARIANT milu) const
{
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, milu, comm, amg );
ISTLUtility::template createAMGPreconditionerPointer<pressureEqnIndex, pressureVarIndex>( *opA, relax, milu, comm, amg );
}
@@ -411,7 +481,7 @@ protected:
const MILU_VARIANT milu ) const
{
ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax,
comm, amg, parameters_ );
comm, amg, parameters_, weights_ );
}
@@ -574,17 +644,208 @@ protected:
}
}
// Weights to make approximate pressure equations.
// Calculated from the storage terms (only) of the
// conservation equations, ignoring all other terms.
Vector getStorageWeights() const
{
Vector weights(rhs_->size());
BlockVector rhs(0.0);
rhs[pressureVarIndex] = 1.0;
int index = 0;
ElementContext elemCtx(simulator_);
const auto& vanguard = simulator_.vanguard();
auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
Dune::FieldVector<Evaluation, numEq> storage;
unsigned threadId = ThreadManager::threadId();
simulator_.model().localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,/*spaceIdx=*/0, /*timeIdx=*/0);
Scalar extrusionFactor = elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor();
Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor;
Scalar storage_scale = scvVolume / elemCtx.simulator().timeStepSize();
MatrixBlockType block;
double pressure_scale = 50e5;
for (int ii = 0; ii < numEq; ++ii) {
for (int jj = 0; jj < numEq; ++jj) {
block[ii][jj] = storage[ii].derivative(jj)/storage_scale;
if (jj == pressureVarIndex) {
block[ii][jj] *= pressure_scale;
}
}
}
BlockVector bweights;
MatrixBlockType block_transpose = Opm::transposeDenseMatrix(block);
block_transpose.solve(bweights, rhs);
bweights /= 1000.0; // given normal densities this scales weights to about 1.
weights[index] = bweights;
++index;
}
return weights;
}
// Interaction between the CPR weights (the function argument 'weights')
// and the variable and equation weights from
// simulator_.model().primaryVarWeight() and
// simulator_.model().eqWeight() is nontrivial and does not work
// at the moment. Possibly refactoring of ewoms weight treatment
// is needed. In the meantime this function shows what needs to be
// done to integrate the weights properly.
void scaleEquationsAndVariables(Vector& weights)
{
// loop over primary variables
const auto endi = matrix_->end();
for (auto i = matrix_->begin(); i != endi; ++i) {
const auto endj = (*i).end();
BlockVector& brhs = (*rhs_)[i.index()];
for (auto j = (*i).begin(); j != endj; ++j) {
MatrixBlockType& block = *j;
for (std::size_t ii = 0; ii < block.rows; ii++ ) {
for (std::size_t jj = 0; jj < block.cols; jj++) {
double var_scale = simulator_.model().primaryVarWeight(i.index(),jj);
block[ii][jj] /= var_scale;
block[ii][jj] *= simulator_.model().eqWeight(i.index(), ii);
}
}
}
for (std::size_t ii = 0; ii < brhs.size(); ii++) {
brhs[ii] *= simulator_.model().eqWeight(i.index(), ii);
}
if (weights.size() == matrix_->N()) {
BlockVector& bw = weights[i.index()];
for (std::size_t ii = 0; ii < brhs.size(); ii++) {
bw[ii] /= simulator_.model().eqWeight(i.index(), ii);
}
double abs_max =
*std::max_element(bw.begin(), bw.end(), [](double a, double b){ return std::abs(a) < std::abs(b); } );
bw /= abs_max;
}
}
}
void scaleSolution(Vector& x)
{
for (std::size_t i = 0; i < x.size(); ++i) {
auto& bx = x[i];
for (std::size_t jj = 0; jj < bx.size(); jj++) {
double var_scale = simulator_.model().primaryVarWeight(i,jj);
bx[jj] /= var_scale;
}
}
}
Vector getQuasiImpesWeights()
{
Matrix& A = *matrix_;
Vector weights(rhs_->size());
BlockVector rhs(0.0);
rhs[pressureVarIndex] = 1;
const auto endi = A.end();
for (auto i = A.begin(); i!=endi; ++i) {
const auto endj = (*i).end();
MatrixBlockType diag_block(0.0);
for (auto j=(*i).begin(); j!=endj; ++j) {
if (i.index() == j.index()) {
diag_block = (*j);
break;
}
}
BlockVector bweights;
auto diag_block_transpose = Opm::transposeDenseMatrix(diag_block);
diag_block_transpose.solve(bweights, rhs);
double abs_max =
*std::max_element(bweights.begin(), bweights.end(), [](double a, double b){ return std::abs(a) < std::abs(b); } );
bweights /= std::abs(abs_max);
weights[i.index()] = bweights;
}
return weights;
}
Vector getSimpleWeights(const BlockVector& rhs)
{
Vector weights(rhs_->size(), 0);
for (auto& bw : weights) {
bw = rhs;
}
return weights;
}
void scaleMatrixAndRhs(const Vector& weights)
{
using Block = typename Matrix::block_type;
const auto endi = matrix_->end();
for (auto i = matrix_->begin(); i !=endi; ++i) {
const BlockVector& bweights = weights[i.index()];
BlockVector& brhs = (*rhs_)[i.index()];
const auto endj = (*i).end();
for (auto j = (*i).begin(); j != endj; ++j) {
// assume it is something on all rows
Block& block = (*j);
BlockVector neweq(0.0);
for (std::size_t ii = 0; ii < block.rows; ii++) {
for (std::size_t jj = 0; jj < block.cols; jj++) {
neweq[jj] += bweights[ii]*block[ii][jj];
}
}
block[pressureEqnIndex] = neweq;
}
Scalar newrhs(0.0);
for (std::size_t ii = 0; ii < brhs.size(); ii++) {
newrhs += bweights[ii]*brhs[ii];
}
brhs[pressureEqnIndex] = newrhs;
}
}
static void multBlocksInMatrix(Matrix& ebosJac, const MatrixBlockType& trans, const bool left = true)
{
const int n = ebosJac.N();
for (int row_index = 0; row_index < n; ++row_index) {
auto& row = ebosJac[row_index];
auto* dataptr = row.getptr();
for (int elem = 0; elem < row.N(); ++elem) {
auto& block = dataptr[elem];
if (left) {
block = block.leftmultiply(trans);
} else {
block = block.rightmultiply(trans);
}
}
}
}
static void multBlocksVector(Vector& ebosResid_cp, const MatrixBlockType& leftTrans)
{
for (auto& bvec : ebosResid_cp) {
auto bvec_new = bvec;
leftTrans.mv(bvec, bvec_new);
bvec = bvec_new;
}
}
static void scaleCPRSystem(Matrix& M_cp, Vector& b_cp, const MatrixBlockType& leftTrans)
{
multBlocksInMatrix(M_cp, leftTrans, true);
multBlocksVector(b_cp, leftTrans);
}
const Simulator& simulator_;
mutable int iterations_;
mutable bool converged_;
boost::any parallelInformation_;
bool isIORank_;
const Matrix *matrix_;
std::unique_ptr<Matrix> matrix_;
Vector *rhs_;
std::unique_ptr<Matrix> matrix_for_preconditioner_;
std::vector<std::pair<int,std::vector<int>>> overlapRowAndColumns_;
FlowLinearSolverParameters parameters_;
Vector weights_;
bool scale_variables_;
}; // end ISTLSolver
} // namespace Opm

3
tests/test_blackoil_amg.cpp
View File

@@ -303,7 +303,8 @@ void runBlackoilAmgLaplace()
smootherArgs.iterations = 1;
Opm::CPRParameter param;
Opm::BlackoilAmg<Operator,ParSmoother,Criterion,Communication,0> amg(param,
Opm::BlackoilAmg<Operator,ParSmoother,Criterion,Communication,0,0> amg(param,
{},
fop, criterion,
smootherArgs,
comm);

3
tests/test_invert.cpp
View File

@@ -21,7 +21,8 @@
#define BOOST_TEST_MODULE InvertSpecializationTest
#include <boost/test/unit_test.hpp>
#include <opm/autodiff/ISTLSolverEbos.hpp>
#include <opm/autodiff/MatrixBlock.hpp>
void checkIdentity(Dune::FieldMatrix<double, 4, 4> M) {
double diag = 0.0;

2
tests/test_multmatrixtransposed.cpp
View File

@@ -21,7 +21,7 @@
#define BOOST_TEST_MODULE MultMatrixTransposed
#include <boost/test/unit_test.hpp>
#include <opm/autodiff/ISTLSolverEbos.hpp>
#include <opm/autodiff/MatrixBlock.hpp>
using namespace Dune;
using namespace Opm::Detail;