From 43c0c3989bec91d1c185d64802c0d63a4f0b66f5 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Fri, 16 Feb 2024 17:35:37 +0100 Subject: [PATCH 1/2] ENH: handle matrix demand-driven internals with unique_ptr - align member order between template and non-template versions - move construct, move assignment --- .../matrices/LduMatrix/LduMatrix/LduMatrix.C | 326 ++++++------- .../matrices/LduMatrix/LduMatrix/LduMatrix.H | 196 ++++---- .../LduMatrix/LduMatrix/LduMatrixOperations.C | 60 ++- .../lduMatrix/lduAddressing/lduAddressing.C | 38 +- .../lduMatrix/lduAddressing/lduAddressing.H | 42 +- .../matrices/lduMatrix/lduMatrix/lduMatrix.C | 445 +++++++++--------- .../matrices/lduMatrix/lduMatrix/lduMatrix.H | 152 +++--- .../lduMatrix/lduMatrix/lduMatrixOperations.C | 77 ++- 8 files changed, 672 insertions(+), 664 deletions(-) diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C index dbdaa0fbc3..c76026bd53 100644 --- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation + Copyright (C) 2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -33,109 +34,86 @@ License template Foam::LduMatrix::LduMatrix(const lduMesh& mesh) : - lduMesh_(mesh), - diagPtr_(nullptr), - upperPtr_(nullptr), - lowerPtr_(nullptr), - sourcePtr_(nullptr), - interfaces_(0), - interfacesUpper_(0), - interfacesLower_(0) + lduMesh_(mesh) {} template Foam::LduMatrix::LduMatrix(const LduMatrix& A) : - lduMesh_(A.lduMesh_), - diagPtr_(nullptr), - upperPtr_(nullptr), - lowerPtr_(nullptr), - sourcePtr_(nullptr), - interfaces_(0), - interfacesUpper_(0), - interfacesLower_(0) + lduMesh_(A.lduMesh_) { if (A.diagPtr_) { - diagPtr_ = new Field(*(A.diagPtr_)); + diagPtr_ = std::make_unique>(*(A.diagPtr_)); } if (A.upperPtr_) { - upperPtr_ = new Field(*(A.upperPtr_)); + upperPtr_ = std::make_unique>(*(A.upperPtr_)); } if (A.lowerPtr_) { - lowerPtr_ = new Field(*(A.lowerPtr_)); + lowerPtr_ = std::make_unique>(*(A.lowerPtr_)); } if (A.sourcePtr_) { - sourcePtr_ = new Field(*(A.sourcePtr_)); + sourcePtr_ = std::make_unique>(*(A.sourcePtr_)); } } +template +Foam::LduMatrix::LduMatrix(LduMatrix&& A) +: + lduMesh_(A.lduMesh_), + diagPtr_(std::move(A.diagPtr_)), + lowerPtr_(std::move(A.lowerPtr_)), + upperPtr_(std::move(A.upperPtr_)), + sourcePtr_(std::move(A.sourcePtr_)) +{ + // Clear the old interfaces? +} + + template Foam::LduMatrix::LduMatrix(LduMatrix& A, bool reuse) : - lduMesh_(A.lduMesh_), - diagPtr_(nullptr), - upperPtr_(nullptr), - lowerPtr_(nullptr), - sourcePtr_(nullptr), - interfaces_(0), - interfacesUpper_(0), - interfacesLower_(0) + lduMesh_(A.lduMesh_) { if (reuse) { - if (A.diagPtr_) - { - diagPtr_ = A.diagPtr_; - A.diagPtr_ = nullptr; - } + // Move assignment + diagPtr_ = std::move(A.diagPtr_); + upperPtr_ = std::move(A.upperPtr_); + lowerPtr_ = std::move(A.lowerPtr_); + sourcePtr_ = std::move(A.sourcePtr_); - if (A.upperPtr_) - { - upperPtr_ = A.upperPtr_; - A.upperPtr_ = nullptr; - } - - if (A.lowerPtr_) - { - lowerPtr_ = A.lowerPtr_; - A.lowerPtr_ = nullptr; - } - - if (A.sourcePtr_) - { - sourcePtr_ = A.sourcePtr_; - A.sourcePtr_ = nullptr; - } + // Clear the old interfaces? } else { + // Copy assignment if (A.diagPtr_) { - diagPtr_ = new Field(*(A.diagPtr_)); + diagPtr_ = std::make_unique>(*(A.diagPtr_)); } if (A.upperPtr_) { - upperPtr_ = new Field(*(A.upperPtr_)); + upperPtr_ = std::make_unique>(*(A.upperPtr_)); } if (A.lowerPtr_) { - lowerPtr_ = new Field(*(A.lowerPtr_)); + lowerPtr_ = std::make_unique>(*(A.lowerPtr_)); } if (A.sourcePtr_) { - sourcePtr_ = new Field(*(A.sourcePtr_)); + sourcePtr_ = std::make_unique>(*(A.sourcePtr_)); } } } @@ -152,109 +130,27 @@ Foam::LduMatrix::LduMatrix diagPtr_(new Field(is)), upperPtr_(new Field(is)), lowerPtr_(new Field(is)), - sourcePtr_(new Field(is)), - interfaces_(0), - interfacesUpper_(0), - interfacesLower_(0) + sourcePtr_(new Field(is)) {} -// * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * // - -template -Foam::LduMatrix::~LduMatrix() -{ - if (diagPtr_) - { - delete diagPtr_; - } - - if (upperPtr_) - { - delete upperPtr_; - } - - if (lowerPtr_) - { - delete lowerPtr_; - } - - if (sourcePtr_) - { - delete sourcePtr_; - } -} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -Foam::Field& Foam::LduMatrix::diag() +Foam::word Foam::LduMatrix::matrixTypeName() const { - if (!diagPtr_) + if (diagPtr_) { - diagPtr_ = new Field(lduAddr().size(), Zero); + return + ( + (!upperPtr_) + ? (!lowerPtr_ ? "diagonal" : "diagonal-lower") + : (!lowerPtr_ ? "symmetric" : "asymmetric") + ); } - return *diagPtr_; -} - - -template -Foam::Field& Foam::LduMatrix::upper() -{ - if (!upperPtr_) - { - if (lowerPtr_) - { - upperPtr_ = new Field(*lowerPtr_); - } - else - { - upperPtr_ = new Field - ( - lduAddr().lowerAddr().size(), - Zero - ); - } - } - - return *upperPtr_; -} - - -template -Foam::Field& Foam::LduMatrix::lower() -{ - if (!lowerPtr_) - { - if (upperPtr_) - { - lowerPtr_ = new Field(*upperPtr_); - } - else - { - lowerPtr_ = new Field - ( - lduAddr().lowerAddr().size(), - Zero - ); - } - } - - return *lowerPtr_; -} - - -template -Foam::Field& Foam::LduMatrix::source() -{ - if (!sourcePtr_) - { - sourcePtr_ = new Field(lduAddr().size(), Zero); - } - - return *sourcePtr_; + // is empty (or just wrong) + return (!upperPtr_ && !lowerPtr_ ? "empty" : "ill-defined"); } @@ -273,47 +169,107 @@ const Foam::Field& Foam::LduMatrix::diag() const template -const Foam::Field& Foam::LduMatrix::upper() const +Foam::Field& Foam::LduMatrix::diag() { - if (!lowerPtr_ && !upperPtr_) + if (!diagPtr_) { - FatalErrorInFunction - << "lowerPtr_ or upperPtr_ unallocated" - << abort(FatalError); + diagPtr_ = + std::make_unique>(lduAddr().size(), Foam::zero{}); } + return *diagPtr_; +} + + +template +const Foam::Field& Foam::LduMatrix::upper() const +{ if (upperPtr_) { return *upperPtr_; } else { + if (!lowerPtr_) + { + FatalErrorInFunction + << "lowerPtr_ and upperPtr_ unallocated" + << abort(FatalError); + } + return *lowerPtr_; } } template -const Foam::Field& Foam::LduMatrix::lower() const +Foam::Field& Foam::LduMatrix::upper() { - if (!lowerPtr_ && !upperPtr_) + if (!upperPtr_) { - FatalErrorInFunction - << "lowerPtr_ or upperPtr_ unallocated" - << abort(FatalError); + if (lowerPtr_) + { + upperPtr_ = std::make_unique>(*lowerPtr_); + } + else + { + upperPtr_ = + std::make_unique> + ( + lduAddr().lowerAddr().size(), + Foam::zero{} + ); + } } + return *upperPtr_; +} + + +template +const Foam::Field& Foam::LduMatrix::lower() const +{ if (lowerPtr_) { return *lowerPtr_; } else { + if (!upperPtr_) + { + FatalErrorInFunction + << "lowerPtr_ and upperPtr_ unallocated" + << abort(FatalError); + } + return *upperPtr_; } } +template +Foam::Field& Foam::LduMatrix::lower() +{ + if (!lowerPtr_) + { + if (upperPtr_) + { + lowerPtr_ = std::make_unique>(*upperPtr_); + } + else + { + lowerPtr_ = std::make_unique> + ( + lduAddr().lowerAddr().size(), + Foam::zero{} + ); + } + } + + return *lowerPtr_; +} + + template const Foam::Field& Foam::LduMatrix::source() const { @@ -328,45 +284,65 @@ const Foam::Field& Foam::LduMatrix::source() const } +template +Foam::Field& Foam::LduMatrix::source() +{ + if (!sourcePtr_) + { + sourcePtr_ = + std::make_unique>(lduAddr().size(), Foam::zero{}); + } + + return *sourcePtr_; +} + + // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // +// template +// Foam::Ostream& Foam::operator<< +// ( +// Ostream& os, +// const InfoProxy& iproxy +// ) +// { +// const auto& mat = *iproxy; +// +// ... +// +// os.check(FUNCTION_NAME); +// return os; +// } + + template Foam::Ostream& Foam::operator<< ( Ostream& os, - const LduMatrix& ldum + const LduMatrix& mat ) { - if (ldum.diagPtr_) + if (mat.hasDiag()) { - os << "Diagonal = " - << *ldum.diagPtr_ - << endl << endl; + os << "Diagonal = " << mat.diag() << nl << nl; } - if (ldum.upperPtr_) + if (mat.hasUpper()) { - os << "Upper triangle = " - << *ldum.upperPtr_ - << endl << endl; + os << "Upper triangle = " << mat.upper() << nl << nl; } - if (ldum.lowerPtr_) + if (mat.hasLower()) { - os << "Lower triangle = " - << *ldum.lowerPtr_ - << endl << endl; + os << "Lower triangle = " << mat.lower() << nl << nl; } - if (ldum.sourcePtr_) + if (mat.hasSource()) { - os << "Source = " - << *ldum.sourcePtr_ - << endl << endl; + os << "Source = " << mat.source() << nl << nl; } os.check(FUNCTION_NAME); - return os; } diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H index 96e9fd14f4..407f7ee322 100644 --- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation - Copyright (C) 2019-2023 OpenCFD Ltd. + Copyright (C) 2019-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -93,25 +93,32 @@ class LduMatrix const lduMesh& lduMesh_; //- Diagonal coefficients - Field *diagPtr_; + std::unique_ptr> diagPtr_; //- Off-diagonal coefficients - Field *upperPtr_, *lowerPtr_; + std::unique_ptr> upperPtr_; + + //- Off-diagonal coefficients + std::unique_ptr> lowerPtr_; //- Source - Field *sourcePtr_; + std::unique_ptr> sourcePtr_; //- Field interfaces (processor patches etc.) LduInterfaceFieldPtrsList interfaces_; //- Off-diagonal coefficients for interfaces - FieldField interfacesUpper_, interfacesLower_; + FieldField interfacesUpper_; + + //- Off-diagonal coefficients for interfaces + FieldField interfacesLower_; public: friend class SolverPerformance; + // ----------------------------------------------------------------------- //- Abstract base-class for LduMatrix solvers class solver { @@ -274,6 +281,7 @@ public: }; + // ----------------------------------------------------------------------- //- Abstract base-class for LduMatrix smoothers class smoother { @@ -371,6 +379,7 @@ public: }; + // ----------------------------------------------------------------------- //- Abstract base-class for LduMatrix preconditioners class preconditioner { @@ -466,6 +475,8 @@ public: }; + // ----------------------------------------------------------------------- + // Static Data // Declare name of the class and its debug switch @@ -476,126 +487,125 @@ public: //- Construct given an LDU addressed mesh. // The coefficients are initially empty for subsequent setting. - LduMatrix(const lduMesh&); + // Not yet 'explicit' (legacy code may rely on implicit construct) + LduMatrix(const lduMesh& mesh); - //- Construct as copy + //- Copy construct LduMatrix(const LduMatrix&); + //- Move construct + LduMatrix(LduMatrix&&); + //- Construct as copy or re-use as specified. LduMatrix(LduMatrix&, bool reuse); //- Construct given an LDU addressed mesh and an Istream - // from which the coefficients are read - LduMatrix(const lduMesh&, Istream&); + //- from which the coefficients are read + LduMatrix(const lduMesh& mesh, Istream& is); //- Destructor - ~LduMatrix(); + ~LduMatrix() = default; // Member Functions - // Access to addressing + // Addressing - //- Return the LDU mesh from which the addressing is obtained - const lduMesh& mesh() const noexcept - { - return lduMesh_; - } + //- Return the LDU mesh from which the addressing is obtained + const lduMesh& mesh() const noexcept + { + return lduMesh_; + } - //- Return the LDU addressing - const lduAddressing& lduAddr() const - { - return lduMesh_.lduAddr(); - } + //- Return the LDU addressing + const lduAddressing& lduAddr() const + { + return lduMesh_.lduAddr(); + } - //- Return the patch evaluation schedule - const lduSchedule& patchSchedule() const - { - return lduAddr().patchSchedule(); - } + //- Return the patch evaluation schedule + const lduSchedule& patchSchedule() const + { + return lduMesh_.lduAddr().patchSchedule(); + } - //- Return interfaces - const LduInterfaceFieldPtrsList& interfaces() const - { - return interfaces_; - } + //- Const access to the interfaces + const LduInterfaceFieldPtrsList& interfaces() const noexcept + { + return interfaces_; + } - //- Return interfaces - LduInterfaceFieldPtrsList& interfaces() - { - return interfaces_; - } + //- Non-const access to the interfaces + LduInterfaceFieldPtrsList& interfaces() noexcept + { + return interfaces_; + } - // Access to coefficients + // Coefficients - Field& diag(); - Field& upper(); - Field& lower(); - Field& source(); + const Field& diag() const; + const Field& upper() const; + const Field& lower() const; + const Field& source() const; - FieldField& interfacesUpper() - { - return interfacesUpper_; - } - - FieldField& interfacesLower() - { - return interfacesLower_; - } + Field& diag(); + Field& upper(); + Field& lower(); + Field& source(); - const Field& diag() const; - const Field& upper() const; - const Field& lower() const; - const Field& source() const; + // Interfaces - const FieldField& interfacesUpper() const - { - return interfacesUpper_; - } + const FieldField& interfacesUpper() const noexcept + { + return interfacesUpper_; + } - const FieldField& interfacesLower() const - { - return interfacesLower_; - } + const FieldField& interfacesLower() const noexcept + { + return interfacesLower_; + } + + FieldField& interfacesUpper() noexcept + { + return interfacesUpper_; + } + + FieldField& interfacesLower() noexcept + { + return interfacesLower_; + } - bool hasDiag() const noexcept - { - return (diagPtr_); - } + // Characteristics - bool hasUpper() const noexcept - { - return (upperPtr_); - } + //- The matrix type (empty, diagonal, symmetric, ...) + word matrixTypeName() const; - bool hasLower() const noexcept - { - return (lowerPtr_); - } + bool hasDiag() const noexcept { return bool(diagPtr_); } + bool hasUpper() const noexcept { return bool(upperPtr_); } + bool hasLower() const noexcept { return bool(lowerPtr_); } + bool hasSource() const noexcept { return bool(sourcePtr_); } - bool hasSource() const noexcept - { - return (sourcePtr_); - } + //- Matrix has diagonal only + bool diagonal() const noexcept + { + return (diagPtr_ && !lowerPtr_ && !upperPtr_); + } - bool diagonal() const noexcept - { - return (diagPtr_ && !lowerPtr_ && !upperPtr_); - } + //- Matrix is symmetric + bool symmetric() const noexcept + { + return (diagPtr_ && !lowerPtr_ && upperPtr_); + } - bool symmetric() const noexcept - { - return (diagPtr_ && (!lowerPtr_ && upperPtr_)); - } - - bool asymmetric() const noexcept - { - return (diagPtr_ && lowerPtr_ && upperPtr_); - } + //- Matrix is asymmetric (ie, full) + bool asymmetric() const noexcept + { + return (diagPtr_ && lowerPtr_ && upperPtr_); + } // Operations @@ -651,8 +661,12 @@ public: // Member Operators + //- Copy assignment void operator=(const LduMatrix&); + //- Move assignment + void operator=(LduMatrix&&); + void negate(); void operator+=(const LduMatrix&); diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C index 7945a8cd2d..46e7ad7b60 100644 --- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -92,7 +92,7 @@ Foam::LduMatrix::H(const Field& psi) const { auto tHpsi = tmp>::New(lduAddr().size(), Foam::zero{}); - if (lowerPtr_ || upperPtr_) + if (hasLower() || hasUpper()) { Type* __restrict__ HpsiPtr = tHpsi.ref().begin(); @@ -169,32 +169,30 @@ void Foam::LduMatrix::operator=(const LduMatrix& A) return; // Self-assignment is a no-op } - if (A.diagPtr_) + if (A.hasDiag()) { diag() = A.diag(); } - if (A.upperPtr_) + if (A.hasUpper()) { upper() = A.upper(); } - else if (upperPtr_) + else { - delete upperPtr_; - upperPtr_ = nullptr; + upperPtr_.reset(nullptr); } - if (A.lowerPtr_) + if (A.hasLower()) { lower() = A.lower(); } - else if (lowerPtr_) + else { - delete lowerPtr_; - lowerPtr_ = nullptr; + lowerPtr_.reset(nullptr); } - if (A.sourcePtr_) + if (A.hasSource()) { source() = A.source(); } @@ -204,6 +202,24 @@ void Foam::LduMatrix::operator=(const LduMatrix& A) } +template +void Foam::LduMatrix::operator=(LduMatrix&& A) +{ + if (this == &A) + { + return; // Self-assignment is a no-op + } + + diagPtr_ = std::move(A.diagPtr_); + upperPtr_ = std::move(A.upperPtr_); + lowerPtr_ = std::move(A.lowerPtr_); + sourcePtr_ = std::move(A.sourcePtr_); + + interfacesUpper_ = std::move(A.interfacesUpper_); + interfacesLower_ = std::move(A.interfacesLower_); +} + + template void Foam::LduMatrix::negate() { @@ -235,12 +251,12 @@ void Foam::LduMatrix::negate() template void Foam::LduMatrix::operator+=(const LduMatrix& A) { - if (A.diagPtr_) + if (A.hasDiag()) { diag() += A.diag(); } - if (A.sourcePtr_) + if (A.hasSource()) { source() += A.source(); } @@ -265,7 +281,7 @@ void Foam::LduMatrix::operator+=(const LduMatrix& A) } else if (asymmetric() && A.symmetric()) { - if (A.upperPtr_) + if (A.hasUpper()) { lower() += A.upper(); upper() += A.upper(); @@ -284,12 +300,12 @@ void Foam::LduMatrix::operator+=(const LduMatrix& A) } else if (diagonal()) { - if (A.upperPtr_) + if (A.hasUpper()) { upper() = A.upper(); } - if (A.lowerPtr_) + if (A.hasLower()) { lower() = A.lower(); } @@ -312,12 +328,12 @@ void Foam::LduMatrix::operator+=(const LduMatrix& A) template void Foam::LduMatrix::operator-=(const LduMatrix& A) { - if (A.diagPtr_) + if (A.hasDiag()) { diag() -= A.diag(); } - if (A.sourcePtr_) + if (A.hasSource()) { source() -= A.source(); } @@ -342,7 +358,7 @@ void Foam::LduMatrix::operator-=(const LduMatrix& A) } else if (asymmetric() && A.symmetric()) { - if (A.upperPtr_) + if (A.hasUpper()) { lower() -= A.upper(); upper() -= A.upper(); @@ -361,12 +377,12 @@ void Foam::LduMatrix::operator-=(const LduMatrix& A) } else if (diagonal()) { - if (A.upperPtr_) + if (A.hasUpper()) { upper() = -A.upper(); } - if (A.lowerPtr_) + if (A.hasLower()) { lower() = -A.lower(); } diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C index 61006d077c..91180115b7 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2016 OpenCFD Ltd. + Copyright (C) 2016-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -27,7 +27,6 @@ License \*---------------------------------------------------------------------------*/ #include "lduAddressing.H" -#include "demandDrivenData.H" #include "scalarField.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -44,7 +43,7 @@ void Foam::lduAddressing::calcLosort() const // Scan the neighbour list to find out how many times the cell // appears as a neighbour of the face. Done this way to avoid guessing // and resizing list - labelList nNbrOfFace(size(), Zero); + labelList nNbrOfFace(size(), Foam::zero{}); const labelUList& nbr = upperAddr(); @@ -73,9 +72,8 @@ void Foam::lduAddressing::calcLosort() const } // Gather the neighbours into the losort array - losortPtr_ = new labelList(nbr.size(), -1); - - labelList& lst = *losortPtr_; + losortPtr_ = std::make_unique(nbr.size(), -1); + auto& lst = *losortPtr_; // Set counter for losort label lstI = 0; @@ -104,9 +102,8 @@ void Foam::lduAddressing::calcOwnerStart() const const labelList& own = lowerAddr(); - ownerStartPtr_ = new labelList(size() + 1, own.size()); - - labelList& ownStart = *ownerStartPtr_; + ownerStartPtr_ = std::make_unique(size() + 1, own.size()); + auto& ownStart = *ownerStartPtr_; // Set up first lookup by hand ownStart[0] = 0; @@ -139,9 +136,8 @@ void Foam::lduAddressing::calcLosortStart() const << abort(FatalError); } - losortStartPtr_ = new labelList(size() + 1, Zero); - - labelList& lsrtStart = *losortStartPtr_; + losortStartPtr_ = std::make_unique(size() + 1, Foam::zero{}); + auto& lsrtStart = *losortStartPtr_; const labelList& nbr = upperAddr(); @@ -173,16 +169,6 @@ void Foam::lduAddressing::calcLosortStart() const } -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::lduAddressing::~lduAddressing() -{ - deleteDemandDrivenData(losortPtr_); - deleteDemandDrivenData(ownerStartPtr_); - deleteDemandDrivenData(losortStartPtr_); -} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // const Foam::labelUList& Foam::lduAddressing::losortAddr() const @@ -220,9 +206,9 @@ const Foam::labelUList& Foam::lduAddressing::losortStartAddr() const void Foam::lduAddressing::clearOut() { - deleteDemandDrivenData(losortPtr_); - deleteDemandDrivenData(ownerStartPtr_); - deleteDemandDrivenData(losortStartPtr_); + losortPtr_.reset(nullptr); + ownerStartPtr_.reset(nullptr); + losortStartPtr_.reset(nullptr); } @@ -262,7 +248,7 @@ Foam::Tuple2 Foam::lduAddressing::band() const const labelUList& owner = lowerAddr(); const labelUList& neighbour = upperAddr(); - labelList cellBandwidth(size(), Zero); + labelList cellBandwidth(size(), Foam::zero{}); forAll(neighbour, facei) { diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H index c37cd61aec..b610a6539d 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2016-2019 OpenCFD Ltd. + Copyright (C) 2016-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -123,23 +123,17 @@ class lduAddressing //- Demand-driven data //- Losort addressing - mutable labelList* losortPtr_; + mutable std::unique_ptr losortPtr_; //- Owner start addressing - mutable labelList* ownerStartPtr_; + mutable std::unique_ptr ownerStartPtr_; //- Losort start addressing - mutable labelList* losortStartPtr_; + mutable std::unique_ptr losortStartPtr_; // Private Member Functions - //- No copy construct - lduAddressing(const lduAddressing&) = delete; - - //- No copy assignment - void operator=(const lduAddressing&) = delete; - //- Calculate losort void calcLosort() const; @@ -152,24 +146,32 @@ class lduAddressing public: - //- Construct with size (number of equations) - explicit lduAddressing(const label nEqns) - : - size_(nEqns), - losortPtr_(nullptr), - ownerStartPtr_(nullptr), - losortStartPtr_(nullptr) - {} + // Generated Methods + + //- No copy construct + lduAddressing(const lduAddressing&) = delete; + + //- No copy assignment + void operator=(const lduAddressing&) = delete; + + + // Constructors + + //- Construct with size (number of equations) + explicit lduAddressing(const label nEqns) noexcept + : + size_(nEqns) + {} //- Destructor - virtual ~lduAddressing(); + virtual ~lduAddressing() = default; // Member Functions //- Return number of equations - label size() const + label size() const noexcept { return size_; } diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C index dd32325419..e2767c6b52 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation - Copyright (C) 2019-2021 OpenCFD Ltd. + Copyright (C) 2019-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -60,79 +60,66 @@ Foam::lduMatrix::normTypesNames_ Foam::lduMatrix::lduMatrix(const lduMesh& mesh) : - lduMesh_(mesh), - lowerPtr_(nullptr), - diagPtr_(nullptr), - upperPtr_(nullptr) + lduMesh_(mesh) {} Foam::lduMatrix::lduMatrix(const lduMatrix& A) : - lduMesh_(A.lduMesh_), - lowerPtr_(nullptr), - diagPtr_(nullptr), - upperPtr_(nullptr) + lduMesh_(A.lduMesh_) { - if (A.lowerPtr_) - { - lowerPtr_ = new scalarField(*(A.lowerPtr_)); - } - if (A.diagPtr_) { - diagPtr_ = new scalarField(*(A.diagPtr_)); + diagPtr_ = std::make_unique(*(A.diagPtr_)); } if (A.upperPtr_) { - upperPtr_ = new scalarField(*(A.upperPtr_)); + upperPtr_ = std::make_unique(*(A.upperPtr_)); + } + + if (A.lowerPtr_) + { + lowerPtr_ = std::make_unique(*(A.lowerPtr_)); } } -Foam::lduMatrix::lduMatrix(lduMatrix& A, bool reuse) +Foam::lduMatrix::lduMatrix(lduMatrix&& A) : lduMesh_(A.lduMesh_), - lowerPtr_(nullptr), - diagPtr_(nullptr), - upperPtr_(nullptr) + diagPtr_(std::move(A.diagPtr_)), + lowerPtr_(std::move(A.lowerPtr_)), + upperPtr_(std::move(A.upperPtr_)) +{} + + +Foam::lduMatrix::lduMatrix(lduMatrix& A, bool reuse) +: + lduMesh_(A.lduMesh_) { if (reuse) { - if (A.lowerPtr_) - { - lowerPtr_ = A.lowerPtr_; - A.lowerPtr_ = nullptr; - } - - if (A.diagPtr_) - { - diagPtr_ = A.diagPtr_; - A.diagPtr_ = nullptr; - } - - if (A.upperPtr_) - { - upperPtr_ = A.upperPtr_; - A.upperPtr_ = nullptr; - } + diagPtr_ = std::move(A.diagPtr_); + upperPtr_ = std::move(A.upperPtr_); + lowerPtr_ = std::move(A.lowerPtr_); } else { - if (A.lowerPtr_) - { - lowerPtr_ = new scalarField(*(A.lowerPtr_)); - } - + // Copy assignment if (A.diagPtr_) { - diagPtr_ = new scalarField(*(A.diagPtr_)); + diagPtr_ = std::make_unique(*(A.diagPtr_)); } if (A.upperPtr_) { - upperPtr_ = new scalarField(*(A.upperPtr_)); + upperPtr_ = std::make_unique(*(A.upperPtr_)); + } + + if (A.lowerPtr_) + { + lowerPtr_ = std::make_unique(*(A.lowerPtr_)); } } } @@ -140,160 +127,43 @@ Foam::lduMatrix::lduMatrix(lduMatrix& A, bool reuse) Foam::lduMatrix::lduMatrix(const lduMesh& mesh, Istream& is) : - lduMesh_(mesh), - lowerPtr_(nullptr), - diagPtr_(nullptr), - upperPtr_(nullptr) + lduMesh_(mesh) { - Switch hasLow(is); - Switch hasDiag(is); - Switch hasUp(is); + bool withLower, withDiag, withUpper; - if (hasLow) + is >> withLower >> withDiag >> withUpper; + + if (withLower) { - lowerPtr_ = new scalarField(is); + lowerPtr_ = std::make_unique(is); } - if (hasDiag) + if (withDiag) { - diagPtr_ = new scalarField(is); + diagPtr_ = std::make_unique(is); } - if (hasUp) + if (withUpper) { - upperPtr_ = new scalarField(is); + upperPtr_ = std::make_unique(is); } } -Foam::lduMatrix::~lduMatrix() -{ - if (lowerPtr_) - { - delete lowerPtr_; - } +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +Foam::word Foam::lduMatrix::matrixTypeName() const +{ if (diagPtr_) { - delete diagPtr_; + return + ( + (!upperPtr_) + ? (!lowerPtr_ ? "diagonal" : "diagonal-lower") + : (!lowerPtr_ ? "symmetric" : "asymmetric") + ); } - if (upperPtr_) - { - delete upperPtr_; - } -} - - -Foam::scalarField& Foam::lduMatrix::lower() -{ - if (!lowerPtr_) - { - if (upperPtr_) - { - lowerPtr_ = new scalarField(*upperPtr_); - } - else - { - lowerPtr_ = new scalarField(lduAddr().lowerAddr().size(), Zero); - } - } - - return *lowerPtr_; -} - - -Foam::scalarField& Foam::lduMatrix::diag() -{ - if (!diagPtr_) - { - diagPtr_ = new scalarField(lduAddr().size(), Zero); - } - - return *diagPtr_; -} - - -Foam::scalarField& Foam::lduMatrix::upper() -{ - if (!upperPtr_) - { - if (lowerPtr_) - { - upperPtr_ = new scalarField(*lowerPtr_); - } - else - { - upperPtr_ = new scalarField(lduAddr().lowerAddr().size(), Zero); - } - } - - return *upperPtr_; -} - - -Foam::scalarField& Foam::lduMatrix::lower(const label nCoeffs) -{ - if (!lowerPtr_) - { - if (upperPtr_) - { - lowerPtr_ = new scalarField(*upperPtr_); - } - else - { - lowerPtr_ = new scalarField(nCoeffs, Zero); - } - } - - return *lowerPtr_; -} - - -Foam::scalarField& Foam::lduMatrix::diag(const label size) -{ - if (!diagPtr_) - { - diagPtr_ = new scalarField(size, Zero); - } - - return *diagPtr_; -} - - -Foam::scalarField& Foam::lduMatrix::upper(const label nCoeffs) -{ - if (!upperPtr_) - { - if (lowerPtr_) - { - upperPtr_ = new scalarField(*lowerPtr_); - } - else - { - upperPtr_ = new scalarField(nCoeffs, Zero); - } - } - - return *upperPtr_; -} - - -const Foam::scalarField& Foam::lduMatrix::lower() const -{ - if (!lowerPtr_ && !upperPtr_) - { - FatalErrorInFunction - << "lowerPtr_ or upperPtr_ unallocated" - << abort(FatalError); - } - - if (lowerPtr_) - { - return *lowerPtr_; - } - else - { - return *upperPtr_; - } + // is empty (or just wrong) + return (!upperPtr_ && !lowerPtr_ ? "empty" : "ill-defined"); } @@ -310,26 +180,164 @@ const Foam::scalarField& Foam::lduMatrix::diag() const } -const Foam::scalarField& Foam::lduMatrix::upper() const +Foam::scalarField& Foam::lduMatrix::diag() { - if (!lowerPtr_ && !upperPtr_) + if (!diagPtr_) { - FatalErrorInFunction - << "lowerPtr_ or upperPtr_ unallocated" - << abort(FatalError); + diagPtr_ = + std::make_unique(lduAddr().size(), Foam::zero{}); } + return *diagPtr_; +} + + +Foam::scalarField& Foam::lduMatrix::diag(label size) +{ + if (!diagPtr_) + { + // if (size < 0) + // { + // size = lduAddr().size(); + // } + diagPtr_ = std::make_unique(size, Foam::zero{}); + } + + return *diagPtr_; +} + + +const Foam::scalarField& Foam::lduMatrix::upper() const +{ if (upperPtr_) { return *upperPtr_; } else { + if (!lowerPtr_) + { + FatalErrorInFunction + << "lowerPtr_ and upperPtr_ unallocated" + << abort(FatalError); + } + return *lowerPtr_; } } +Foam::scalarField& Foam::lduMatrix::upper() +{ + if (!upperPtr_) + { + if (lowerPtr_) + { + upperPtr_ = std::make_unique(*lowerPtr_); + } + else + { + upperPtr_ = + std::make_unique + ( + lduAddr().lowerAddr().size(), + Foam::zero{} + ); + } + } + + return *upperPtr_; +} + + +Foam::scalarField& Foam::lduMatrix::upper(label nCoeffs) +{ + if (!upperPtr_) + { + if (lowerPtr_) + { + upperPtr_ = std::make_unique(*lowerPtr_); + } + else + { + // if (nCoeffs < 0) + // { + // nCoeffs = lduAddr().lowerAddr().size(); + // } + upperPtr_ = std::make_unique(nCoeffs, Foam::zero{}); + } + } + + return *upperPtr_; +} + + +const Foam::scalarField& Foam::lduMatrix::lower() const +{ + if (lowerPtr_) + { + return *lowerPtr_; + } + else + { + if (!upperPtr_) + { + FatalErrorInFunction + << "lowerPtr_ and upperPtr_ unallocated" + << abort(FatalError); + } + + return *upperPtr_; + } +} + + +Foam::scalarField& Foam::lduMatrix::lower() +{ + if (!lowerPtr_) + { + if (upperPtr_) + { + lowerPtr_ = std::make_unique(*upperPtr_); + } + else + { + lowerPtr_ = + std::make_unique + ( + lduAddr().lowerAddr().size(), + Foam::zero{} + ); + } + } + + return *lowerPtr_; +} + + +Foam::scalarField& Foam::lduMatrix::lower(label nCoeffs) +{ + if (!lowerPtr_) + { + if (upperPtr_) + { + lowerPtr_ = std::make_unique(*upperPtr_); + } + else + { + // if (nCoeffs < 0) + // { + // nCoeffs = lduAddr().lowerAddr().size(); + // } + lowerPtr_ = + std::make_unique(nCoeffs, Foam::zero{}); + } + } + + return *lowerPtr_; +} + + void Foam::lduMatrix::setResidualField ( const scalarField& residual, @@ -377,28 +385,25 @@ void Foam::lduMatrix::setResidualField // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // -Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& ldum) +Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& mat) { - Switch hasLow = ldum.hasLower(); - Switch hasDiag = ldum.hasDiag(); - Switch hasUp = ldum.hasUpper(); + os << mat.hasLower() << token::SPACE + << mat.hasDiag() << token::SPACE + << mat.hasUpper() << token::SPACE; - os << hasLow << token::SPACE << hasDiag << token::SPACE - << hasUp << token::SPACE; - - if (hasLow) + if (mat.hasLower()) { - os << ldum.lower(); + os << mat.lower(); } - if (hasDiag) + if (mat.hasDiag()) { - os << ldum.diag(); + os << mat.diag(); } - if (hasUp) + if (mat.hasUpper()) { - os << ldum.upper(); + os << mat.upper(); } os.check(FUNCTION_NAME); @@ -413,54 +418,50 @@ Foam::Ostream& Foam::operator<< const InfoProxy& iproxy ) { - const auto& ldum = *iproxy; + const auto& mat = *iproxy; - Switch hasLow = ldum.hasLower(); - Switch hasDiag = ldum.hasDiag(); - Switch hasUp = ldum.hasUpper(); + os << "Lower:" << Switch::name(mat.hasLower()) + << " Diag:" << Switch::name(mat.hasDiag()) + << " Upper:" << Switch::name(mat.hasUpper()) << endl; - os << "Lower:" << hasLow - << " Diag:" << hasDiag - << " Upper:" << hasUp << endl; - - if (hasLow) + if (mat.hasLower()) { - os << "lower:" << ldum.lower().size() << endl; + os << "lower:" << mat.lower().size() << endl; } - if (hasDiag) + if (mat.hasDiag()) { - os << "diag :" << ldum.diag().size() << endl; + os << "diag :" << mat.diag().size() << endl; } - if (hasUp) + if (mat.hasUpper()) { - os << "upper:" << ldum.upper().size() << endl; + os << "upper:" << mat.upper().size() << endl; } - //if (hasLow) + //if (hasLower) //{ // os << "lower contents:" << endl; - // forAll(ldum.lower(), i) + // forAll(mat.lower(), i) // { - // os << "i:" << i << "\t" << ldum.lower()[i] << endl; + // os << "i:" << i << "\t" << mat.lower()[i] << endl; // } // os << endl; //} //if (hasDiag) //{ // os << "diag contents:" << endl; - // forAll(ldum.diag(), i) + // forAll(mat.diag(), i) // { - // os << "i:" << i << "\t" << ldum.diag()[i] << endl; + // os << "i:" << i << "\t" << mat.diag()[i] << endl; // } // os << endl; //} - //if (hasUp) + //if (hasUpper) //{ // os << "upper contents:" << endl; - // forAll(ldum.upper(), i) + // forAll(mat.upper(), i) // { - // os << "i:" << i << "\t" << ldum.upper()[i] << endl; + // os << "i:" << i << "\t" << mat.upper()[i] << endl; // } // os << endl; //} diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H index e13bb4cd05..a4dd1b0601 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation - Copyright (C) 2016-2023 OpenCFD Ltd. + Copyright (C) 2016-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -57,13 +57,14 @@ SourceFiles #include "primitiveFieldsFwd.H" #include "FieldField.H" #include "lduInterfaceFieldPtrsList.H" +#include "solverPerformance.H" #include "typeInfo.H" #include "autoPtr.H" #include "runTimeSelectionTables.H" -#include "solverPerformance.H" #include "InfoProxy.H" #include "Enum.H" #include "profilingTrigger.H" +#include // For reference_wrapper // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -89,8 +90,14 @@ class lduMatrix //const lduMesh& lduMesh_; std::reference_wrapper lduMesh_; - //- Coefficients (not including interfaces) - scalarField *lowerPtr_, *diagPtr_, *upperPtr_; + //- Diagonal coefficients + std::unique_ptr diagPtr_; + + //- Off-diagonal coefficients (not including interfaces) + std::unique_ptr lowerPtr_; + + //- Off-diagonal coefficients (not including interfaces) + std::unique_ptr upperPtr_; public: @@ -115,6 +122,7 @@ public: static const scalar defaultTolerance; + // ----------------------------------------------------------------------- //- Abstract base-class for lduMatrix solvers class solver { @@ -331,6 +339,7 @@ public: }; + // ----------------------------------------------------------------------- //- Abstract base-class for lduMatrix smoothers class smoother { @@ -478,6 +487,7 @@ public: }; + // ----------------------------------------------------------------------- //- Abstract base-class for lduMatrix preconditioners class preconditioner { @@ -582,6 +592,8 @@ public: }; + // ----------------------------------------------------------------------- + // Static Data // Declare name of the class and its debug switch @@ -590,101 +602,101 @@ public: // Constructors - //- Construct given an LDU addressed mesh. - // The coefficients are initially empty for subsequent setting. - lduMatrix(const lduMesh&); + //- Construct (without coefficients) for an LDU addressed mesh. + // Not yet 'explicit' (legacy code may rely on implicit construct) + lduMatrix(const lduMesh& mesh); - //- Construct as copy + //- Copy construct lduMatrix(const lduMatrix&); + //- Move construct + lduMatrix(lduMatrix&&); + //- Construct as copy or re-use as specified. lduMatrix(lduMatrix&, bool reuse); //- Construct given an LDU addressed mesh and an Istream //- from which the coefficients are read - lduMatrix(const lduMesh&, Istream&); + lduMatrix(const lduMesh& mesh, Istream& is); //- Destructor - ~lduMatrix(); + ~lduMatrix() = default; // Member Functions - // Access to addressing + // Addressing - //- Return the LDU mesh from which the addressing is obtained - const lduMesh& mesh() const noexcept - { - return lduMesh_; - } + //- Return the LDU mesh from which the addressing is obtained + const lduMesh& mesh() const noexcept + { + return lduMesh_; + } - //- Set the LDU mesh containing the addressing is obtained - void setLduMesh(const lduMesh& m) - { - lduMesh_ = m; - } + //- Set the LDU mesh containing the addressing + void setLduMesh(const lduMesh& m) + { + lduMesh_ = m; + } - //- Return the LDU addressing - const lduAddressing& lduAddr() const - { - return mesh().lduAddr(); - } + //- Return the LDU addressing + const lduAddressing& lduAddr() const + { + return mesh().lduAddr(); + } - //- Return the patch evaluation schedule - const lduSchedule& patchSchedule() const - { - return lduAddr().patchSchedule(); - } + //- Return the patch evaluation schedule + const lduSchedule& patchSchedule() const + { + return mesh().lduAddr().patchSchedule(); + } - // Access to coefficients + // Coefficients - scalarField& lower(); - scalarField& diag(); - scalarField& upper(); + const scalarField& diag() const; + const scalarField& upper() const; + const scalarField& lower() const; - // Size with externally provided sizes (for constructing with 'fake' - // mesh in GAMG) + scalarField& diag(); + scalarField& upper(); + scalarField& lower(); - scalarField& lower(const label size); - scalarField& diag(const label nCoeffs); - scalarField& upper(const label nCoeffs); + // Size with externally provided sizes + // (for constructing with 'fake' mesh in GAMG) + + scalarField& diag(label size); + scalarField& upper(label nCoeffs); + scalarField& lower(label nCoeffs); - const scalarField& lower() const; - const scalarField& diag() const; - const scalarField& upper() const; + // Characteristics - bool hasDiag() const noexcept - { - return (diagPtr_); - } + //- The matrix type (empty, diagonal, symmetric, ...) + word matrixTypeName() const; - bool hasUpper() const noexcept - { - return (upperPtr_); - } + bool hasDiag() const noexcept { return bool(diagPtr_); } + bool hasUpper() const noexcept { return bool(upperPtr_); } + bool hasLower() const noexcept { return bool(lowerPtr_); } - bool hasLower() const noexcept - { - return (lowerPtr_); - } + //- Matrix has diagonal only + bool diagonal() const noexcept + { + return (diagPtr_ && !lowerPtr_ && !upperPtr_); + } - bool diagonal() const noexcept - { - return (diagPtr_ && !lowerPtr_ && !upperPtr_); - } + //- Matrix is symmetric + bool symmetric() const noexcept + { + return (diagPtr_ && !lowerPtr_ && upperPtr_); + } - bool symmetric() const noexcept - { - return (diagPtr_ && (!lowerPtr_ && upperPtr_)); - } - - bool asymmetric() const noexcept - { - return (diagPtr_ && lowerPtr_ && upperPtr_); - } + //- Matrix is asymmetric (ie, full) + bool asymmetric() const noexcept + { + return (diagPtr_ && lowerPtr_ && upperPtr_); + } // Operations @@ -801,8 +813,12 @@ public: // Member Operators + //- Copy assignment void operator=(const lduMatrix&); + //- Move assignment + void operator=(lduMatrix&&); + void negate(); void operator+=(const lduMatrix&); diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C index 4d307f5710..1e83a03c75 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -86,7 +86,7 @@ void Foam::lduMatrix::sumMagOffDiag } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // void Foam::lduMatrix::operator=(const lduMatrix& A) { @@ -95,38 +95,49 @@ void Foam::lduMatrix::operator=(const lduMatrix& A) return; // Self-assignment is a no-op } - if (A.lowerPtr_) + if (A.hasLower()) { lower() = A.lower(); } - else if (lowerPtr_) + else { - delete lowerPtr_; - lowerPtr_ = nullptr; + lowerPtr_.reset(nullptr); } - if (A.upperPtr_) + if (A.hasUpper()) { upper() = A.upper(); } - else if (upperPtr_) + else { - delete upperPtr_; - upperPtr_ = nullptr; + upperPtr_.reset(nullptr); } - if (A.diagPtr_) + if (A.hasDiag()) { diag() = A.diag(); } } +void Foam::lduMatrix::operator=(lduMatrix&& A) +{ + if (this == &A) + { + return; // Self-assignment is a no-op + } + + diagPtr_ = std::move(A.diagPtr_); + upperPtr_ = std::move(A.upperPtr_); + lowerPtr_ = std::move(A.lowerPtr_); +} + + void Foam::lduMatrix::negate() { - if (lowerPtr_) + if (diagPtr_) { - lowerPtr_->negate(); + diagPtr_->negate(); } if (upperPtr_) @@ -134,16 +145,16 @@ void Foam::lduMatrix::negate() upperPtr_->negate(); } - if (diagPtr_) + if (lowerPtr_) { - diagPtr_->negate(); + lowerPtr_->negate(); } } void Foam::lduMatrix::operator+=(const lduMatrix& A) { - if (A.diagPtr_) + if (A.hasDiag()) { diag() += A.diag(); } @@ -168,7 +179,7 @@ void Foam::lduMatrix::operator+=(const lduMatrix& A) } else if (asymmetric() && A.symmetric()) { - if (A.upperPtr_) + if (A.hasUpper()) { lower() += A.upper(); upper() += A.upper(); @@ -187,12 +198,12 @@ void Foam::lduMatrix::operator+=(const lduMatrix& A) } else if (diagonal()) { - if (A.upperPtr_) + if (A.hasUpper()) { upper() = A.upper(); } - if (A.lowerPtr_) + if (A.hasLower()) { lower() = A.lower(); } @@ -206,15 +217,8 @@ void Foam::lduMatrix::operator+=(const lduMatrix& A) { WarningInFunction << "Unknown matrix type combination" << nl - << " this :" - << " diagonal:" << diagonal() - << " symmetric:" << symmetric() - << " asymmetric:" << asymmetric() << nl - << " A :" - << " diagonal:" << A.diagonal() - << " symmetric:" << A.symmetric() - << " asymmetric:" << A.asymmetric() - << endl; + << " this : " << this->matrixTypeName() + << " A : " << A.matrixTypeName() << endl; } } } @@ -247,7 +251,7 @@ void Foam::lduMatrix::operator-=(const lduMatrix& A) } else if (asymmetric() && A.symmetric()) { - if (A.upperPtr_) + if (A.hasUpper()) { lower() -= A.upper(); upper() -= A.upper(); @@ -266,12 +270,12 @@ void Foam::lduMatrix::operator-=(const lduMatrix& A) } else if (diagonal()) { - if (A.upperPtr_) + if (A.hasUpper()) { upper() = -A.upper(); } - if (A.lowerPtr_) + if (A.hasLower()) { lower() = -A.lower(); } @@ -285,15 +289,8 @@ void Foam::lduMatrix::operator-=(const lduMatrix& A) { WarningInFunction << "Unknown matrix type combination" << nl - << " this :" - << " diagonal:" << diagonal() - << " symmetric:" << symmetric() - << " asymmetric:" << asymmetric() << nl - << " A :" - << " diagonal:" << A.diagonal() - << " symmetric:" << A.symmetric() - << " asymmetric:" << A.asymmetric() - << endl; + << " this : " << this->matrixTypeName() + << " A : " << A.matrixTypeName() << endl; } } } From 27aa7e4e91b1a9fd18f0ea0e5ad2758845cee84d Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Thu, 8 Feb 2024 11:20:08 +0100 Subject: [PATCH 2/2] ENH: improve procLduMatrix streaming and assembly of LUscalarMatrix - reduce local overhead prior to sending for master assembly - non-blocking mode when assembling solution vector - automatic resizing of pivots --- .../matrices/LUscalarMatrix/LUscalarMatrix.C | 168 ++++++++---------- .../matrices/LUscalarMatrix/LUscalarMatrix.H | 28 +-- .../LUscalarMatrix/LUscalarMatrixTemplates.C | 168 ++++++++++++------ .../LUscalarMatrix/procLduInterface.C | 51 ++++-- .../LUscalarMatrix/procLduInterface.H | 14 +- .../matrices/LUscalarMatrix/procLduMatrix.C | 53 ++++-- .../matrices/LUscalarMatrix/procLduMatrix.H | 48 +++-- .../matrices/SquareMatrix/SquareMatrix.C | 6 +- .../matrices/scalarMatrices/scalarMatrices.C | 26 +-- .../matrices/scalarMatrices/scalarMatrices.H | 6 +- .../scalarMatrices/scalarMatricesTemplates.C | 2 +- 11 files changed, 340 insertions(+), 230 deletions(-) diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C index 2d03159abe..1c401c232b 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation - Copyright (C) 2020 OpenCFD Ltd. + Copyright (C) 2020-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -42,17 +42,16 @@ namespace Foam // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::LUscalarMatrix::LUscalarMatrix() +Foam::LUscalarMatrix::LUscalarMatrix() noexcept : comm_(UPstream::worldComm) {} -Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix) +Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& mat) : - scalarSquareMatrix(matrix), - comm_(UPstream::worldComm), - pivotIndices_(m()) + scalarSquareMatrix(mat), + comm_(UPstream::worldComm) { LUDecompose(*this, pivotIndices_); } @@ -69,13 +68,14 @@ Foam::LUscalarMatrix::LUscalarMatrix { if (UPstream::parRun()) { - PtrList lduMatrices(UPstream::nProcs(comm_)); - - label lduMatrixi = 0; + PtrList lduMatrices + ( + UPstream::master(comm_) ? UPstream::nProcs(comm_) : 1 + ); lduMatrices.set ( - lduMatrixi++, + 0, // rank-local matrix (and/or master) new procLduMatrix ( ldum, @@ -88,101 +88,66 @@ Foam::LUscalarMatrix::LUscalarMatrix { for (const int proci : UPstream::subProcs(comm_)) { - lduMatrices.set - ( - lduMatrixi++, - new procLduMatrix - ( - IPstream - ( - UPstream::commsTypes::scheduled, - proci, - 0, // bufSize - UPstream::msgType(), - comm_ - )() - ) - ); + auto& mat = lduMatrices.emplace_set(proci); + + IPstream::recv(mat, proci, UPstream::msgType(), comm_); } + + convert(lduMatrices); } else { - OPstream toMaster + OPstream::send ( - UPstream::commsTypes::scheduled, + lduMatrices[0], // rank-local matrix UPstream::masterNo(), - 0, // bufSize UPstream::msgType(), comm_ ); - procLduMatrix cldum - ( - ldum, - interfaceCoeffs, - interfaces - ); - toMaster<< cldum; - - } - - if (UPstream::master(comm_)) - { - label nCells = 0; - forAll(lduMatrices, i) - { - nCells += lduMatrices[i].size(); - } - - scalarSquareMatrix m(nCells, 0.0); - transfer(m); - convert(lduMatrices); } } else { - label nCells = ldum.lduAddr().size(); - scalarSquareMatrix m(nCells, Zero); - transfer(m); convert(ldum, interfaceCoeffs, interfaces); } + + if (debug && UPstream::master(comm_)) + { + const label numRows = nRows(); + const label numCols = nCols(); + + Pout<< "LUscalarMatrix : size:" << numRows << endl; + for (label rowi = 0; rowi < numRows; ++rowi) + { + const scalar* row = operator[](rowi); + + Pout<< "cell:" << rowi << " diagCoeff:" << row[rowi] << nl; + + Pout<< " connects to upper cells :"; + for (label coli = rowi+1; coli < numCols; ++coli) + { + if (mag(row[coli]) > SMALL) + { + Pout<< ' ' << coli << " (coeff:" << row[coli] << ')'; + } + } + Pout<< nl; + Pout<< " connects to lower cells :"; + for (label coli = 0; coli < rowi; ++coli) + { + if (mag(row[coli]) > SMALL) + { + Pout<< ' ' << coli << " (coeff:" << row[coli] << ')'; + } + } + Pout<< nl; + } + Pout<< endl; + } + if (UPstream::master(comm_)) { - if (debug) - { - const label numRows = m(); - const label numCols = n(); - - Pout<< "LUscalarMatrix : size:" << numRows << endl; - for (label rowi = 0; rowi < numRows; ++rowi) - { - const scalar* row = operator[](rowi); - - Pout<< "cell:" << rowi << " diagCoeff:" << row[rowi] << endl; - - Pout<< " connects to upper cells :"; - for (label coli = rowi+1; coli < numCols; ++coli) - { - if (mag(row[coli]) > SMALL) - { - Pout<< ' ' << coli << " (coeff:" << row[coli] << ')'; - } - } - Pout<< endl; - Pout<< " connects to lower cells :"; - for (label coli = 0; coli < rowi; ++coli) - { - if (mag(row[coli]) > SMALL) - { - Pout<< ' ' << coli << " (coeff:" << row[coli] << ')'; - } - } - Pout<< nl; - } - Pout<< nl; - } - - pivotIndices_.setSize(m()); LUDecompose(*this, pivotIndices_); } } @@ -197,6 +162,10 @@ void Foam::LUscalarMatrix::convert const lduInterfaceFieldPtrsList& interfaces ) { + // Resize and fill with zero + scalarSquareMatrix::resize_nocopy(ldum.lduAddr().size()); + scalarSquareMatrix::operator=(Foam::zero{}); + const label* __restrict__ uPtr = ldum.lduAddr().upperAddr().begin(); const label* __restrict__ lPtr = ldum.lduAddr().lowerAddr().begin(); @@ -259,14 +228,26 @@ void Foam::LUscalarMatrix::convert const PtrList& lduMatrices ) { - procOffsets_.setSize(lduMatrices.size() + 1); - procOffsets_[0] = 0; + procOffsets_.resize_nocopy(lduMatrices.size() + 1); - forAll(lduMatrices, ldumi) { - procOffsets_[ldumi+1] = procOffsets_[ldumi] + lduMatrices[ldumi].size(); + auto iter = procOffsets_.begin(); + + label nCellsTotal = 0; + *iter++ = nCellsTotal; + + for (const auto& mat : lduMatrices) + { + nCellsTotal += mat.size(); + *iter++ = nCellsTotal; + } + + // Resize and fill with zero + scalarSquareMatrix::resize_nocopy(nCellsTotal); + scalarSquareMatrix::operator=(Foam::zero{}); } + forAll(lduMatrices, ldumi) { const procLduMatrix& lduMatrixi = lduMatrices[ldumi]; @@ -400,10 +381,9 @@ void Foam::LUscalarMatrix::printDiagonalDominance() const } -void Foam::LUscalarMatrix::decompose(const scalarSquareMatrix& M) +void Foam::LUscalarMatrix::decompose(const scalarSquareMatrix& mat) { - scalarSquareMatrix::operator=(M); - pivotIndices_.setSize(m()); + scalarSquareMatrix::operator=(mat); LUDecompose(*this, pivotIndices_); } diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H index ada33620be..682ec7e6cd 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H @@ -34,8 +34,8 @@ SourceFiles \*---------------------------------------------------------------------------*/ -#ifndef LUscalarMatrix_H -#define LUscalarMatrix_H +#ifndef Foam_LUscalarMatrix_H +#define Foam_LUscalarMatrix_H #include "scalarMatrices.H" #include "labelList.H" @@ -47,6 +47,7 @@ SourceFiles namespace Foam { +// Forward Declarations class lduMatrix; class procLduMatrix; @@ -58,7 +59,7 @@ class LUscalarMatrix : public scalarSquareMatrix { - // Private data + // Private Data //- Communicator to use const label comm_; @@ -70,7 +71,7 @@ class LUscalarMatrix labelList pivotIndices_; - // Private member functions + // Private Member Functions //- Convert the given lduMatrix into this LUscalarMatrix void convert @@ -81,12 +82,12 @@ class LUscalarMatrix ); //- Convert the given list of procLduMatrix into this LUscalarMatrix - // on the master processor + //- on the master processor void convert(const PtrList& lduMatrices); //- Print the ratio of the mag-sum of the off-diagonal coefficients - // to the mag-diagonal + //- to the mag-diagonal void printDiagonalDominance() const; @@ -98,13 +99,14 @@ public: // Constructors - //- Construct null - LUscalarMatrix(); + //- Default construct + LUscalarMatrix() noexcept; - //- Construct from and perform LU decomposition of the matrix M - LUscalarMatrix(const scalarSquareMatrix& M); + //- Construct from and perform LU decomposition of the given matrix + explicit LUscalarMatrix(const scalarSquareMatrix& mat); - //- Construct from lduMatrix and perform LU decomposition + //- Construct from lduMatrix and perform LU decomposition. + //- In parallel it assembles the matrix on the master. LUscalarMatrix ( const lduMatrix& ldum, @@ -115,8 +117,8 @@ public: // Member Functions - //- Perform the LU decomposition of the matrix M - void decompose(const scalarSquareMatrix& M); + //- Perform the LU decomposition of the matrix + void decompose(const scalarSquareMatrix& mat); //- Solve the linear system with the given source // and returning the solution in the Field argument x. diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C index 568275e54c..9113a8a124 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation - Copyright (C) 2019-2020 OpenCFD Ltd. + Copyright (C) 2019-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -44,79 +44,141 @@ void Foam::LUscalarMatrix::solve x = source; } - if (Pstream::parRun()) + const auto tag = UPstream::msgType(); + + if (UPstream::parRun()) { - List X; // scratch space (on master) + List allx; // scratch space (on master) - if (Pstream::master(comm_)) + const label startOfRequests = UPstream::nRequests(); + + // Like globalIndex::gather() + if (UPstream::master(comm_)) { - X.resize(m()); + allx.resize(m()); - SubList(X, x.size()) = x; + SubList(allx, x.size()) = x; - for (const int proci : Pstream::subProcs(comm_)) + for (const int proci : UPstream::subProcs(comm_)) { - UIPstream::read + SubList procSlot ( - Pstream::commsTypes::scheduled, - proci, - reinterpret_cast - ( - &(X[procOffsets_[proci]]) - ), - (procOffsets_[proci+1]-procOffsets_[proci])*sizeof(Type), - Pstream::msgType(), - comm_ + allx, + procOffsets_[proci+1]-procOffsets_[proci], + procOffsets_[proci] ); + + if (procSlot.empty()) + { + // Nothing to do + } + else if (is_contiguous::value) + { + UIPstream::read + ( + UPstream::commsTypes::nonBlocking, + proci, + procSlot.data_bytes(), + procSlot.size_bytes(), + tag, + comm_ + ); + } + else + { + IPstream::recv(procSlot, proci, tag, comm_); + } } } else { - UOPstream::write - ( - Pstream::commsTypes::scheduled, - Pstream::masterNo(), - x.cdata_bytes(), - x.byteSize(), - Pstream::msgType(), - comm_ - ); - } - - if (Pstream::master(comm_)) - { - LUBacksubstitute(*this, pivotIndices_, X); - - x = SubList(X, x.size()); - - for (const int proci : Pstream::subProcs(comm_)) + if (x.empty()) + { + // Nothing to do + } + else if (is_contiguous::value) { UOPstream::write ( - Pstream::commsTypes::scheduled, - proci, - reinterpret_cast - ( - &(X[procOffsets_[proci]]) - ), - (procOffsets_[proci+1]-procOffsets_[proci])*sizeof(Type), - Pstream::msgType(), + UPstream::commsTypes::nonBlocking, + UPstream::masterNo(), + x.cdata_bytes(), + x.size_bytes(), + tag, comm_ ); } + else + { + OPstream::send(x, UPstream::masterNo(), tag, comm_); + } + } + + UPstream::waitRequests(startOfRequests); + + // LUBacksubstitute and then like globalIndex::scatter() + if (UPstream::master(comm_)) + { + LUBacksubstitute(*this, pivotIndices_, allx); + + x = SubList(allx, x.size()); + + for (const int proci : UPstream::subProcs(comm_)) + { + SubList procSlot + ( + allx, + procOffsets_[proci+1]-procOffsets_[proci], + procOffsets_[proci] + ); + + if (procSlot.empty()) + { + // Nothing to do + } + else if (is_contiguous::value) + { + UOPstream::write + ( + UPstream::commsTypes::nonBlocking, + proci, + procSlot.cdata_bytes(), + procSlot.size_bytes(), + tag, + comm_ + ); + } + else + { + OPstream::send(procSlot, proci, tag, comm_); + } + } } else { - UIPstream::read - ( - Pstream::commsTypes::scheduled, - Pstream::masterNo(), - x.data_bytes(), - x.byteSize(), - Pstream::msgType(), - comm_ - ); + if (x.empty()) + { + // Nothing to do + } + else if (is_contiguous::value) + { + UIPstream::read + ( + UPstream::commsTypes::nonBlocking, + UPstream::masterNo(), + x.data_bytes(), + x.size_bytes(), + tag, + comm_ + ); + } + else + { + IPstream::recv(x, UPstream::masterNo(), tag, comm_); + } } + + UPstream::waitRequests(startOfRequests); } else { @@ -131,7 +193,7 @@ Foam::tmp> Foam::LUscalarMatrix::solve const UList& source ) const { - auto tx(tmp>::New(m())); + auto tx = tmp>::New(m()); solve(tx.ref(), source); diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.C b/src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.C index 02309fbab5..6680dca7f1 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.C @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2015 OpenFOAM Foundation + Copyright (C) 2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -69,27 +70,47 @@ Foam::procLduInterface::procLduInterface Foam::procLduInterface::procLduInterface(Istream& is) -: - faceCells_(is), - coeffs_(is), - myProcNo_(readLabel(is)), - neighbProcNo_(readLabel(is)), - tag_(readLabel(is)), - comm_(readLabel(is)) -{} +{ + read(is); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::procLduInterface::read(Istream& is) +{ + is >> faceCells_ + >> coeffs_ + >> myProcNo_ + >> neighbProcNo_ + >> tag_ + >> comm_; +} + + +void Foam::procLduInterface::write(Ostream& os) const +{ + os << faceCells_ + << coeffs_ + << myProcNo_ + << neighbProcNo_ + << tag_ + << comm_; +} // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // -Foam::Ostream& Foam::operator<<(Ostream& os, const procLduInterface& cldui) +Foam::Istream& Foam::operator>>(Istream& is, procLduInterface& intf) { - os << cldui.faceCells_ - << cldui.coeffs_ - << cldui.myProcNo_ - << cldui.neighbProcNo_ - << cldui.tag_ - << cldui.comm_; + intf.read(is); + return is; +} + +Foam::Ostream& Foam::operator<<(Ostream& os, const procLduInterface& intf) +{ + intf.write(os); return os; } diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.H b/src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.H index c0e74aeb97..22f1582bc2 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.H +++ b/src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.H @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation + Copyright (C) 2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -49,10 +50,11 @@ namespace Foam class lduInterfaceField; class procLduInterface; +Istream& operator>>(Istream&, procLduInterface&); Ostream& operator<<(Ostream&, const procLduInterface&); /*---------------------------------------------------------------------------*\ - Class procLduInterface Declaration + Class procLduInterface Declaration \*---------------------------------------------------------------------------*/ class procLduInterface @@ -68,6 +70,7 @@ class procLduInterface public: + //- Friendship friend class LUscalarMatrix; @@ -85,7 +88,8 @@ public: const scalarField& coeffs ); - procLduInterface(Istream& is); + //- Read construct from Istream + explicit procLduInterface(Istream& is); autoPtr clone() { @@ -99,8 +103,12 @@ public: } - // Ostream Operator + // IO Operations + void read(Istream& is); + void write(Ostream& os) const; + + friend Istream& operator>>(Istream&, procLduInterface&); friend Ostream& operator<<(Ostream&, const procLduInterface&); }; diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/procLduMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/procLduMatrix.C index 9670e99e14..ff2e032d97 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/procLduMatrix.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/procLduMatrix.C @@ -50,7 +50,7 @@ Foam::procLduMatrix::procLduMatrix forAll(interfaces, i) { - if (interfaces.set(i)) + if (interfaces.test(i)) { interfaces_.set ( @@ -67,27 +67,48 @@ Foam::procLduMatrix::procLduMatrix Foam::procLduMatrix::procLduMatrix(Istream& is) -: - upperAddr_(is), - lowerAddr_(is), - diag_(is), - upper_(is), - lower_(is), - interfaces_(is) -{} +{ + read(is); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +// void Foam::procLduMatrix::clear() +// { +// upperAddr_.clear(); +// lowerAddr_.clear(); +// diag_.clear(); +// upper_.clear(); +// lower_.clear(); +// interfaces_.clear(); +// } + + +void Foam::procLduMatrix::read(Istream& is) +{ + is >> upperAddr_ >> lowerAddr_ >> diag_ >> upper_ >> lower_ >> interfaces_; +} + + +void Foam::procLduMatrix::write(Ostream& os) const +{ + os << upperAddr_ << lowerAddr_ << diag_ << upper_ << lower_ << interfaces_; +} // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // -Foam::Ostream& Foam::operator<<(Ostream& os, const procLduMatrix& cldum) +Foam::Istream& Foam::operator>>(Istream& is, procLduMatrix& mat) { - os << cldum.upperAddr_ - << cldum.lowerAddr_ - << cldum.diag_ - << cldum.upper_ - << cldum.lower_ - << cldum.interfaces_; + mat.read(is); + return is; +} + +Foam::Ostream& Foam::operator<<(Ostream& os, const procLduMatrix& mat) +{ + mat.write(os); return os; } diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/procLduMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/procLduMatrix.H index 179b6a8860..71c63e4ef9 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/procLduMatrix.H +++ b/src/OpenFOAM/matrices/LUscalarMatrix/procLduMatrix.H @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2013 OpenFOAM Foundation + Copyright (C) 2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -34,8 +35,8 @@ SourceFiles \*---------------------------------------------------------------------------*/ -#ifndef procLduMatrix_H -#define procLduMatrix_H +#ifndef Foam_procLduMatrix_H +#define Foam_procLduMatrix_H #include "labelList.H" #include "scalarField.H" @@ -47,23 +48,21 @@ SourceFiles namespace Foam { +// Forward Declarations +class procLduMatrix; class procLduInterface; class lduMatrix; -// Forward declaration of friend functions and operators - -class procLduMatrix; - +Istream& operator>>(Istream&, procLduMatrix&); Ostream& operator<<(Ostream&, const procLduMatrix&); - /*---------------------------------------------------------------------------*\ - Class procLduMatrix Declaration + Class procLduMatrix Declaration \*---------------------------------------------------------------------------*/ class procLduMatrix { - // Private data + // Private Data labelList upperAddr_; labelList lowerAddr_; @@ -72,20 +71,24 @@ class procLduMatrix scalarField lower_; PtrList interfaces_; +public: - // Private Member Functions + //- Friendship + friend class LUscalarMatrix; + + + // Generated Methods + + //- Default construct + procLduMatrix() = default; //- No copy construct procLduMatrix(const procLduMatrix&) = delete; -public: - - friend class LUscalarMatrix; - - // Constructors + //- Construct from components. Extracts active interfaces procLduMatrix ( const lduMatrix& ldum, @@ -93,19 +96,26 @@ public: const lduInterfaceFieldPtrsList& interfaces ); - procLduMatrix(Istream& is); + //- Read construct from Istream + explicit procLduMatrix(Istream& is); - // Member functions + // Member Functions - label size() const + label size() const noexcept { return diag_.size(); } + // void clear(); - // Ostream operator + // IO Operations + + void read(Istream& is); + void write(Ostream& os) const; + + friend Istream& operator>>(Istream&, procLduMatrix&); friend Ostream& operator<<(Ostream&, const procLduMatrix&); }; diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C index 2c64598fd6..fc5e124184 100644 --- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2013-2016 OpenFOAM Foundation - Copyright (C) 2019-2023 OpenCFD Ltd. + Copyright (C) 2019-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -124,7 +124,7 @@ scalar det(const SquareMatrix& matrix) { SquareMatrix matrixTmp = matrix; - labelList pivotIndices(matrix.m()); + labelList pivotIndices; label sign; LUDecompose(matrixTmp, pivotIndices, sign); @@ -136,7 +136,7 @@ scalar det(const SquareMatrix& matrix) template scalar det(SquareMatrix& matrix) { - labelList pivotIndices(matrix.m()); + labelList pivotIndices; label sign; LUDecompose(matrix, pivotIndices, sign); diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C index 39d3cc37b5..1bd42a3608 100644 --- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation + Copyright (C) 2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -48,17 +49,20 @@ void Foam::LUDecompose label& sign ) { - label m = matrix.m(); - scalar vv[m]; + const label size = matrix.m(); + + pivotIndices.resize_nocopy(size); + + List vv(size); sign = 1; - for (label i = 0; i < m; ++i) + for (label i = 0; i < size; ++i) { scalar largestCoeff = 0.0; scalar temp; const scalar* __restrict__ matrixi = matrix[i]; - for (label j = 0; j < m; ++j) + for (label j = 0; j < size; ++j) { if ((temp = mag(matrixi[j])) > largestCoeff) { @@ -75,7 +79,7 @@ void Foam::LUDecompose vv[i] = 1.0/largestCoeff; } - for (label j = 0; j < m; ++j) + for (label j = 0; j < size; ++j) { scalar* __restrict__ matrixj = matrix[j]; @@ -94,7 +98,7 @@ void Foam::LUDecompose label iMax = 0; scalar largestCoeff = 0.0; - for (label i = j; i < m; ++i) + for (label i = j; i < size; ++i) { scalar* __restrict__ matrixi = matrix[i]; scalar sum = matrixi[j]; @@ -120,7 +124,7 @@ void Foam::LUDecompose { scalar* __restrict__ matrixiMax = matrix[iMax]; - for (label k = 0; k < m; ++k) + for (label k = 0; k < size; ++k) { std::swap(matrixj[k], matrixiMax[k]); } @@ -134,11 +138,11 @@ void Foam::LUDecompose matrixj[j] = SMALL; } - if (j != m-1) + if (j != size-1) { scalar rDiag = 1.0/matrixj[j]; - for (label i = j + 1; i < m; ++i) + for (label i = j + 1; i < size; ++i) { matrix(i, j) *= rDiag; } @@ -150,7 +154,7 @@ void Foam::LUDecompose void Foam::LUDecompose(scalarSymmetricSquareMatrix& matrix) { // Store result in upper triangular part of matrix - label size = matrix.m(); + const label size = matrix.m(); // Set upper triangular parts to zero. for (label j = 0; j < size; ++j) @@ -223,7 +227,7 @@ void Foam::multiply << abort(FatalError); } - ans = scalarRectangularMatrix(A.m(), C.n(), Zero); + ans = scalarRectangularMatrix(A.m(), C.n(), Foam::zero{}); for (label i = 0; i < A.m(); ++i) { diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H index d6ba156a2b..8035cdfebf 100644 --- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H @@ -74,19 +74,21 @@ void solve const List& source ); -//- LU decompose the matrix with pivoting +//- LU decompose the matrix with pivoting. void LUDecompose ( scalarSquareMatrix& matrix, + //! [out] size is adjusted as required labelList& pivotIndices ); //- LU decompose the matrix with pivoting. -//- sign is -1 for odd number of row interchanges and 1 for even number. void LUDecompose ( scalarSquareMatrix& matrix, + //! [out] size is adjusted as required labelList& pivotIndices, + //! [out] is -1 for odd number of row interchanges and 1 for even number label& sign ); diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C index 0c3ae68228..7a66cbe6d3 100644 --- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C @@ -216,7 +216,7 @@ void Foam::LUsolve List& sourceSol ) { - labelList pivotIndices(matrix.m()); + labelList pivotIndices; LUDecompose(matrix, pivotIndices); LUBacksubstitute(matrix, pivotIndices, sourceSol); }