Make FlexibleSolver feature-complete.

In particular:
 - Add parallel solvers and preconditioners.
 - Add the update() interface to preconditioners, and use it with CPR.

CMakeLists_files.cmake

@@ -169,12 +169,14 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/linalg/ISTLSolverEbosFlexible.hpp
   opm/simulators/linalg/MatrixBlock.hpp
   opm/simulators/linalg/MatrixMarketUtils.hpp
+  opm/simulators/linalg/OwningBlockPreconditioner.hpp
   opm/simulators/linalg/OwningTwolevelPreconditioner.hpp
   opm/simulators/linalg/ParallelOverlappingILU0.hpp
   opm/simulators/linalg/ParallelRestrictedAdditiveSchwarz.hpp
   opm/simulators/linalg/ParallelIstlInformation.hpp
   opm/simulators/linalg/PressureSolverPolicy.hpp
   opm/simulators/linalg/PressureTransferPolicy.hpp
+  opm/simulators/linalg/PreconditionerWithUpdate.hpp
   opm/simulators/linalg/makePreconditioner.hpp
   opm/simulators/linalg/setupPropertyTree.hpp
   opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp

opm/simulators/linalg/BlackoilAmg.hpp

@@ -395,7 +395,7 @@ private:
             }
         }
 
-        void updateAmgPreconditioner(typename AMGType::Operator& op)
+        void updatePreconditioner(typename AMGType::Operator& op)
         {
             amg_->updateSolver(crit_, op, comm_);
         }

opm/simulators/linalg/FlexibleSolver.hpp

@@ -34,17 +34,32 @@ namespace Dune
 {
 
 
+template <class X, class Y>
+class SolverWithUpdate : public Dune::InverseOperator<X, Y>
+{
+public:
+    virtual void updatePreconditioner() = 0;
+};
+
 
 template <class MatrixTypeT, class VectorTypeT>
-class FlexibleSolver : Dune::InverseOperator<VectorTypeT, VectorTypeT>
+class FlexibleSolver : public Dune::SolverWithUpdate<VectorTypeT, VectorTypeT>
 {
 public:
     using MatrixType = MatrixTypeT;
     using VectorType = VectorTypeT;
 
+    /// Create a sequential solver.
     FlexibleSolver(const boost::property_tree::ptree& prm, const MatrixType& matrix)
     {
-        makeSolver(prm, matrix);
+        init(prm, matrix, Dune::Amg::SequentialInformation());
+    }
+
+    /// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
+    template <class Comm>
+    FlexibleSolver(const boost::property_tree::ptree& prm, const MatrixType& matrix, const Comm& comm)
+    {
+        init(prm, matrix, comm);
     }
 
     virtual void apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res) override
@@ -57,28 +72,61 @@ public:
         linsolver_->apply(x, rhs, reduction, res);
     }
 
+    virtual void updatePreconditioner() override
+    {
+        preconditioner_->update();
+    }
+
     virtual Dune::SolverCategory::Category category() const override
     {
-        return Dune::SolverCategory::sequential;
+        return linearoperator_->category();
     }
 
 private:
-    void makeSolver(const boost::property_tree::ptree& prm, const MatrixType& matrix)
+
+    using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
+    using AbstractPrecondType = Dune::PreconditionerWithUpdate<VectorType, VectorType>;
+    using AbstractScalarProductType = Dune::ScalarProduct<VectorType>;
+    using AbstractSolverType = Dune::InverseOperator<VectorType, VectorType>;
+
+    // Machinery for making sequential or parallel operators/preconditioners/scalar products.
+    template <class Comm>
+    void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm, const Comm& comm)
+    {
+        // Parallel case.
+        using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
+        auto linop = std::make_shared<ParOperatorType>(matrix, comm);
+        linearoperator_ = linop;
+        preconditioner_ = Dune::makePreconditioner<ParOperatorType, VectorType, Comm>(*linop, prm, comm);
+        scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, linearoperator_->category());
+    }
+    template <>
+    void initOpPrecSp<Dune::Amg::SequentialInformation>(const MatrixType& matrix, const boost::property_tree::ptree& prm, const Dune::Amg::SequentialInformation&)
+    {
+        // Sequential case.
+        using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
+        auto linop = std::make_shared<SeqOperatorType>(matrix);
+        linearoperator_ = linop;
+        preconditioner_ = Dune::makePreconditioner<SeqOperatorType, VectorType>(*linop, prm);
+        scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
+    }
+
+    void initSolver(const boost::property_tree::ptree& prm)
     {
         const double tol = prm.get<double>("tol");
         const int maxiter = prm.get<int>("maxiter");
-        linearoperator_.reset(new Dune::MatrixAdapter<MatrixType, VectorType, VectorType>(matrix));
+        const int verbosity = prm.get<int>("verbosity");
-        preconditioner_ = Dune::makePreconditioner<MatrixType, VectorType>(*linearoperator_, prm);
+        const std::string solver_type = prm.get<std::string>("solver");
-        int verbosity = prm.get<int>("verbosity");
-        std::string solver_type = prm.get<std::string>("solver");
         if (solver_type == "bicgstab") {
             linsolver_.reset(new Dune::BiCGSTABSolver<VectorType>(*linearoperator_,
+                                                                  *scalarproduct_,
                                                                   *preconditioner_,
                                                                   tol, // desired residual reduction factor
                                                                   maxiter, // maximum number of iterations
                                                                   verbosity));
         } else if (solver_type == "loopsolver") {
             linsolver_.reset(new Dune::LoopSolver<VectorType>(*linearoperator_,
+                                                              *scalarproduct_,
                                                               *preconditioner_,
                                                               tol, // desired residual reduction factor
                                                               maxiter, // maximum number of iterations
@@ -86,6 +134,7 @@ private:
         } else if (solver_type == "gmres") {
             int restart = prm.get<int>("restart");
             linsolver_.reset(new Dune::RestartedGMResSolver<VectorType>(*linearoperator_,
+                                                                        *scalarproduct_,
                                                                         *preconditioner_,
                                                                         tol,
                                                                         restart, // desired residual reduction factor
@@ -103,9 +152,20 @@ private:
         }
     }
 
-    std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> preconditioner_;
+
-    std::shared_ptr<Dune::MatrixAdapter<MatrixType, VectorType, VectorType>> linearoperator_;
+    // Main initialization routine.
-    std::shared_ptr<Dune::InverseOperator<VectorType, VectorType>> linsolver_;
+    // Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
+    template <class Comm>
+    void init(const boost::property_tree::ptree& prm, const MatrixType& matrix, const Comm& comm)
+    {
+        initOpPrecSp(matrix, prm, comm);
+        initSolver(prm);
+    }
+
+    std::shared_ptr<AbstractOperatorType> linearoperator_;
+    std::shared_ptr<AbstractPrecondType> preconditioner_;
+    std::shared_ptr<AbstractScalarProductType> scalarproduct_;
+    std::shared_ptr<AbstractSolverType> linsolver_;
 };
 
 } // namespace Dune

opm/simulators/linalg/ISTLSolverEbosFlexible.hpp

@@ -58,10 +58,12 @@ class ISTLSolverEbosFlexible
     using VectorType = typename GET_PROP_TYPE(TypeTag, GlobalEqVector);
     using Simulator = typename GET_PROP_TYPE(TypeTag, Simulator);
     using MatrixType = typename SparseMatrixAdapter::IstlMatrix;
-    using POrComm = Dune::Amg::SequentialInformation;
+#if HAVE_MPI
-    using MatrixAdapterType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
+    using Communication = Dune::OwnerOverlapCopyCommunication<int, int>;
+#endif
     using SolverType = Dune::FlexibleSolver<MatrixType, VectorType>;
 
+
 public:
     static void registerParameters()
     {
@@ -70,8 +72,19 @@ public:
 
     explicit ISTLSolverEbosFlexible(const Simulator& simulator)
         : simulator_(simulator)
+        , comm_(nullptr)
     {
         parameters_.template init<TypeTag>();
+        prm_ = setupPropertyTree(parameters_);
+        extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
+#if HAVE_MPI
+        if (parallelInformation_.type() == typeid(ParallelISTLInformation)) {
+            // Parallel case.
+            const ParallelISTLInformation* parinfo = boost::any_cast<ParallelISTLInformation>(&parallelInformation_);
+            assert(parinfo);
+            comm_.reset(new Communication(parinfo->communicator()));
+        }
+#endif
     }
 
     void eraseMatrix()
@@ -80,9 +93,47 @@ public:
 
     void prepare(const SparseMatrixAdapter& mat, VectorType& b)
     {
-        boost::property_tree::ptree prm = setupPropertyTree(parameters_);
+#if HAVE_MPI
-        solver_.reset(new SolverType(prm, mat.istlMatrix()));
+        static bool firstcall = true;
-        rhs_ = b;
+        if (firstcall && parallelInformation_.type() == typeid(ParallelISTLInformation)) {
+            // Parallel case.
+            const ParallelISTLInformation* parinfo = boost::any_cast<ParallelISTLInformation>(&parallelInformation_);
+            assert(parinfo);
+            const size_t size = mat.istlMatrix().N();
+            parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1);
+            firstcall = false;
+        }
+#endif
+        // Decide if we should recreate the solver or just do
+        // a minimal preconditioner update.
+        const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
+        bool recreate_solver = false;
+        if (this->parameters_.cpr_reuse_setup_ == 0) {
+            recreate_solver = true;
+        } else if (this->parameters_.cpr_reuse_setup_ == 1) {
+            if (newton_iteration == 0) {
+                recreate_solver = true;
+            }
+        } else if (this->parameters_.cpr_reuse_setup_ == 2) {
+            if (this->iterations() > 10) {
+                recreate_solver = true;
+            }
+        } else {
+            assert(this->parameters_.cpr_reuse_setup_ == 3);
+            assert(recreate_solver == false);
+        }
+
+        if (recreate_solver || !solver_) {
+            if (isParallel()) {
+                solver_.reset(new SolverType(prm_, mat.istlMatrix(), *comm_));
+            } else {
+                solver_.reset(new SolverType(prm_, mat.istlMatrix()));
+            }
+            rhs_ = b;
+        } else {
+            solver_->updatePreconditioner();
+            rhs_ = b;
+        }
     }
 
     bool solve(VectorType& x)
@@ -91,9 +142,13 @@ public:
         return res_.converged;
     }
 
-    bool isParallel()
+    bool isParallel() const
     {
+#if HAVE_MPI
+        return parallelInformation_.type() == typeid(ParallelISTLInformation);
+#else
         return false;
+#endif
     }
 
     int iterations() const
@@ -116,8 +171,11 @@ protected:
 
     std::unique_ptr<SolverType> solver_;
     FlowLinearSolverParameters parameters_;
+    boost::property_tree::ptree prm_;
     VectorType rhs_;
     Dune::InverseOperatorResult res_;
+    boost::any parallelInformation_;
+    std::unique_ptr<Communication> comm_;
 }; // end ISTLSolverEbosFlexible
 
 } // namespace Opm

opm/simulators/linalg/MatrixMarketUtils.hpp

@@ -22,6 +22,7 @@
 #define OPM_MATRIXMARKETUTILS_HEADER_INCLUDED
 
 #include <dune/istl/matrixmarket.hh>
+
 #include <cstddef>
 
 namespace Dune

opm/simulators/linalg/OwningBlockPreconditioner.hpp

@@ -0,0 +1,86 @@
+/*
+  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_OWNINGBLOCKPRECONDITIONER_HEADER_INCLUDED
+#define OPM_OWNINGBLOCKPRECONDITIONER_HEADER_INCLUDED
+
+#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
+
+#include <dune/istl/schwarz.hh>
+
+namespace Dune
+{
+
+template <class OriginalPreconditioner, class Comm>
+class OwningBlockPreconditioner : public PreconditionerWithUpdate<typename OriginalPreconditioner::domain_type,
+                                                                  typename OriginalPreconditioner::range_type>
+{
+public:
+    template <class... Args>
+    OwningBlockPreconditioner(const Comm& comm, Args&&... args)
+        : orig_precond_(std::forward<Args>(args)...)
+        , block_precond_(orig_precond_, comm)
+    {
+    }
+
+    using X = typename OriginalPreconditioner::domain_type;
+    using Y = typename OriginalPreconditioner::range_type;
+
+    virtual void pre(X& x, Y& b) override
+    {
+        block_precond_.pre(x, b);
+    }
+
+    virtual void apply(X& v, const Y& d) override
+    {
+        block_precond_.apply(v, d);
+    }
+
+    virtual void post(X& x) override
+    {
+        block_precond_.post(x);
+    }
+
+    virtual SolverCategory::Category category() const override
+    {
+        return block_precond_.category();
+    }
+
+    // The update() function does nothing for a wrapped preconditioner.
+    virtual void update() override
+    {
+        orig_precond_.update();
+    }
+
+private:
+    OriginalPreconditioner orig_precond_;
+    BlockPreconditioner<X, Y, Comm, OriginalPreconditioner> block_precond_;
+};
+
+
+template <class OriginalPreconditioner, class Comm, class... Args>
+std::shared_ptr<OwningBlockPreconditioner<OriginalPreconditioner, Comm>>
+wrapBlockPreconditioner(const Comm& comm, Args&&... args)
+{
+    return std::make_shared<OwningBlockPreconditioner<OriginalPreconditioner, Comm>>(comm, std::forward<Args>(args)...);
+}
+
+} // namespace Dune
+
+#endif // OPM_OWNINGBLOCKPRECONDITIONER_HEADER_INCLUDED

opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp

@@ -20,9 +20,11 @@
 #ifndef OPM_OWNINGTWOLEVELPRECONDITIONER_HEADER_INCLUDED
 #define OPM_OWNINGTWOLEVELPRECONDITIONER_HEADER_INCLUDED
 
+#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
 #include <opm/simulators/linalg/PressureSolverPolicy.hpp>
 #include <opm/simulators/linalg/PressureTransferPolicy.hpp>
 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
+#include <opm/simulators/linalg/twolevelmethodcpr.hh>
 
 #include <dune/common/fmatrix.hh>
 #include <dune/istl/bcrsmatrix.hh>
@@ -31,6 +33,7 @@
 #include <boost/property_tree/ptree.hpp>
 
 #include <fstream>
+#include <type_traits>
 
 namespace Dune
 {
@@ -38,30 +41,36 @@ namespace Dune
 // Circular dependency between makePreconditioner() [which can make an OwningTwoLevelPreconditioner]
 // and OwningTwoLevelPreconditioner [which uses makePreconditioner() to choose the fine-level smoother]
 // must be broken, accomplished by forward-declaration here.
-template <class MatrixType, class VectorType>
+template <class OperatorType, class VectorType>
-std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>>
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
-makePreconditioner(const Dune::MatrixAdapter<MatrixType, VectorType, VectorType>& linearoperator,
+makePreconditioner(OperatorType& linearoperator, const boost::property_tree::ptree& prm);
-                   const boost::property_tree::ptree& prm);
+
+template <class OperatorType, class VectorType, class Comm>
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
+makePreconditioner(OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm);
+
 
 // Must forward-declare FlexibleSolver as we want to use it as solver for the pressure system.
 template <class MatrixTypeT, class VectorTypeT>
 class FlexibleSolver;
 
-template <class MatrixTypeT, class VectorTypeT, bool transpose = false>
+template <class OperatorType,
-class OwningTwoLevelPreconditioner : public Dune::Preconditioner<VectorTypeT, VectorTypeT>
+          class VectorType,
+          bool transpose = false,
+          class Communication = Dune::Amg::SequentialInformation>
+class OwningTwoLevelPreconditioner : public Dune::PreconditionerWithUpdate<VectorType, VectorType>
 {
 public:
     using pt = boost::property_tree::ptree;
-    using MatrixType = MatrixTypeT;
+    using MatrixType = typename OperatorType::matrix_type;
-    using VectorType = VectorTypeT;
-    using OperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
 
     OwningTwoLevelPreconditioner(OperatorType& linearoperator, pt& prm)
-        : finesmoother_(makePreconditioner<MatrixType, VectorType>(linearoperator, prm.get_child("finesmoother")))
+        : linear_operator_(linearoperator)
-        , comm_()
+        , finesmoother_(makePreconditioner<OperatorType, VectorType>(linearoperator, prm.get_child("finesmoother")))
+        , comm_(nullptr)
         , weights_(Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
               linearoperator.getmat(), prm.get<int>("pressure_var_index"), transpose))
-        , levelTransferPolicy_(comm_, weights_, prm.get<int>("pressure_var_index"))
+        , levelTransferPolicy_(*comm_, weights_, prm.get<int>("pressure_var_index"))
         , coarseSolverPolicy_(prm.get_child("coarsesolver"))
         , twolevel_method_(linearoperator,
                            finesmoother_,
@@ -69,6 +78,33 @@ public:
                            coarseSolverPolicy_,
                            transpose ? 1 : 0,
                            transpose ? 0 : 1)
+        , prm_(prm)
+    {
+        if (prm.get<int>("verbosity") > 10) {
+            std::ofstream outfile(prm.get<std::string>("weights_filename"));
+            if (!outfile) {
+                throw std::runtime_error("Could not write weights");
+            }
+            Dune::writeMatrixMarket(weights_, outfile);
+        }
+    }
+
+    OwningTwoLevelPreconditioner(OperatorType& linearoperator, pt& prm, const Communication& comm)
+        : linear_operator_(linearoperator)
+        , finesmoother_(makePreconditioner<OperatorType, VectorType, Communication>(
+              linearoperator, prm.get_child("finesmoother"), comm))
+        , comm_(&comm)
+        , weights_(Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
+              linearoperator.getmat(), prm.get<int>("pressure_var_index"), transpose))
+        , levelTransferPolicy_(*comm_, weights_, prm.get<int>("pressure_var_index"))
+        , coarseSolverPolicy_(prm.get_child("coarsesolver"))
+        , twolevel_method_(linearoperator,
+                           finesmoother_,
+                           levelTransferPolicy_,
+                           coarseSolverPolicy_,
+                           transpose ? 1 : 0,
+                           transpose ? 0 : 1)
+        , prm_(prm)
     {
         if (prm.get<int>("verbosity") > 10) {
             std::ofstream outfile(prm.get<std::string>("weights_filename"));
@@ -83,36 +119,71 @@ public:
     {
         twolevel_method_.pre(x, b);
     }
+
     virtual void apply(VectorType& v, const VectorType& d) override
     {
         twolevel_method_.apply(v, d);
     }
+
     virtual void post(VectorType& x) override
     {
         twolevel_method_.post(x);
     }
+
+    virtual void update() override
+    {
+        Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
+            linear_operator_.getmat(), prm_.get<int>("pressure_var_index"), transpose, weights_);
+        updateImpl<Communication>();
+    }
+
     virtual Dune::SolverCategory::Category category() const override
     {
-        return Dune::SolverCategory::sequential;
+        return linear_operator_.category();
     }
 
 private:
     using PressureMatrixType = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
     using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
-    using CoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType, PressureVectorType, PressureVectorType>;
+    using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType, PressureVectorType, PressureVectorType>;
-    using Communication = Dune::Amg::SequentialInformation;
+    using ParCoarseOperatorType
+        = Dune::OverlappingSchwarzOperator<PressureMatrixType, PressureVectorType, PressureVectorType, Communication>;
+    using CoarseOperatorType = std::conditional_t<std::is_same<Communication, Dune::Amg::SequentialInformation>::value,
+                                                  SeqCoarseOperatorType,
+                                                  ParCoarseOperatorType>;
     using LevelTransferPolicy = Opm::PressureTransferPolicy<OperatorType, CoarseOperatorType, Communication, transpose>;
-    using CoarseSolverPolicy
+    using CoarseSolverPolicy = Dune::Amg::PressureSolverPolicy<CoarseOperatorType,
-        = Dune::Amg::PressureSolverPolicy<CoarseOperatorType, FlexibleSolver<PressureMatrixType, PressureVectorType>>;
+                                                               FlexibleSolver<PressureMatrixType, PressureVectorType>,
+                                                               LevelTransferPolicy>;
     using TwoLevelMethod
-        = Dune::Amg::TwoLevelMethod<OperatorType, CoarseSolverPolicy, Dune::Preconditioner<VectorType, VectorType>>;
+        = Dune::Amg::TwoLevelMethodCpr<OperatorType, CoarseSolverPolicy, Dune::Preconditioner<VectorType, VectorType>>;
 
+    // Handling parallel vs serial instantiation of makePreconditioner().
+    template <class Comm>
+    void updateImpl()
+    {
+        // Parallel case.
+        finesmoother_ = makePreconditioner<OperatorType, VectorType, Communication>(
+            linear_operator_, prm_.get_child("finesmoother"), *comm_);
+        twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
+    }
+
+    template <>
+    void updateImpl<Dune::Amg::SequentialInformation>()
+    {
+        // Serial case.
+        finesmoother_ = makePreconditioner<OperatorType, VectorType>(linear_operator_, prm_.get_child("finesmoother"));
+        twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
+    }
+
+    OperatorType& linear_operator_;
     std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> finesmoother_;
-    Communication comm_;
+    const Communication* comm_;
     VectorType weights_;
     LevelTransferPolicy levelTransferPolicy_;
     CoarseSolverPolicy coarseSolverPolicy_;
     TwoLevelMethod twolevel_method_;
+    boost::property_tree::ptree prm_;
 };
 
 } // namespace Dune

opm/simulators/linalg/PreconditionerWithUpdate.hpp

@@ -0,0 +1,91 @@
+/*
+  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_PRECONDITIONERWITHUPDATE_HEADER_INCLUDED
+#define OPM_PRECONDITIONERWITHUPDATE_HEADER_INCLUDED
+
+#include <dune/istl/preconditioner.hh>
+#include <memory>
+
+namespace Dune
+{
+
+/// Interface class adding the update() method to the preconditioner interface.
+template <class X, class Y>
+class PreconditionerWithUpdate : public Preconditioner<X, Y>
+{
+public:
+    virtual void update() = 0;
+};
+
+template <class OriginalPreconditioner>
+class DummyUpdatePreconditioner : public PreconditionerWithUpdate<typename OriginalPreconditioner::domain_type,
+                                                                  typename OriginalPreconditioner::range_type>
+{
+public:
+    template <class... Args>
+    DummyUpdatePreconditioner(Args&&... args)
+        : orig_precond_(std::forward<Args>(args)...)
+    {
+    }
+
+    using X = typename OriginalPreconditioner::domain_type;
+    using Y = typename OriginalPreconditioner::range_type;
+
+    virtual void pre (X& x, Y& b) override
+    {
+        orig_precond_.pre(x, b);
+    }
+
+    virtual void apply (X& v, const Y& d) override
+    {
+        orig_precond_.apply(v, d);
+    }
+
+    virtual void post (X& x) override
+    {
+        orig_precond_.post(x);
+    }
+
+    virtual SolverCategory::Category category() const override
+    {
+        return orig_precond_.category();
+    }
+
+    // The update() function does nothing for a wrapped preconditioner.
+    virtual void update() override
+    {
+    }
+
+private:
+    OriginalPreconditioner orig_precond_;
+};
+
+
+template <class OriginalPreconditioner, class... Args>
+std::shared_ptr<DummyUpdatePreconditioner<OriginalPreconditioner>>
+wrapPreconditioner(Args&&... args)
+{
+    return std::make_shared<DummyUpdatePreconditioner<OriginalPreconditioner>>(std::forward<Args>(args)...);
+}
+
+
+} // namespace Dune
+
+#endif // OPM_PRECONDITIONERWITHUPDATE_HEADER_INCLUDED

opm/simulators/linalg/PressureSolverPolicy.hpp

@@ -20,6 +20,8 @@
 #ifndef OPM_PRESSURE_SOLVER_POLICY_HEADER_INCLUDED
 #define OPM_PRESSURE_SOLVER_POLICY_HEADER_INCLUDED
 
+#include <opm/simulators/linalg/PressureTransferPolicy.hpp>
+
 #include <boost/property_tree/ptree.hpp>
 
 #include <dune/istl/solver.hh>
@@ -30,7 +32,7 @@ namespace Amg
 {
     namespace pt = boost::property_tree;
 
-    template <class OperatorType, class Solver>
+    template <class OperatorType, class Solver, class LevelTransferPolicy>
     class PressureSolverPolicy
     {
     public:
@@ -54,10 +56,15 @@ namespace Amg
          * the coarse level system.
          */
         struct PressureInverseOperator : public Dune::InverseOperator<X, X> {
-            PressureInverseOperator(Operator& op, const boost::property_tree::ptree& prm)
+            template <class Comm>
+            PressureInverseOperator(Operator& op, const boost::property_tree::ptree& prm, const Comm& comm)
                 : linsolver_()
             {
-                linsolver_.reset(new Solver(prm, op.getmat()));
+                if (op.category() == Dune::SolverCategory::overlapping) {
+                    linsolver_.reset(new Solver(prm, op.getmat(), comm));
+                } else {
+                    linsolver_.reset(new Solver(prm, op.getmat()));
+                }
             }
 
             Dune::SolverCategory::Category category() const override
@@ -75,6 +82,11 @@ namespace Amg
                 linsolver_->apply(x, b, res);
             }
 
+            void updatePreconditioner()
+            {
+                linsolver_->updatePreconditioner();
+            }
+
         private:
             std::shared_ptr<Solver> linsolver_;
         };
@@ -98,7 +110,9 @@ namespace Amg
         CoarseLevelSolver* createCoarseLevelSolver(LTP& transferPolicy)
         {
             coarseOperator_ = transferPolicy.getCoarseLevelOperator();
-            PressureInverseOperator* inv = new PressureInverseOperator(*coarseOperator_, prm_);
+            auto& tp = dynamic_cast<LevelTransferPolicy&>(transferPolicy); // TODO: make this unnecessary.
+            PressureInverseOperator* inv
+                = new PressureInverseOperator(*coarseOperator_, prm_, tp.getCoarseLevelCommunication());
             return inv;
         }
 

opm/simulators/linalg/PressureTransferPolicy.hpp

@@ -22,17 +22,17 @@
 #define OPM_PRESSURE_TRANSFER_POLICY_HEADER_INCLUDED
 
 
-#include <dune/istl/paamg/twolevelmethod.hh>
+#include <opm/simulators/linalg/twolevelmethodcpr.hh>
 
 
 namespace Opm
 {
 
 template <class FineOperator, class CoarseOperator, class Communication, bool transpose = false>
-class PressureTransferPolicy : public Dune::Amg::LevelTransferPolicy<FineOperator, CoarseOperator>
+class PressureTransferPolicy : public Dune::Amg::LevelTransferPolicyCpr<FineOperator, CoarseOperator>
 {
 public:
-    typedef Dune::Amg::LevelTransferPolicy<FineOperator, CoarseOperator> FatherType;
+    typedef Dune::Amg::LevelTransferPolicyCpr<FineOperator, CoarseOperator> ParentType;
     typedef Communication ParallelInformation;
     typedef typename FineOperator::domain_type FineVectorType;
 
@@ -44,7 +44,7 @@ public:
     {
     }
 
-    void createCoarseLevelSystem(const FineOperator& fineOperator)
+    virtual void createCoarseLevelSystem(const FineOperator& fineOperator) override
     {
         using CoarseMatrix = typename CoarseOperator::matrix_type;
         const auto& fineLevelMatrix = fineOperator.getmat();
@@ -68,7 +68,7 @@ public:
         this->operator_.reset(Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs));
     }
 
-    void calculateCoarseEntries(const FineOperator& fineOperator)
+    virtual void calculateCoarseEntries(const FineOperator& fineOperator) override
     {
         const auto& fineMatrix = fineOperator.getmat();
         *coarseLevelMatrix_ = 0;
@@ -96,7 +96,7 @@ public:
         assert(rowCoarse == coarseLevelMatrix_->end());
     }
 
-    void moveToCoarseLevel(const typename FatherType::FineRangeType& fine)
+    virtual void moveToCoarseLevel(const typename ParentType::FineRangeType& fine) override
     {
         // Set coarse vector to zero
         this->rhs_ = 0;
@@ -119,7 +119,7 @@ public:
         this->lhs_ = 0;
     }
 
-    void moveToFineLevel(typename FatherType::FineDomainType& fine)
+    virtual void moveToFineLevel(typename ParentType::FineDomainType& fine) override
     {
         auto end = fine.end(), begin = fine.begin();
 
@@ -135,7 +135,7 @@ public:
         }
     }
 
-    PressureTransferPolicy* clone() const
+    virtual PressureTransferPolicy* clone() const override
     {
         return new PressureTransferPolicy(*this);
     }

opm/simulators/linalg/amgcpr.hh

@@ -6,7 +6,8 @@
 // NOTE: This file is a modified version of dune/istl/paamg/amg.hh from
 // dune-istl release 2.6.0. Modifications have been kept as minimal as possible.
 
-#include <memory>
+#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
+
 #include <dune/common/exceptions.hh>
 #include <dune/istl/paamg/smoother.hh>
 #include <dune/istl/paamg/transfer.hh>
@@ -19,6 +20,8 @@
 #include <dune/common/typetraits.hh>
 #include <dune/common/exceptions.hh>
 
+#include <memory>
+
 namespace Dune
 {
   namespace Amg
@@ -130,7 +133,7 @@ namespace Dune
      */
     template<class M, class X, class S, class PI=SequentialInformation,
         class A=std::allocator<X> >
-    class AMGCPR : public Preconditioner<X,X>
+    class AMGCPR : public PreconditionerWithUpdate<X,X>
     {
       template<class M1, class X1, class S1, class P1, class K1, class A1>
       friend class KAMG;
@@ -293,6 +296,11 @@ namespace Dune
       template<class C>
       void updateSolver(C& criterion, Operator& /* matrix */, const PI& pinfo);
 
+      /**
+       * @brief Update the coarse solver and the hierarchies.
+       */
+      virtual void update();
+
       /**
        * @brief Check whether the coarse solver used is a direct solver.
        * @return True if the coarse level solver is a direct solver.
@@ -526,6 +534,12 @@ namespace Dune
     template<class M, class X, class S, class PI, class A>
     template<class C>
     void AMGCPR<M,X,S,PI,A>::updateSolver(C& /* criterion */, Operator& /* matrix */, const PI& /* pinfo */)
+    {
+      update();
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    void AMGCPR<M,X,S,PI,A>::update()
     {
       Timer watch;
       smoothers_.reset(new Hierarchy<Smoother,A>);

opm/simulators/linalg/getQuasiImpesWeights.hpp

@@ -43,7 +43,7 @@ namespace Details
 namespace Amg
 {
     template <class Matrix, class Vector>
-    void getQuasiImpesWeights(const Matrix& matrix, const int pressureVarIndex, Vector& weights, bool transpose = false)
+    void getQuasiImpesWeights(const Matrix& matrix, const int pressureVarIndex, const bool transpose, Vector& weights)
     {
         using VectorBlockType = typename Vector::block_type;
         using MatrixBlockType = typename Matrix::block_type;
@@ -76,10 +76,10 @@ namespace Amg
     }
 
     template <class Matrix, class Vector>
-    Vector getQuasiImpesWeights(const Matrix& matrix, const int pressureVarIndex, bool transpose = false)
+    Vector getQuasiImpesWeights(const Matrix& matrix, const int pressureVarIndex, const bool transpose)
     {
         Vector weights(matrix.N());
-        getQuasiImpesWeights(matrix, pressureVarIndex, weights, transpose);
+        getQuasiImpesWeights(matrix, pressureVarIndex, transpose, weights);
         return weights;
     }
 

opm/simulators/linalg/makePreconditioner.hpp

@@ -20,8 +20,11 @@
 #ifndef OPM_MAKEPRECONDITIONER_HEADER_INCLUDED
 #define OPM_MAKEPRECONDITIONER_HEADER_INCLUDED
 
+#include <opm/simulators/linalg/OwningBlockPreconditioner.hpp>
 #include <opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp>
 #include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
+#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
+#include <opm/simulators/linalg/amgcpr.hh>
 
 #include <dune/istl/paamg/amg.hh>
 #include <dune/istl/paamg/fastamg.hh>
@@ -32,49 +35,85 @@
 namespace Dune
 {
 
-template <class MatrixType, class VectorType>
+template <class OperatorType, class VectorType>
-std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>>
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
-makeSeqPreconditioner(const Dune::MatrixAdapter<MatrixType, VectorType, VectorType>& linearoperator,
+makeSeqPreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm)
-                      const boost::property_tree::ptree& prm)
 {
+    using MatrixType = typename OperatorType::matrix_type;
     auto& matrix = linearoperator.getmat();
-    std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> preconditioner;
     double w = prm.get<double>("w");
     int n = prm.get<int>("n");
     std::string precond(prm.get<std::string>("preconditioner"));
     if (precond == "ILU0") {
-        preconditioner.reset(new Dune::SeqILU0<MatrixType, VectorType, VectorType>(matrix, w));
+        return wrapPreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>>(matrix, w);
     } else if (precond == "ParOverILU0") {
-        preconditioner.reset(new Opm::ParallelOverlappingILU0<MatrixType, VectorType, VectorType>(matrix, n, w, Opm::MILU_VARIANT::ILU));
+        return wrapPreconditioner<Opm::ParallelOverlappingILU0<MatrixType, VectorType, VectorType>>(
+            matrix, n, w, Opm::MILU_VARIANT::ILU);
     } else if (precond == "Jac") {
-        preconditioner.reset(new Dune::SeqJac<MatrixType, VectorType, VectorType>(matrix, n, w));
+        return wrapPreconditioner<Dune::SeqJac<MatrixType, VectorType, VectorType>>(matrix, n, w);
     } else if (precond == "GS") {
-        preconditioner.reset(new Dune::SeqGS<MatrixType, VectorType, VectorType>(matrix, n, w));
+        return wrapPreconditioner<Dune::SeqGS<MatrixType, VectorType, VectorType>>(matrix, n, w);
     } else if (precond == "SOR") {
-        preconditioner.reset(new Dune::SeqSOR<MatrixType, VectorType, VectorType>(matrix, n, w));
+        return wrapPreconditioner<Dune::SeqSOR<MatrixType, VectorType, VectorType>>(matrix, n, w);
     } else if (precond == "SSOR") {
-        preconditioner.reset(new Dune::SeqSSOR<MatrixType, VectorType, VectorType>(matrix, n, w));
+        return wrapPreconditioner<Dune::SeqSSOR<MatrixType, VectorType, VectorType>>(matrix, n, w);
     } else if (precond == "ILUn") {
-        preconditioner.reset(new Dune::SeqILUn<MatrixType, VectorType, VectorType>(matrix, n, w));
+        return wrapPreconditioner<Dune::SeqILUn<MatrixType, VectorType, VectorType>>(matrix, n, w);
     } else {
-        std::string msg("No such seq preconditioner : ");
+        std::string msg("No such preconditioner : ");
         msg += precond;
         throw std::runtime_error(msg);
     }
-    return preconditioner;
 }
 
-template <class Smoother, class MatrixType, class VectorType>
+template <class OperatorType, class VectorType, class Comm>
-std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>>
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
-makeAmgPreconditioner(Dune::MatrixAdapter<MatrixType, VectorType, VectorType>& linearoperator,
+makeParPreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm)
-                      const boost::property_tree::ptree& global_prm)
+{
+    auto& matrix = linearoperator.getmat();
+    using MatrixType = typename OperatorType::matrix_type;
+    double w = prm.get<double>("w");
+    int n = prm.get<int>("n");
+    std::string precond(prm.get<std::string>("preconditioner"));
+    if (precond == "ILU0") {
+        return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>>>(
+            comm, matrix, w);
+    } else if (precond == "ParOverILU0") {
+        // Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner.
+        return wrapPreconditioner<
+            DummyUpdatePreconditioner<Opm::ParallelOverlappingILU0<MatrixType, VectorType, VectorType, Comm>>>(
+            matrix, comm, n, w, Opm::MILU_VARIANT::ILU);
+    } else if (precond == "Jac") {
+        return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqJac<MatrixType, VectorType, VectorType>>>(
+            comm, matrix, n, w);
+    } else if (precond == "GS") {
+        return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqGS<MatrixType, VectorType, VectorType>>>(
+            comm, matrix, n, w);
+    } else if (precond == "SOR") {
+        return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqSOR<MatrixType, VectorType, VectorType>>>(
+            comm, matrix, n, w);
+    } else if (precond == "SSOR") {
+        return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqSSOR<MatrixType, VectorType, VectorType>>>(
+            comm, matrix, n, w);
+    } else if (precond == "ILUn") {
+        return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqILUn<MatrixType, VectorType, VectorType>>>(
+            comm, matrix, n, w);
+    } else {
+        std::string msg("No such preconditioner : ");
+        msg += precond;
+        throw std::runtime_error(msg);
+    }
+}
+
+template <class Smoother, class OperatorType, class VectorType>
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
+makeAmgPreconditioner(OperatorType& linearoperator, const boost::property_tree::ptree& global_prm)
 {
     boost::property_tree::ptree prm = global_prm.get_child("amg");
-    typedef Dune::MatrixAdapter<MatrixType, VectorType, VectorType> OperatorType;
+    using MatrixType = typename OperatorType::matrix_type;
-    std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> preconditioner;
+    using CriterionBase
-    typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<MatrixType, Dune::Amg::FirstDiagonal>>
+        = Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<MatrixType, Dune::Amg::FirstDiagonal>>;
-        CriterionBase;
+    using Criterion = Dune::Amg::CoarsenCriterion<CriterionBase>;
-    typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
     int coarsenTarget = prm.get<int>("coarsenTarget");
     int ml = prm.get<int>("maxlevel");
     Criterion criterion(15, coarsenTarget);
@@ -88,7 +127,7 @@ makeAmgPreconditioner(Dune::MatrixAdapter<MatrixType, VectorType, VectorType>& l
         Dune::Amg::Parameters parms;
         parms.setNoPreSmoothSteps(1);
         parms.setNoPostSmoothSteps(1);
-        preconditioner.reset(new Dune::Amg::FastAMG<OperatorType, VectorType>(linearoperator, criterion, parms));
+        return wrapPreconditioner<Dune::Amg::FastAMG<OperatorType, VectorType>>(linearoperator, criterion, parms);
     } else {
         typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
         SmootherArgs smootherArgs;
@@ -97,92 +136,187 @@ makeAmgPreconditioner(Dune::MatrixAdapter<MatrixType, VectorType, VectorType>& l
         // smootherArgs.overlap=SmootherArgs::none;
         // smootherArgs.overlap=SmootherArgs::aggregate;
         smootherArgs.relaxationFactor = prm.get<double>("w");
-        preconditioner.reset(
+        return std::make_shared<Dune::Amg::AMGCPR<OperatorType, VectorType, Smoother>>(
-            new Dune::Amg::AMG<OperatorType, VectorType, Smoother>(linearoperator, criterion, smootherArgs));
+            linearoperator, criterion, smootherArgs);
+    }
+    return std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>(nullptr);
+}
+
+template <class Smoother, class OperatorType, class VectorType, class Comm>
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
+makeParAmgPreconditioner(OperatorType& linearoperator, const boost::property_tree::ptree& global_prm, const Comm& comm)
+{
+    boost::property_tree::ptree prm = global_prm.get_child("amg");
+    using MatrixType = typename OperatorType::matrix_type;
+    using CriterionBase
+        = Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<MatrixType, Dune::Amg::FirstDiagonal>>;
+    using Criterion = Dune::Amg::CoarsenCriterion<CriterionBase>;
+    int coarsenTarget = prm.get<int>("coarsenTarget");
+    int ml = prm.get<int>("maxlevel");
+    Criterion criterion(15, coarsenTarget);
+    criterion.setDefaultValuesIsotropic(2);
+    criterion.setAlpha(prm.get<double>("alpha"));
+    criterion.setBeta(prm.get<double>("beta"));
+    criterion.setMaxLevel(ml);
+    criterion.setSkipIsolated(false);
+    criterion.setDebugLevel(prm.get<int>("verbosity"));
+    if (global_prm.get<std::string>("preconditioner") == "famg") {
+        throw std::runtime_error("The FastAMG preconditioner cannot be used in parallel.");
+    } else {
+        typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
+        SmootherArgs smootherArgs;
+        smootherArgs.iterations = prm.get<int>("n");
+        // smootherArgs.overlap=SmootherArgs::vertex;
+        // smootherArgs.overlap=SmootherArgs::none;
+        // smootherArgs.overlap=SmootherArgs::aggregate;
+        smootherArgs.relaxationFactor = prm.get<double>("w");
+        return std::make_shared<Dune::Amg::AMGCPR<OperatorType, VectorType, Smoother, Comm>>(
+            linearoperator, criterion, smootherArgs, comm);
     }
-    return preconditioner;
 }
 
 
 
 
-template <class MatrixType, class VectorType>
+
-std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>>
+template <class OperatorType, class VectorType>
-makeAmgPreconditioners(Dune::MatrixAdapter<MatrixType, VectorType, VectorType>& linearoperator,
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
-                       const boost::property_tree::ptree& prm)
+makeAmgPreconditioners(OperatorType& linearoperator, const boost::property_tree::ptree& prm)
 {
-    std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> preconditioner;
+    using MatrixType = typename OperatorType::matrix_type;
     if (prm.get<std::string>("preconditioner") == "famg") {
         // smoother type should not be used
-        preconditioner
+        return makeAmgPreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
-            = makeAmgPreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
+            linearoperator, prm);
-                linearoperator, prm);
-        return preconditioner;
     }
 
     std::string precond = prm.get<std::string>("amg.smoother");
     if (precond == "ILU0") {
-        preconditioner
+        return makeAmgPreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
-            = makeAmgPreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
+            linearoperator, prm);
-                linearoperator, prm);
     } else if (precond == "Jac") {
-        preconditioner
+        return makeAmgPreconditioner<Dune::SeqJac<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
-            = makeAmgPreconditioner<Dune::SeqJac<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
+            linearoperator, prm);
-                linearoperator, prm);
+        // } else if (precond == "GS") {
-        // }else if(precond == "GS"){
+        //   return makeAmgPreconditioner<
-        //   preconditioner = makeAmgPreconditioner<
+        // 	Dune::SeqGS<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(linearoperator, prm);
-        // 	Dune::SeqGS<MatrixType, VectorType, VectorType>,
-        // 	MatrixType, VectorType>(linearoperator,prm);
     } else if (precond == "SOR") {
-        preconditioner
+        return makeAmgPreconditioner<Dune::SeqSOR<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
-            = makeAmgPreconditioner<Dune::SeqSOR<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
+            linearoperator, prm);
-                linearoperator, prm);
     } else if (precond == "SSOR") {
-        preconditioner
+        return makeAmgPreconditioner<Dune::SeqSSOR<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
-            = makeAmgPreconditioner<Dune::SeqSSOR<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
+            linearoperator, prm);
-                linearoperator, prm);
     } else if (precond == "ILUn") {
-        preconditioner
+        return makeAmgPreconditioner<Dune::SeqILUn<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
-            = makeAmgPreconditioner<Dune::SeqILUn<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
+            linearoperator, prm);
-                linearoperator, prm);
     } else {
-        std::string msg("No such seq preconditioner : ");
+        std::string msg("No such sequential preconditioner : ");
         msg += precond;
         throw std::runtime_error(msg);
     }
-    return preconditioner;
 }
 
-template <class MatrixType, class VectorType>
+
-std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>>
+template <class OperatorType, class VectorType, class Comm>
-makeTwoLevelPreconditioner(Dune::MatrixAdapter<MatrixType, VectorType, VectorType>& linearoperator,
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
-                           const boost::property_tree::ptree& global_prm)
+makeParAmgPreconditioners(OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm)
+{
+    if (prm.get<std::string>("preconditioner") == "famg") {
+        throw std::runtime_error("The FastAMG preconditioner cannot be used in parallel.");
+    }
+
+    using MatrixType = typename OperatorType::matrix_type;
+    std::string precond = prm.get<std::string>("amg.smoother");
+    if (precond == "ILU0") {
+        using SmootherType = Opm::ParallelOverlappingILU0<MatrixType, VectorType, VectorType, Comm>;
+        return makeParAmgPreconditioner<SmootherType, OperatorType, VectorType, Comm>(linearoperator, prm, comm);
+        /*
+        return makeParAmgPreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
+                linearoperator, prm, comm);
+    } else if (precond == "Jac") {
+        return makeParAmgPreconditioner<Dune::SeqJac<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
+                linearoperator, prm, comm);
+        // } else if (precond == "GS") {
+        //   return makeParAmgPreconditioner<
+        // 	Dune::SeqGS<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(linearoperator, prm, comm);
+    } else if (precond == "SOR") {
+        return makeParAmgPreconditioner<Dune::SeqSOR<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
+                linearoperator, prm, comm);
+    } else if (precond == "SSOR") {
+        return makeParAmgPreconditioner<Dune::SeqSSOR<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
+                linearoperator, prm, comm);
+    } else if (precond == "ILUn") {
+        return makeParAmgPreconditioner<Dune::SeqILUn<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
+                linearoperator, prm, comm);
+        */
+    } else {
+        std::string msg("No such parallel preconditioner : ");
+        msg += precond;
+        throw std::runtime_error(msg);
+    }
+}
+
+
+
+template <class OperatorType, class VectorType>
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
+makeTwoLevelPreconditioner(OperatorType& linearoperator, const boost::property_tree::ptree& global_prm)
 {
     boost::property_tree::ptree prm = global_prm.get_child("cpr");
-    std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> preconditioner;
     if (global_prm.get<std::string>("preconditioner") == "cpr") {
-        preconditioner.reset(new OwningTwoLevelPreconditioner<MatrixType, VectorType, false>(linearoperator, prm));
+        return std::make_shared<OwningTwoLevelPreconditioner<OperatorType, VectorType, false>>(linearoperator, prm);
     } else if (global_prm.get<std::string>("preconditioner") == "cprt") {
-        preconditioner.reset(new OwningTwoLevelPreconditioner<MatrixType, VectorType, true>(linearoperator, prm));
+        return std::make_shared<OwningTwoLevelPreconditioner<OperatorType, VectorType, true>>(linearoperator, prm);
     } else {
         std::string msg("Wrong cpr Should not happen");
         throw std::runtime_error(msg);
     }
-    return preconditioner;
 }
 
-template <class MatrixType, class VectorType>
+template <class OperatorType, class VectorType, class Comm>
-std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>>
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
-makePreconditioner(Dune::MatrixAdapter<MatrixType, VectorType, VectorType>& linearoperator,
+makeParTwoLevelPreconditioner(OperatorType& linearoperator,
-                   const boost::property_tree::ptree& prm)
+                              const boost::property_tree::ptree& global_prm,
+                              const Comm& comm)
+{
+    boost::property_tree::ptree prm = global_prm.get_child("cpr");
+    if (global_prm.get<std::string>("preconditioner") == "cpr") {
+        return std::make_shared<OwningTwoLevelPreconditioner<OperatorType, VectorType, false, Comm>>(
+            linearoperator, prm, comm);
+    } else if (global_prm.get<std::string>("preconditioner") == "cprt") {
+        return std::make_shared<OwningTwoLevelPreconditioner<OperatorType, VectorType, true, Comm>>(
+            linearoperator, prm, comm);
+    } else {
+        std::string msg("Wrong cpr Should not happen");
+        throw std::runtime_error(msg);
+    }
+}
+
+template <class OperatorType, class VectorType>
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
+makePreconditioner(OperatorType& linearoperator, const boost::property_tree::ptree& prm)
 {
     if ((prm.get<std::string>("preconditioner") == "famg") or (prm.get<std::string>("preconditioner") == "amg")) {
-        return makeAmgPreconditioners<MatrixType, VectorType>(linearoperator, prm);
+        return makeAmgPreconditioners<OperatorType, VectorType>(linearoperator, prm);
     } else if ((prm.get<std::string>("preconditioner") == "cpr")
                or (prm.get<std::string>("preconditioner") == "cprt")) {
-        return makeTwoLevelPreconditioner<MatrixType, VectorType>(linearoperator, prm);
+        return makeTwoLevelPreconditioner<OperatorType, VectorType>(linearoperator, prm);
     } else {
-        return makeSeqPreconditioner<MatrixType, VectorType>(linearoperator, prm);
+        return makeSeqPreconditioner<OperatorType, VectorType>(linearoperator, prm);
+    }
+}
+
+template <class OperatorType, class VectorType, class Comm>
+std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
+makePreconditioner(OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm)
+{
+    if ((prm.get<std::string>("preconditioner") == "famg") or (prm.get<std::string>("preconditioner") == "amg")) {
+        return makeParAmgPreconditioners<OperatorType, VectorType, Comm>(linearoperator, prm, comm);
+    } else if ((prm.get<std::string>("preconditioner") == "cpr")
+               or (prm.get<std::string>("preconditioner") == "cprt")) {
+        return makeParTwoLevelPreconditioner<OperatorType, VectorType, Comm>(linearoperator, prm, comm);
+    } else {
+        return makeParPreconditioner<OperatorType, VectorType, Comm>(linearoperator, prm, comm);
     }
 }
 

opm/simulators/linalg/setupPropertyTree.cpp

@@ -24,7 +24,8 @@
 namespace Opm
 {
 
-boost::property_tree::ptree setupPropertyTree(const FlowLinearSolverParameters& p)
+boost::property_tree::ptree
+setupPropertyTree(const FlowLinearSolverParameters& p)
 {
     boost::property_tree::ptree prm;
     if (p.linear_solver_configuration_json_file_ != "none") {

opm/simulators/linalg/twolevelmethodcpr.hh

@@ -16,6 +16,7 @@
 #include<dune/istl/solver.hh>
 
 #include<dune/common/unused.hh>
+#include<dune/common/version.hh>
 
 /**
  * @addtogroup ISTL_PAAMG
@@ -454,13 +455,19 @@ public:
   void updatePreconditioner(FineOperatorType& /* op */,
                             std::shared_ptr<SmootherType> smoother,
                             CoarseLevelSolverPolicy& coarsePolicy)
+  {
+    updatePreconditioner(smoother, coarsePolicy);
+  }
+
+  void updatePreconditioner(std::shared_ptr<SmootherType> smoother,
+                            CoarseLevelSolverPolicy& coarsePolicy)
   {
     //assume new matrix is not reallocated the new precondition should anyway be made
     smoother_ = smoother;
     if (coarseSolver_) {
       policy_->calculateCoarseEntries(*operator_);
       coarsePolicy.setCoarseOperator(*policy_);
-      coarseSolver_->updateAmgPreconditioner(*(policy_->getCoarseLevelOperator()));
+      coarseSolver_->updatePreconditioner();
     } else {
       // we should probably not be heere
       policy_->createCoarseLevelSystem(*operator_);

tests/test_flexiblesolver.cpp

@@ -22,11 +22,12 @@
 #define BOOST_TEST_MODULE OPM_test_FlexibleSolver
 #include <boost/test/unit_test.hpp>
 
+#include <opm/simulators/linalg/FlexibleSolver.hpp>
+
 #include <boost/property_tree/json_parser.hpp>
 #include <boost/property_tree/ptree.hpp>
 #include <fstream>
 #include <iostream>
-#include <opm/simulators/linalg/FlexibleSolver.hpp>
 
 
 template <int bz>