TpsaMatrix.hpp

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  Copyright 2025 NORCE 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 2 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.
*/
#ifndef OPM_TPSA_MATRIX_HPP
#define OPM_TPSA_MATRIX_HPP
 
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/linalg/matrixblock.hh>
#include <opm/simulators/linalg/tpsa/TpsaTypes.hpp>
 
#include <array>
#include <cstddef>
#include <numeric>
#include <ranges>
#include <tuple>
#include <utility>
#include <vector>
 
namespace Opm::Linear
{
 
template <class Scalar>
class TpsaMatrix;
 
template <class Scalar>
class TpsaBlockRef
{
public:
    using MatrixBlock = Opm::MatrixBlock<Scalar, numTpsaEq, numTpsaEq>;
 
    TpsaBlockRef() = default;
 
    TpsaBlockRef(const TpsaMatrix<Scalar>& matrix, std::size_t flatIdx)
        : matrix_(&matrix)
        , k_(flatIdx)
    {
    }
 
    TpsaBlockRef& operator*()
    {
        return *this;
    }
 
    TpsaBlockRef& operator+=(const MatrixBlock& b);
 
    TpsaBlockRef& operator=(const MatrixBlock& b);
 
    TpsaBlockRef& operator=(Scalar value);
 
    void gather(MatrixBlock& b) const;
 
private:
    template <class Op, class Block>
    void apply_(Op op, Block& b) const;
 
    const TpsaMatrix<Scalar>* matrix_{nullptr};
    std::size_t k_{0};
};
 
template <class Scalar>
class TpsaMatrix
{
    friend class TpsaBlockRef<Scalar>;
 
public:
    using IstlMatrix = TpsaMatrixView<Scalar>;
 
    using MatrixBlock = Opm::MatrixBlock<Scalar, numTpsaEq, numTpsaEq>;
 
    using BlockAddress = TpsaBlockRef<Scalar>;
 
    using field_type = Scalar;
 
    TpsaMatrix(std::size_t rows, std::size_t columns)
        : rows_(rows)
        , columns_(columns)
    {
    }
 
    template <class Simulator>
    explicit TpsaMatrix(const Simulator& simulator)
        : TpsaMatrix(simulator.model().numTotalDof(), simulator.model().numTotalDof())
    {
    }
 
    // The cached value-array base pointers and the sub-matrix view point into
    // this object, so it must not be copied or moved.
    TpsaMatrix(const TpsaMatrix&) = delete;
 
    TpsaMatrix(TpsaMatrix&&) = delete;
 
    TpsaMatrix& operator=(const TpsaMatrix&) = delete;
 
    TpsaMatrix& operator=(TpsaMatrix&&) = delete;
 
    template <class Set>
    void reserve(const std::vector<Set>& sparsityPattern)
    {
        if (sparsityPattern.size() != rows_) {
            OPM_THROW(std::logic_error,
                      "TPSA: sparsity pattern does not match the number of matrix rows");
        }
 
        // Flatten the pattern once; the column indices of a std::set are
        // already ascending, which is the order BCRSMatrix stores them in.
        rowStart_.resize(rows_ + 1);
        rowStart_[0] = 0;
        const auto sizes = sparsityPattern
            | std::views::transform([](const auto& a) { return a.size(); });
        std::partial_sum(sizes.begin(), sizes.end(), rowStart_.begin() + 1);
        nnz_ = rowStart_[rows_];
 
        colIdx_.clear();
        colIdx_.reserve(nnz_);
        for (std::size_t row = 0; row < rows_; ++row) {
            colIdx_.insert(colIdx_.end(),
                           sparsityPattern[row].begin(),
                           sparsityPattern[row].end());
        }
 
        forEachSubMatrix_([&](auto& subMatrix) {
            reserveSubMatrix_(subMatrix, sparsityPattern);
        });
 
        // Initialize base_ pointer to sub-matrices and set up TpsaMatrixView
        cacheValueArrays_();
        setMatrixView_();
    }
 
    BlockAddress blockAddress(const std::size_t rowIdx, const std::size_t colIdx) const
    {
        return BlockAddress(*this, flatIndex_(rowIdx, colIdx));
    }
 
    void clear();
 
    void clearRow(const std::size_t row, const Scalar diag = 1.0);
 
    void makeOverlapRowsInvalid(const std::vector<int>& overlapRows);
 
    void scaleFields(const std::array<Scalar, numTpsaFields>& rowFac,
                     const std::array<Scalar, numTpsaFields>& colFac);
 
    void block(const std::size_t rowIdx, const std::size_t colIdx, MatrixBlock& value) const
    {
        blockAddress(rowIdx, colIdx).gather(value);
    }
 
    void setBlock(const std::size_t rowIdx, const std::size_t colIdx, const MatrixBlock& value)
    {
        blockAddress(rowIdx, colIdx) = value;
    }
 
    void addToBlock(const std::size_t rowIdx, const std::size_t colIdx, const MatrixBlock& value)
    {
        blockAddress(rowIdx, colIdx) += value;
    }
 
    void commit()
    {
    }
 
    void finalize()
    {
    }
 
    IstlMatrix& istlMatrix()
    {
        return view_;
    }
 
    const IstlMatrix& istlMatrix() const
    {
        return view_;
    }
 
    std::size_t rows() const
    {
        return rows_;
    }
 
    std::size_t cols() const
    {
        return columns_;
    }
 
    std::size_t N() const
    {
        return rows_;
    }
 
    std::size_t M() const
    {
        return columns_;
    }
 
    std::size_t nonzeroes() const
    {
        return nnz_;
    }
 
    // Sub-matrix accessors.  dd00/dd11/dd22 and spsp are the scalar blocks Hypre
    // can precondition.
 
    DispDispMatrix00T<Scalar>& dd00()
    {
        return dd00_;
    }
 
    DispDispMatrix11T<Scalar>& dd11()
    {
        return dd11_;
    }
 
    DispDispMatrix22T<Scalar>& dd22()
    {
        return dd22_;
    }
 
    RotRotMatrixT<Scalar>& rr()
    {
        return rr_;
    }
 
    SPresSPresMatrixT<Scalar>& spsp()
    {
        return spsp_;
    }
 
    const DispDispMatrix00T<Scalar>& dd00() const
    {
        return dd00_;
    }
 
    const DispDispMatrix11T<Scalar>& dd11() const
    {
        return dd11_;
    }
 
    const DispDispMatrix22T<Scalar>& dd22() const
    {
        return dd22_;
    }
 
    const RotRotMatrixT<Scalar>& rr() const
    {
        return rr_;
    }
 
    const SPresSPresMatrixT<Scalar>& spsp() const
    {
        return spsp_;
    }
 
private:
    enum SubMatrixIdx : std::size_t
    {
        DD00 = 0, DD11, DD22,
        DR0, DR1, DR2,
        DSP0, DSP1, DSP2,
        RD0, RD1, RD2,
        RR, RSP,
        SPD0, SPD1, SPD2,
        SPR, SPSP,
        numSubMatrices
    };
 
    static constexpr std::array<std::pair<std::size_t, std::size_t>, numSubMatrices>
    subMatrixFields_ {{
        {0, 0}, {1, 1}, {2, 2},   // DD00, DD11, DD22
        {0, 3}, {1, 3}, {2, 3},   // DR0,  DR1,  DR2
        {0, 4}, {1, 4}, {2, 4},   // DSP0, DSP1, DSP2
        {3, 0}, {3, 1}, {3, 2},   // RD0,  RD1,  RD2
        {3, 3}, {3, 4},           // RR,   RSP
        {4, 0}, {4, 1}, {4, 2},   // SPD0, SPD1, SPD2
        {4, 3}, {4, 4}            // SPR,  SPSP
    }};
 
    template <class SubMatrix>
    static constexpr std::size_t blockScalars_(const SubMatrix&)
    {
        using Block = typename SubMatrix::block_type;
 
        return static_cast<std::size_t>(Block::rows) * static_cast<std::size_t>(Block::cols);
    }
 
    auto subMatrices_()
    {
        return std::tie(dd00_,
                        dd11_,
                        dd22_,
                        dr0_,
                        dr1_,
                        dr2_,
                        dsp0_,
                        dsp1_,
                        dsp2_,
                        rd0_,
                        rd1_,
                        rd2_,
                        rr_,
                        rsp_,
                        spd0_,
                        spd1_,
                        spd2_,
                        spr_,
                        spsp_);
    }
 
    template <class Op>
    void forEachSubMatrix_(Op op)
    {
        std::apply([&op](auto&... subMatrix) {
                       (op(subMatrix), ...);
                   },
                   subMatrices_());
    }
 
    template <class Op>
    void forEachSubMatrixWithIndex_(Op op)
    {
        std::apply([&op, subIdx = std::size_t{0}](auto&... subMatrix) mutable {
                       (op(subMatrix, subIdx++), ...);
                   },
                   subMatrices_());
    }
 
    template <class SubMatrix, class Set>
    void reserveSubMatrix_(SubMatrix& subMatrix, const std::vector<Set>& sparsityPattern)
    {
        subMatrix.setBuildMode(SubMatrix::random);
        subMatrix.setSize(rows_, columns_);
 
        for (std::size_t row = 0; row < rows_; ++row) {
            subMatrix.setrowsize(row, sparsityPattern[row].size());
        }
        subMatrix.endrowsizes();
 
        for (std::size_t row = 0; row < rows_; ++row) {
            for (const auto& col : sparsityPattern[row]) {
                subMatrix.addindex(row, col);
            }
        }
        // Note: all entries in subMatrix are Scalar(0.0) by default construction in endindices()
        subMatrix.endindices();
    }
 
    void cacheValueArrays_();
 
    void setMatrixView_();
 
    std::size_t flatIndex_(std::size_t rowIdx, std::size_t colIdx) const;
 
    std::size_t rows_{0};
    std::size_t columns_{0};
    std::size_t nnz_{0};
 
    // Flattened sparsity pattern, shared by all sub-matrices.
    std::vector<std::size_t> rowStart_{};
    std::vector<unsigned> colIdx_{};
 
    // Base pointers into the sub-matrices' contiguous value arrays.
    std::array<Scalar*, numSubMatrices> base_{};
 
    DispDispMatrix00T<Scalar> dd00_{};
    DispDispMatrix11T<Scalar> dd11_{};
    DispDispMatrix22T<Scalar> dd22_{};
 
    DispRotMatrix0T<Scalar> dr0_{};
    DispRotMatrix1T<Scalar> dr1_{};
    DispRotMatrix2T<Scalar> dr2_{};
 
    DispSPresMatrix0T<Scalar> dsp0_{};
    DispSPresMatrix1T<Scalar> dsp1_{};
    DispSPresMatrix2T<Scalar> dsp2_{};
 
    RotDispMatrix0T<Scalar> rd0_{};
    RotDispMatrix1T<Scalar> rd1_{};
    RotDispMatrix2T<Scalar> rd2_{};
 
    RotRotMatrixT<Scalar> rr_{};
    RotSPresMatrixT<Scalar> rsp_{};
 
    SPresDispMatrix0T<Scalar> spd0_{};
    SPresDispMatrix1T<Scalar> spd1_{};
    SPresDispMatrix2T<Scalar> spd2_{};
 
    SPresRotMatrixT<Scalar> spr_{};
 
    SPresSPresMatrixT<Scalar> spsp_{};
 
    IstlMatrix view_{};
};
 
//
// TpsaBlockRef implementation
//
template <class Scalar>
template <class Op, class Block>
void
TpsaBlockRef<Scalar>::apply_(Op op, Block& b) const
{
    using Matrix = TpsaMatrix<Scalar>;
    const auto& base = matrix_->base_;
    const std::size_t k = k_;
 
    // Displacement-displacement (diagonal components only)
    op(base[Matrix::DD00][k], b[0][0]);
    op(base[Matrix::DD11][k], b[1][1]);
    op(base[Matrix::DD22][k], b[2][2]);
 
    // Displacement-rotation: three 1x3 blocks
    for (int j = 0; j < numRotDofs; ++j) {
        op(base[Matrix::DR0][3 * k + j], b[0][3 + j]);
        op(base[Matrix::DR1][3 * k + j], b[1][3 + j]);
        op(base[Matrix::DR2][3 * k + j], b[2][3 + j]);
    }
 
    // Displacement-solid pressure: three 1x1 blocks
    op(base[Matrix::DSP0][k], b[0][6]);
    op(base[Matrix::DSP1][k], b[1][6]);
    op(base[Matrix::DSP2][k], b[2][6]);
 
    // Rotation-displacement: three 3x1 blocks
    for (int i = 0; i < numRotDofs; ++i) {
        op(base[Matrix::RD0][3 * k + i], b[3 + i][0]);
        op(base[Matrix::RD1][3 * k + i], b[3 + i][1]);
        op(base[Matrix::RD2][3 * k + i], b[3 + i][2]);
    }
 
    // Rotation-rotation: one 3x3 block
    for (int i = 0; i < numRotDofs; ++i) {
        for (int j = 0; j < numRotDofs; ++j) {
            op(base[Matrix::RR][9 * k + 3 * i + j], b[3 + i][3 + j]);
        }
    }
 
    // Rotation-solid pressure: one 3x1 block
    for (int i = 0; i < numRotDofs; ++i) {
        op(base[Matrix::RSP][3 * k + i], b[3 + i][6]);
    }
 
    // Solid pressure-displacement: three 1x1 blocks
    op(base[Matrix::SPD0][k], b[6][0]);
    op(base[Matrix::SPD1][k], b[6][1]);
    op(base[Matrix::SPD2][k], b[6][2]);
 
    // Solid pressure-rotation: one 1x3 block
    for (int j = 0; j < numRotDofs; ++j) {
        op(base[Matrix::SPR][3 * k + j], b[6][3 + j]);
    }
 
    // Solid pressure-solid pressure: one 1x1 block
    op(base[Matrix::SPSP][k], b[6][6]);
}
 
} // namespace Opm::Linear
 
#endif // OPM_TPSA_MATRIX_HPP