Merge pull request #1878 from atgeirr/flexible-solver-revision

Small fixes for the FlexibleSolver system.

CMakeLists_files.cmake

@@ -169,7 +169,6 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/linalg/ISTLSolverEbosCpr.hpp
   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/BlackoilAmg.hpp

@@ -838,10 +838,11 @@ public:
         this->lhs_.resize(this->coarseLevelMatrix_->M());
         this->rhs_.resize(this->coarseLevelMatrix_->N());
         using OperatorArgs = typename Dune::Amg::ConstructionTraits<CoarseOperator>::Arguments;
-        OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_);
 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
+        OperatorArgs oargs(coarseLevelMatrix_, *coarseLevelCommunication_);
         this->operator_ = Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs);
 #else
+        OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_);
         this->operator_.reset(Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs));
 #endif
     }

opm/simulators/linalg/FlexibleSolver.hpp

@@ -34,16 +34,12 @@ namespace Dune
 {
 
 
-template <class X, class Y>
+/// A solver class that encapsulates all needed objects for a linear solver
-class SolverWithUpdate : public Dune::InverseOperator<X, Y>
+/// (operator, scalar product, iterative solver and preconditioner) and sets
-{
+/// them up based on runtime parameters, using the PreconditionerFactory for
-public:
+/// setting up preconditioners.
-    virtual void updatePreconditioner() = 0;
-};
-
-
 template <class MatrixTypeT, class VectorTypeT>
-class FlexibleSolver : public Dune::SolverWithUpdate<VectorTypeT, VectorTypeT>
+class FlexibleSolver : public Dune::InverseOperator<VectorTypeT, VectorTypeT>
 {
 public:
     using MatrixType = MatrixTypeT;
@@ -72,9 +68,13 @@ public:
         linsolver_->apply(x, rhs, reduction, res);
     }
 
-    virtual void updatePreconditioner() override
+    /// Type of the contained preconditioner.
+    using AbstractPrecondType = Dune::PreconditionerWithUpdate<VectorType, VectorType>;
+
+    /// Access the contained preconditioner.
+    AbstractPrecondType& preconditioner()
     {
-        preconditioner_->update();
+        return *preconditioner_;
     }
 
     virtual Dune::SolverCategory::Category category() const override
@@ -84,7 +84,6 @@ public:
 
 private:
     using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
-    using AbstractPrecondType = Dune::PreconditionerWithUpdate<VectorType, VectorType>;
     using AbstractScalarProductType = Dune::ScalarProduct<VectorType>;
     using AbstractSolverType = Dune::InverseOperator<VectorType, VectorType>;
 

opm/simulators/linalg/ISTLSolverEbosFlexible.hpp

@@ -113,18 +113,22 @@ public:
         const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
         bool recreate_solver = false;
         if (this->parameters_.cpr_reuse_setup_ == 0) {
+            // Always recreate solver.
             recreate_solver = true;
         } else if (this->parameters_.cpr_reuse_setup_ == 1) {
+            // Recreate solver on the first iteration of every timestep.
             if (newton_iteration == 0) {
                 recreate_solver = true;
             }
         } else if (this->parameters_.cpr_reuse_setup_ == 2) {
+            // Recreate solver if the last solve used more than 10 iterations.
             if (this->iterations() > 10) {
                 recreate_solver = true;
             }
         } else {
             assert(this->parameters_.cpr_reuse_setup_ == 3);
             assert(recreate_solver == false);
+            // Never recreate solver.
         }
 
         if (recreate_solver || !solver_) {
@@ -135,7 +139,7 @@ public:
             }
             rhs_ = b;
         } else {
-            solver_->updatePreconditioner();
+            solver_->preconditioner().update();
             rhs_ = b;
         }
     }

opm/simulators/linalg/MatrixMarketUtils.hpp

@@ -1,101 +0,0 @@
-/*
-  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
-  Copyright 2019 Dr. Blatt - HPC-Simulation-Software & Services.
-
-  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_MATRIXMARKETUTILS_HEADER_INCLUDED
-#define OPM_MATRIXMARKETUTILS_HEADER_INCLUDED
-
-#include <dune/istl/matrixmarket.hh>
-
-#include <cstddef>
-
-namespace Dune
-{
-
-namespace MatrixMarketImpl
-{
-    template <typename T, typename A, int brows, int bcols, typename D>
-    void makeSparseEntries(Dune::BCRSMatrix<Dune::FieldMatrix<T, brows, bcols>, A>& matrix,
-                           std::vector<int>& i,
-                           std::vector<int>& j,
-                           std::vector<double>& val,
-                           const D&)
-    {
-        // addapted from readSparseEntries
-        typedef Dune::BCRSMatrix<Dune::FieldMatrix<T, brows, bcols>, A> Matrix;
-        std::vector<std::set<IndexData<D>>> rows(matrix.N() * brows);
-        for (std::size_t kk = 0; kk < i.size(); ++kk) {
-            std::size_t row;
-            IndexData<D> data;
-            row = i[kk];
-            assert(row / bcols < matrix.N());
-            data.number = val[kk];
-            data.index = j[kk];
-            assert(data.index / bcols < matrix.M());
-            rows[row].insert(data);
-        }
-
-        // Setup the matrix sparsity pattern
-        int nnz = 0;
-        for (typename Matrix::CreateIterator iter = matrix.createbegin(); iter != matrix.createend(); ++iter) {
-            for (std::size_t brow = iter.index() * brows, browend = iter.index() * brows + brows; brow < browend;
-                 ++brow) {
-                typedef typename std::set<IndexData<D>>::const_iterator Siter;
-                for (Siter siter = rows[brow].begin(), send = rows[brow].end(); siter != send; ++siter, ++nnz)
-                    iter.insert(siter->index / bcols);
-            }
-        }
-
-        // Set the matrix values
-        matrix = 0;
-
-        MatrixValuesSetter<D, brows, bcols> Setter;
-
-        Setter(rows, matrix);
-    }
-} // end namespace MatrixMarketImpl
-
-
-template <typename T, typename A, int brows, int bcols>
-void
-makeMatrixMarket(Dune::BCRSMatrix<Dune::FieldMatrix<T, brows, bcols>, A>& matrix,
-                 std::vector<int> i,
-                 std::vector<int> j,
-                 std::vector<T> val,
-                 size_t rows,
-                 size_t cols)
-{
-    // addapted from readMatrixMarket
-    using namespace MatrixMarketImpl;
-    // std::size_t rows, cols, entries;
-    // std::size_t nnz, blockrows, blockcols;
-    // std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
-    std::size_t blockrows = rows / brows;
-    std::size_t blockcols = cols / bcols;
-    matrix.setSize(blockrows, blockcols);
-    matrix.setBuildMode(Dune::BCRSMatrix<Dune::FieldMatrix<T, brows, bcols>, A>::row_wise);
-    makeSparseEntries(matrix, i, j, val, NumericWrapper<T>());
-}
-
-
-} // end namespace Dune
-
-
-
-#endif // OPM_MATRIXMARKETUTILS_HEADER_INCLUDED

opm/simulators/linalg/PressureSolverPolicy.hpp

@@ -69,7 +69,7 @@ namespace Amg
 
             Dune::SolverCategory::Category category() const override
             {
-                return Dune::SolverCategory::sequential;
+                return linsolver_->category();
             }
 
             void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override
@@ -84,7 +84,7 @@ namespace Amg
 
             void updatePreconditioner()
             {
-                linsolver_->updatePreconditioner();
+                linsolver_->preconditioner().update();
             }
 
         private:

opm/simulators/linalg/PressureTransferPolicy.hpp

@@ -64,8 +64,13 @@ public:
         this->lhs_.resize(this->coarseLevelMatrix_->M());
         this->rhs_.resize(this->coarseLevelMatrix_->N());
         using OperatorArgs = typename Dune::Amg::ConstructionTraits<CoarseOperator>::Arguments;
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
+        OperatorArgs oargs(coarseLevelMatrix_, *coarseLevelCommunication_);
+        this->operator_ = Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs);
+#else
         OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_);
         this->operator_.reset(Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs));
+#endif
     }
 
     virtual void calculateCoarseEntries(const FineOperator& fineOperator) override