zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'update-matrix-streaming' into 'develop'

Browse Source 
Improve streaming matrix and internal memory management

See merge request Development/openfoam!671
This commit is contained in:
Kutalmış Berçin
2024-03-07 10:18:27 +00:00
parent 05194796ef 27aa7e4e91
commit d89ecc74be
19 changed files with 1012 additions and 894 deletions
Download Patch File Download Diff File Expand all files Collapse all files

168
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
View File

@ -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 (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<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_);
}

28
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
View File

@ -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.

168
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
View File

@ -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)
if (Pstream::master(comm_))
const label startOfRequests = UPstream::nRequests();
// 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 (Pstream::master(comm_))
{
LUBacksubstitute(*this, pivotIndices_, X);
x = SubList<Type>(X, x.size());
for (const int proci : Pstream::subProcs(comm_))
if (x.empty())
{
// Nothing to do
}
else if (is_contiguous<Type>::value)
{
UOPstream::write
(
Pstream::commsTypes::scheduled,
proci,
reinterpret_cast<const char*>
(
&(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<Type>(allx, x.size());
for (const int proci : UPstream::subProcs(comm_))
{
SubList<Type> procSlot
(
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);

51
src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.C
View File

@ -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;
}

14
src/OpenFOAM/matrices/LUscalarMatrix/procLduInterface.H
View File

@ -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&);
};

53
src/OpenFOAM/matrices/LUscalarMatrix/procLduMatrix.C
View File

@ -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;
}

48
src/OpenFOAM/matrices/LUscalarMatrix/procLduMatrix.H
View File

@ -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;
// 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&);
};

326
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C
View File

@ -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)
:
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<class Type, class DType, class LUType>
Foam::LduMatrix<Type, DType, LUType>::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<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,109 +130,27 @@ 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 * * * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::LduMatrix<Type, DType, LUType>::~LduMatrix()
{
if (diagPtr_)
{
delete diagPtr_;
}
if (upperPtr_)
{
delete upperPtr_;
}
if (lowerPtr_)
{
delete lowerPtr_;
}
if (sourcePtr_)
{
delete sourcePtr_;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag()
Foam::word Foam::LduMatrix<Type, DType, LUType>::matrixTypeName() const
{
if (!diagPtr_)
if (diagPtr_)
{
diagPtr_ = new Field<DType>(lduAddr().size(), Zero);
return
(
(!upperPtr_)
? (!lowerPtr_ ? "diagonal" : "diagonal-lower")
: (!lowerPtr_ ? "symmetric" : "asymmetric")
);
}
return *diagPtr_;
}
template<class Type, class DType, class LUType>
Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper()
{
if (!upperPtr_)
{
if (lowerPtr_)
{
upperPtr_ = new Field<LUType>(*lowerPtr_);
}
else
{
upperPtr_ = new Field<LUType>
(
lduAddr().lowerAddr().size(),
Zero
);
}
}
return *upperPtr_;
}
template<class Type, class DType, class LUType>
Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower()
{
if (!lowerPtr_)
{
if (upperPtr_)
{
lowerPtr_ = new Field<LUType>(*upperPtr_);
}
else
{
lowerPtr_ = new Field<LUType>
(
lduAddr().lowerAddr().size(),
Zero
);
}
}
return *lowerPtr_;
}
template<class Type, class DType, class LUType>
Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source()
{
if (!sourcePtr_)
{
sourcePtr_ = new Field<Type>(lduAddr().size(), Zero);
}
return *sourcePtr_;
// is empty (or just wrong)
return (!upperPtr_ && !lowerPtr_ ? "empty" : "ill-defined");
}
@ -273,47 +169,107 @@ const Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag() const
template<class Type, class DType, class LUType>
const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper() const
Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag()
{
if (!lowerPtr_ && !upperPtr_)
if (!diagPtr_)
{
FatalErrorInFunction
<< "lowerPtr_ or upperPtr_ unallocated"
<< abort(FatalError);
diagPtr_ =
std::make_unique<Field<DType>>(lduAddr().size(), Foam::zero{});
}
return *diagPtr_;
}
template<class Type, class DType, class LUType>
const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper() const
{
if (upperPtr_)
{
return *upperPtr_;
}
else
{
if (!lowerPtr_)
{
FatalErrorInFunction
<< "lowerPtr_ and upperPtr_ unallocated"
<< abort(FatalError);
}
return *lowerPtr_;
}
}
template<class Type, class DType, class LUType>
const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower() const
Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper()
{
if (!lowerPtr_ && !upperPtr_)
if (!upperPtr_)
{
FatalErrorInFunction
<< "lowerPtr_ or upperPtr_ unallocated"
<< abort(FatalError);
if (lowerPtr_)
{
upperPtr_ = std::make_unique<Field<LUType>>(*lowerPtr_);
}
else
{
upperPtr_ =
std::make_unique<Field<LUType>>
(
lduAddr().lowerAddr().size(),
Foam::zero{}
);
}
}
return *upperPtr_;
}
template<class Type, class DType, class LUType>
const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower() const
{
if (lowerPtr_)
{
return *lowerPtr_;
}
else
{
if (!upperPtr_)
{
FatalErrorInFunction
<< "lowerPtr_ and upperPtr_ unallocated"
<< abort(FatalError);
}
return *upperPtr_;
}
}
template<class Type, class DType, class LUType>
Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower()
{
if (!lowerPtr_)
{
if (upperPtr_)
{
lowerPtr_ = std::make_unique<Field<LUType>>(*upperPtr_);
}
else
{
lowerPtr_ = std::make_unique<Field<LUType>>
(
lduAddr().lowerAddr().size(),
Foam::zero{}
);
}
}
return *lowerPtr_;
}
template<class Type, class DType, class LUType>
const Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source() const
{
@ -328,45 +284,65 @@ const Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source() const
}
template<class Type, class DType, class LUType>
Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source()
{
if (!sourcePtr_)
{
sourcePtr_ =
std::make_unique<Field<Type>>(lduAddr().size(), Foam::zero{});
}
return *sourcePtr_;
}
// * * * * * * * * * * * * * * * 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;
}

196
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
View File

@ -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_;
}
FieldField<Field, LUType>& interfacesLower()
{
return interfacesLower_;
}
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;
// Interfaces
const FieldField<Field, LUType>& interfacesUpper() const
{
return interfacesUpper_;
}
const FieldField<Field, LUType>& interfacesUpper() const noexcept
{
return interfacesUpper_;
}
const FieldField<Field, LUType>& interfacesLower() const
{
return interfacesLower_;
}
const FieldField<Field, LUType>& interfacesLower() const noexcept
{
return interfacesLower_;
}
FieldField<Field, LUType>& interfacesUpper() noexcept
{
return interfacesUpper_;
}
FieldField<Field, LUType>& 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<Type, DType, LUType>&);
//- Move assignment
void operator=(LduMatrix<Type, DType, LUType>&&);
void negate();
void operator+=(const LduMatrix<Type, DType, LUType>&);

60
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C
View File

@ -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();
}

6
src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
View File

@ -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);

38
src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C
View File

@ -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)
{

42
src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H
View File

@ -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_;
}

445
src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
View File

@ -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<scalarField>(*(A.diagPtr_));
}
if (A.upperPtr_)
{
upperPtr_ = new scalarField(*(A.upperPtr_));
upperPtr_ = std::make_unique<scalarField>(*(A.upperPtr_));
}
if (A.lowerPtr_)
{
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_)
{
lowerPtr_ = new scalarField(*(A.lowerPtr_));
}
// Copy assignment
if (A.diagPtr_)
{
diagPtr_ = new scalarField(*(A.diagPtr_));
diagPtr_ = std::make_unique<scalarField>(*(A.diagPtr_));
}
if (A.upperPtr_)
{
upperPtr_ = new scalarField(*(A.upperPtr_));
upperPtr_ = std::make_unique<scalarField>(*(A.upperPtr_));
}
if (A.lowerPtr_)
{
lowerPtr_ = std::make_unique<scalarField>(*(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<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_;
}
}
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<scalarField>(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<scalarField>(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<scalarField>(*lowerPtr_);
}
else
{
upperPtr_ =
std::make_unique<scalarField>
(
lduAddr().lowerAddr().size(),
Foam::zero{}
);
}
}
return *upperPtr_;
}
Foam::scalarField& Foam::lduMatrix::upper(label nCoeffs)
{
if (!upperPtr_)
{
if (lowerPtr_)
{
upperPtr_ = std::make_unique<scalarField>(*lowerPtr_);
}
else
{
// if (nCoeffs < 0)
// {
// nCoeffs = lduAddr().lowerAddr().size();
// }
upperPtr_ = std::make_unique<scalarField>(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<scalarField>(*upperPtr_);
}
else
{
lowerPtr_ =
std::make_unique<scalarField>
(
lduAddr().lowerAddr().size(),
Foam::zero{}
);
}
}
return *lowerPtr_;
}
Foam::scalarField& Foam::lduMatrix::lower(label nCoeffs)
{
if (!lowerPtr_)
{
if (upperPtr_)
{
lowerPtr_ = std::make_unique<scalarField>(*upperPtr_);
}
else
{
// if (nCoeffs < 0)
// {
// nCoeffs = lduAddr().lowerAddr().size();
// }
lowerPtr_ =
std::make_unique<scalarField>(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<lduMatrix>& 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;
//}

152
src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
View File

@ -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;
// 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&);

77
src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C
View File

@ -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;
}
}
}

26
src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
View File

@ -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)
{

6
src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
View File

@ -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
);

2
src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
View File

@ -216,7 +216,7 @@ void Foam::LUsolve
List<Type>& sourceSol
)
{
labelList pivotIndices(matrix.m());
labelList pivotIndices;
LUDecompose(matrix, pivotIndices);
LUBacksubstitute(matrix, pivotIndices, sourceSol);
}