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 #342 from dr-robertk/PR/some-improvemnts-on-solvers

Browse Source 
More parameter for solvers.
This commit is contained in:
Atgeirr Flø Rasmussen 2015-03-31 18:00:40 +02:00
commit 62d3f37b1d
5 changed files with 113 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

119
opm/autodiff/CPRPreconditioner.hpp
View File

@ -27,6 +27,7 @@
#include <type_traits>
#include <opm/core/utility/platform_dependent/disable_warnings.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
@ -39,6 +40,7 @@
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/kamg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/paamg/fastamg.hh>
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
@ -93,6 +95,12 @@ struct CPRSelector
typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
/// \brief The type of the unique pointer to the preconditioner of the elliptic part.
typedef std::unique_ptr<EllipticPreconditioner> EllipticPreconditionerPointer;
/// \brief type of AMG used to precondition the elliptic system.
typedef EllipticPreconditioner Smoother;
typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
//typedef Dune::Amg::FastAMG<Operator, X> AMG;
/// \brief creates an Operator from the matrix
/// \param M The matrix to use.
/// \param p The parallel information to use.
@ -119,6 +127,9 @@ struct CPRSelector<M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
ParallelPreconditionerDeleter<Dune::SeqILU0<M,X,X> > >
EllipticPreconditionerPointer;
typedef EllipticPreconditioner Smoother;
typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
/// \brief creates an Operator from the matrix
/// \param M The matrix to use.
/// \param p The parallel information to use.
@ -256,6 +267,44 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
#endif
} // end namespace
struct CPRParameter
{
double cpr_relax_;
double cpr_solver_tol_;
int cpr_ilu_n_;
int cpr_max_ell_iter_;
bool cpr_use_amg_;
bool cpr_use_bicgstab_;
bool cpr_solver_verbose_;
CPRParameter() { reset(); }
CPRParameter( const parameter::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_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_);
}
void reset()
{
cpr_relax_ = 1.0;
cpr_solver_tol_ = 1e-4;
cpr_ilu_n_ = 0;
cpr_max_ell_iter_ = 5000;
cpr_use_amg_ = false;
cpr_use_bicgstab_ = true;
cpr_solver_verbose_ = false;
}
};
/*!
\brief CPR preconditioner.
@ -303,18 +352,13 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
//! \brief preconditioner for the whole system (here either ILU(0) or ILU(n)
typedef Dune::Preconditioner<X,X> WholeSystemPreconditioner;
//! \brief the ilu-0 preconditioner used the for the elliptic system
typedef typename CPRSelector<M,X,X,P>::EllipticPreconditioner
EllipticPreconditioner;
//! \brief type of the unique pointer to the ilu-0 preconditioner
//! used the for the elliptic system
typedef typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
EllipticPreconditionerPointer;
//! \brief amg preconditioner for the elliptic system
typedef EllipticPreconditioner Smoother;
typedef Dune::Amg::AMG<Operator, X, Smoother, P> AMG;
typedef typename CPRSelector<M,X,X,P>::AMG AMG;
/*! \brief Constructor.
@ -327,12 +371,10 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
\param paralleInformation The information about the parallelization, if this is a
parallel run
*/
CPRPreconditioner (const M& A, const M& Ae, const field_type relax,
const unsigned int ilu_n,
const bool useAMG,
const bool useBiCG,
CPRPreconditioner (const CPRParameter& param, const M& A, const M& Ae,
const ParallelInformation& comm=ParallelInformation())
: A_(A),
: param_( param ),
A_(A),
Ae_(Ae),
de_( Ae_.N() ),
ve_( Ae_.M() ),
@ -342,18 +384,16 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
amg_(), // amg preconditioner for elliptic system
pre_(), // copy A will be made be the preconditioner
vilu_( A_.N() ),
relax_(relax),
use_bicg_solver_( useBiCG ),
comm_(comm)
{
// create appropriate preconditioner for elliptic system
createPreconditioner( useAMG, comm );
createPreconditioner( param_.cpr_use_amg_, comm );
if( ilu_n == 0 ) {
pre_ = createILU0Ptr<M,X>( A_, relax_, comm );
if( param_.cpr_ilu_n_ == 0 ) {
pre_ = createILU0Ptr<M,X>( A_, param_.cpr_relax_, comm );
}
else {
pre_ = createILUnPtr<M,X>( A_, ilu_n, relax_, comm );
pre_ = createILUnPtr<M,X>( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm );
}
}
@ -396,11 +436,11 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
pre_->apply( vilu_, dmodified_);
// don't apply relaxation if relax_ == 1
if( std::abs( relax_ - 1.0 ) < 1e-12 ) {
if( std::abs( param_.cpr_relax_ - 1.0 ) < 1e-12 ) {
v += vilu_;
}
else {
v *= relax_;
v *= param_.cpr_relax_;
v += vilu_;
}
}
@ -418,9 +458,9 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
void solveElliptic(Y& x, Y& de)
{
// Linear solver parameters
const double tolerance = 1e-4;
const int maxit = 5000;
const int verbosity = 0;
const double tolerance = param_.cpr_solver_tol_;
const int maxit = param_.cpr_max_ell_iter_;
const int verbosity = param_.cpr_solver_verbose_ ? 1 : 0;
// operator result containing iterations etc.
Dune::InverseOperatorResult result;
@ -435,7 +475,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
if( amg_ )
{
// Solve system with AMG
if( use_bicg_solver_ ) {
if( param_.cpr_use_bicgstab_ ) {
Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
linsolve.apply(x, de, result);
}
@ -448,7 +488,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
{
assert( precond_ );
// Solve system with ILU-0
if( use_bicg_solver_ ) {
if( param_.cpr_use_bicgstab_ ) {
Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
linsolve.apply(x, de, result);
}
@ -464,6 +504,9 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
}
}
//! \brief Parameter collection for CPR
const CPRParameter& param_;
//! \brief The matrix for the full linear problem.
const matrix_type& A_;
//! \brief The elliptic part of the matrix.
@ -491,12 +534,6 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
//! \brief temporary variables for ILU solve
Y vilu_;
//! \brief The relaxation factor to use.
field_type relax_;
//! \brief true if ISTL BiCGSTABSolver is used, otherwise ISTL CGSolver is used
const bool use_bicg_solver_;
//! \brief The information about the parallelization
const P& comm_;
protected:
@ -504,30 +541,32 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
{
if( amg )
{
typedef Dune::Amg::CoarsenCriterion< Dune::Amg::SymmetricCriterion<M, Dune::Amg::FirstDiagonal> > Criterion;
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
// The coupling metric used in the AMG
typedef Dune::Amg::FirstDiagonal CouplingMetric;
SmootherArgs smootherArgs;
// The coupling criterion used in the AMG
typedef Dune::Amg::SymmetricCriterion<M, CouplingMetric> CritBase;
smootherArgs.iterations = 1;
smootherArgs.relaxationFactor = relax_;
// The coarsening criterion used in the AMG
typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
// TODO: revise choice of parameters
int coarsenTarget=1200;
Criterion criterion(15,coarsenTarget);
criterion.setDebugLevel( 0 ); // no debug information, 1 for printing hierarchy information
criterion.setDefaultValuesIsotropic(2);
criterion.setAlpha(.67);
criterion.setBeta(1.0e-6);
criterion.setMaxLevel(10);
amg_ = std::unique_ptr< AMG > (new AMG(*opAe_, criterion, smootherArgs));
criterion.setNoPostSmoothSteps( 1 );
criterion.setNoPreSmoothSteps( 1 );
amg_ = std::unique_ptr< AMG > (new AMG(*opAe_, criterion));//, smootherArgs));
}
else
{
precond_ = createEllipticPreconditionerPointer<M,X>( Ae_, relax_, comm);
precond_ = createEllipticPreconditionerPointer<M,X>( Ae_, param_.cpr_relax_, comm);
}
}
};
} // namespace Opm
#endif // OPM_CPRPRECONDITIONER_HEADER_INCLUDED

4
opm/autodiff/FullyImplicitBlackoilSolver.hpp
View File

@ -73,7 +73,8 @@ namespace Opm {
double tolerance_mb_;
double tolerance_cnv_;
double tolerance_wells_;
int max_iter_;
int max_iter_; // max newton iterations
int min_iter_; // min newton iterations
SolverParameter( const parameter::ParameterGroup& param );
SolverParameter();
@ -418,6 +419,7 @@ namespace Opm {
double relaxIncrement() const { return param_.relax_increment_; };
double relaxRelTol() const { return param_.relax_rel_tol_; };
double maxIter() const { return param_.max_iter_; }
double minIter() const { return param_.min_iter_; }
double maxResidualAllowed() const { return param_.max_residual_allowed_; }
};

12
opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
View File

@ -2,6 +2,7 @@
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
Copyright 2015 IRIS AS
This file is part of the Open Porous Media project (OPM).
@ -143,7 +144,8 @@ namespace detail {
relax_max_ = 0.5;
relax_increment_ = 0.1;
relax_rel_tol_ = 0.2;
max_iter_ = 15;
max_iter_ = 15; // not more then 15 its by default
min_iter_ = 0; // keep the default as it was
max_residual_allowed_ = std::numeric_limits< double >::max();
tolerance_mb_ = 1.0e-7;
tolerance_cnv_ = 1.0e-3;
@ -171,6 +173,7 @@ namespace detail {
dr_max_rel_ = param.getDefault("dr_max_rel", dr_max_rel_);
relax_max_ = param.getDefault("relax_max", relax_max_);
max_iter_ = param.getDefault("max_iter", max_iter_);
min_iter_ = param.getDefault("min_iter", min_iter_);
max_residual_allowed_ = param.getDefault("max_residual_allowed", max_residual_allowed_);
tolerance_mb_ = param.getDefault("tolerance_mb", tolerance_mb_);
@ -294,7 +297,7 @@ namespace detail {
const enum RelaxType relaxtype = relaxType();
int linearIterations = 0;
while ((!converged) && (it < maxIter())) {
while ( (!converged && (it < maxIter())) || (minIter() > it)) {
V dx = solveJacobianSystem();
// store number of linear iterations used
@ -326,9 +329,8 @@ namespace detail {
}
if (!converged) {
// the runtime_error is caught by the AdaptiveTimeStepping
OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
return -1;
std::cerr << "WARNING: Failed to compute converged solution in " << it << " iterations." << std::endl;
return -1; // -1 indicates that the solver has to be restarted
}
linearIterations_ += linearIterations;

9
opm/autodiff/NewtonIterationBlackoilCPR.cpp
View File

@ -99,13 +99,10 @@ namespace Opm
/// Construct a system solver.
NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param,
const boost::any& parallelInformation)
: iterations_( 0 ), parallelInformation_(parallelInformation)
: cpr_param_( param ),
iterations_( 0 ),
parallelInformation_(parallelInformation)
{
cpr_relax_ = param.getDefault("cpr_relax", 1.0);
cpr_ilu_n_ = param.getDefault("cpr_ilu_n", 0);
cpr_use_amg_ = param.getDefault("cpr_use_amg", false);
cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", true);
linear_solver_reduction_ = param.getDefault("linear_solver_reduction", 1e-3 );
linear_solver_maxiter_ = param.getDefault("linear_solver_maxiter", 150 );
linear_solver_restart_ = param.getDefault("linear_solver_restart", 40 );

32
opm/autodiff/NewtonIterationBlackoilCPR.hpp
View File

@ -1,5 +1,6 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2015 IRIS AS
This file is part of the Open Porous Media project (OPM).
@ -45,6 +46,7 @@ namespace Opm
typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> Vector;
public:
/// Construct a system solver.
@ -90,23 +92,31 @@ namespace Opm
sp(ScalarProductChooser::construct(parallelInformation));
// Construct preconditioner.
// typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
typedef Opm::CPRPreconditioner<Mat,Vector,Vector,P> Preconditioner;
typedef Opm::CPRPreconditioner<Mat,Vector,Vector,P> Preconditioner;
parallelInformation.copyOwnerToAll(istlb, istlb);
Preconditioner precond(opA.getmat(), istlAe, cpr_relax_, cpr_ilu_n_, cpr_use_amg_, cpr_use_bicgstab_, parallelInformation);
Preconditioner precond(cpr_param_, opA.getmat(), istlAe, parallelInformation);
// TODO: Revise when linear solvers interface opm-core is done
// Construct linear solver.
Dune::RestartedGMResSolver<Vector> linsolve(opA, *sp, precond,
linear_solver_reduction_, linear_solver_restart_, linear_solver_maxiter_, linear_solver_verbosity_);
// Solve system.
linsolve.apply(x, istlb, result);
// GMRes solver
if ( newton_use_gmres_ ) {
Dune::RestartedGMResSolver<Vector> linsolve(opA, *sp, precond,
linear_solver_reduction_, linear_solver_restart_, linear_solver_maxiter_, linear_solver_verbosity_);
// Solve system.
linsolve.apply(x, istlb, result);
}
else { // BiCGstab solver
Dune::BiCGSTABSolver<Vector> linsolve(opA, *sp, precond,
linear_solver_reduction_, linear_solver_maxiter_, linear_solver_verbosity_);
// Solve system.
linsolve.apply(x, istlb, result);
}
}
CPRParameter cpr_param_;
mutable int iterations_;
double cpr_relax_;
unsigned int cpr_ilu_n_;
bool cpr_use_amg_;
bool cpr_use_bicgstab_;
bool newton_use_gmres_;
boost::any parallelInformation_;
double linear_solver_reduction_;