Add utilities that will be needed for NLDD solvers.

Includes:
 - SubDomain struct,
 - simple partitioning utilities for testing,
 - some matrix and vector manipulation functions,
 - tests for the above.

CMakeLists_files.cmake

@@ -50,6 +50,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/flow/NonlinearSolverEbos.cpp
   opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.cpp
   opm/simulators/flow/ValidationFunctions.cpp
+  opm/simulators/flow/partitionCells.cpp
   opm/simulators/linalg/bda/WellContributions.cpp
   opm/simulators/linalg/bda/MultisegmentWellContribution.cpp
   opm/simulators/linalg/ExtractParallelGridInformationToISTL.cpp
@@ -226,12 +227,14 @@ endif()
 # find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort
 list (APPEND TEST_SOURCE_FILES
   tests/test_ALQState.cpp
+  tests/test_partitionCells.cpp
   tests/test_blackoil_amg.cpp
   tests/test_convergenceoutputconfiguration.cpp
   tests/test_convergencereport.cpp
   tests/test_deferredlogger.cpp
   tests/test_eclinterregflows.cpp
   tests/test_equil.cc
+  tests/test_extractMatrix.cpp
   tests/test_flexiblesolver.cpp
   tests/test_glift1.cpp
   tests/test_graphcoloring.cpp
@@ -343,6 +346,7 @@ list (APPEND TEST_DATA_FILES
   tests/include/summary.inc
   tests/include/test1_20x30x10.grdecl
   tests/include/well_vfp.ecl
+  tests/test10.partition
   )
 
 
@@ -362,6 +366,8 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp
   opm/simulators/flow/KeywordValidation.hpp
   opm/simulators/flow/ValidationFunctions.hpp
+  opm/simulators/flow/partitionCells.hpp
+  opm/simulators/flow/SubDomain.hpp
   opm/core/props/BlackoilPhases.hpp
   opm/core/props/phaseUsageFromDeck.hpp
   opm/core/props/satfunc/RelpermDiagnostics.hpp
@@ -421,6 +427,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/linalg/SmallDenseMatrixUtils.hpp
   opm/simulators/linalg/WellOperators.hpp
   opm/simulators/linalg/WriteSystemMatrixHelper.hpp
+  opm/simulators/linalg/extractMatrix.hpp
   opm/simulators/linalg/findOverlapRowsAndColumns.hpp
   opm/simulators/linalg/getQuasiImpesWeights.hpp
   opm/simulators/linalg/setupPropertyTree.hpp

opm/simulators/flow/SubDomain.hpp

@@ -0,0 +1,58 @@
+/*
+  Copyright 2021 Total SE
+
+  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_SUBDOMAIN_HEADER_INCLUDED
+#define OPM_SUBDOMAIN_HEADER_INCLUDED
+
+#include <opm/grid/common/SubGridPart.hpp>
+
+#include <utility>
+#include <vector>
+
+namespace Opm
+{
+
+    /// Representing a part of a grid, in a way suitable for performing
+    /// local solves.
+    template <class Grid>
+    struct SubDomain
+    {
+        // The index of a subdomain is arbitrary, but can be used by the
+        // solvers to keep track of well locations etc.
+        int index;
+        // The set of cells that make up a SubDomain, stored as cell indices
+        // in the local numbering of the current MPI rank.
+        std::vector<int> cells;
+        // Flag for each cell of the current MPI rank, true if the cell is part
+        // of the subdomain. If empty, assumed to be all true. Not required for
+        // all nonlinear solver algorithms.
+        std::vector<bool> interior;
+        // Enables subdomain solves and linearization using the generic linearization
+        // approach (i.e. FvBaseLinearizer as opposed to TpfaLinearizer).
+        Dune::SubGridPart<Grid> view;
+        // Constructor that moves from its argument.
+        SubDomain(const int i, std::vector<int>&& c, std::vector<bool>&& in, Dune::SubGridPart<Grid>&& v)
+            : index(i), cells(std::move(c)), interior(std::move(in)), view(std::move(v))
+        {}
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_SUBDOMAIN_HEADER_INCLUDED

opm/simulators/flow/partitionCells.cpp

@@ -0,0 +1,74 @@
+/*
+  Copyright 2021 Total SE
+
+  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/>.
+*/
+
+#include <opm/simulators/flow/partitionCells.hpp>
+
+#include <fmt/format.h>
+
+#include <algorithm>
+#include <cassert>
+#include <fstream>
+#include <iterator>
+#include <numeric>
+#include <stdexcept>
+#include <string>
+
+namespace Opm
+{
+
+
+    // Read from file, containing one number per cell.
+    std::pair<std::vector<int>, int> partitionCellsFromFile(const std::string& filename, const int num_cells)
+    {
+        // Read file into single vector.
+        std::ifstream is(filename);
+        const std::vector<int> cellpart{std::istream_iterator<int>(is), std::istream_iterator<int>()};
+        if (cellpart.size() != size_t(num_cells)) {
+            auto msg = fmt::format("Partition file contains {} entries, but there are {} cells.",
+                                   cellpart.size(), num_cells);
+            throw std::runtime_error(msg);
+        }
+
+        // Create and return the output domain vector.
+        const int num_domains = (*std::max_element(cellpart.begin(), cellpart.end())) + 1;
+        return { cellpart, num_domains };
+    }
+
+
+    // Trivially simple partitioner
+    std::pair<std::vector<int>, int> partitionCellsSimple(const int num_cells, const int num_domains)
+    {
+        // Build the partitions.
+        const int dom_sz = num_cells / num_domains;
+        std::vector<int> bounds(num_domains + 1, dom_sz);
+        bounds[0] = 0;
+        for (int i = 0; i < num_cells % num_domains; ++i) {
+            ++bounds[i + 1];
+        }
+        std::partial_sum(bounds.begin(), bounds.end(), bounds.begin());
+        assert(bounds[num_domains] == num_cells);
+        std::vector<int> part(num_cells);
+        for (int i = 0; i < num_domains; ++i) {
+            std::fill(part.begin() + bounds[i], part.begin() + bounds[i + 1], i);
+        }
+        return { part, num_domains };
+    }
+
+
+} // namespace Opm

opm/simulators/flow/partitionCells.hpp

@@ -0,0 +1,41 @@
+/*
+  Copyright 2021 Total SE
+
+  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_ASPINPARTITION_HEADER_INCLUDED
+#define OPM_ASPINPARTITION_HEADER_INCLUDED
+
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace Opm
+{
+
+    /// Read a partitioning from file, assumed to contain one number per cell, its partition number.
+    /// \return pair containing a partition vector (partition number for each cell), and the number of partitions.
+    std::pair<std::vector<int>, int> partitionCellsFromFile(const std::string& filename, const int num_cells);
+
+    /// Trivially simple partitioner assigning partitions en bloc, consecutively by cell index.
+    /// \return pair containing a partition vector (partition number for each cell), and the number of partitions.
+    std::pair<std::vector<int>, int> partitionCellsSimple(const int num_cells, const int num_domains);
+
+} // namespace Opm
+
+
+#endif // OPM_ASPINPARTITION_HEADER_INCLUDED

opm/simulators/linalg/extractMatrix.hpp

@@ -0,0 +1,154 @@
+/*
+  Copyright 2021 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_EXTRACT_MATRIX_HEADER_INCLUDED
+#define OPM_EXTRACT_MATRIX_HEADER_INCLUDED
+
+#include <algorithm>
+#include <cassert>
+#include <cstddef>
+#include <vector>
+
+namespace Opm
+{
+
+namespace Details
+{
+
+    template <class Matrix>
+    void copySubMatrix(const Matrix& A, const std::vector<int>& indices, Matrix& B)
+    {
+        // Copy elements so that B_{i, j} = A_{indices[i], indices[j]}.
+        for (auto row = B.begin(); row != B.end(); ++row) {
+            for (auto col = row->begin(); col != row->end(); ++col) {
+                *col = A[indices[row.index()]][indices[col.index()]];
+            }
+        }
+    }
+
+
+    template <class Matrix>
+    Matrix extractMatrix(const Matrix& m, const std::vector<int>& indices)
+    {
+        assert(std::is_sorted(indices.begin(), indices.end()));
+
+        // Set up reverse index map.
+        const std::size_t n = indices.size();
+        std::size_t newrow = 0;
+        enum { NotIncluded = -1 };
+        std::vector<int> index_map(m.N(), NotIncluded);
+        for (auto row = m.begin(); row != m.end(); ++row) {
+            const int row_index = row.index();
+            if (row_index == indices[newrow]) {
+                index_map[row_index] = newrow;
+                ++newrow;
+                if (newrow == n) {
+                    break;
+                }
+            }
+        }
+        assert(newrow == n);
+
+        // Count nonzeroes.
+        std::size_t nnz = 0;
+        for (auto row = m.begin(); row != m.end(); ++row) {
+            if (index_map[row.index()] != NotIncluded) {
+                for (auto col = row->begin(); col != row->end(); ++col) {
+                    if (index_map[col.index()] != NotIncluded) {
+                        ++nnz;
+                    }
+                }
+            }
+        }
+
+        // Create the matrix structure.
+        Matrix res(n, n, nnz, Matrix::row_wise);
+        auto from_row = m.begin();
+        for (auto row = res.createbegin(); row != res.createend(); ++row) {
+            // Move from_row to point to the row to be extracted.
+            while (static_cast<int>(from_row.index()) < indices[row.index()]) {
+                ++from_row;
+            }
+            assert(static_cast<int>(from_row.index()) == indices[row.index()]);
+            // Insert nonzeros for row.
+            for (auto from_col = from_row->begin(); from_col != from_row->end(); ++from_col) {
+                const int new_col = index_map[from_col.index()];
+                if (new_col != NotIncluded) {
+                    row.insert(new_col);
+                }
+            }
+        }
+
+        copySubMatrix(m, indices, res);
+        return res;
+    }
+
+
+    template <class Vector>
+    Vector extractVector(const Vector& x, const std::vector<int>& indices)
+    {
+        Vector res(indices.size());
+        for (std::size_t ii = 0; ii < indices.size(); ++ii) {
+            res[ii] = x[indices[ii]];
+        }
+        return res;
+    }
+
+
+    template <class Vector>
+    void setGlobal(const Vector& x, const std::vector<int>& indices, Vector& global_x)
+    {
+        for (std::size_t ii = 0; ii < indices.size(); ++ii) {
+            global_x[indices[ii]] = x[ii];
+        }
+    }
+
+
+
+    template <class Matrix>
+    bool matrixEqual(const Matrix& m1, const Matrix& m2)
+    {
+        // Compare size and nonzeroes.
+        if (m1.N() != m2.N()) return false;
+        if (m1.M() != m2.M()) return false;
+        if (m1.nonzeroes() != m2.nonzeroes()) return false;
+
+        auto row1 = m1.begin();
+        auto row2 = m2.begin();
+        for (; row1 != m1.end(); ++row1, ++row2) {
+            if (row2 == m2.end()) return false;
+            if (row1.index() != row2.index()) return false;
+            auto col1 = row1->begin();
+            auto col2 = row2->begin();
+            for (; col1 != row1->end(); ++col1, ++col2) {
+                if (col2 == row2->end()) return false;
+                if (col1.index() != col2.index()) return false;
+                if (*col1 != *col2) return false;
+            }
+        }
+        return true;
+    }
+
+
+} // namespace Details
+
+
+} // namespace Opm
+
+#endif // OPM_EXTRACT_MATRIX_HEADER_INCLUDED

tests/test10.partition

@@ -0,0 +1,10 @@
+0
+0
+1
+1
+2
+2
+1
+1
+0
+0

tests/test_extractMatrix.cpp

@@ -0,0 +1,112 @@
+/*
+  Copyright 2023 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/>.
+*/
+
+#include <config.h>
+
+#define BOOST_TEST_MODULE OPM_test_extractMatrix
+#include <boost/test/unit_test.hpp>
+
+#include <opm/simulators/linalg/extractMatrix.hpp>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+using B = Dune::FieldMatrix<double, 2, 2>;
+using M = Dune::BCRSMatrix<B>;
+using V = Dune::BlockVector<Dune::FieldVector<double, 2>>;
+
+M build3x3BlockMatrix()
+{
+    // Build matrix with this pattern:
+    //   ( x     ) 
+    //   (   x   )
+    //   (   x x )
+    M res(3, 3, 4, M::row_wise);
+    for (auto row = res.createbegin(); row != res.createend(); ++row) {
+        row.insert(row.index());
+        if (row.index() == 2) {
+            row.insert(1);
+        }
+    }
+    // Set all entries to the block (1, 2; 3, x), where x starts at 4
+    // but gets incremented for each entry.
+    B b;
+    b[0][0] = 1.0;
+    b[0][1] = 2.0;
+    b[1][0] = 3.0;
+    b[1][1] = 4.0;
+    res[0][0] = b;
+    b[1][1] = b[1][1] + 1.0;
+    res[1][1] = b;
+    b[1][1] = b[1][1] + 1.0;
+    res[2][2] = b;
+    b[1][1] = b[1][1] + 1.0;
+    res[2][1] = b;
+
+    return res;
+}
+
+BOOST_AUTO_TEST_CASE(extractMatrix)
+{
+    auto m1 = build3x3BlockMatrix();
+    BOOST_CHECK_EQUAL(m1[2][1][1][1], 7.0);
+    std::vector<int> indices = {1, 2};
+    auto m2 = Opm::Details::extractMatrix(m1, indices);
+    BOOST_CHECK_EQUAL(m2[0][0], m1[1][1]);
+    BOOST_CHECK_EQUAL(m2[1][1], m1[2][2]);
+    BOOST_CHECK_EQUAL(m2[1][0], m1[2][1]);
+}
+
+BOOST_AUTO_TEST_CASE(copySubMatrixAndMatrixEqual)
+{
+    auto m1 = build3x3BlockMatrix();
+    std::vector<int> indices = {1, 2};
+    auto m2 = Opm::Details::extractMatrix(m1, indices);
+    m1[2][1][0][0] = 0.1234;
+    auto m3 = Opm::Details::extractMatrix(m1, indices);
+    BOOST_CHECK(m2[1][0] != m1[2][1]);
+    BOOST_CHECK(m3[1][0] == m1[2][1]);
+    BOOST_CHECK(!Opm::Details::matrixEqual(m2, m3));
+    Opm::Details::copySubMatrix(m1, indices, m2);
+    BOOST_CHECK(Opm::Details::matrixEqual(m2, m3));
+}
+
+BOOST_AUTO_TEST_CASE(vectorOps)
+{
+    V v1(4);
+    v1[0] = { 0.1, 0.2 };
+    v1[1] = { 0.3, 0.4 };
+    v1[2] = { 0.5, 0.6 };
+    v1[3] = { 0.6, 0.8 };
+    std::vector<int> indices = { 0, 2, 3 };
+
+    V vref(3);
+    vref[0] = { 0.1, 0.2 };
+    vref[1] = { 0.5, 0.6 };
+    vref[2] = { 0.6, 0.8 };
+
+    auto v2 = Opm::Details::extractVector(v1, indices);
+    BOOST_CHECK_EQUAL_COLLECTIONS(vref.begin(), vref.end(), v2.begin(), v2.end());
+
+    V v3 = v1;
+    v3[2] = { 1.1, 1.2 };
+    v2[1] = { 1.1, 1.2 };
+    Opm::Details::setGlobal(v2, indices, v1);
+    BOOST_CHECK_EQUAL_COLLECTIONS(v1.begin(), v1.end(), v3.begin(), v3.end());
+}

tests/test_partitionCells.cpp

@@ -0,0 +1,59 @@
+/*
+  Copyright 2021 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/>.
+*/
+
+#include <config.h>
+
+#define BOOST_TEST_MODULE OPM_test_aspinPartition
+#include <boost/test/unit_test.hpp>
+
+#include <opm/simulators/flow/partitionCells.hpp>
+
+
+BOOST_AUTO_TEST_CASE(FileBased)
+{
+    auto [part, num_part] = Opm::partitionCellsFromFile("test10.partition", 10);
+    BOOST_CHECK_EQUAL(num_part, 3);
+    std::vector<int> expected = { 0, 0, 1, 1, 2, 2, 1, 1, 0, 0 };
+    BOOST_CHECK_EQUAL_COLLECTIONS(expected.begin(), expected.end(), part.begin(), part.end());
+}
+
+BOOST_AUTO_TEST_CASE(FileBasedWrongNumberOfCells)
+{
+    auto func = []() { auto [part, num_part] = Opm::partitionCellsFromFile("test10.partition", 11); };
+    BOOST_CHECK_THROW(func(), std::runtime_error);
+}
+
+BOOST_AUTO_TEST_CASE(Simple1)
+{
+    auto [part, num_part]  = Opm::partitionCellsSimple(10, 3);
+    BOOST_CHECK_EQUAL(num_part, 3);
+    std::vector<int> expected = { 0, 0, 0, 0, 1, 1, 1, 2, 2, 2 };
+    BOOST_CHECK_EQUAL_COLLECTIONS(expected.begin(), expected.end(), part.begin(), part.end());
+
+}
+
+BOOST_AUTO_TEST_CASE(Simple2)
+{
+    auto [part, num_part]  = Opm::partitionCellsSimple(10, 7);
+    BOOST_CHECK_EQUAL(num_part, 7);
+    std::vector<int> expected = { 0, 0, 1, 1, 2, 2, 3, 4, 5, 6 };
+    BOOST_CHECK_EQUAL_COLLECTIONS(expected.begin(), expected.end(), part.begin(), part.end());
+
+}
+