24 #ifndef OPM_CPRPRECONDITIONER_HEADER_INCLUDED
25 #define OPM_CPRPRECONDITIONER_HEADER_INCLUDED
28 #include <type_traits>
30 #include <opm/common/utility/platform_dependent/disable_warnings.h>
31 #include <opm/core/utility/parameters/ParameterGroup.hpp>
33 #include <dune/istl/bvector.hh>
34 #include <dune/istl/bcrsmatrix.hh>
35 #include <dune/istl/operators.hh>
36 #include <dune/istl/io.hh>
37 #include <dune/istl/owneroverlapcopy.hh>
38 #include <dune/istl/preconditioners.hh>
39 #include <dune/istl/schwarz.hh>
40 #include <dune/istl/solvers.hh>
41 #include <dune/istl/paamg/amg.hh>
42 #include <dune/istl/paamg/kamg.hh>
43 #include <dune/istl/paamg/pinfo.hh>
45 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
47 #include <opm/common/ErrorMacros.hpp>
48 #include <opm/common/Exceptions.hpp>
50 #include <opm/autodiff/AdditionalObjectDeleter.hpp>
51 #include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
52 #include <opm/autodiff/ParallelOverlappingILU0.hpp>
66 template<
class M,
class X,
class Y,
class P>
80 typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
93 template<
class M,
class X,
class Y,
class I1,
class I2>
95 struct CPRSelector<
M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
100 typedef Dune::OverlappingSchwarzOperator<M,X,X,ParallelInformation>
Operator;
114 typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
131 template<
class ILU,
class I1,
class I2>
132 AdditionalObjectDeleter<ILU>
133 createParallelDeleter(ILU& ilu,
const Dune::OwnerOverlapCopyCommunication<I1,I2>& p)
136 return AdditionalObjectDeleter<ILU>(ilu);
143 template<
class M,
class X>
144 std::shared_ptr<Dune::SeqILU0<M,X,X> >
145 createILU0Ptr(
const M& A,
double relax,
const Dune::Amg::SequentialInformation&)
147 return std::shared_ptr<Dune::SeqILU0<M,X,X> >(
new Dune::SeqILU0<M,X,X>( A, relax) );
153 template<
class M,
class X>
154 std::shared_ptr<Dune::SeqILUn<M,X,X> >
155 createILUnPtr(
const M& A,
int ilu_n,
double relax,
const Dune::Amg::SequentialInformation&)
157 return std::shared_ptr<Dune::SeqILUn<M,X,X> >(
new Dune::SeqILUn<M,X,X>( A, ilu_n, relax) );
161 template<
class ILU,
class I1,
class I2>
162 struct SelectParallelILUSharedPtr
164 typedef std::shared_ptr<
165 Dune::BlockPreconditioner<
166 typename ILU::range_type,
167 typename ILU::domain_type,
168 Dune::OwnerOverlapCopyCommunication<I1,I2>,
178 template<
class M,
class X,
class I1,
class I2>
179 typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>::type
180 createILU0Ptr(
const M& A,
double relax,
181 const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
183 typedef Dune::BlockPreconditioner<
186 Dune::OwnerOverlapCopyCommunication<I1,I2>,
189 Dune::SeqILU0<M,X,X>* ilu =
nullptr;
190 int ilu_setup_successful = 1;
193 ilu =
new Dune::SeqILU0<M,X,X>(A, relax);
195 catch ( Dune::MatrixBlockError error )
197 message = error.what();
198 std::cerr<<
"Exception occured on process " <<
199 comm.communicator().rank() <<
" during " <<
200 "setup of ILU0 preconditioner with message: " <<
202 ilu_setup_successful = 0;
205 if ( comm.communicator().min(ilu_setup_successful) == 0 )
212 throw Dune::MatrixBlockError();
214 return typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>
215 ::type (
new PointerType(*ilu, comm), createParallelDeleter(*ilu, comm));
223 template<
class M,
class X,
class I1,
class I2>
224 typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
225 createILUnPtr(
const M& A,
int ilu_n,
double relax,
226 const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
228 typedef Dune::BlockPreconditioner<
231 Dune::OwnerOverlapCopyCommunication<I1,I2>,
234 Dune::SeqILUn<M,X,X>* ilu =
new Dune::SeqILUn<M,X,X>( A, ilu_n, relax);
236 return typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
237 (
new PointerType(*ilu, comm),createParallelDeleter(*ilu, comm));
244 template<
class M,
class X=
typename M::range_type>
245 std::unique_ptr<Dune::SeqILU0<M,X,X> >
246 createEllipticPreconditionerPointer(
const M& Ae,
double relax,
247 const Dune::Amg::SequentialInformation&)
249 return std::unique_ptr<Dune::SeqILU0<M,X,X> >(
new Dune::SeqILU0<M,X,X>(Ae, relax));
253 template<
class M,
class X=
typename M::range_type,
class I1,
class I2>
258 typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
259 ::EllipticPreconditionerPointer
260 createEllipticPreconditionerPointer(
const M& Ae,
double relax,
261 const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
263 typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
264 ::EllipticPreconditioner ParallelPreconditioner;
267 typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
268 ::EllipticPreconditionerPointer EllipticPreconditionerPointer;
269 return EllipticPreconditionerPointer(
new ParallelPreconditioner(Ae, comm, relax));
272 #endif // #if HAVE_MPI
279 template <
int pressureIndex=0,
class Op,
class P,
class AMG >
281 createAMGPreconditionerPointer( Op& opA,
const double relax,
const P& comm, std::unique_ptr< AMG >& amgPtr )
284 typedef typename Op::matrix_type M;
287 typedef Dune::Amg::Diagonal<pressureIndex> CouplingMetric;
290 typedef Dune::Amg::SymmetricCriterion<M, CouplingMetric> CritBase;
293 typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
296 int coarsenTarget=1200;
297 Criterion criterion(15,coarsenTarget);
298 criterion.setDebugLevel( 0 );
299 criterion.setDefaultValuesIsotropic(2);
300 criterion.setNoPostSmoothSteps( 1 );
301 criterion.setNoPreSmoothSteps( 1 );
304 typedef typename AMG::Smoother Smoother;
305 typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
306 SmootherArgs smootherArgs;
307 smootherArgs.iterations = 1;
308 smootherArgs.relaxationFactor = relax;
310 amgPtr.reset(
new AMG(opA, criterion, smootherArgs, comm ) );
318 double cpr_solver_tol_;
320 int cpr_max_ell_iter_;
322 bool cpr_use_bicgstab_;
323 bool cpr_solver_verbose_;
332 cpr_relax_ = param.getDefault(
"cpr_relax", cpr_relax_);
333 cpr_solver_tol_ = param.getDefault(
"cpr_solver_tol", cpr_solver_tol_);
334 cpr_ilu_n_ = param.getDefault(
"cpr_ilu_n", cpr_ilu_n_);
335 cpr_max_ell_iter_ = param.getDefault(
"cpr_max_elliptic_iter",cpr_max_ell_iter_);
336 cpr_use_amg_ = param.getDefault(
"cpr_use_amg", cpr_use_amg_);
337 cpr_use_bicgstab_ = param.getDefault(
"cpr_use_bicgstab", cpr_use_bicgstab_);
338 cpr_solver_verbose_ = param.getDefault(
"cpr_solver_verbose", cpr_solver_verbose_);
344 cpr_solver_tol_ = 1e-2;
346 cpr_max_ell_iter_ = 25;
348 cpr_use_bicgstab_ =
true;
349 cpr_solver_verbose_ =
false;
368 template<
class M,
class X,
class Y,
369 class P=Dune::Amg::SequentialInformation>
390 category = std::is_same<P,Dune::Amg::SequentialInformation>::value?
391 Dune::SolverCategory::sequential:Dune::SolverCategory::overlapping
408 typedef typename CPRSelectorType::AMG
AMG;
429 dmodified_(
A_.N() ),
439 createEllipticPreconditioner(
param_.cpr_use_amg_,
commAe_ );
442 if(
param_.cpr_ilu_n_ == 0 ) {
455 virtual void pre (X& , Y& )
464 virtual void apply (X& v,
const Y& d)
468 std::copy_n(d.begin(),
de_.size(),
de_.begin());
473 solveElliptic( ve_,
de_ );
478 std::copy(ve_.begin(), ve_.end(), v.begin());
483 A_.mmv(v, dmodified_);
485 comm_.copyOwnerToAll(dmodified_, dmodified_);
491 if( std::abs(
param_.cpr_relax_ - 1.0 ) < 1e-12 ) {
510 void solveElliptic(Y& x, Y& de)
513 const double tolerance =
param_.cpr_solver_tol_;
514 const int maxit =
param_.cpr_max_ell_iter_;
515 const int verbosity = (
param_.cpr_solver_verbose_ &&
516 comm_.communicator().rank()==0 ) ? 1 : 0;
519 Dune::InverseOperatorResult result;
522 typedef Dune::ScalarProductChooser<X,ParallelInformation,category>
523 ScalarProductChooser;
525 std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
526 sp(ScalarProductChooser::construct(
commAe_));
531 if(
param_.cpr_use_bicgstab_ ) {
532 Dune::BiCGSTABSolver<X> linsolve(*
opAe_, *sp, (*
amg_), tolerance, maxit, verbosity);
533 linsolve.apply(x, de, result);
536 Dune::CGSolver<X> linsolve(*
opAe_, *sp, (*
amg_), tolerance, maxit, verbosity);
537 linsolve.apply(x, de, result);
544 if(
param_.cpr_use_bicgstab_ ) {
545 Dune::BiCGSTABSolver<X> linsolve(*
opAe_, *sp, (*
precond_), tolerance, maxit, verbosity);
546 linsolve.apply(x, de, result);
549 Dune::CGSolver<X> linsolve(*
opAe_, *sp, (*
precond_), tolerance, maxit, verbosity);
550 linsolve.apply(x, de, result);
555 if (!result.converged) {
556 OPM_THROW(LinearSolverProblem,
"CPRPreconditioner failed to solve elliptic subsystem.");
585 std::shared_ptr< WholeSystemPreconditioner >
pre_;
596 void createEllipticPreconditioner(
const bool amg,
const P& comm )
600 ISTLUtility::createAMGPreconditionerPointer( *
opAe_ ,
param_.cpr_relax_, comm,
amg_ );
604 precond_ = ISTLUtility::createEllipticPreconditionerPointer<M,X>(
Ae_,
param_.cpr_relax_, comm);
612 #endif // OPM_CPRPRECONDITIONER_HEADER_INCLUDED
Y vilu_
temporary variables for ILU solve
Definition: CPRPreconditioner.hpp:588
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition: ParallelOverlappingILU0.hpp:35
std::unique_ptr< EllipticPreconditioner > EllipticPreconditionerPointer
The type of the unique pointer to the preconditioner of the elliptic part.
Definition: CPRPreconditioner.hpp:76
static Operator * makeOperator(const M &m, const P &)
creates an Operator from the matrix
Definition: CPRPreconditioner.hpp:85
AutoDiffMatrix is a wrapper class that optimizes matrix operations.
Definition: AutoDiffMatrix.hpp:43
CPRSelectorType::EllipticPreconditionerPointer EllipticPreconditionerPointer
type of the unique pointer to the ilu-0 preconditioner used the for the elliptic system ...
Definition: CPRPreconditioner.hpp:405
const P & commAe_
The information about the parallelization of the elliptic part of the system.
Definition: CPRPreconditioner.hpp:594
std::unique_ptr< AMG > amg_
AMG preconditioner with ILU0 smoother.
Definition: CPRPreconditioner.hpp:577
virtual void pre(X &, Y &)
Prepare the preconditioner.
Definition: CPRPreconditioner.hpp:455
CPR preconditioner.
Definition: CPRPreconditioner.hpp:370
CPRSelectorType::Operator Operator
Elliptic Operator.
Definition: CPRPreconditioner.hpp:397
A custom deleter that will delete an additional object, too.
Definition: AdditionalObjectDeleter.hpp:36
CPRSelectorType::AMG AMG
amg preconditioner for the elliptic system
Definition: CPRPreconditioner.hpp:408
X domain_type
The domain type of the preconditioner.
Definition: CPRPreconditioner.hpp:381
Dune::SeqILU0< M, X, X > EllipticPreconditioner
The type of the preconditioner used for the elliptic part.
Definition: CPRPreconditioner.hpp:74
The category the preconditioner is part of.
Definition: CPRPreconditioner.hpp:390
const P & comm_
The information about the parallelization of the whole system.
Definition: CPRPreconditioner.hpp:591
const CPRParameter & param_
Parameter collection for CPR.
Definition: CPRPreconditioner.hpp:561
EllipticPreconditioner Smoother
type of AMG used to precondition the elliptic system.
Definition: CPRPreconditioner.hpp:79
P ParallelInformation
The type describing the parallel information.
Definition: CPRPreconditioner.hpp:377
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: CPRPreconditioner.hpp:379
Dune::MatrixAdapter< M, X, Y > Operator
The operator type;.
Definition: CPRPreconditioner.hpp:72
EllipticPreconditionerPointer precond_
ILU0 preconditioner for the elliptic system.
Definition: CPRPreconditioner.hpp:575
std::unique_ptr< Operator > opAe_
elliptic operator
Definition: CPRPreconditioner.hpp:572
virtual void apply(X &v, const Y &d)
Apply the preconditoner.
Definition: CPRPreconditioner.hpp:464
const matrix_type & A_
The matrix for the full linear problem.
Definition: CPRPreconditioner.hpp:564
Dune::Preconditioner< X, X > WholeSystemPreconditioner
preconditioner for the whole system (here either ILU(0) or ILU(n)
Definition: CPRPreconditioner.hpp:400
Y range_type
The range type of the preconditioner.
Definition: CPRPreconditioner.hpp:383
Y de_
temporary variables for elliptic solve
Definition: CPRPreconditioner.hpp:569
const matrix_type & Ae_
The elliptic part of the matrix.
Definition: CPRPreconditioner.hpp:566
CPRPreconditioner(const CPRParameter ¶m, const M &A, const M &Ae, const ParallelInformation &comm=ParallelInformation(), const ParallelInformation &commAe=ParallelInformation())
Constructor.
Definition: CPRPreconditioner.hpp:421
A traits class for selecting the types of the preconditioner.
Definition: CPRPreconditioner.hpp:67
Dune::Amg::SequentialInformation ParallelInformation
The information about the parallelization and communication.
Definition: CPRPreconditioner.hpp:70
Definition: CPRPreconditioner.hpp:315
X::field_type field_type
The field type of the preconditioner.
Definition: CPRPreconditioner.hpp:385
std::shared_ptr< WholeSystemPreconditioner > pre_
The preconditioner for the whole system.
Definition: CPRPreconditioner.hpp:585
virtual void post(X &)
Clean up.
Definition: CPRPreconditioner.hpp:505