Setup parallel solvers for tracers.

CMakeLists_files.cmake

@@ -310,6 +310,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/utils/ParallelEclipseState.hpp
   opm/simulators/utils/ParallelRestart.hpp
   opm/simulators/utils/PropsCentroidsDataHandle.hpp
+  opm/simulators/utils/VectorVectorDataHandle.hpp
   opm/simulators/wells/PerfData.hpp
   opm/simulators/wells/PerforationData.hpp
   opm/simulators/wells/RateConverter.hpp

ebos/eclgenerictracermodel.cc

@@ -24,6 +24,8 @@
 #include <config.h>
 #include <ebos/eclgenerictracermodel.hh>
 
+#include <opm/simulators/linalg/PropertyTree.hpp>
+#include <opm/simulators/linalg/FlexibleSolver.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
 #include <opm/grid/CpGrid.hpp>
 #include <opm/grid/polyhedralgrid.hh>
@@ -31,10 +33,13 @@
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/EclipseState/Runspec.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/TracerVdTable.hpp>
+#include <opm/simulators/linalg/FlexibleSolver.hpp>
 
 #include <dune/istl/operators.hh>
 #include <dune/istl/solvers.hh>
+#include <dune/istl/schwarz.hh>
 #include <dune/istl/preconditioners.hh>
+#include <dune/istl/schwarz.hh>
 
 #if HAVE_DUNE_FEM
 #include <dune/fem/gridpart/adaptiveleafgridpart.hh>
@@ -47,9 +52,44 @@
 #include <stdexcept>
 #include <functional>
 #include <array>
+#include <string>
 
 namespace Opm {
 
+#if HAVE_MPI
+template<class M, class V>
+struct TracerSolverSelector
+{
+    using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
+    using TracerOperator = Dune::OverlappingSchwarzOperator<M, V, V, Comm>;
+    using type = Dune::FlexibleSolver<M, V>;
+};
+template<class Vector, class Grid, class Matrix>
+std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
+                                                            Dune::OwnerOverlapCopyCommunication<int,int>>>,
+           std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
+createParallelFlexibleSolver(const Grid&, const Matrix&, const PropertyTree&)
+{
+    OPM_THROW(std::logic_error, "Grid not supported for parallel Tracers.");
+    return {nullptr, nullptr};
+}
+
+template<class Vector, class Matrix>
+std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
+                                                            Dune::OwnerOverlapCopyCommunication<int,int>>>,
+           std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
+createParallelFlexibleSolver(const Dune::CpGrid& grid, const Matrix& M, const PropertyTree& prm)
+{
+        using TracerOperator = Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
+                                                                Dune::OwnerOverlapCopyCommunication<int,int>>;
+        using TracerSolver = Dune::FlexibleSolver<Matrix, Vector>;
+        const auto& cellComm = grid.cellCommunication();
+        auto op = std::make_unique<TracerOperator>(M, cellComm);
+        auto dummyWeights = [](){ return Vector();};
+        return {std::move(op), std::make_unique<TracerSolver>(*op, cellComm, prm, dummyWeights, 0)};
+}
+#endif
+
 template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar>
 EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
 EclGenericTracerModel(const GridView& gridView,
@@ -217,24 +257,50 @@ linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
     int maxIter = 100;
 
     int verbosity = 0;
-    using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
-    using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
-    using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
-    using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
+    PropertyTree prm;
+    prm.put("maxiter", maxIter);
+    prm.put("tol", tolerance);
+    prm.put("verbosity", verbosity);
+    prm.put("solver", std::string("bicgstab"));
+    prm.put("preconditioner.type", std::string("ParOverILU0"));
 
-    TracerOperator tracerOperator(M);
-    TracerScalarProduct tracerScalarProduct;
-    TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
+#if HAVE_MPI
+    if(gridView_.grid().comm().size() > 1)
+    {
+        auto [tracerOperator, solver] =
+            createParallelFlexibleSolver<TracerVector>(gridView_.grid(), M, prm);
+        (void) tracerOperator;
 
-    TracerSolver solver (tracerOperator, tracerScalarProduct,
-                         tracerPreconditioner, tolerance, maxIter,
-                         verbosity);
+        Dune::InverseOperatorResult result;
+        solver->apply(x, b, result);
 
-    Dune::InverseOperatorResult result;
-    solver.apply(x, b, result);
+        // return the result of the solver
+        return result.converged;
+    }
+    else
+    {
+#endif
+        using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
+        using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
+        using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
+        using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
 
-    // return the result of the solver
-    return result.converged;
+        TracerOperator tracerOperator(M);
+        TracerScalarProduct tracerScalarProduct;
+        TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
+
+        TracerSolver solver (tracerOperator, tracerScalarProduct,
+                             tracerPreconditioner, tolerance, maxIter,
+                             verbosity);
+
+        Dune::InverseOperatorResult result;
+        solver.apply(x, b, result);
+
+        // return the result of the solver
+        return result.converged;
+#if HAVE_MPI
+    }
+#endif
 }
 
 template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
@@ -249,29 +315,57 @@ linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::
     int maxIter = 100;
 
     int verbosity = 0;
-    using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
-    using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
-    using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
-    using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
+    PropertyTree prm;
+    prm.put("maxiter", maxIter);
+    prm.put("tol", tolerance);
+    prm.put("verbosity", verbosity);
+    prm.put("solver", std::string("bicgstab"));
+    prm.put("preconditioner.type", std::string("ParOverILU0"));
 
-    TracerOperator tracerOperator(M);
-    TracerScalarProduct tracerScalarProduct;
-    TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
-
-    TracerSolver solver (tracerOperator, tracerScalarProduct,
-                         tracerPreconditioner, tolerance, maxIter,
-                         verbosity);
-
-    bool converged = true;
-    for (size_t nrhs =0; nrhs < b.size(); ++nrhs) {
-        x[nrhs] = 0.0;
-        Dune::InverseOperatorResult result;
-        solver.apply(x[nrhs], b[nrhs], result);
-        converged = (converged && result.converged);
+#if HAVE_MPI
+    if(gridView_.grid().comm().size() > 1)
+    {
+        auto [tracerOperator, solver] =
+            createParallelFlexibleSolver<TracerVector>(gridView_.grid(), M, prm);
+        (void) tracerOperator;
+        bool converged = true;
+        for (size_t nrhs =0; nrhs < b.size(); ++nrhs) {
+            x[nrhs] = 0.0;
+            Dune::InverseOperatorResult result;
+            solver->apply(x[nrhs], b[nrhs], result);
+            converged = (converged && result.converged);
+        }
+        return converged;
     }
+    else
+    {
+#endif
+        using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
+        using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
+        using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
+        using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
 
-    // return the result of the solver
-    return converged;
+        TracerOperator tracerOperator(M);
+        TracerScalarProduct tracerScalarProduct;
+        TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
+
+        TracerSolver solver (tracerOperator, tracerScalarProduct,
+                             tracerPreconditioner, tolerance, maxIter,
+                             verbosity);
+
+        bool converged = true;
+        for (size_t nrhs =0; nrhs < b.size(); ++nrhs) {
+            x[nrhs] = 0.0;
+            Dune::InverseOperatorResult result;
+            solver.apply(x[nrhs], b[nrhs], result);
+            converged = (converged && result.converged);
+        }
+
+        // return the result of the solver
+        return converged;
+#if HAVE_MPI
+    }
+#endif
 }
 
 #if HAVE_DUNE_FEM

ebos/ecltracermodel.hh

@@ -31,6 +31,7 @@
 #include <ebos/eclgenerictracermodel.hh>
 
 #include <opm/models/utils/propertysystem.hh>
+#include <opm/simulators/utils/VectorVectorDataHandle.hpp>
 
 #include <string>
 #include <vector>
@@ -235,6 +236,15 @@ protected:
         for (; elemIt != elemEndIt; ++ elemIt) {
             elemCtx.updateAll(*elemIt);
 
+            size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timIdx=*/0);
+
+            if (elemIt->partitionType() != Dune::InteriorEntity)
+            {
+                // Dirichlet boundary conditions needed for the parallel matrix
+                (*this->tracerMatrix_)[I][I][0][0] = 1.;
+                continue;
+            }
+
             Scalar extrusionFactor =
                     elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
             Valgrind::CheckDefined(extrusionFactor);
@@ -245,7 +255,6 @@ protected:
                     * extrusionFactor;
             Scalar dt = elemCtx.simulator().timeStepSize();
 
-            size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timIdx=*/0);
             size_t I1 = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timIdx=*/1);
 
             std::vector<Scalar> storageOfTimeIndex1(tr.numTracer());
@@ -323,6 +332,12 @@ protected:
                 }
             }
         }
+
+        // Communicate overlap using grid Communication
+        auto handle = VectorVectorDataHandle<GridView, std::vector<TracerVector>>(tr.residual_,
+                                                                                  simulator_.gridView());
+        simulator_.gridView().communicate(handle, Dune::InteriorBorder_All_Interface,
+                                          Dune::ForwardCommunication);
     }
 
     template <class TrRe>

opm/simulators/utils/VectorVectorDataHandle.hpp

@@ -0,0 +1,102 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=4:
+/*
+  Copyright 2021 Equinor AS.
+
+  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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/**
+ * \file
+ * \brief A datahandle sending data located in multiple vectors
+ * \author Markus Blatt, OPM-OP AS
+ */
+
+namespace Opm
+{
+
+/// \brief A data handle sending multiple data store in vectors attached
+///        to cells.
+///
+/// Each data is assumed to a container with operator[] and the
+/// class operates on a vector of these.
+/// \tparam GridView the type of the grid view the data associated with
+/// \tparam The type of the vector of vectors.
+template<class GridView, class Vector>
+class VectorVectorDataHandle
+  : public Dune::CommDataHandleIF<VectorVectorDataHandle<GridView,Vector>,
+                                  std::decay_t<decltype(Vector()[0][0])>>
+{
+public:
+
+  /// \brief the data type we send
+  using DataType = std::decay_t<decltype(Vector()[0][0])>;
+
+  /// \brief Constructor
+  /// \param data The vector of data vectors
+  /// \param gridView The gridview the data is attached to.
+  VectorVectorDataHandle(Vector& data, const GridView& gridView)
+    : data_(data), gridView_(gridView)
+  {}
+
+  bool contains(int /* dim */, int codim) const
+  {
+    return codim == 0;
+  }
+
+  bool fixedsize(int /* dim */, int /* codim */) const
+  {
+    return true;
+  }
+  bool fixedSize(int /* dim */, int /* codim */) const
+  {
+    return true;
+  }
+
+  template<class EntityType>
+  std::size_t size(const EntityType /* entity */) const
+  {
+    return data_.size();
+  }
+
+
+  template<class BufferType, class EntityType>
+  void gather(BufferType& buffer, const EntityType& e) const
+  {
+    for(const auto& vec: data_)
+    {
+      buffer.write(vec[gridView_.indexSet().index(e)]);
+    }
+  }
+
+  template<class BufferType, class EntityType>
+  void scatter(BufferType& buffer, const EntityType& e, std::size_t n)
+  {
+    assert(n == data_.size());
+    for(auto& vec: data_)
+    {
+      buffer.read(vec[gridView_.indexSet().index(e)]);
+    }
+  }
+private:
+  Vector& data_;
+  const GridView& gridView_;
+};
+
+} // end namespace Opm