Merge pull request #4699 from atgeirr/nldd-utilities

Add utilities that will be needed for NLDD solvers.

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