From 267bd1af366774d146f034b93254c1db3f78fb12 Mon Sep 17 00:00:00 2001 From: mattijs Date: Sat, 23 Nov 2024 10:44:34 +0000 Subject: [PATCH] ENH: CSR: cell-cell addressing --- .../lduMatrix/lduAddressing/lduAddressing.C | 26 ++++ .../lduMatrix/lduAddressing/lduAddressing.H | 27 ++++ .../lduAddressing/lduAddressingTemplates.C | 60 ++++++++ .../matrices/lduMatrix/lduMatrix/lduMatrix.C | 141 +++++++++++++++++- .../matrices/lduMatrix/lduMatrix/lduMatrix.H | 27 +++- .../lduMatrix/lduMatrix/lduMatrixATmul.C | 83 +++++++++-- .../GaussSeidel/GaussSeidelSmoother.C | 4 +- .../symGaussSeidel/symGaussSeidelSmoother.C | 4 +- 8 files changed, 350 insertions(+), 22 deletions(-) create mode 100644 src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressingTemplates.C diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C index 91180115b7..81b80f0efb 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C @@ -169,6 +169,20 @@ void Foam::lduAddressing::calcLosortStart() const } +void Foam::lduAddressing::calcLoCSR() const +{ + if (lowerCSRAddrPtr_) + { + FatalErrorInFunction + << "lowerCSRAddr already calculated" + << abort(FatalError); + } + + lowerCSRAddrPtr_ = std::make_unique(lowerAddr().size()); + map(lowerAddr(), *lowerCSRAddrPtr_); +} + + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // const Foam::labelUList& Foam::lduAddressing::losortAddr() const @@ -204,11 +218,23 @@ const Foam::labelUList& Foam::lduAddressing::losortStartAddr() const } +const Foam::labelUList& Foam::lduAddressing::lowerCSRAddr() const +{ + if (!lowerCSRAddrPtr_) + { + calcLoCSR(); + } + + return *lowerCSRAddrPtr_; +} + + void Foam::lduAddressing::clearOut() { losortPtr_.reset(nullptr); ownerStartPtr_.reset(nullptr); losortStartPtr_.reset(nullptr); + lowerCSRAddrPtr_.reset(nullptr); } diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H index b610a6539d..3c57108897 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H @@ -91,6 +91,10 @@ Description list. Thus, for every point the losort start gives the address of the first face to neighbour this point. + Instead of using losort to lookup the face and then using the lowerAddr + to find the neighbour cell one can also directly lookup the neighbour cell + using the lowerCSRAddr (upperAddr is already in CSR order). + SourceFiles lduAddressing.C @@ -131,6 +135,9 @@ class lduAddressing //- Losort start addressing mutable std::unique_ptr losortStartPtr_; + //- Lower addressing + mutable std::unique_ptr lowerCSRAddrPtr_; + // Private Member Functions @@ -143,6 +150,9 @@ class lduAddressing //- Calculate losort start void calcLosortStart() const; + //- Calculate CSR lower addressing + void calcLoCSR() const; + public: @@ -203,11 +213,22 @@ public: //- Return losort start addressing const labelUList& losortStartAddr() const; + //- Return CSR addressing + const labelUList& lowerCSRAddr() const; + //- Return off-diagonal index given owner and neighbour label label triIndex(const label a, const label b) const; //- Calculate bandwidth and profile of addressing Tuple2 band() const; + + //- Helper to convert lower addressing & data into CSR format + template + void map + ( + const UList& faceVals, + List& vals + ) const; }; @@ -217,6 +238,12 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#ifdef NoRepository + #include "lduAddressingTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + #endif // ************************************************************************* // diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressingTemplates.C b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressingTemplates.C new file mode 100644 index 0000000000..93cd848534 --- /dev/null +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressingTemplates.C @@ -0,0 +1,60 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2024 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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. + + OpenFOAM 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 OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "lduAddressing.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::lduAddressing::map +( + const UList& faceVals, + List& vals +) const +{ + const auto& offsets = losortStartAddr(); + const auto& indexToFace = losortAddr(); + + const label n = offsets.size()-1; + + vals.resize_nocopy(faceVals.size()); + + for (label celli = 0; celli < n; celli++) + { + const label start = offsets[celli]; + const label end = offsets[celli+1]; + + for (label i = start; i < end; i++) + { + const label facei = indexToFace[i]; + vals[i] = faceVals[facei]; + } + } +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C index e2767c6b52..8862b3a368 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C @@ -82,6 +82,13 @@ Foam::lduMatrix::lduMatrix(const lduMatrix& A) { lowerPtr_ = std::make_unique(*(A.lowerPtr_)); } + + if (A.lowerCSRPtr_) + { + lowerCSRPtr_ = std::make_unique(*(A.lowerCSRPtr_)); + } + + // No need to copy work } @@ -90,7 +97,9 @@ Foam::lduMatrix::lduMatrix(lduMatrix&& A) lduMesh_(A.lduMesh_), diagPtr_(std::move(A.diagPtr_)), lowerPtr_(std::move(A.lowerPtr_)), - upperPtr_(std::move(A.upperPtr_)) + upperPtr_(std::move(A.upperPtr_)), + lowerCSRPtr_(std::move(A.lowerCSRPtr_)), + workPtr_(std::move(A.workPtr_)) {} @@ -103,6 +112,8 @@ Foam::lduMatrix::lduMatrix(lduMatrix& A, bool reuse) diagPtr_ = std::move(A.diagPtr_); upperPtr_ = std::move(A.upperPtr_); lowerPtr_ = std::move(A.lowerPtr_); + lowerCSRPtr_ = std::move(A.lowerCSRPtr_); + workPtr_ = std::move(A.workPtr_); } else { @@ -121,6 +132,13 @@ Foam::lduMatrix::lduMatrix(lduMatrix& A, bool reuse) { lowerPtr_ = std::make_unique(*(A.lowerPtr_)); } + + // Note: no real need to keep lowerCSR except we use it (hasLowerCSR()) + // to trigger certain actions + if (A.lowerCSRPtr_) + { + lowerCSRPtr_ = std::make_unique(*(A.lowerCSRPtr_)); + } } } @@ -129,9 +147,9 @@ Foam::lduMatrix::lduMatrix(const lduMesh& mesh, Istream& is) : lduMesh_(mesh) { - bool withLower, withDiag, withUpper; + bool withLower, withDiag, withUpper, withLowerCSR; - is >> withLower >> withDiag >> withUpper; + is >> withLower >> withDiag >> withUpper >> withLowerCSR; if (withLower) { @@ -145,6 +163,20 @@ Foam::lduMatrix::lduMatrix(const lduMesh& mesh, Istream& is) { upperPtr_ = std::make_unique(is); } + if (withLowerCSR) + { + // Check if it has been sent over or we need to construct it. We need + // to set the lowerCSRPtr since it determines behaviour. + if (withLower) + { + // Construct from lower + (void)lowerCSR(); + } + else + { + lowerCSRPtr_ = std::make_unique(is); + } + } } @@ -237,6 +269,9 @@ Foam::scalarField& Foam::lduMatrix::upper() } else { + // no lowerPtr so any lowerCSR was constructed from upper + lowerCSRPtr_.reset(nullptr); + upperPtr_ = std::make_unique ( @@ -260,6 +295,9 @@ Foam::scalarField& Foam::lduMatrix::upper(label nCoeffs) } else { + // no lowerPtr so any lowerCSR was constructed from upper + lowerCSRPtr_.reset(nullptr); + // if (nCoeffs < 0) // { // nCoeffs = lduAddr().lowerAddr().size(); @@ -296,6 +334,8 @@ Foam::scalarField& Foam::lduMatrix::lower() { if (!lowerPtr_) { + lowerCSRPtr_.reset(nullptr); + if (upperPtr_) { lowerPtr_ = std::make_unique(*upperPtr_); @@ -319,6 +359,8 @@ Foam::scalarField& Foam::lduMatrix::lower(label nCoeffs) { if (!lowerPtr_) { + lowerCSRPtr_.reset(nullptr); + if (upperPtr_) { lowerPtr_ = std::make_unique(*upperPtr_); @@ -338,6 +380,86 @@ Foam::scalarField& Foam::lduMatrix::lower(label nCoeffs) } +const Foam::scalarField& Foam::lduMatrix::lowerCSR() const +{ + if (!lowerCSRPtr_) + { + const label nLower = lduAddr().losortAddr().size(); + + lowerCSRPtr_ = std::make_unique(nLower); + + if (lowerPtr_) + { + lduAddr().map(*lowerPtr_, *lowerCSRPtr_); + } + else if (upperPtr_) + { + lduAddr().map(*upperPtr_, *lowerCSRPtr_); + } + else + { + FatalErrorInFunction + << "lowerPtr_ and upperPtr_ unallocated" + << abort(FatalError); + } + } + + return *lowerCSRPtr_; +} + + +Foam::scalarField& Foam::lduMatrix::lowerCSR() +{ + if (!lowerCSRPtr_) + { + const label nLower = lduAddr().losortAddr().size(); + + lowerCSRPtr_ = std::make_unique(nLower); + + if (lowerPtr_) + { + lduAddr().map(*lowerPtr_, *lowerCSRPtr_); + } + else if (upperPtr_) + { + lduAddr().map(*upperPtr_, *lowerCSRPtr_); + } + else + { + FatalErrorInFunction + << "lowerPtr_ and upperPtr_ unallocated" + << abort(FatalError); + } + } + + return *lowerCSRPtr_; +} + + +Foam::solveScalarField& Foam::lduMatrix::work(label size) const +{ + if (!workPtr_ || workPtr_->size() != size) + { + workPtr_ = std::make_unique(size, Foam::zero{}); + } + + return *workPtr_; +} + + +const Foam::solveScalarField& Foam::lduMatrix::work() const +{ + if (!workPtr_) + { + FatalErrorInFunction + << "workPtr_ unallocated" + << abort(FatalError); + } + + return *workPtr_; +} + + void Foam::lduMatrix::setResidualField ( const scalarField& residual, @@ -389,7 +511,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& mat) { os << mat.hasLower() << token::SPACE << mat.hasDiag() << token::SPACE - << mat.hasUpper() << token::SPACE; + << mat.hasUpper() << token::SPACE + << mat.hasLowerCSR() << token::SPACE; if (mat.hasLower()) { @@ -406,6 +529,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& mat) os << mat.upper(); } + if (mat.hasLowerCSR() && !mat.hasLower()) + { + // Only send over if can not be reconstructed locally + os << mat.lowerCSR(); + } + os.check(FUNCTION_NAME); return os; @@ -422,7 +551,9 @@ Foam::Ostream& Foam::operator<< os << "Lower:" << Switch::name(mat.hasLower()) << " Diag:" << Switch::name(mat.hasDiag()) - << " Upper:" << Switch::name(mat.hasUpper()) << endl; + << " Upper:" << Switch::name(mat.hasUpper()) + << " lowerCSR:" << Switch::name(mat.hasLowerCSR()) + << endl; if (mat.hasLower()) { diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H index a4dd1b0601..269c7caec6 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H @@ -99,6 +99,12 @@ class lduMatrix //- Off-diagonal coefficients (not including interfaces) std::unique_ptr upperPtr_; + //- Off-diagonal coefficients (not including interfaces) in CSR ordering + mutable std::unique_ptr lowerCSRPtr_; + + //- Work space + mutable std::unique_ptr workPtr_; + public: @@ -671,6 +677,21 @@ public: scalarField& lower(label nCoeffs); + // CSR + + bool hasLowerCSR() const noexcept { return bool(lowerCSRPtr_); } + + const scalarField& lowerCSR() const; + + scalarField& lowerCSR(); + + //- Work array + const solveScalarField& work() const; + + //- Work array + solveScalarField& work(label size) const; + + // Characteristics //- The matrix type (empty, diagonal, symmetric, ...) @@ -683,19 +704,19 @@ public: //- Matrix has diagonal only bool diagonal() const noexcept { - return (diagPtr_ && !lowerPtr_ && !upperPtr_); + return (diagPtr_ && !lowerPtr_ && !lowerCSRPtr_ && !upperPtr_); } //- Matrix is symmetric bool symmetric() const noexcept { - return (diagPtr_ && !lowerPtr_ && upperPtr_); + return (diagPtr_ && !(lowerPtr_ || lowerCSRPtr_) && upperPtr_); } //- Matrix is asymmetric (ie, full) bool asymmetric() const noexcept { - return (diagPtr_ && lowerPtr_ && upperPtr_); + return (diagPtr_ && (lowerPtr_ || lowerCSRPtr_) && upperPtr_); } diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C index 91f6f45027..59ceeec19a 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2017-2019 OpenCFD Ltd. + Copyright (C) 2017-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -43,6 +43,8 @@ void Foam::lduMatrix::Amul const direction cmpt ) const { + const auto& addr = lduAddr(); + solveScalar* __restrict__ ApsiPtr = Apsi.begin(); const solveScalarField& psi = tpsi(); @@ -50,8 +52,8 @@ void Foam::lduMatrix::Amul const scalar* const __restrict__ diagPtr = diag().begin(); - const label* const __restrict__ uPtr = lduAddr().upperAddr().begin(); - const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin(); + const label* const __restrict__ uPtr = addr.upperAddr().begin(); + const label* const __restrict__ lPtr = addr.lowerAddr().begin(); const scalar* const __restrict__ upperPtr = upper().begin(); const scalar* const __restrict__ lowerPtr = lower().begin(); @@ -70,18 +72,75 @@ void Foam::lduMatrix::Amul ); const label nCells = diag().size(); - for (label cell=0; cell