diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
index 2d03159abe84e52446523856553a3511829a29c6..1c401c232b7aa646b212c272dbd3bf54e49f9e66 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<procLduMatrix> lduMatrices(UPstream::nProcs(comm_));
-
- label lduMatrixi = 0;
+ PtrList<procLduMatrix> 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 (UPstream::master(comm_))
+
+ if (debug && UPstream::master(comm_))
{
- if (debug)
- {
- const label numRows = m();
- const label numCols = n();
+ 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<< "LUscalarMatrix : size:" << numRows << endl;
+ for (label rowi = 0; rowi < numRows; ++rowi)
+ {
+ const scalar* row = operator[](rowi);
- Pout<< "cell:" << rowi << " diagCoeff:" << row[rowi] << endl;
+ Pout<< "cell:" << rowi << " diagCoeff:" << row[rowi] << nl;
- Pout<< " connects to upper cells :";
- for (label coli = rowi+1; coli < numCols; ++coli)
+ Pout<< " connects to upper cells :";
+ for (label coli = rowi+1; coli < numCols; ++coli)
+ {
+ if (mag(row[coli]) > SMALL)
{
- if (mag(row[coli]) > SMALL)
- {
- Pout<< ' ' << coli << " (coeff:" << row[coli] << ')';
- }
+ Pout<< ' ' << coli << " (coeff:" << row[coli] << ')';
}
- Pout<< endl;
- Pout<< " connects to lower cells :";
- for (label coli = 0; coli < rowi; ++coli)
+ }
+ Pout<< nl;
+ Pout<< " connects to lower cells :";
+ for (label coli = 0; coli < rowi; ++coli)
+ {
+ if (mag(row[coli]) > SMALL)
{
- if (mag(row[coli]) > SMALL)
- {
- Pout<< ' ' << coli << " (coeff:" << row[coli] << ')';
- }
+ Pout<< ' ' << coli << " (coeff:" << row[coli] << ')';
}
- Pout<< nl;
}
Pout<< nl;
}
+ Pout<< endl;
+ }
- pivotIndices_.setSize(m());
+ if (UPstream::master(comm_))
+ {
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<procLduMatrix>& 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 ada33620bee9e3da7fc7bf0c9308dec474e25c2f..682ec7e6cdfe4869867bd373264f236144c28f8a 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<procLduMatrix>& 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 568275e54c8ae398a3ef56dddff9ed1a336fbb75..9113a8a1243255d3822b52c091c8f71f0ea70b8c 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<Type> X; // scratch space (on master)
+ List<Type> allx; // scratch space (on master)
+
+ const label startOfRequests = UPstream::nRequests();
- if (Pstream::master(comm_))
+ // Like globalIndex::gather()
+ if (UPstream::master(comm_))
{
- X.resize(m());
+ allx.resize(m());
- SubList<Type>(X, x.size()) = x;
+ SubList<Type>(allx, x.size()) = x;
- for (const int proci : Pstream::subProcs(comm_))
+ for (const int proci : UPstream::subProcs(comm_))
{
- UIPstream::read
+ SubList<Type> procSlot
(
- Pstream::commsTypes::scheduled,
- proci,
- reinterpret_cast<char*>
- (
- &(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<Type>::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 (x.empty())
+ {
+ // Nothing to do
+ }
+ else if (is_contiguous<Type>::value)
+ {
+ UOPstream::write
+ (
+ UPstream::commsTypes::nonBlocking,
+ UPstream::masterNo(),
+ x.cdata_bytes(),
+ x.size_bytes(),
+ tag,
+ comm_
+ );
+ }
+ else
+ {
+ OPstream::send(x, UPstream::masterNo(), tag, comm_);
+ }
}
- if (Pstream::master(comm_))
+ UPstream::waitRequests(startOfRequests);
+
+ // LUBacksubstitute and then like globalIndex::scatter()
+ if (UPstream::master(comm_))
{
- LUBacksubstitute(*this, pivotIndices_, X);
+ LUBacksubstitute(*this, pivotIndices_, allx);
- x = SubList<Type>(X, x.size());
+ x = SubList<Type>(allx, x.size());
- for (const int proci : Pstream::subProcs(comm_))
+ for (const int proci : UPstream::subProcs(comm_))
{
- UOPstream::write
+ SubList<Type> procSlot
(
- Pstream::commsTypes::scheduled,
- proci,
- reinterpret_cast<const char*>
- (
- &(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<Type>::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<Type>::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::Field<Type>> Foam::LUscalarMatrix::solve
const UList<Type>& source
) const
{
- auto tx(tmp<Field<Type>>::New(m()));
+ auto tx = tmp<Field<Type>>::New(m());
solve(tx.ref(), source);
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.C b/src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.C
index 02309fbab5851ba5091b9c1d7a8d0ed7ca2f5532..6680dca7f1d0a54cc200b6a2b3cd77e2fc5dd3b3 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 c0e74aeb978aea8361822a01ccaba75b77ab9203..22f1582bc21e9e09ca128ea5c8ab7072ae1b5547 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<procLduInterface> 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 9670e99e14657811473c9c9de61a37b2d2880458..ff2e032d97dec24bc9d7854756646e3c9b7e5075 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 179b6a88602e60b443fc600ec7b48dd2ee1a0831..71c63e4ef97df153b3e991f4157401dcc128a771 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<procLduInterface> interfaces_;
+public:
- // Private Member Functions
+ //- Friendship
+ friend class LUscalarMatrix;
- //- No copy construct
- procLduMatrix(const procLduMatrix&) = delete;
+ // Generated Methods
-public:
+ //- Default construct
+ procLduMatrix() = default;
- friend class LUscalarMatrix;
+ //- No copy construct
+ procLduMatrix(const procLduMatrix&) = delete;
// 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();
+
+
+ // IO Operations
- // Ostream operator
+ 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/LduMatrix/LduMatrix/LduMatrix.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C
index dbdaa0fbc3f050bcc54a3eae48887d773536225a..c76026bd5353c776dccd68516803b636184e37e8 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<class Type, class DType, class LUType>
Foam::LduMatrix<Type, DType, LUType>::LduMatrix(const lduMesh& mesh)
:
- lduMesh_(mesh),
- diagPtr_(nullptr),
- upperPtr_(nullptr),
- lowerPtr_(nullptr),
- sourcePtr_(nullptr),
- interfaces_(0),
- interfacesUpper_(0),
- interfacesLower_(0)
+ lduMesh_(mesh)
{}
template<class Type, class DType, class LUType>
Foam::LduMatrix<Type, DType, LUType>::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<DType>(*(A.diagPtr_));
+ diagPtr_ = std::make_unique<Field<DType>>(*(A.diagPtr_));
}
if (A.upperPtr_)
{
- upperPtr_ = new Field<LUType>(*(A.upperPtr_));
+ upperPtr_ = std::make_unique<Field<LUType>>(*(A.upperPtr_));
}
if (A.lowerPtr_)
{
- lowerPtr_ = new Field<LUType>(*(A.lowerPtr_));
+ lowerPtr_ = std::make_unique<Field<LUType>>(*(A.lowerPtr_));
}
if (A.sourcePtr_)
{
- sourcePtr_ = new Field<Type>(*(A.sourcePtr_));
+ sourcePtr_ = std::make_unique<Field<Type>>(*(A.sourcePtr_));
}
}
template<class Type, class DType, class LUType>
-Foam::LduMatrix<Type, DType, LUType>::LduMatrix(LduMatrix& A, bool reuse)
+Foam::LduMatrix<Type, DType, LUType>::LduMatrix(LduMatrix&& A)
:
lduMesh_(A.lduMesh_),
- diagPtr_(nullptr),
- upperPtr_(nullptr),
- lowerPtr_(nullptr),
- sourcePtr_(nullptr),
- interfaces_(0),
- interfacesUpper_(0),
- interfacesLower_(0)
+ diagPtr_(std::move(A.diagPtr_)),
+ lowerPtr_(std::move(A.lowerPtr_)),
+ upperPtr_(std::move(A.upperPtr_)),
+ sourcePtr_(std::move(A.sourcePtr_))
{
- if (reuse)
- {
- if (A.diagPtr_)
- {
- diagPtr_ = A.diagPtr_;
- A.diagPtr_ = nullptr;
- }
+ // Clear the old interfaces?
+}
- if (A.upperPtr_)
- {
- upperPtr_ = A.upperPtr_;
- A.upperPtr_ = nullptr;
- }
- if (A.lowerPtr_)
- {
- lowerPtr_ = A.lowerPtr_;
- A.lowerPtr_ = nullptr;
- }
+template<class Type, class DType, class LUType>
+Foam::LduMatrix<Type, DType, LUType>::LduMatrix(LduMatrix& A, bool reuse)
+:
+ lduMesh_(A.lduMesh_)
+{
+ if (reuse)
+ {
+ // Move assignment
+ diagPtr_ = std::move(A.diagPtr_);
+ upperPtr_ = std::move(A.upperPtr_);
+ lowerPtr_ = std::move(A.lowerPtr_);
+ sourcePtr_ = std::move(A.sourcePtr_);
- if (A.sourcePtr_)
- {
- sourcePtr_ = A.sourcePtr_;
- A.sourcePtr_ = nullptr;
- }
+ // Clear the old interfaces?
}
else
{
+ // Copy assignment
if (A.diagPtr_)
{
- diagPtr_ = new Field<DType>(*(A.diagPtr_));
+ diagPtr_ = std::make_unique<Field<DType>>(*(A.diagPtr_));
}
if (A.upperPtr_)
{
- upperPtr_ = new Field<LUType>(*(A.upperPtr_));
+ upperPtr_ = std::make_unique<Field<LUType>>(*(A.upperPtr_));
}
if (A.lowerPtr_)
{
- lowerPtr_ = new Field<LUType>(*(A.lowerPtr_));
+ lowerPtr_ = std::make_unique<Field<LUType>>(*(A.lowerPtr_));
}
if (A.sourcePtr_)
{
- sourcePtr_ = new Field<Type>(*(A.sourcePtr_));
+ sourcePtr_ = std::make_unique<Field<Type>>(*(A.sourcePtr_));
}
}
}
@@ -152,48 +130,51 @@ Foam::LduMatrix<Type, DType, LUType>::LduMatrix
diagPtr_(new Field<DType>(is)),
upperPtr_(new Field<LUType>(is)),
lowerPtr_(new Field<LUType>(is)),
- sourcePtr_(new Field<Type>(is)),
- interfaces_(0),
- interfacesUpper_(0),
- interfacesLower_(0)
+ sourcePtr_(new Field<Type>(is))
{}
-// * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
-Foam::LduMatrix<Type, DType, LUType>::~LduMatrix()
+Foam::word Foam::LduMatrix<Type, DType, LUType>::matrixTypeName() const
{
if (diagPtr_)
{
- delete diagPtr_;
+ return
+ (
+ (!upperPtr_)
+ ? (!lowerPtr_ ? "diagonal" : "diagonal-lower")
+ : (!lowerPtr_ ? "symmetric" : "asymmetric")
+ );
}
- if (upperPtr_)
- {
- delete upperPtr_;
- }
+ // is empty (or just wrong)
+ return (!upperPtr_ && !lowerPtr_ ? "empty" : "ill-defined");
+}
- if (lowerPtr_)
- {
- delete lowerPtr_;
- }
- if (sourcePtr_)
+template<class Type, class DType, class LUType>
+const Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag() const
+{
+ if (!diagPtr_)
{
- delete sourcePtr_;
+ FatalErrorInFunction
+ << "diagPtr_ unallocated"
+ << abort(FatalError);
}
-}
+ return *diagPtr_;
+}
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag()
{
if (!diagPtr_)
{
- diagPtr_ = new Field<DType>(lduAddr().size(), Zero);
+ diagPtr_ =
+ std::make_unique<Field<DType>>(lduAddr().size(), Foam::zero{});
}
return *diagPtr_;
@@ -201,127 +182,115 @@ Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag()
template<class Type, class DType, class LUType>
-Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper()
+const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper() const
{
- if (!upperPtr_)
+ if (upperPtr_)
{
- if (lowerPtr_)
- {
- upperPtr_ = new Field<LUType>(*lowerPtr_);
- }
- else
+ return *upperPtr_;
+ }
+ else
+ {
+ if (!lowerPtr_)
{
- upperPtr_ = new Field<LUType>
- (
- lduAddr().lowerAddr().size(),
- Zero
- );
+ FatalErrorInFunction
+ << "lowerPtr_ and upperPtr_ unallocated"
+ << abort(FatalError);
}
- }
- return *upperPtr_;
+ return *lowerPtr_;
+ }
}
template<class Type, class DType, class LUType>
-Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower()
+Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper()
{
- if (!lowerPtr_)
+ if (!upperPtr_)
{
- if (upperPtr_)
+ if (lowerPtr_)
{
- lowerPtr_ = new Field<LUType>(*upperPtr_);
+ upperPtr_ = std::make_unique<Field<LUType>>(*lowerPtr_);
}
else
{
- lowerPtr_ = new Field<LUType>
- (
- lduAddr().lowerAddr().size(),
- Zero
- );
+ upperPtr_ =
+ std::make_unique<Field<LUType>>
+ (
+ lduAddr().lowerAddr().size(),
+ Foam::zero{}
+ );
}
}
- return *lowerPtr_;
+ return *upperPtr_;
}
template<class Type, class DType, class LUType>
-Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source()
+const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower() const
{
- if (!sourcePtr_)
+ if (lowerPtr_)
{
- sourcePtr_ = new Field<Type>(lduAddr().size(), Zero);
+ return *lowerPtr_;
}
-
- return *sourcePtr_;
-}
-
-
-template<class Type, class DType, class LUType>
-const Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag() const
-{
- if (!diagPtr_)
+ else
{
- FatalErrorInFunction
- << "diagPtr_ unallocated"
- << abort(FatalError);
- }
+ if (!upperPtr_)
+ {
+ FatalErrorInFunction
+ << "lowerPtr_ and upperPtr_ unallocated"
+ << abort(FatalError);
+ }
- return *diagPtr_;
+ return *upperPtr_;
+ }
}
template<class Type, class DType, class LUType>
-const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper() const
+Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower()
{
- if (!lowerPtr_ && !upperPtr_)
+ if (!lowerPtr_)
{
- FatalErrorInFunction
- << "lowerPtr_ or upperPtr_ unallocated"
- << abort(FatalError);
+ if (upperPtr_)
+ {
+ lowerPtr_ = std::make_unique<Field<LUType>>(*upperPtr_);
+ }
+ else
+ {
+ lowerPtr_ = std::make_unique<Field<LUType>>
+ (
+ lduAddr().lowerAddr().size(),
+ Foam::zero{}
+ );
+ }
}
- if (upperPtr_)
- {
- return *upperPtr_;
- }
- else
- {
- return *lowerPtr_;
- }
+ return *lowerPtr_;
}
template<class Type, class DType, class LUType>
-const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower() const
+const Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source() const
{
- if (!lowerPtr_ && !upperPtr_)
+ if (!sourcePtr_)
{
FatalErrorInFunction
- << "lowerPtr_ or upperPtr_ unallocated"
+ << "sourcePtr_ unallocated"
<< abort(FatalError);
}
- if (lowerPtr_)
- {
- return *lowerPtr_;
- }
- else
- {
- return *upperPtr_;
- }
+ return *sourcePtr_;
}
template<class Type, class DType, class LUType>
-const Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source() const
+Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source()
{
if (!sourcePtr_)
{
- FatalErrorInFunction
- << "sourcePtr_ unallocated"
- << abort(FatalError);
+ sourcePtr_ =
+ std::make_unique<Field<Type>>(lduAddr().size(), Foam::zero{});
}
return *sourcePtr_;
@@ -330,43 +299,50 @@ const Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source() const
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
+// template<class Type, class DType, class LUType>
+// Foam::Ostream& Foam::operator<<
+// (
+// Ostream& os,
+// const InfoProxy<Type, DType, LUType>& iproxy
+// )
+// {
+// const auto& mat = *iproxy;
+//
+// ...
+//
+// os.check(FUNCTION_NAME);
+// return os;
+// }
+
+
template<class Type, class DType, class LUType>
Foam::Ostream& Foam::operator<<
(
Ostream& os,
- const LduMatrix<Type, DType, LUType>& ldum
+ const LduMatrix<Type, DType, LUType>& 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 96e9fd14f4e303bf5ca51a51a0f739db6ed20f95..407f7ee3226666147dcf8029533ce8f5f04d988b 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<DType> *diagPtr_;
+ std::unique_ptr<Field<DType>> diagPtr_;
//- Off-diagonal coefficients
- Field<LUType> *upperPtr_, *lowerPtr_;
+ std::unique_ptr<Field<LUType>> upperPtr_;
+
+ //- Off-diagonal coefficients
+ std::unique_ptr<Field<LUType>> lowerPtr_;
//- Source
- Field<Type> *sourcePtr_;
+ std::unique_ptr<Field<Type>> sourcePtr_;
//- Field interfaces (processor patches etc.)
LduInterfaceFieldPtrsList<Type> interfaces_;
//- Off-diagonal coefficients for interfaces
- FieldField<Field, LUType> interfacesUpper_, interfacesLower_;
+ FieldField<Field, LUType> interfacesUpper_;
+
+ //- Off-diagonal coefficients for interfaces
+ FieldField<Field, LUType> interfacesLower_;
public:
friend class SolverPerformance<Type>;
+ // -----------------------------------------------------------------------
//- 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<Type, DType, LUType>&);
+ //- Move construct
+ LduMatrix(LduMatrix<Type, DType, LUType>&&);
+
//- Construct as copy or re-use as specified.
LduMatrix(LduMatrix<Type, DType, LUType>&, 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<Type>& interfaces() const
- {
- return interfaces_;
- }
+ //- Const access to the interfaces
+ const LduInterfaceFieldPtrsList<Type>& interfaces() const noexcept
+ {
+ return interfaces_;
+ }
- //- Return interfaces
- LduInterfaceFieldPtrsList<Type>& interfaces()
- {
- return interfaces_;
- }
+ //- Non-const access to the interfaces
+ LduInterfaceFieldPtrsList<Type>& interfaces() noexcept
+ {
+ return interfaces_;
+ }
- // Access to coefficients
+ // Coefficients
- Field<DType>& diag();
- Field<LUType>& upper();
- Field<LUType>& lower();
- Field<Type>& source();
+ const Field<DType>& diag() const;
+ const Field<LUType>& upper() const;
+ const Field<LUType>& lower() const;
+ const Field<Type>& source() const;
- FieldField<Field, LUType>& interfacesUpper()
- {
- return interfacesUpper_;
- }
+ Field<DType>& diag();
+ Field<LUType>& upper();
+ Field<LUType>& lower();
+ Field<Type>& source();
- FieldField<Field, LUType>& interfacesLower()
- {
- return interfacesLower_;
- }
+ // Interfaces
- const Field<DType>& diag() const;
- const Field<LUType>& upper() const;
- const Field<LUType>& lower() const;
- const Field<Type>& source() const;
+ const FieldField<Field, LUType>& interfacesUpper() const noexcept
+ {
+ return interfacesUpper_;
+ }
- const FieldField<Field, LUType>& interfacesUpper() const
- {
- return interfacesUpper_;
- }
+ const FieldField<Field, LUType>& interfacesLower() const noexcept
+ {
+ return interfacesLower_;
+ }
- const FieldField<Field, LUType>& interfacesLower() const
- {
- return interfacesLower_;
- }
+ FieldField<Field, LUType>& interfacesUpper() noexcept
+ {
+ return interfacesUpper_;
+ }
+ FieldField<Field, LUType>& interfacesLower() noexcept
+ {
+ return interfacesLower_;
+ }
- bool hasDiag() const noexcept
- {
- return (diagPtr_);
- }
- bool hasUpper() const noexcept
- {
- return (upperPtr_);
- }
+ // Characteristics
- bool hasLower() const noexcept
- {
- return (lowerPtr_);
- }
+ //- The matrix type (empty, diagonal, symmetric, ...)
+ word matrixTypeName() const;
- bool hasSource() const noexcept
- {
- return (sourcePtr_);
- }
+ 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 diagonal() const noexcept
- {
- return (diagPtr_ && !lowerPtr_ && !upperPtr_);
- }
+ //- Matrix has diagonal only
+ bool diagonal() const noexcept
+ {
+ return (diagPtr_ && !lowerPtr_ && !upperPtr_);
+ }
- bool symmetric() const noexcept
- {
- return (diagPtr_ && (!lowerPtr_ && upperPtr_));
- }
+ //- Matrix is symmetric
+ 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<Type, DType, LUType>&);
+ //- Move assignment
+ void operator=(LduMatrix<Type, DType, LUType>&&);
+
void negate();
void operator+=(const LduMatrix<Type, DType, LUType>&);
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C
index 7945a8cd2da82752cbb4740102af3fd788d0f911..46e7ad7b605cd058fdcb2898aca7e1eb58bb2d4c 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<Type, DType, LUType>::H(const Field<Type>& psi) const
{
auto tHpsi = tmp<Field<Type>>::New(lduAddr().size(), Foam::zero{});
- if (lowerPtr_ || upperPtr_)
+ if (hasLower() || hasUpper())
{
Type* __restrict__ HpsiPtr = tHpsi.ref().begin();
@@ -169,32 +169,30 @@ void Foam::LduMatrix<Type, DType, LUType>::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<Type, DType, LUType>::operator=(const LduMatrix& A)
}
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::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<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::negate()
{
@@ -235,12 +251,12 @@ void Foam::LduMatrix<Type, DType, LUType>::negate()
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::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<Type, DType, LUType>::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<Type, DType, LUType>::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<Type, DType, LUType>::operator+=(const LduMatrix& A)
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::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<Type, DType, LUType>::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<Type, DType, LUType>::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/SquareMatrix/SquareMatrix.C b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
index 2c64598fd6918851e3a7c9bf7379388c6a018229..fc5e12418479ce116e72c9f6e03bbf2b87a30953 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<Type>& matrix)
{
SquareMatrix<Type> matrixTmp = matrix;
- labelList pivotIndices(matrix.m());
+ labelList pivotIndices;
label sign;
LUDecompose(matrixTmp, pivotIndices, sign);
@@ -136,7 +136,7 @@ scalar det(const SquareMatrix<Type>& matrix)
template<class Type>
scalar det(SquareMatrix<Type>& matrix)
{
- labelList pivotIndices(matrix.m());
+ labelList pivotIndices;
label sign;
LUDecompose(matrix, pivotIndices, sign);
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C
index 61006d077c3ffc57cddbc9741cf296d11992774f..91180115b7f2c568dc0fb042fbb525349dfce927 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<labelList>(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<labelList>(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<labelList>(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::label, Foam::scalar> 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 c37cd61aec852a29cf1a44e4f87e0118f9bc5ff3..b610a6539d4b8805519015220ea80aa4dc573e4d 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<labelList> losortPtr_;
//- Owner start addressing
- mutable labelList* ownerStartPtr_;
+ mutable std::unique_ptr<labelList> ownerStartPtr_;
//- Losort start addressing
- mutable labelList* losortStartPtr_;
+ mutable std::unique_ptr<labelList> 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 dd3232541991f1077ccf41d0808c4da12282b712..e2767c6b5262da32109a81c1fd087ccee9b71d3e 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_)
+ if (A.diagPtr_)
{
- lowerPtr_ = new scalarField(*(A.lowerPtr_));
+ diagPtr_ = std::make_unique<scalarField>(*(A.diagPtr_));
}
- if (A.diagPtr_)
+ if (A.upperPtr_)
{
- diagPtr_ = new scalarField(*(A.diagPtr_));
+ upperPtr_ = std::make_unique<scalarField>(*(A.upperPtr_));
}
- if (A.upperPtr_)
+ if (A.lowerPtr_)
{
- upperPtr_ = new scalarField(*(A.upperPtr_));
+ lowerPtr_ = std::make_unique<scalarField>(*(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_)
+ // Copy assignment
+ if (A.diagPtr_)
{
- lowerPtr_ = new scalarField(*(A.lowerPtr_));
+ diagPtr_ = std::make_unique<scalarField>(*(A.diagPtr_));
}
- if (A.diagPtr_)
+ if (A.upperPtr_)
{
- diagPtr_ = new scalarField(*(A.diagPtr_));
+ upperPtr_ = std::make_unique<scalarField>(*(A.upperPtr_));
}
- if (A.upperPtr_)
+ if (A.lowerPtr_)
{
- upperPtr_ = new scalarField(*(A.upperPtr_));
+ lowerPtr_ = std::make_unique<scalarField>(*(A.lowerPtr_));
}
}
}
@@ -140,64 +127,56 @@ 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;
+
+ is >> withLower >> withDiag >> withUpper;
- if (hasLow)
+ if (withLower)
{
- lowerPtr_ = new scalarField(is);
+ lowerPtr_ = std::make_unique<scalarField>(is);
}
- if (hasDiag)
+ if (withDiag)
{
- diagPtr_ = new scalarField(is);
+ diagPtr_ = std::make_unique<scalarField>(is);
}
- if (hasUp)
+ if (withUpper)
{
- upperPtr_ = new scalarField(is);
+ upperPtr_ = std::make_unique<scalarField>(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_;
- }
+ // is empty (or just wrong)
+ return (!upperPtr_ && !lowerPtr_ ? "empty" : "ill-defined");
}
-Foam::scalarField& Foam::lduMatrix::lower()
+const Foam::scalarField& Foam::lduMatrix::diag() const
{
- if (!lowerPtr_)
+ if (!diagPtr_)
{
- if (upperPtr_)
- {
- lowerPtr_ = new scalarField(*upperPtr_);
- }
- else
- {
- lowerPtr_ = new scalarField(lduAddr().lowerAddr().size(), Zero);
- }
+ FatalErrorInFunction
+ << "diagPtr_ unallocated"
+ << abort(FatalError);
}
- return *lowerPtr_;
+ return *diagPtr_;
}
@@ -205,71 +184,87 @@ Foam::scalarField& Foam::lduMatrix::diag()
{
if (!diagPtr_)
{
- diagPtr_ = new scalarField(lduAddr().size(), Zero);
+ diagPtr_ =
+ std::make_unique<scalarField>(lduAddr().size(), Foam::zero{});
}
return *diagPtr_;
}
-Foam::scalarField& Foam::lduMatrix::upper()
+Foam::scalarField& Foam::lduMatrix::diag(label size)
{
- if (!upperPtr_)
+ if (!diagPtr_)
{
- if (lowerPtr_)
- {
- upperPtr_ = new scalarField(*lowerPtr_);
- }
- else
- {
- upperPtr_ = new scalarField(lduAddr().lowerAddr().size(), Zero);
- }
+ // if (size < 0)
+ // {
+ // size = lduAddr().size();
+ // }
+ diagPtr_ = std::make_unique<scalarField>(size, Foam::zero{});
}
- return *upperPtr_;
+ return *diagPtr_;
}
-Foam::scalarField& Foam::lduMatrix::lower(const label nCoeffs)
+const Foam::scalarField& Foam::lduMatrix::upper() const
{
- if (!lowerPtr_)
+ if (upperPtr_)
{
- if (upperPtr_)
- {
- lowerPtr_ = new scalarField(*upperPtr_);
- }
- else
+ return *upperPtr_;
+ }
+ else
+ {
+ if (!lowerPtr_)
{
- lowerPtr_ = new scalarField(nCoeffs, Zero);
+ FatalErrorInFunction
+ << "lowerPtr_ and upperPtr_ unallocated"
+ << abort(FatalError);
}
- }
- return *lowerPtr_;
+ return *lowerPtr_;
+ }
}
-Foam::scalarField& Foam::lduMatrix::diag(const label size)
+Foam::scalarField& Foam::lduMatrix::upper()
{
- if (!diagPtr_)
+ if (!upperPtr_)
{
- diagPtr_ = new scalarField(size, Zero);
+ if (lowerPtr_)
+ {
+ upperPtr_ = std::make_unique<scalarField>(*lowerPtr_);
+ }
+ else
+ {
+ upperPtr_ =
+ std::make_unique<scalarField>
+ (
+ lduAddr().lowerAddr().size(),
+ Foam::zero{}
+ );
+ }
}
- return *diagPtr_;
+ return *upperPtr_;
}
-Foam::scalarField& Foam::lduMatrix::upper(const label nCoeffs)
+Foam::scalarField& Foam::lduMatrix::upper(label nCoeffs)
{
if (!upperPtr_)
{
if (lowerPtr_)
{
- upperPtr_ = new scalarField(*lowerPtr_);
+ upperPtr_ = std::make_unique<scalarField>(*lowerPtr_);
}
else
{
- upperPtr_ = new scalarField(nCoeffs, Zero);
+ // if (nCoeffs < 0)
+ // {
+ // nCoeffs = lduAddr().lowerAddr().size();
+ // }
+ upperPtr_ = std::make_unique<scalarField>(nCoeffs, Foam::zero{});
}
}
@@ -279,54 +274,67 @@ Foam::scalarField& Foam::lduMatrix::upper(const label nCoeffs)
const Foam::scalarField& Foam::lduMatrix::lower() const
{
- if (!lowerPtr_ && !upperPtr_)
- {
- FatalErrorInFunction
- << "lowerPtr_ or upperPtr_ unallocated"
- << abort(FatalError);
- }
-
if (lowerPtr_)
{
return *lowerPtr_;
}
else
{
+ if (!upperPtr_)
+ {
+ FatalErrorInFunction
+ << "lowerPtr_ and upperPtr_ unallocated"
+ << abort(FatalError);
+ }
+
return *upperPtr_;
}
}
-const Foam::scalarField& Foam::lduMatrix::diag() const
+Foam::scalarField& Foam::lduMatrix::lower()
{
- if (!diagPtr_)
+ if (!lowerPtr_)
{
- FatalErrorInFunction
- << "diagPtr_ unallocated"
- << abort(FatalError);
+ if (upperPtr_)
+ {
+ lowerPtr_ = std::make_unique<scalarField>(*upperPtr_);
+ }
+ else
+ {
+ lowerPtr_ =
+ std::make_unique<scalarField>
+ (
+ lduAddr().lowerAddr().size(),
+ Foam::zero{}
+ );
+ }
}
- return *diagPtr_;
+ return *lowerPtr_;
}
-const Foam::scalarField& Foam::lduMatrix::upper() const
+Foam::scalarField& Foam::lduMatrix::lower(label nCoeffs)
{
- if (!lowerPtr_ && !upperPtr_)
+ if (!lowerPtr_)
{
- FatalErrorInFunction
- << "lowerPtr_ or upperPtr_ unallocated"
- << abort(FatalError);
+ if (upperPtr_)
+ {
+ lowerPtr_ = std::make_unique<scalarField>(*upperPtr_);
+ }
+ else
+ {
+ // if (nCoeffs < 0)
+ // {
+ // nCoeffs = lduAddr().lowerAddr().size();
+ // }
+ lowerPtr_ =
+ std::make_unique<scalarField>(nCoeffs, Foam::zero{});
+ }
}
- if (upperPtr_)
- {
- return *upperPtr_;
- }
- else
- {
- return *lowerPtr_;
- }
+ return *lowerPtr_;
}
@@ -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<lduMatrix>& iproxy
)
{
- const auto& ldum = *iproxy;
-
- Switch hasLow = ldum.hasLower();
- Switch hasDiag = ldum.hasDiag();
- Switch hasUp = ldum.hasUpper();
+ const auto& mat = *iproxy;
- os << "Lower:" << hasLow
- << " Diag:" << hasDiag
- << " Upper:" << hasUp << endl;
+ os << "Lower:" << Switch::name(mat.hasLower())
+ << " Diag:" << Switch::name(mat.hasDiag())
+ << " Upper:" << Switch::name(mat.hasUpper()) << 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 e13bb4cd05764a1561598805afbaa48c99372b00..a4dd1b0601732f8457b896258cc32025dcb2e88a 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 <functional> // For reference_wrapper
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -89,8 +90,14 @@ class lduMatrix
//const lduMesh& lduMesh_;
std::reference_wrapper<const lduMesh> lduMesh_;
- //- Coefficients (not including interfaces)
- scalarField *lowerPtr_, *diagPtr_, *upperPtr_;
+ //- Diagonal coefficients
+ std::unique_ptr<scalarField> diagPtr_;
+
+ //- Off-diagonal coefficients (not including interfaces)
+ std::unique_ptr<scalarField> lowerPtr_;
+
+ //- Off-diagonal coefficients (not including interfaces)
+ std::unique_ptr<scalarField> 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;
- 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 diagonal() const noexcept
- {
- return (diagPtr_ && !lowerPtr_ && !upperPtr_);
- }
+ //- Matrix has diagonal only
+ bool diagonal() const noexcept
+ {
+ return (diagPtr_ && !lowerPtr_ && !upperPtr_);
+ }
- bool symmetric() const noexcept
- {
- return (diagPtr_ && (!lowerPtr_ && upperPtr_));
- }
+ //- Matrix is symmetric
+ 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 4d307f57100f385789528a2a34001fa70f0f4217..1e83a03c75b85e8d277ed75820bfc21bca41861a 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;
}
}
}
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
index 39d3cc37b50ce5b88f5362489809e55162e5bed0..1bd42a36080dad51a0be66547ee16e168659252a 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<scalar> 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 d6ba156a2bcaa6479fd82f12d9c12415ca39f3a8..8035cdfebfe01854b606a7c68fec852a3e5c53bd 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
@@ -74,19 +74,21 @@ void solve
const List<Type>& 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 0c3ae682287d35e9666a09d32407975a65780643..7a66cbe6d3459f5bdeef62d2cb0649dd655053f2 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
@@ -216,7 +216,7 @@ void Foam::LUsolve
List<Type>& sourceSol
)
{
- labelList pivotIndices(matrix.m());
+ labelList pivotIndices;
LUDecompose(matrix, pivotIndices);
LUBacksubstitute(matrix, pivotIndices, sourceSol);
}