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_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED
00025 #define OPM_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED
00026
00027 #include <opm/autodiff/DuneMatrix.hpp>
00028 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
00029 #include <opm/autodiff/CPRPreconditioner.hpp>
00030 #include <opm/core/utility/parameters/ParameterGroup.hpp>
00031 #include <opm/core/linalg/LinearSolverInterface.hpp>
00032 #include <dune/istl/scalarproducts.hh>
00033 #include <dune/istl/operators.hh>
00034 #include <dune/istl/bvector.hh>
00035 #include <memory>
00036
00037 namespace Opm
00038 {
00039
00046 class NewtonIterationBlackoilCPR : public NewtonIterationBlackoilInterface
00047 {
00048 typedef Dune::FieldVector<double, 1 > VectorBlockType;
00049 typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
00050 typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
00051 typedef Dune::BlockVector<VectorBlockType> Vector;
00052
00053 public:
00054
00066 NewtonIterationBlackoilCPR(const ParameterGroup& param,
00067 const boost::any& parallelInformation=boost::any());
00068
00074 virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const;
00075
00077 virtual int iterations () const { return iterations_; }
00078
00080 virtual const boost::any& parallelInformation() const;
00081
00082 private:
00083
00087 template<int category=Dune::SolverCategory::sequential, class O, class P>
00088 void constructPreconditionerAndSolve(O& opA, DuneMatrix& istlAe,
00089 Vector& x, Vector& istlb,
00090 const P& parallelInformation_arg,
00091 const P& parallelInformationAe,
00092 Dune::InverseOperatorResult& result) const
00093 {
00094 typedef Dune::ScalarProductChooser<Vector,P,category> ScalarProductChooser;
00095 std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
00096 sp(ScalarProductChooser::construct(parallelInformation_arg));
00097
00098
00099 typedef Opm::CPRPreconditioner<Mat,Vector,Vector,P> Preconditioner;
00100 parallelInformation_arg.copyOwnerToAll(istlb, istlb);
00101 Preconditioner precond(cpr_param_, opA.getmat(), istlAe, parallelInformation_arg,
00102 parallelInformationAe);
00103
00104
00105
00106
00107 if ( newton_use_gmres_ ) {
00108 Dune::RestartedGMResSolver<Vector> linsolve(opA, *sp, precond,
00109 linear_solver_reduction_, linear_solver_restart_, linear_solver_maxiter_, linear_solver_verbosity_);
00110
00111 linsolve.apply(x, istlb, result);
00112 }
00113 else {
00114 Dune::BiCGSTABSolver<Vector> linsolve(opA, *sp, precond,
00115 linear_solver_reduction_, linear_solver_maxiter_, linear_solver_verbosity_);
00116
00117 linsolve.apply(x, istlb, result);
00118 }
00119 }
00120
00121 CPRParameter cpr_param_;
00122
00123 mutable int iterations_;
00124 boost::any parallelInformation_;
00125
00126 const bool newton_use_gmres_;
00127 const double linear_solver_reduction_;
00128 const int linear_solver_maxiter_;
00129 const int linear_solver_restart_;
00130 const int linear_solver_verbosity_;
00131 const bool linear_solver_ignoreconvergencefailure_;
00132 };
00133
00134 }
00135
00136 #endif // OPM_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED