All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
CPRPreconditioner.hpp
1 /*
2  Copyright 2014 SINTEF ICT, Applied Mathematics.
3  Copyright 2014 IRIS AS.
4  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
5  Copyright 2015 NTNU
6  Copyright 2015 Statoil AS
7 
8  This file is part of the Open Porous Media project (OPM).
9 
10  OPM is free software: you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  OPM is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with OPM. If not, see <http://www.gnu.org/licenses/>.
22 */
23 
24 #ifndef OPM_CPRPRECONDITIONER_HEADER_INCLUDED
25 #define OPM_CPRPRECONDITIONER_HEADER_INCLUDED
26 
27 #include <memory>
28 #include <type_traits>
29 
30 #include <opm/common/utility/platform_dependent/disable_warnings.h>
31 #include <opm/core/utility/parameters/ParameterGroup.hpp>
32 
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>
44 
45 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
46 
47 #include <opm/common/ErrorMacros.hpp>
48 #include <opm/common/Exceptions.hpp>
49 
50 #include <opm/autodiff/AdditionalObjectDeleter.hpp>
51 #include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
52 #include <opm/autodiff/ParallelOverlappingILU0.hpp>
53 namespace Opm
54 {
55 
56 namespace ISTLUtility
57 {
66 template<class M, class X, class Y, class P>
68 {
70  typedef Dune::Amg::SequentialInformation ParallelInformation;
72  typedef Dune::MatrixAdapter<M, X, Y> Operator;
74  typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
76  typedef std::unique_ptr<EllipticPreconditioner> EllipticPreconditionerPointer;
77 
80  typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
81 
85  static Operator* makeOperator(const M& m, const P&)
86  {
87  return new Operator(m);
88  }
89 
90 };
91 
92 #if HAVE_MPI
93 template<class M, class X, class Y, class I1, class I2>
95 struct CPRSelector<M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
96 {
98  typedef Dune::OwnerOverlapCopyCommunication<I1,I2> ParallelInformation;
100  typedef Dune::OverlappingSchwarzOperator<M,X,X,ParallelInformation> Operator;
102  //typedef Dune::BlockPreconditioner<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
103  //EllipticPreconditioner;
104  //typedef ParallelRestrictedOverlappingSchwarz<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
105  //EllipticPreconditioner;
109  typedef std::unique_ptr<EllipticPreconditioner,
112 
114  typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
115 
119  static Operator* makeOperator(const M& m, const ParallelInformation& p)
120  {
121  return new Operator(m, p);
122  }
123 };
124 
131 template<class ILU, class I1, class I2>
132 AdditionalObjectDeleter<ILU>
133 createParallelDeleter(ILU& ilu, const Dune::OwnerOverlapCopyCommunication<I1,I2>& p)
134  {
135  (void) p;
136  return AdditionalObjectDeleter<ILU>(ilu);
137  }
138 #endif
139 
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&)
146 {
147  return std::shared_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>( A, relax) );
148 }
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&)
156 {
157  return std::shared_ptr<Dune::SeqILUn<M,X,X> >(new Dune::SeqILUn<M,X,X>( A, ilu_n, relax) );
158 }
159 
160 #if HAVE_MPI
161 template<class ILU, class I1, class I2>
162 struct SelectParallelILUSharedPtr
163 {
164  typedef std::shared_ptr<
165  Dune::BlockPreconditioner<
166  typename ILU::range_type,
167  typename ILU::domain_type,
168  Dune::OwnerOverlapCopyCommunication<I1,I2>,
169  ILU
170  >
171  > type;
172 };
173 
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)
182 {
183  typedef Dune::BlockPreconditioner<
184  X,
185  X,
186  Dune::OwnerOverlapCopyCommunication<I1,I2>,
187  Dune::SeqILU0<M,X,X>
188  > PointerType;
189  Dune::SeqILU0<M,X,X>* ilu = nullptr;
190  int ilu_setup_successful = 1;
191  std::string message;
192  try {
193  ilu = new Dune::SeqILU0<M,X,X>(A, relax);
194  }
195  catch ( Dune::MatrixBlockError error )
196  {
197  message = error.what();
198  std::cerr<<"Exception occured on process " <<
199  comm.communicator().rank() << " during " <<
200  "setup of ILU0 preconditioner with message: " <<
201  message<<std::endl;
202  ilu_setup_successful = 0;
203  }
204  // Check whether there was a problem on some process
205  if ( comm.communicator().min(ilu_setup_successful) == 0 )
206  {
207  if ( ilu ) // not null if constructor succeeded
208  {
209  // prevent memory leak
210  delete ilu;
211  }
212  throw Dune::MatrixBlockError();
213  }
214  return typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>
215  ::type ( new PointerType(*ilu, comm), createParallelDeleter(*ilu, comm));
216 }
217 
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)
227 {
228  typedef Dune::BlockPreconditioner<
229  X,
230  X,
231  Dune::OwnerOverlapCopyCommunication<I1,I2>,
232  Dune::SeqILUn<M,X,X>
233  > PointerType;
234  Dune::SeqILUn<M,X,X>* ilu = new Dune::SeqILUn<M,X,X>( A, ilu_n, relax);
235 
236  return typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
237  (new PointerType(*ilu, comm),createParallelDeleter(*ilu, comm));
238 }
239 #endif
240 
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&)
248 {
249  return std::unique_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>(Ae, relax));
250 }
251 
252 #if HAVE_MPI
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)
262 {
263  typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
264  ::EllipticPreconditioner ParallelPreconditioner;
265 
266  //Dune::SeqILU0<M,X,X>* ilu=new Dune::SeqILU0<M,X,X>(Ae, relax);
267  typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
268  ::EllipticPreconditionerPointer EllipticPreconditionerPointer;
269  return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax));
270 }
271 
272 #endif // #if HAVE_MPI
273 
278 // \param amgPtr The unique_ptr to be filled (return)
279 template < int pressureIndex=0, class Op, class P, class AMG >
280 inline void
281 createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std::unique_ptr< AMG >& amgPtr )
282 {
283  // type of matrix
284  typedef typename Op::matrix_type M;
285 
286  // The coupling metric used in the AMG
287  typedef Dune::Amg::Diagonal<pressureIndex> CouplingMetric;
288 
289  // The coupling criterion used in the AMG
290  typedef Dune::Amg::SymmetricCriterion<M, CouplingMetric> CritBase;
291 
292  // The coarsening criterion used in the AMG
293  typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
294 
295  // TODO: revise choice of parameters
296  int coarsenTarget=1200;
297  Criterion criterion(15,coarsenTarget);
298  criterion.setDebugLevel( 0 ); // no debug information, 1 for printing hierarchy information
299  criterion.setDefaultValuesIsotropic(2);
300  criterion.setNoPostSmoothSteps( 1 );
301  criterion.setNoPreSmoothSteps( 1 );
302 
303  // for DUNE 2.2 we also need to pass the smoother args
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;
309 
310  amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) );
311 }
312 
313 } // end namespace ISTLUtility
314 
316  {
317  double cpr_relax_;
318  double cpr_solver_tol_;
319  int cpr_ilu_n_;
320  int cpr_max_ell_iter_;
321  bool cpr_use_amg_;
322  bool cpr_use_bicgstab_;
323  bool cpr_solver_verbose_;
324 
325  CPRParameter() { reset(); }
326 
327  CPRParameter( const ParameterGroup& param)
328  {
329  // reset values to default
330  reset();
331 
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_);
339  }
340 
341  void reset()
342  {
343  cpr_relax_ = 1.0;
344  cpr_solver_tol_ = 1e-2;
345  cpr_ilu_n_ = 0;
346  cpr_max_ell_iter_ = 25;
347  cpr_use_amg_ = true;
348  cpr_use_bicgstab_ = true;
349  cpr_solver_verbose_ = false;
350  }
351  };
352 
353 
368  template<class M, class X, class Y,
369  class P=Dune::Amg::SequentialInformation>
370  class CPRPreconditioner : public Dune::Preconditioner<X,Y>
371  {
372  // prohibit copying for now
374 
375  public:
379  typedef typename std::remove_const<M>::type matrix_type;
381  typedef X domain_type;
383  typedef Y range_type;
385  typedef typename X::field_type field_type;
386 
387  // define the category
388  enum {
390  category = std::is_same<P,Dune::Amg::SequentialInformation>::value?
391  Dune::SolverCategory::sequential:Dune::SolverCategory::overlapping
392  };
393 
394  typedef ISTLUtility::CPRSelector<M,X,X,P> CPRSelectorType ;
395 
398 
400  typedef Dune::Preconditioner<X,X> WholeSystemPreconditioner;
401 
406 
408  typedef typename CPRSelectorType::AMG AMG;
409 
421  CPRPreconditioner (const CPRParameter& param, const M& A, const M& Ae,
424  : param_( param ),
425  A_(A),
426  Ae_(Ae),
427  de_( Ae_.N() ),
428  ve_( Ae_.M() ),
429  dmodified_( A_.N() ),
430  opAe_(CPRSelectorType::makeOperator(Ae_, commAe)),
431  precond_(), // ilu0 preconditioner for elliptic system
432  amg_(), // amg preconditioner for elliptic system
433  pre_(), // copy A will be made be the preconditioner
434  vilu_( A_.N() ),
435  comm_(comm),
436  commAe_(commAe)
437  {
438  // create appropriate preconditioner for elliptic system
439  createEllipticPreconditioner( param_.cpr_use_amg_, commAe_ );
440 
441  // create the preconditioner for the whole system.
442  if( param_.cpr_ilu_n_ == 0 ) {
443  pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, param_.cpr_relax_, comm_ );
444  }
445  else {
446  pre_ = ISTLUtility::createILUnPtr<M,X>( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm_ );
447  }
448  }
449 
455  virtual void pre (X& /*x*/, Y& /*b*/)
456  {
457  }
458 
464  virtual void apply (X& v, const Y& d)
465  {
466  // Extract part of d corresponding to elliptic part.
467  // Note: Assumes that the elliptic part comes first.
468  std::copy_n(d.begin(), de_.size(), de_.begin());
469 
470  // Solve elliptic part, extend solution to full.
471  // reset result
472  ve_ = 0;
473  solveElliptic( ve_, de_ );
474 
475  //reset return value
476  v = 0.0;
477  // Again assuming that the elliptic part comes first.
478  std::copy(ve_.begin(), ve_.end(), v.begin());
479 
480  // Subtract elliptic residual from initial residual.
481  // dmodified = d - A * vfull
482  dmodified_ = d;
483  A_.mmv(v, dmodified_);
484  // A is not parallel, do communication manually.
485  comm_.copyOwnerToAll(dmodified_, dmodified_);
486 
487  // Apply Preconditioner for whole system (relax will be applied already)
488  pre_->apply( vilu_, dmodified_);
489 
490  // don't apply relaxation if relax_ == 1
491  if( std::abs( param_.cpr_relax_ - 1.0 ) < 1e-12 ) {
492  v += vilu_;
493  }
494  else {
495  v *= param_.cpr_relax_;
496  v += vilu_;
497  }
498  }
499 
505  virtual void post (X& /*x*/)
506  {
507  }
508 
509  protected:
510  void solveElliptic(Y& x, Y& de)
511  {
512  // Linear solver parameters
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;
517 
518  // operator result containing iterations etc.
519  Dune::InverseOperatorResult result;
520 
521  // the scalar product chooser
522  typedef Dune::ScalarProductChooser<X,ParallelInformation,category>
523  ScalarProductChooser;
524  // the scalar product.
525  std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
526  sp(ScalarProductChooser::construct(commAe_));
527 
528  if( amg_ )
529  {
530  // Solve system with AMG
531  if( param_.cpr_use_bicgstab_ ) {
532  Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
533  linsolve.apply(x, de, result);
534  }
535  else {
536  Dune::CGSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
537  linsolve.apply(x, de, result);
538  }
539  }
540  else
541  {
542  assert( precond_ );
543  // Solve system with ILU-0
544  if( param_.cpr_use_bicgstab_ ) {
545  Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
546  linsolve.apply(x, de, result);
547  }
548  else {
549  Dune::CGSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
550  linsolve.apply(x, de, result);
551  }
552 
553  }
554 
555  if (!result.converged) {
556  OPM_THROW(LinearSolverProblem, "CPRPreconditioner failed to solve elliptic subsystem.");
557  }
558  }
559 
562 
564  const matrix_type& A_;
566  const matrix_type& Ae_;
567 
569  Y de_, ve_, dmodified_;
570 
572  std::unique_ptr<Operator> opAe_;
573 
575  EllipticPreconditionerPointer precond_;
577  std::unique_ptr< AMG > amg_;
578 
585  std::shared_ptr< WholeSystemPreconditioner > pre_;
586 
588  Y vilu_;
589 
591  const P& comm_;
594  const P& commAe_;
595  protected:
596  void createEllipticPreconditioner( const bool amg, const P& comm )
597  {
598  if( amg )
599  {
600  ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, comm, amg_ );
601  }
602  else
603  {
604  precond_ = ISTLUtility::createEllipticPreconditionerPointer<M,X>( Ae_, param_.cpr_relax_, comm);
605  }
606  }
607  };
608 
609 
610 } // namespace Opm
611 
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 &param, 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