00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #ifndef OPM_CPRPRECONDITIONER_HEADER_INCLUDED
00025 #define OPM_CPRPRECONDITIONER_HEADER_INCLUDED
00026
00027 #include <memory>
00028 #include <type_traits>
00029
00030 #include <opm/common/utility/platform_dependent/disable_warnings.h>
00031 #include <opm/core/utility/parameters/ParameterGroup.hpp>
00032
00033 #include <dune/istl/bvector.hh>
00034 #include <dune/istl/bcrsmatrix.hh>
00035 #include <dune/istl/operators.hh>
00036 #include <dune/istl/io.hh>
00037 #include <dune/istl/owneroverlapcopy.hh>
00038 #include <dune/istl/preconditioners.hh>
00039 #include <dune/istl/schwarz.hh>
00040 #include <dune/istl/solvers.hh>
00041 #include <dune/istl/paamg/amg.hh>
00042 #include <dune/istl/paamg/kamg.hh>
00043 #include <dune/istl/paamg/pinfo.hh>
00044
00045 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
00046
00047 #include <opm/common/ErrorMacros.hpp>
00048 #include <opm/common/Exceptions.hpp>
00049
00050 #include <opm/autodiff/AdditionalObjectDeleter.hpp>
00051 #include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
00052 #include <opm/autodiff/ParallelOverlappingILU0.hpp>
00053 namespace Opm
00054 {
00055
00056 namespace ISTLUtility
00057 {
00066 template<class M, class X, class Y, class P>
00067 struct CPRSelector
00068 {
00070 typedef Dune::Amg::SequentialInformation ParallelInformation;
00072 typedef Dune::MatrixAdapter<M, X, Y> Operator;
00074 typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
00076 typedef std::unique_ptr<EllipticPreconditioner> EllipticPreconditionerPointer;
00077
00079 typedef EllipticPreconditioner Smoother;
00080 typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
00081
00085 static Operator* makeOperator(const M& m, const P&)
00086 {
00087 return new Operator(m);
00088 }
00089
00090 };
00091
00092 #if HAVE_MPI
00094 template<class M, class X, class Y, class I1, class I2>
00095 struct CPRSelector<M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
00096 {
00098 typedef Dune::OwnerOverlapCopyCommunication<I1,I2> ParallelInformation;
00100 typedef Dune::OverlappingSchwarzOperator<M,X,X,ParallelInformation> Operator;
00102
00103
00104
00105
00106 typedef ParallelOverlappingILU0<M,X, X, ParallelInformation>
00107 EllipticPreconditioner;
00109 typedef std::unique_ptr<EllipticPreconditioner,
00110 AdditionalObjectDeleter<Dune::SeqILU0<M,X,X> > >
00111 EllipticPreconditionerPointer;
00112
00113 typedef EllipticPreconditioner Smoother;
00114 typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
00115
00119 static Operator* makeOperator(const M& m, const ParallelInformation& p)
00120 {
00121 return new Operator(m, p);
00122 }
00123 };
00124
00131 template<class ILU, class I1, class I2>
00132 AdditionalObjectDeleter<ILU>
00133 createParallelDeleter(ILU& ilu, const Dune::OwnerOverlapCopyCommunication<I1,I2>& p)
00134 {
00135 (void) p;
00136 return AdditionalObjectDeleter<ILU>(ilu);
00137 }
00138 #endif
00139
00143 template<class M, class X>
00144 std::shared_ptr<Dune::SeqILU0<M,X,X> >
00145 createILU0Ptr(const M& A, double relax, const Dune::Amg::SequentialInformation&)
00146 {
00147 return std::shared_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>( A, relax) );
00148 }
00153 template<class M, class X>
00154 std::shared_ptr<Dune::SeqILUn<M,X,X> >
00155 createILUnPtr(const M& A, int ilu_n, double relax, const Dune::Amg::SequentialInformation&)
00156 {
00157 return std::shared_ptr<Dune::SeqILUn<M,X,X> >(new Dune::SeqILUn<M,X,X>( A, ilu_n, relax) );
00158 }
00159
00160 #if HAVE_MPI
00161 template<class ILU, class I1, class I2>
00162 struct SelectParallelILUSharedPtr
00163 {
00164 typedef std::shared_ptr<
00165 Dune::BlockPreconditioner<
00166 typename ILU::range_type,
00167 typename ILU::domain_type,
00168 Dune::OwnerOverlapCopyCommunication<I1,I2>,
00169 ILU
00170 >
00171 > type;
00172 };
00173
00178 template<class M, class X, class I1, class I2>
00179 typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>::type
00180 createILU0Ptr(const M& A, double relax,
00181 const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
00182 {
00183 typedef Dune::BlockPreconditioner<
00184 X,
00185 X,
00186 Dune::OwnerOverlapCopyCommunication<I1,I2>,
00187 Dune::SeqILU0<M,X,X>
00188 > PointerType;
00189 Dune::SeqILU0<M,X,X>* ilu = nullptr;
00190 int ilu_setup_successful = 1;
00191 std::string message;
00192 try {
00193 ilu = new Dune::SeqILU0<M,X,X>(A, relax);
00194 }
00195 catch ( Dune::MatrixBlockError error )
00196 {
00197 message = error.what();
00198 std::cerr<<"Exception occured on process " <<
00199 comm.communicator().rank() << " during " <<
00200 "setup of ILU0 preconditioner with message: " <<
00201 message<<std::endl;
00202 ilu_setup_successful = 0;
00203 }
00204
00205 if ( comm.communicator().min(ilu_setup_successful) == 0 )
00206 {
00207 if ( ilu )
00208 {
00209
00210 delete ilu;
00211 }
00212 throw Dune::MatrixBlockError();
00213 }
00214 return typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>
00215 ::type ( new PointerType(*ilu, comm), createParallelDeleter(*ilu, comm));
00216 }
00217
00223 template<class M, class X, class I1, class I2>
00224 typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
00225 createILUnPtr(const M& A, int ilu_n, double relax,
00226 const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
00227 {
00228 typedef Dune::BlockPreconditioner<
00229 X,
00230 X,
00231 Dune::OwnerOverlapCopyCommunication<I1,I2>,
00232 Dune::SeqILUn<M,X,X>
00233 > PointerType;
00234 Dune::SeqILUn<M,X,X>* ilu = new Dune::SeqILUn<M,X,X>( A, ilu_n, relax);
00235
00236 return typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
00237 (new PointerType(*ilu, comm),createParallelDeleter(*ilu, comm));
00238 }
00239 #endif
00240
00244 template<class M, class X=typename M::range_type>
00245 std::unique_ptr<Dune::SeqILU0<M,X,X> >
00246 createEllipticPreconditionerPointer(const M& Ae, double relax,
00247 const Dune::Amg::SequentialInformation&)
00248 {
00249 return std::unique_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>(Ae, relax));
00250 }
00251
00252 #if HAVE_MPI
00257 template<class M, class X=typename M::range_type, class I1, class I2>
00258 typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
00259 ::EllipticPreconditionerPointer
00260 createEllipticPreconditionerPointer(const M& Ae, double relax,
00261 const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
00262 {
00263 typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
00264 ::EllipticPreconditioner ParallelPreconditioner;
00265
00266
00267 typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
00268 ::EllipticPreconditionerPointer EllipticPreconditionerPointer;
00269 return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax));
00270 }
00271
00272 #endif // #if HAVE_MPI
00273
00278
00279 template < int pressureIndex=0, class Op, class P, class AMG >
00280 inline void
00281 createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std::unique_ptr< AMG >& amgPtr )
00282 {
00283
00284 typedef typename Op::matrix_type M;
00285
00286
00287 typedef Dune::Amg::Diagonal<pressureIndex> CouplingMetric;
00288
00289
00290 typedef Dune::Amg::SymmetricCriterion<M, CouplingMetric> CritBase;
00291
00292
00293 typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
00294
00295
00296 int coarsenTarget=1200;
00297 Criterion criterion(15,coarsenTarget);
00298 criterion.setDebugLevel( 0 );
00299 criterion.setDefaultValuesIsotropic(2);
00300 criterion.setNoPostSmoothSteps( 1 );
00301 criterion.setNoPreSmoothSteps( 1 );
00302
00303
00304 typedef typename AMG::Smoother Smoother;
00305 typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
00306 SmootherArgs smootherArgs;
00307 smootherArgs.iterations = 1;
00308 smootherArgs.relaxationFactor = relax;
00309
00310 amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) );
00311 }
00312
00313 }
00314
00315 struct CPRParameter
00316 {
00317 double cpr_relax_;
00318 double cpr_solver_tol_;
00319 int cpr_ilu_n_;
00320 int cpr_max_ell_iter_;
00321 bool cpr_use_amg_;
00322 bool cpr_use_bicgstab_;
00323 bool cpr_solver_verbose_;
00324
00325 CPRParameter() { reset(); }
00326
00327 CPRParameter( const ParameterGroup& param)
00328 {
00329
00330 reset();
00331
00332 cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_);
00333 cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_);
00334 cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_);
00335 cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_);
00336 cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_);
00337 cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_);
00338 cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_);
00339 }
00340
00341 void reset()
00342 {
00343 cpr_relax_ = 1.0;
00344 cpr_solver_tol_ = 1e-2;
00345 cpr_ilu_n_ = 0;
00346 cpr_max_ell_iter_ = 25;
00347 cpr_use_amg_ = true;
00348 cpr_use_bicgstab_ = true;
00349 cpr_solver_verbose_ = false;
00350 }
00351 };
00352
00353
00368 template<class M, class X, class Y,
00369 class P=Dune::Amg::SequentialInformation>
00370 class CPRPreconditioner : public Dune::Preconditioner<X,Y>
00371 {
00372
00373 CPRPreconditioner( const CPRPreconditioner& );
00374
00375 public:
00377 typedef P ParallelInformation;
00379 typedef typename std::remove_const<M>::type matrix_type;
00381 typedef X domain_type;
00383 typedef Y range_type;
00385 typedef typename X::field_type field_type;
00386
00387
00388 enum {
00390 category = std::is_same<P,Dune::Amg::SequentialInformation>::value?
00391 Dune::SolverCategory::sequential:Dune::SolverCategory::overlapping
00392 };
00393
00394 typedef ISTLUtility::CPRSelector<M,X,X,P> CPRSelectorType ;
00395
00397 typedef typename CPRSelectorType::Operator Operator;
00398
00400 typedef Dune::Preconditioner<X,X> WholeSystemPreconditioner;
00401
00404 typedef typename CPRSelectorType::EllipticPreconditionerPointer
00405 EllipticPreconditionerPointer;
00406
00408 typedef typename CPRSelectorType::AMG AMG;
00409
00421 CPRPreconditioner (const CPRParameter& param, const M& A, const M& Ae,
00422 const ParallelInformation& comm=ParallelInformation(),
00423 const ParallelInformation& commAe=ParallelInformation())
00424 : param_( param ),
00425 A_(A),
00426 Ae_(Ae),
00427 de_( Ae_.N() ),
00428 ve_( Ae_.M() ),
00429 dmodified_( A_.N() ),
00430 opAe_(CPRSelectorType::makeOperator(Ae_, commAe)),
00431 precond_(),
00432 amg_(),
00433 pre_(),
00434 vilu_( A_.N() ),
00435 comm_(comm),
00436 commAe_(commAe)
00437 {
00438
00439 createEllipticPreconditioner( param_.cpr_use_amg_, commAe_ );
00440
00441
00442 if( param_.cpr_ilu_n_ == 0 ) {
00443 pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, param_.cpr_relax_, comm_ );
00444 }
00445 else {
00446 pre_ = ISTLUtility::createILUnPtr<M,X>( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm_ );
00447 }
00448 }
00449
00455 virtual void pre (X& , Y& )
00456 {
00457 }
00458
00464 virtual void apply (X& v, const Y& d)
00465 {
00466
00467
00468 std::copy_n(d.begin(), de_.size(), de_.begin());
00469
00470
00471
00472 ve_ = 0;
00473 solveElliptic( ve_, de_ );
00474
00475
00476 v = 0.0;
00477
00478 std::copy(ve_.begin(), ve_.end(), v.begin());
00479
00480
00481
00482 dmodified_ = d;
00483 A_.mmv(v, dmodified_);
00484
00485 comm_.copyOwnerToAll(dmodified_, dmodified_);
00486
00487
00488 pre_->apply( vilu_, dmodified_);
00489
00490
00491 if( std::abs( param_.cpr_relax_ - 1.0 ) < 1e-12 ) {
00492 v += vilu_;
00493 }
00494 else {
00495 v *= param_.cpr_relax_;
00496 v += vilu_;
00497 }
00498 }
00499
00505 virtual void post (X& )
00506 {
00507 }
00508
00509 protected:
00510 void solveElliptic(Y& x, Y& de)
00511 {
00512
00513 const double tolerance = param_.cpr_solver_tol_;
00514 const int maxit = param_.cpr_max_ell_iter_;
00515 const int verbosity = ( param_.cpr_solver_verbose_ &&
00516 comm_.communicator().rank()==0 ) ? 1 : 0;
00517
00518
00519 Dune::InverseOperatorResult result;
00520
00521
00522 typedef Dune::ScalarProductChooser<X,ParallelInformation,category>
00523 ScalarProductChooser;
00524
00525 std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
00526 sp(ScalarProductChooser::construct(commAe_));
00527
00528 if( amg_ )
00529 {
00530
00531 if( param_.cpr_use_bicgstab_ ) {
00532 Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
00533 linsolve.apply(x, de, result);
00534 }
00535 else {
00536 Dune::CGSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
00537 linsolve.apply(x, de, result);
00538 }
00539 }
00540 else
00541 {
00542 assert( precond_ );
00543
00544 if( param_.cpr_use_bicgstab_ ) {
00545 Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
00546 linsolve.apply(x, de, result);
00547 }
00548 else {
00549 Dune::CGSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
00550 linsolve.apply(x, de, result);
00551 }
00552
00553 }
00554
00555 if (!result.converged) {
00556 OPM_THROW(LinearSolverProblem, "CPRPreconditioner failed to solve elliptic subsystem.");
00557 }
00558 }
00559
00561 const CPRParameter& param_;
00562
00564 const matrix_type& A_;
00566 const matrix_type& Ae_;
00567
00569 Y de_, ve_, dmodified_;
00570
00572 std::unique_ptr<Operator> opAe_;
00573
00575 EllipticPreconditionerPointer precond_;
00577 std::unique_ptr< AMG > amg_;
00578
00585 std::shared_ptr< WholeSystemPreconditioner > pre_;
00586
00588 Y vilu_;
00589
00591 const P& comm_;
00594 const P& commAe_;
00595 protected:
00596 void createEllipticPreconditioner( const bool amg, const P& comm )
00597 {
00598 if( amg )
00599 {
00600 ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, comm, amg_ );
00601 }
00602 else
00603 {
00604 precond_ = ISTLUtility::createEllipticPreconditionerPointer<M,X>( Ae_, param_.cpr_relax_, comm);
00605 }
00606 }
00607 };
00608
00609
00610 }
00611
00612 #endif // OPM_CPRPRECONDITIONER_HEADER_INCLUDED