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>
63 template<
class M,
class X,
class Y,
class P>
67 typedef Dune::Amg::SequentialInformation ParallelInformation;
69 typedef Dune::MatrixAdapter<M, X, Y> Operator;
71 typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
73 typedef std::unique_ptr<EllipticPreconditioner> EllipticPreconditionerPointer;
76 typedef EllipticPreconditioner Smoother;
77 typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
82 static Operator* makeOperator(
const M& m,
const P&)
84 return new Operator(m);
89 template<
class M,
class X,
class Y,
class I1,
class I2>
91 struct CPRSelector<
M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
94 typedef Dune::OwnerOverlapCopyCommunication<I1,I2> ParallelInformation;
96 typedef Dune::OverlappingSchwarzOperator<M,X,X,ParallelInformation> Operator;
98 typedef Dune::BlockPreconditioner<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
99 EllipticPreconditioner;
101 typedef std::unique_ptr<EllipticPreconditioner,
102 AdditionalObjectDeleter<Dune::SeqILU0<M,X,X> > >
103 EllipticPreconditionerPointer;
105 typedef EllipticPreconditioner Smoother;
106 typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
111 static Operator* makeOperator(
const M& m,
const ParallelInformation& p)
113 return new Operator(m, p);
123 template<
class ILU,
class I1,
class I2>
124 AdditionalObjectDeleter<ILU>
125 createParallelDeleter(ILU& ilu,
const Dune::OwnerOverlapCopyCommunication<I1,I2>& p)
128 return AdditionalObjectDeleter<ILU>(ilu);
135 template<
class M,
class X>
136 std::shared_ptr<Dune::SeqILU0<M,X,X> >
137 createILU0Ptr(
const M& A,
double relax,
const Dune::Amg::SequentialInformation&)
139 return std::shared_ptr<Dune::SeqILU0<M,X,X> >(
new Dune::SeqILU0<M,X,X>( A, relax) );
145 template<
class M,
class X>
146 std::shared_ptr<Dune::SeqILUn<M,X,X> >
147 createILUnPtr(
const M& A,
int ilu_n,
double relax,
const Dune::Amg::SequentialInformation&)
149 return std::shared_ptr<Dune::SeqILUn<M,X,X> >(
new Dune::SeqILUn<M,X,X>( A, ilu_n, relax) );
153 template<
class ILU,
class I1,
class I2>
154 struct SelectParallelILUSharedPtr
156 typedef std::shared_ptr<
157 Dune::BlockPreconditioner<
158 typename ILU::range_type,
159 typename ILU::domain_type,
160 Dune::OwnerOverlapCopyCommunication<I1,I2>,
170 template<
class M,
class X,
class I1,
class I2>
171 typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>::type
172 createILU0Ptr(
const M& A,
double relax,
173 const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
175 typedef Dune::BlockPreconditioner<
178 Dune::OwnerOverlapCopyCommunication<I1,I2>,
181 Dune::SeqILU0<M,X,X>* ilu =
nullptr;
182 int ilu_setup_successful = 1;
185 ilu =
new Dune::SeqILU0<M,X,X>(A, relax);
187 catch ( Dune::MatrixBlockError error )
189 message = error.what();
190 std::cerr<<
"Exception occured on process " <<
191 comm.communicator().rank() <<
" during " <<
192 "setup of ILU0 preconditioner with message: " <<
194 ilu_setup_successful = 0;
197 if ( comm.communicator().min(ilu_setup_successful) == 0 )
204 throw Dune::MatrixBlockError();
206 return typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>
207 ::type (
new PointerType(*ilu, comm), createParallelDeleter(*ilu, comm));
215 template<
class M,
class X,
class I1,
class I2>
216 typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
217 createILUnPtr(
const M& A,
int ilu_n,
double relax,
218 const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
220 typedef Dune::BlockPreconditioner<
223 Dune::OwnerOverlapCopyCommunication<I1,I2>,
226 Dune::SeqILUn<M,X,X>* ilu =
new Dune::SeqILUn<M,X,X>( A, ilu_n, relax);
228 return typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
229 (
new PointerType(*ilu, comm),createParallelDeleter(*ilu, comm));
236 template<
class M,
class X=
typename M::range_type>
237 std::unique_ptr<Dune::SeqILU0<M,X,X> >
238 createEllipticPreconditionerPointer(
const M& Ae,
double relax,
239 const Dune::Amg::SequentialInformation&)
241 return std::unique_ptr<Dune::SeqILU0<M,X,X> >(
new Dune::SeqILU0<M,X,X>(Ae, relax));
245 template<
class M,
class X=
typename M::range_type,
class I1,
class I2>
250 typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
251 ::EllipticPreconditionerPointer
252 createEllipticPreconditionerPointer(
const M& Ae,
double relax,
253 const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
255 typedef Dune::BlockPreconditioner<X, X,
256 Dune::OwnerOverlapCopyCommunication<I1,I2>,
257 Dune::SeqILU0<M,X,X> >
258 ParallelPreconditioner;
260 Dune::SeqILU0<M,X,X>* ilu=
new Dune::SeqILU0<M,X,X>(Ae, relax);
261 typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
262 ::EllipticPreconditionerPointer EllipticPreconditionerPointer;
263 return EllipticPreconditionerPointer(
new ParallelPreconditioner(*ilu, comm),
264 createParallelDeleter(*ilu, comm));
286 cpr_relax_ = param.getDefault(
"cpr_relax", cpr_relax_);
287 cpr_solver_tol_ = param.getDefault(
"cpr_solver_tol", cpr_solver_tol_);
288 cpr_ilu_n_ = param.getDefault(
"cpr_ilu_n", cpr_ilu_n_);
289 cpr_max_ell_iter_ = param.getDefault(
"cpr_max_elliptic_iter",cpr_max_ell_iter_);
290 cpr_use_amg_ = param.getDefault(
"cpr_use_amg", cpr_use_amg_);
291 cpr_use_bicgstab_ = param.getDefault(
"cpr_use_bicgstab", cpr_use_bicgstab_);
292 cpr_solver_verbose_ = param.getDefault(
"cpr_solver_verbose", cpr_solver_verbose_);
298 cpr_solver_tol_ = 1e-2;
300 cpr_max_ell_iter_ = 25;
302 cpr_use_bicgstab_ =
true;
303 cpr_solver_verbose_ =
false;
322 template<
class M,
class X,
class Y,
323 class P=Dune::Amg::SequentialInformation>
344 category = std::is_same<P,Dune::Amg::SequentialInformation>::value?
345 Dune::SolverCategory::sequential:Dune::SolverCategory::overlapping
349 typedef typename CPRSelector<M,X,X,P>::Operator
Operator;
356 typedef typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
360 typedef typename CPRSelector<M,X,X,P>::AMG
AMG;
382 opAe_(CPRSelector<
M,X,Y,P>::makeOperator(
Ae_, commAe)),
407 virtual void pre (X& , Y& )
416 virtual void apply (X& v,
const Y& d)
420 std::copy_n(d.begin(),
de_.size(),
de_.begin());
430 std::copy(ve_.begin(), ve_.end(), v.begin());
435 A_.mmv(v, dmodified_);
437 comm_.copyOwnerToAll(dmodified_, dmodified_);
468 comm_.communicator().rank()==0 ) ? 1 : 0;
471 Dune::InverseOperatorResult result;
474 typedef Dune::ScalarProductChooser<X,ParallelInformation,category>
475 ScalarProductChooser;
477 std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
478 sp(ScalarProductChooser::construct(
commAe_));
484 Dune::BiCGSTABSolver<X> linsolve(*
opAe_, *sp, (*
amg_), tolerance, maxit, verbosity);
485 linsolve.apply(x, de, result);
488 Dune::CGSolver<X> linsolve(*
opAe_, *sp, (*
amg_), tolerance, maxit, verbosity);
489 linsolve.apply(x, de, result);
497 Dune::BiCGSTABSolver<X> linsolve(*
opAe_, *sp, (*
precond_), tolerance, maxit, verbosity);
498 linsolve.apply(x, de, result);
501 Dune::CGSolver<X> linsolve(*
opAe_, *sp, (*
precond_), tolerance, maxit, verbosity);
502 linsolve.apply(x, de, result);
507 if (!result.converged) {
508 OPM_THROW(LinearSolverProblem,
"CPRPreconditioner failed to solve elliptic subsystem.");
516 const matrix_type&
A_;
537 std::shared_ptr< WholeSystemPreconditioner >
pre_;
553 typedef Dune::Amg::FirstDiagonal CouplingMetric;
556 typedef Dune::Amg::SymmetricCriterion<M, CouplingMetric> CritBase;
559 typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
562 int coarsenTarget=1200;
563 Criterion criterion(15,coarsenTarget);
564 criterion.setDebugLevel( 0 );
565 criterion.setDefaultValuesIsotropic(2);
566 criterion.setNoPostSmoothSteps( 1 );
567 criterion.setNoPreSmoothSteps( 1 );
570 typedef typename AMG::Smoother Smoother;
571 typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
572 SmootherArgs smootherArgs;
573 smootherArgs.iterations = 1;
574 smootherArgs.relaxationFactor = param_.
cpr_relax_;
576 amg_ = std::unique_ptr< AMG > (
new AMG(*opAe_, criterion, smootherArgs, comm ));
580 precond_ = createEllipticPreconditionerPointer<M,X>(
Ae_, param_.
cpr_relax_, comm);
588 #endif // OPM_CPRPRECONDITIONER_HEADER_INCLUDED
Dune::Preconditioner< X, X > WholeSystemPreconditioner
preconditioner for the whole system (here either ILU(0) or ILU(n)
Definition: CPRPreconditioner.hpp:352
virtual void apply(X &v, const Y &d)
Apply the preconditoner.
Definition: CPRPreconditioner.hpp:416
int cpr_ilu_n_
Definition: CPRPreconditioner.hpp:273
bool cpr_use_amg_
Definition: CPRPreconditioner.hpp:275
Definition: CPRPreconditioner.hpp:269
Dune::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: CPRPreconditioner.hpp:333
Definition: AdditionalObjectDeleter.hpp:22
virtual void post(X &)
Clean up.
Definition: CPRPreconditioner.hpp:457
CPRSelector< M, X, X, P >::Operator Operator
Elliptic Operator.
Definition: CPRPreconditioner.hpp:349
void solveElliptic(Y &x, Y &de)
Definition: CPRPreconditioner.hpp:462
CPRParameter()
Definition: CPRPreconditioner.hpp:279
const matrix_type & Ae_
The elliptic part of the matrix.
Definition: CPRPreconditioner.hpp:518
bool cpr_solver_verbose_
Definition: CPRPreconditioner.hpp:277
std::unique_ptr< Operator > opAe_
elliptic operator
Definition: CPRPreconditioner.hpp:524
The category the preconditioner is part of.
Definition: CPRPreconditioner.hpp:344
CPRPreconditioner(const CPRParameter ¶m, const M &A, const M &Ae, const ParallelInformation &comm=ParallelInformation(), const ParallelInformation &commAe=ParallelInformation())
Constructor.
Definition: CPRPreconditioner.hpp:373
void reset()
Definition: CPRPreconditioner.hpp:295
Y de_
temporary variables for elliptic solve
Definition: CPRPreconditioner.hpp:521
const P & comm_
The information about the parallelization of the whole system.
Definition: CPRPreconditioner.hpp:543
Y ve_
Definition: CPRPreconditioner.hpp:521
Y vilu_
temporary variables for ILU solve
Definition: CPRPreconditioner.hpp:540
const P & commAe_
The information about the parallelization of the elliptic part of the system.
Definition: CPRPreconditioner.hpp:546
X domain_type
The domain type of the preconditioner.
Definition: CPRPreconditioner.hpp:335
CPRParameter(const parameter::ParameterGroup ¶m)
Definition: CPRPreconditioner.hpp:281
double cpr_relax_
Definition: CPRPreconditioner.hpp:271
std::unique_ptr< AMG > amg_
AMG preconditioner with ILU0 smoother.
Definition: CPRPreconditioner.hpp:529
CPRSelector< M, X, X, P >::EllipticPreconditionerPointer EllipticPreconditionerPointer
type of the unique pointer to the ilu-0 preconditioner used the for the elliptic system ...
Definition: CPRPreconditioner.hpp:357
double cpr_solver_tol_
Definition: CPRPreconditioner.hpp:272
const CPRParameter & param_
Parameter collection for CPR.
Definition: CPRPreconditioner.hpp:513
Definition: AutoDiffMatrix.hpp:43
P ParallelInformation
The type describing the parallel information.
Definition: CPRPreconditioner.hpp:331
bool cpr_use_bicgstab_
Definition: CPRPreconditioner.hpp:276
Y range_type
The range type of the preconditioner.
Definition: CPRPreconditioner.hpp:337
ADB::M M
Definition: BlackoilModelBase_impl.hpp:85
CPRSelector< M, X, X, P >::AMG AMG
amg preconditioner for the elliptic system
Definition: CPRPreconditioner.hpp:360
Y dmodified_
Definition: CPRPreconditioner.hpp:521
X::field_type field_type
The field type of the preconditioner.
Definition: CPRPreconditioner.hpp:339
virtual void pre(X &, Y &)
Prepare the preconditioner.
Definition: CPRPreconditioner.hpp:407
void createEllipticPreconditioner(const bool amg, const P &comm)
Definition: CPRPreconditioner.hpp:548
int cpr_max_ell_iter_
Definition: CPRPreconditioner.hpp:274
std::shared_ptr< WholeSystemPreconditioner > pre_
The preconditioner for the whole system.
Definition: CPRPreconditioner.hpp:537
const matrix_type & A_
The matrix for the full linear problem.
Definition: CPRPreconditioner.hpp:516
CPR preconditioner.
Definition: CPRPreconditioner.hpp:324
EllipticPreconditionerPointer precond_
ILU0 preconditioner for the elliptic system.
Definition: CPRPreconditioner.hpp:527