UpwindTransportSolver.hpp

Go to the documentation of this file.

/*
  Copyright 2026 Equinor ASA.
 
  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_UPWIND_TRANSPORT_SOLVER_HEADER_INCLUDED
#define OPM_UPWIND_TRANSPORT_SOLVER_HEADER_INCLUDED
 
#include <opm/common/utility/SparseTable.hpp>
 
#include <opm/simulators/utils/Tarjan.hpp>
 
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <span>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
 
namespace Opm
{
 
struct UpwindTransportSettings
{
    int maxDirectSize = 200;
 
    double relTolerance = 1e-4;
 
    double absTolerance = 1e-20;
 
    int maxIterations = 100;
};
 
struct UpwindTransportStatistics
{
    int numComponents = 0;
 
    int numLoops = 0;
 
    int largestComponent = 0;
 
    int maxIterations = 0;
};
 
namespace detail
{
 
    template <class Scalar>
    void solveDenseSystem(std::vector<Scalar>& a,
                          std::vector<Scalar>& b,
                          std::size_t m,
                          std::size_t numRhs)
    {
        for (std::size_t col = 0; col < m; ++col) {
            // Gaussian elimination with partial pivoting.
            std::size_t pivot = col;
            for (std::size_t row = col + 1; row < m; ++row) {
                if (std::abs(a[row * m + col]) > std::abs(a[pivot * m + col])) {
                    pivot = row;
                }
            }
            if (a[pivot * m + col] == Scalar {0}) {
                throw std::runtime_error("Opm::solveReordered: singular matrix");
            }
            if (pivot != col) {
                for (std::size_t k = 0; k < m; ++k) {
                    std::swap(a[pivot * m + k], a[col * m + k]);
                }
                for (std::size_t k = 0; k < numRhs; ++k) {
                    std::swap(b[pivot * numRhs + k], b[col * numRhs + k]);
                }
            }
            for (std::size_t row = col + 1; row < m; ++row) {
                const Scalar factor = a[row * m + col] / a[col * m + col];
                if (factor == Scalar {0}) {
                    continue;
                }
                for (std::size_t k = col; k < m; ++k) {
                    a[row * m + k] -= factor * a[col * m + k];
                }
                for (std::size_t k = 0; k < numRhs; ++k) {
                    b[row * numRhs + k] -= factor * b[col * numRhs + k];
                }
            }
        }
        for (std::size_t col = m; col-- > 0;) {
            for (std::size_t k = 0; k < numRhs; ++k) {
                Scalar sum = b[col * numRhs + k];
                for (std::size_t j = col + 1; j < m; ++j) {
                    sum -= a[col * m + j] * b[j * numRhs + k];
                }
                b[col * numRhs + k] = sum / a[col * m + col];
            }
        }
    }
 
    template <class Block>
    bool isNonzeroBlock(const Block& block)
    {
        for (std::size_t a = 0; a < Block::rows; ++a) {
            for (std::size_t b = 0; b < Block::cols; ++b) {
                if (block[a][b] != 0) {
                    return true;
                }
            }
        }
        return false;
    }
 
    template <class Block>
    bool isFiniteBlock(const Block& block)
    {
        for (std::size_t a = 0; a < Block::rows; ++a) {
            for (std::size_t b = 0; b < Block::cols; ++b) {
                if (!std::isfinite(block[a][b])) {
                    return false;
                }
            }
        }
        return true;
    }
 
    template <class Matrix, class Function>
    void forEachNonzeroCoupling(const Matrix& matrix, Function&& f)
    {
        for (auto row = matrix.begin(); row != matrix.end(); ++row) {
            for (auto entry = row->begin(); entry != row->end(); ++entry) {
                if (entry.index() != row.index() && isNonzeroBlock(*entry)) {
                    f(row.index(), entry.index(), *entry);
                }
            }
        }
    }
 
    template <class Matrix>
    class MatrixCoefficients
    {
        using Block = typename Matrix::block_type;
        static_assert(Block::rows == Block::cols, "The matrix blocks must be square");
 
    public:
        using Scalar = typename Block::field_type;
        static constexpr std::size_t blockSize = Block::rows;
 
        explicit MatrixCoefficients(const Matrix& matrix)
        {
            if (matrix.N() != matrix.M()) {
                throw std::invalid_argument("Opm::solveReordered: the matrix must be square");
            }
            diagonal_.reserve(matrix.N());
            for (std::size_t r = 0; r < matrix.N(); ++r) {
                const auto entry = matrix[r].find(r);
                if (entry == matrix[r].end()) {
                    throw std::invalid_argument("Opm::solveReordered: row " + std::to_string(r)
                                                + " has no diagonal block");
                }
                // Checked up front, since an infinite entry can give a finite solution, e.g.
                // b / inf = 0 or after elimination in a directly solved loop
                if (!isFiniteBlock(*entry)) {
                    throw std::runtime_error("Opm::solveReordered: NaN or infinite entry in the "
                                             "diagonal block of row " + std::to_string(r));
                }
                diagonal_.push_back(&*entry);
            }
 
            std::vector<int> rowStart(matrix.N() + 1, 0);
            std::vector<int> columns;
            std::vector<const Block*> blocks;
            forEachNonzeroCoupling(
                matrix, [&](std::size_t r, std::size_t c, const Block& block) {
                    if (!isFiniteBlock(block)) {
                        throw std::runtime_error("Opm::solveReordered: NaN or infinite entry in "
                                                 "the block of row " + std::to_string(r)
                                                 + ", column " + std::to_string(c));
                    }
                    ++rowStart[r + 1];
                    columns.push_back(static_cast<int>(c));
                    blocks.push_back(&block);
                });
            for (std::size_t r = 0; r < matrix.N(); ++r) {
                rowStart[r + 1] += rowStart[r];
            }
            graph_ = SparseTable<int>(std::move(columns), std::vector<int>(rowStart));
            couplings_ = SparseTable<const Block*>(std::move(blocks), std::move(rowStart));
        }
 
        std::size_t numCells() const
        {
            return diagonal_.size();
        }
 
        const SparseTable<int>& graph() const
        {
            return graph_;
        }
 
        const Block& diagonal(std::size_t c) const
        {
            return *diagonal_[c];
        }
 
        template <class Function>
        void forEachCoupling(std::size_t c, Function&& f) const
        {
            auto block = couplings_[static_cast<int>(c)].begin();
            for (const int u : graph_[static_cast<int>(c)]) {
                f(static_cast<std::size_t>(u), **block++);
            }
        }
 
    private:
        std::vector<const Block*> diagonal_; 
        SparseTable<int> graph_; 
        SparseTable<const Block*> couplings_; 
    };
 
 
    template <class Coefficients, class Vectors>
    class ReorderedSystem
    {
        using Scalar = typename Coefficients::Scalar;
        static constexpr std::size_t N = Coefficients::blockSize;
 
    public:
        ReorderedSystem(const Coefficients& coefficients,
                        const Vectors& rhs,
                        Vectors& x,
                        const UpwindTransportSettings& settings,
                        const std::string& name)
            : coefficients_(coefficients)
            , rhs_(rhs)
            , x_(x)
            , settings_(settings)
            , name_(name)
            , numCells_(coefficients.numCells())
            , numVectors_(rhs.size())
            , localIndex_(numCells_, -1)
            , cellSolution_(N * numVectors_)
            , change_(N * numVectors_)
            , scale_(N * numVectors_)
        {
            validate_();
        }
 
        UpwindTransportStatistics solve()
        {
            const auto components = tarjan(coefficients_.graph());
            UpwindTransportStatistics statistics;
            statistics.numComponents = components.numComponents();
            for (const auto cells : components.components()) {
                const int size = static_cast<int>(cells.size());
                statistics.largestComponent = std::max(statistics.largestComponent, size);
                if (size == 1) {
                    solveCell_(cells);
                } else {
                    ++statistics.numLoops;
                    if (size * static_cast<int>(N) <= settings_.maxDirectSize) {
                        solveDirectly_(cells);
                    } else {
                        const int sweeps = solveByGaussSeidel_(cells);
                        statistics.maxIterations = std::max(statistics.maxIterations, sweeps);
                    }
                }
            }
            return statistics;
        }
 
    private:
        void validate_() const
        {
            if (x_.size() != numVectors_) {
                throw std::invalid_argument(name_ + ": sizes do not match");
            }
            for (std::size_t q = 0; q < numVectors_; ++q) {
                if (static_cast<std::size_t>(rhs_[q].size()) != numCells_
                    || static_cast<std::size_t>(x_[q].size()) != numCells_) {
                    throw std::invalid_argument(name_ + ": wrong vector length");
                }
            }
        }
 
        void checkFinite_(Scalar value,
                          std::size_t q,
                          std::span<const int> cells,
                          bool useInLoop) const
        {
            if (!std::isfinite(value)) [[unlikely]] {
                throw std::runtime_error(nonFiniteMessage_(q, cells, useInLoop));
            }
        }
 
        std::string nonFiniteMessage_(std::size_t q,
                                      std::span<const int> cells,
                                      bool useInLoop) const
        {
            // A sorted copy, so that looking up the cells of a large loop stays cheap
            std::vector<int> sorted(cells.begin(), cells.end());
            std::ranges::sort(sorted);
            const auto inCells = [&sorted](std::size_t u) {
                return std::ranges::binary_search(sorted, static_cast<int>(u));
            };
            std::string cause;
            for (const int cell : cells) {
                const auto c = static_cast<std::size_t>(cell);
                for (std::size_t a = 0; a < N && cause.empty(); ++a) {
                    if (!std::isfinite(rhs_[q][c][a])) {
                        cause = "the right-hand side in row " + std::to_string(c);
                    }
                }
                coefficients_.forEachCoupling(c, [&](std::size_t u, const auto&) {
                    if (!cause.empty()) {
                        return;
                    }
                    if (useInLoop || !inCells(u)) {
                        for (std::size_t b = 0; b < N; ++b) {
                            if (!std::isfinite(x_[q][u][b])) {
                                cause = "the value in row " + std::to_string(u)
                                    + (inCells(u) ? ", an initial value of Gauss-Seidel" : "");
                                return;
                            }
                        }
                    }
                });
                if (!cause.empty()) {
                    break;
                }
            }
            if (cause.empty()) {
                cause = "an overflow in the solve";
            }
            return name_ + ": NaN or infinite solution for right-hand side " + std::to_string(q)
                + " in a component with row " + std::to_string(cells.front()) + ", caused by "
                + cause;
        }
 
        std::array<Scalar, N> rightHandSide_(std::size_t q, std::size_t c) const
        {
            std::array<Scalar, N> sum;
            for (std::size_t a = 0; a < N; ++a) {
                sum[a] = rhs_[q][c][a];
            }
            coefficients_.forEachCoupling(c, [&](std::size_t u, const auto& block) {
                if (localIndex_[u] < 0) {
                    for (std::size_t a = 0; a < N; ++a) {
                        for (std::size_t b = 0; b < N; ++b) {
                            sum[a] -= block[a][b] * x_[q][u][b];
                        }
                    }
                }
            });
            return sum;
        }
 
        void solveDiagonal_(std::size_t c)
        {
            if constexpr (N == 1) {
                // One division
                const Scalar diagonal = coefficients_.diagonal(c)[0][0];
                if (diagonal == Scalar {0}) {
                    throw std::runtime_error(name_ + ": singular diagonal in row "
                                             + std::to_string(c));
                }
                for (std::size_t q = 0; q < numVectors_; ++q) {
                    cellSolution_[q] = rightHandSide_(q, c)[0] / diagonal;
                }
            } else {
                matrix_.resize(N * N);
                const auto& diagonal = coefficients_.diagonal(c);
                for (std::size_t a = 0; a < N; ++a) {
                    for (std::size_t b = 0; b < N; ++b) {
                        matrix_[a * N + b] = diagonal[a][b];
                    }
                }
                for (std::size_t q = 0; q < numVectors_; ++q) {
                    const auto sum = rightHandSide_(q, c);
                    for (std::size_t a = 0; a < N; ++a) {
                        cellSolution_[a * numVectors_ + q] = sum[a];
                    }
                }
                try {
                    solveDenseSystem(matrix_, cellSolution_, N, numVectors_);
                } catch (const std::runtime_error&) {
                    throw std::runtime_error(name_ + ": singular diagonal block in row "
                                             + std::to_string(c));
                }
            }
        }
 
        void solveCell_(std::span<const int> cells)
        {
            const auto c = static_cast<std::size_t>(cells[0]);
            solveDiagonal_(c);
            for (std::size_t q = 0; q < numVectors_; ++q) {
                for (std::size_t a = 0; a < N; ++a) {
                    const Scalar value = cellSolution_[a * numVectors_ + q];
                    checkFinite_(value, q, cells, false);
                    x_[q][c][a] = value;
                }
            }
        }
 
        void solveDirectly_(std::span<const int> cells)
        {
            const std::size_t m = cells.size();
            const std::size_t size = m * N;
            matrix_.assign(size * size, Scalar {0});
            solution_.assign(size * numVectors_, Scalar {0});
            for (std::size_t i = 0; i < m; ++i) {
                localIndex_[cells[i]] = static_cast<int>(i);
            }
            for (std::size_t i = 0; i < m; ++i) {
                const auto c = static_cast<std::size_t>(cells[i]);
                const auto& diagonal = coefficients_.diagonal(c);
                for (std::size_t a = 0; a < N; ++a) {
                    for (std::size_t b = 0; b < N; ++b) {
                        matrix_[(i * N + a) * size + i * N + b] = diagonal[a][b];
                    }
                }
                coefficients_.forEachCoupling(c, [&](std::size_t u, const auto& block) {
                    if (localIndex_[u] >= 0) {
                        const auto j = static_cast<std::size_t>(localIndex_[u]);
                        for (std::size_t a = 0; a < N; ++a) {
                            for (std::size_t b = 0; b < N; ++b) {
                                matrix_[(i * N + a) * size + j * N + b] += block[a][b];
                            }
                        }
                    }
                });
                for (std::size_t q = 0; q < numVectors_; ++q) {
                    // Cells outside the loop are already solved.
                    const auto sum = rightHandSide_(q, c);
                    for (std::size_t a = 0; a < N; ++a) {
                        solution_[(i * N + a) * numVectors_ + q] = sum[a];
                    }
                }
            }
            for (std::size_t i = 0; i < m; ++i) {
                localIndex_[cells[i]] = -1;
            }
            try {
                solveDenseSystem(matrix_, solution_, size, numVectors_);
            } catch (const std::runtime_error&) {
                throw std::runtime_error(name_ + ": singular matrix for a flow loop of "
                                         + std::to_string(m) + " cells");
            }
            for (std::size_t i = 0; i < m; ++i) {
                for (std::size_t q = 0; q < numVectors_; ++q) {
                    for (std::size_t a = 0; a < N; ++a) {
                        const Scalar value = solution_[(i * N + a) * numVectors_ + q];
                        checkFinite_(value, q, cells, false);
                        x_[q][static_cast<std::size_t>(cells[i])][a] = value;
                    }
                }
            }
        }
 
        int solveByGaussSeidel_(std::span<const int> cells)
        {
            for (int sweep = 1; sweep <= settings_.maxIterations; ++sweep) {
                std::fill(change_.begin(), change_.end(), Scalar {0});
                std::fill(scale_.begin(), scale_.end(), Scalar {0});
                for (const int cell : cells) {
                    const auto c = static_cast<std::size_t>(cell);
                    solveDiagonal_(c);
                    for (std::size_t q = 0; q < numVectors_; ++q) {
                        for (std::size_t a = 0; a < N; ++a) {
                            auto& value = x_[q][c][a];
                            const Scalar updated = cellSolution_[a * numVectors_ + q];
                            checkFinite_(updated, q, cells, true);
                            const std::size_t k = q * N + a;
                            change_[k] = std::max(change_[k], std::abs(updated - value));
                            scale_[k] = std::max(scale_[k], std::abs(updated));
                            value = updated;
                        }
                    }
                }
                // Each block component converges on its own scale, so small components are not
                // only converged relative to large ones. The absolute tolerance covers zero
                bool converged = true;
                for (std::size_t k = 0; k < change_.size(); ++k) {
                    converged = converged
                        && (change_[k] <= settings_.relTolerance * scale_[k]
                            || change_[k] <= settings_.absTolerance);
                }
                if (converged) {
                    return sweep;
                }
            }
            throw std::runtime_error(name_ + ": no convergence after "
                                     + std::to_string(settings_.maxIterations)
                                     + " sweeps for a flow loop of " + std::to_string(cells.size())
                                     + " cells");
        }
 
        const Coefficients& coefficients_; 
        const Vectors& rhs_; 
        Vectors& x_; 
        const UpwindTransportSettings& settings_; 
        std::string name_; 
        std::size_t numCells_; 
        std::size_t numVectors_; 
 
        // Work arrays, reused for every cell and flow loop
        std::vector<int> localIndex_; 
        std::vector<Scalar> matrix_; 
        std::vector<Scalar> solution_; 
        std::vector<Scalar> cellSolution_; 
        std::vector<Scalar> change_; 
        std::vector<Scalar> scale_; 
    };
 
} // namespace detail
 
template <class Matrix, class Vectors>
UpwindTransportStatistics
solveReordered(const Matrix& matrix,
               const Vectors& rhs,
               Vectors& x,
               const UpwindTransportSettings& settings = {})
{
    const detail::MatrixCoefficients<Matrix> coefficients(matrix);
    return detail::ReorderedSystem<detail::MatrixCoefficients<Matrix>, Vectors>(
               coefficients, rhs, x, settings, "Opm::solveReordered")
        .solve();
}
 
} // namespace Opm
 
#endif // OPM_UPWIND_TRANSPORT_SOLVER_HEADER_INCLUDED