This commits allows for flexible choice of either ILU(0) or ILU(n) where n is a

dynamical parameter given in the parameter file. The default is 0 (as before).
In addition the relaxation parameter has been added to the parameter with the
default preserving the state from before.
Also, the default parameter for use_amg and use_bicgstab in the constructor of
CPRPrecondition have been removed.
This commit is contained in:
Robert K 2014-12-04 13:55:56 +01:00
parent dacfbcc583
commit 9a2a95c6eb
3 changed files with 41 additions and 24 deletions

View File

@ -78,11 +78,14 @@ namespace Opm
//! \brief Elliptic Operator
typedef Dune::MatrixAdapter<M,X,X> Operator;
//! \brief preconditioner for the whole system (here either ILU(0) or ILU(n)
typedef Dune::Preconditioner<X,X> WholeSystemPreconditioner;
//! \brief ilu-0 preconditioner for the elliptic system
typedef Dune::SeqILU0<M,X,X> Preconditioner;
typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
//! \brief amg preconditioner for the elliptic system
typedef Preconditioner Smoother;
typedef EllipticPreconditioner Smoother;
typedef Dune::Amg::AMG<Operator, X, Smoother> AMG;
/*! \brief Constructor.
@ -95,8 +98,9 @@ namespace Opm
\param useBiCG if true, BiCG solver is used (default), otherwise CG solver
*/
CPRPreconditioner (const M& A, const M& Ae, const field_type relax,
const bool useAMG = false,
const bool useBiCG = true )
const unsigned int ilu_n,
const bool useAMG,
const bool useBiCG )
: A_(A),
Ae_(Ae),
de_( Ae_.N() ),
@ -105,15 +109,20 @@ namespace Opm
opAe_( Ae_ ),
precond_(), // ilu0 preconditioner for elliptic system
amg_(), // amg preconditioner for elliptic system
ILU_(A), // copy A (will be overwritten by ILU decomp)
vilu_( ILU_.N() ),
pre_(), // copy A will be made be the preconditioner
vilu_( A_.N() ),
relax_(relax),
use_bicg_solver_( useBiCG )
{
// create appropriate preconditioner for elliptic system
createPreconditioner( useAMG );
Dune::bilu0_decomposition(ILU_);
if( ilu_n == 0 ) {
pre_.reset( new Dune::SeqILU0<M,X,X>( A_, relax_) );
}
else {
pre_.reset( new Dune::SeqILUn<M,X,X>( A_, relax_, ilu_n) );
}
}
/*!
@ -151,14 +160,17 @@ namespace Opm
dmodified_ = d;
A_.mmv(v, dmodified_);
// Apply ILU0.
Dune::bilu_backsolve(ILU_, vilu_, dmodified_);
v += vilu_;
// Apply Preconditioner for whole system (relax will be applied already)
pre_->apply( vilu_, dmodified_);
// don't apply relaxation if relax_ == 1
if( std::abs( relax_ - 1.0 ) < 1e-12 ) return;
v *= relax_;
// don't apply relaxation if relax_ == 1
if( std::abs( relax_ - 1.0 ) < 1e-12 ) {
v += vilu_;
}
else {
v *= relax_;
v += vilu_;
}
}
/*!
@ -227,12 +239,12 @@ namespace Opm
Operator opAe_;
//! \brief ILU0 preconditioner for the elliptic system
std::unique_ptr< Preconditioner > precond_;
std::unique_ptr< EllipticPreconditioner > precond_;
//! \brief AMG preconditioner with ILU0 smoother
std::unique_ptr< AMG > amg_;
//! \brief The ILU0 decomposition of the matrix.
matrix_type ILU_;
//! \brief The preconditioner for the whole system
std::unique_ptr< WholeSystemPreconditioner > pre_;
//! \brief temporary variables for ILU solve
Y vilu_;
@ -266,7 +278,7 @@ namespace Opm
amg_ = std::unique_ptr< AMG > (new AMG(opAe_, criterion, smootherArgs));
}
else
precond_ = std::unique_ptr< Preconditioner > (new Preconditioner( Ae_, relax_ ));
precond_ = std::unique_ptr< EllipticPreconditioner > (new EllipticPreconditioner( Ae_, relax_ ));
}
};

View File

@ -117,8 +117,10 @@ namespace Opm
NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param)
: iterations_( 0 )
{
use_amg_ = param.getDefault("cpr_use_amg", false);
use_bicgstab_ = param.getDefault("cpr_use_bicgstab", true);
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);
}
@ -193,8 +195,7 @@ namespace Opm
// Construct preconditioner.
// typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
typedef Opm::CPRPreconditioner<Mat,Vector,Vector> Preconditioner;
const double relax = 1.0;
Preconditioner precond(istlA, istlAe, relax, use_amg_, use_bicgstab_);
Preconditioner precond(istlA, istlAe, cpr_relax_, cpr_ilu_n_, cpr_use_amg_, cpr_use_bicgstab_);
// Construct linear solver.
const double tolerance = 1e-3;

View File

@ -43,6 +43,8 @@ namespace Opm
/// the preconditioning and choice of
/// linear solvers.
/// Parameters:
/// cpr_relax (default 1.0) relaxation for the preconditioner
/// cpr_ilu_n (default 0) use ILU(n) for preconditioning of the linear system
/// cpr_use_amg (default false) if true, use AMG preconditioner for elliptic part
/// cpr_use_bicgstab (default true) if true, use BiCGStab (else use CG) for elliptic part
NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param);
@ -58,8 +60,10 @@ namespace Opm
virtual int iterations () const { return iterations_; }
private:
mutable int iterations_;
bool use_amg_;
bool use_bicgstab_;
double cpr_relax_;
unsigned int cpr_ilu_n_;
bool cpr_use_amg_;
bool cpr_use_bicgstab_;
};
} // namespace Opm