diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 6e3cc752f..d8851bc26 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -164,7 +164,7 @@ namespace Opm Dune::InverseOperatorResult result; // Parallel version is deactivated until we figure out how to do it properly. -#if 0 // HAVE_MPI +#if HAVE_MPI if (parallelInformation_.type() == typeid(ParallelISTLInformation)) { typedef Dune::OwnerOverlapCopyCommunication Comm; diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp index eb46ca461..f6aacfb69 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -27,6 +27,7 @@ #include #include +#include #include @@ -34,6 +35,7 @@ #include #include #include +#include #include #include @@ -87,33 +89,66 @@ namespace Opm /// \brief construct the CPR preconditioner and the solver. /// \tparam P The type of the parallel information. /// \param parallelInformation the information about the parallelization. - template + template void constructPreconditionerAndSolve(O& opA, Vector& x, Vector& istlb, - const P& parallelInformation, + const POrComm& parallelInformation, Dune::InverseOperatorResult& result) const { - typedef Dune::ScalarProductChooser ScalarProductChooser; - std::unique_ptr - sp(ScalarProductChooser::construct(parallelInformation)); - // Construct preconditioner. - typedef Dune::SeqILU0 Preconditioner; - // typedef Opm::CPRPreconditioner Preconditioner; - parallelInformation.copyOwnerToAll(istlb, istlb); - const double relax = 1.0; - Preconditioner precond(opA.getmat(), relax); + // Construct scalar product. + typedef Dune::ScalarProductChooser ScalarProductChooser; + auto sp = ScalarProductChooser::construct(parallelInformation); + // Construct preconditioner. + auto precond = constructPrecond(opA, parallelInformation); + + // Communicate if parallel. + parallelInformation.copyOwnerToAll(istlb, istlb); + + // Solve. + solve(opA, x, istlb, *sp, precond, result); + } + + + typedef Dune::SeqILU0 SeqPreconditioner; + typedef Dune::OwnerOverlapCopyCommunication Comm; + typedef Dune::BlockPreconditioner ParPreconditioner; + + + template + SeqPreconditioner constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const + { + const double relax = 1.0; + SeqPreconditioner precond(opA.getmat(), relax); + return precond; + } + + + template + ParPreconditioner constructPrecond(Operator& opA, const Comm& comm) const + { + const double relax = 1.0; + SeqPreconditioner seq_precond(opA.getmat(), relax); + ParPreconditioner precond(seq_precond, comm); + return precond; + } + + + /// \brief Solve the system using the given preconditioner and scalar product. + template + void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const + { // TODO: Revise when linear solvers interface opm-core is done // Construct linear solver. // GMRes solver if ( newton_use_gmres_ ) { - Dune::RestartedGMResSolver linsolve(opA, *sp, precond, + Dune::RestartedGMResSolver linsolve(opA, sp, precond, linear_solver_reduction_, linear_solver_restart_, linear_solver_maxiter_, linear_solver_verbosity_); // Solve system. linsolve.apply(x, istlb, result); } else { // BiCGstab solver - Dune::BiCGSTABSolver linsolve(opA, *sp, precond, + Dune::BiCGSTABSolver linsolve(opA, sp, precond, linear_solver_reduction_, linear_solver_maxiter_, linear_solver_verbosity_); // Solve system. linsolve.apply(x, istlb, result);