songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'feature-csr-MERGE' into 'develop'

Browse Source 
ENH: CSR: cell-cell addressing

See merge request Development/openfoam!713
This commit is contained in:
Andrew Heather 2024-12-10 09:30:32 +00:00
commit ff6ef5bdf8
8 changed files with 350 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -169,6 +169,20 @@ void Foam::lduAddressing::calcLosortStart() const
}
void Foam::lduAddressing::calcLoCSR() const
{
if (lowerCSRAddrPtr_)
{
FatalErrorInFunction
<< "lowerCSRAddr already calculated"
<< abort(FatalError);
}
lowerCSRAddrPtr_ = std::make_unique<labelList>(lowerAddr().size());
map(lowerAddr(), *lowerCSRAddrPtr_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::labelUList& Foam::lduAddressing::losortAddr() const
@ -204,11 +218,23 @@ const Foam::labelUList& Foam::lduAddressing::losortStartAddr() const
}
const Foam::labelUList& Foam::lduAddressing::lowerCSRAddr() const
{
if (!lowerCSRAddrPtr_)
{
calcLoCSR();
}
return *lowerCSRAddrPtr_;
}
void Foam::lduAddressing::clearOut()
{
losortPtr_.reset(nullptr);
ownerStartPtr_.reset(nullptr);
losortStartPtr_.reset(nullptr);
lowerCSRAddrPtr_.reset(nullptr);
}

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

@ -91,6 +91,10 @@ Description
list. Thus, for every point the losort start gives the address of the
first face to neighbour this point.
Instead of using losort to lookup the face and then using the lowerAddr
to find the neighbour cell one can also directly lookup the neighbour cell
using the lowerCSRAddr (upperAddr is already in CSR order).
SourceFiles
lduAddressing.C
@ -131,6 +135,9 @@ class lduAddressing
//- Losort start addressing
mutable std::unique_ptr<labelList> losortStartPtr_;
//- Lower addressing
mutable std::unique_ptr<labelList> lowerCSRAddrPtr_;
// Private Member Functions
@ -143,6 +150,9 @@ class lduAddressing
//- Calculate losort start
void calcLosortStart() const;
//- Calculate CSR lower addressing
void calcLoCSR() const;
public:
@ -203,11 +213,22 @@ public:
//- Return losort start addressing
const labelUList& losortStartAddr() const;
//- Return CSR addressing
const labelUList& lowerCSRAddr() const;
//- Return off-diagonal index given owner and neighbour label
label triIndex(const label a, const label b) const;
//- Calculate bandwidth and profile of addressing
Tuple2<label, scalar> band() const;
//- Helper to convert lower addressing & data into CSR format
template<class Type>
void map
(
const UList<Type>& faceVals,
List<Type>& vals
) const;
};
@ -217,6 +238,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "lduAddressingTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

60
src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressingTemplates.C Normal file
View File

@ -0,0 +1,60 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "lduAddressing.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::lduAddressing::map
(
const UList<Type>& faceVals,
List<Type>& vals
) const
{
const auto& offsets = losortStartAddr();
const auto& indexToFace = losortAddr();
const label n = offsets.size()-1;
vals.resize_nocopy(faceVals.size());
for (label celli = 0; celli < n; celli++)
{
const label start = offsets[celli];
const label end = offsets[celli+1];
for (label i = start; i < end; i++)
{
const label facei = indexToFace[i];
vals[i] = faceVals[facei];
}
}
}
// ************************************************************************* //

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

@ -82,6 +82,13 @@ Foam::lduMatrix::lduMatrix(const lduMatrix& A)
{
lowerPtr_ = std::make_unique<scalarField>(*(A.lowerPtr_));
}
if (A.lowerCSRPtr_)
{
lowerCSRPtr_ = std::make_unique<scalarField>(*(A.lowerCSRPtr_));
}
// No need to copy work
}
@ -90,7 +97,9 @@ Foam::lduMatrix::lduMatrix(lduMatrix&& A)
lduMesh_(A.lduMesh_),
diagPtr_(std::move(A.diagPtr_)),
lowerPtr_(std::move(A.lowerPtr_)),
upperPtr_(std::move(A.upperPtr_))
upperPtr_(std::move(A.upperPtr_)),
lowerCSRPtr_(std::move(A.lowerCSRPtr_)),
workPtr_(std::move(A.workPtr_))
{}
@ -103,6 +112,8 @@ Foam::lduMatrix::lduMatrix(lduMatrix& A, bool reuse)
diagPtr_ = std::move(A.diagPtr_);
upperPtr_ = std::move(A.upperPtr_);
lowerPtr_ = std::move(A.lowerPtr_);
lowerCSRPtr_ = std::move(A.lowerCSRPtr_);
workPtr_ = std::move(A.workPtr_);
}
else
{
@ -121,6 +132,13 @@ Foam::lduMatrix::lduMatrix(lduMatrix& A, bool reuse)
{
lowerPtr_ = std::make_unique<scalarField>(*(A.lowerPtr_));
}
// Note: no real need to keep lowerCSR except we use it (hasLowerCSR())
// to trigger certain actions
if (A.lowerCSRPtr_)
{
lowerCSRPtr_ = std::make_unique<scalarField>(*(A.lowerCSRPtr_));
}
}
}
@ -129,9 +147,9 @@ Foam::lduMatrix::lduMatrix(const lduMesh& mesh, Istream& is)
:
lduMesh_(mesh)
{
bool withLower, withDiag, withUpper;
bool withLower, withDiag, withUpper, withLowerCSR;
is >> withLower >> withDiag >> withUpper;
is >> withLower >> withDiag >> withUpper >> withLowerCSR;
if (withLower)
{
@ -145,6 +163,20 @@ Foam::lduMatrix::lduMatrix(const lduMesh& mesh, Istream& is)
{
upperPtr_ = std::make_unique<scalarField>(is);
}
if (withLowerCSR)
{
// Check if it has been sent over or we need to construct it. We need
// to set the lowerCSRPtr since it determines behaviour.
if (withLower)
{
// Construct from lower
(void)lowerCSR();
}
else
{
lowerCSRPtr_ = std::make_unique<scalarField>(is);
}
}
}
@ -237,6 +269,9 @@ Foam::scalarField& Foam::lduMatrix::upper()
}
else
{
// no lowerPtr so any lowerCSR was constructed from upper
lowerCSRPtr_.reset(nullptr);
upperPtr_ =
std::make_unique<scalarField>
(
@ -260,6 +295,9 @@ Foam::scalarField& Foam::lduMatrix::upper(label nCoeffs)
}
else
{
// no lowerPtr so any lowerCSR was constructed from upper
lowerCSRPtr_.reset(nullptr);
// if (nCoeffs < 0)
// {
// nCoeffs = lduAddr().lowerAddr().size();
@ -296,6 +334,8 @@ Foam::scalarField& Foam::lduMatrix::lower()
{
if (!lowerPtr_)
{
lowerCSRPtr_.reset(nullptr);
if (upperPtr_)
{
lowerPtr_ = std::make_unique<scalarField>(*upperPtr_);
@ -319,6 +359,8 @@ Foam::scalarField& Foam::lduMatrix::lower(label nCoeffs)
{
if (!lowerPtr_)
{
lowerCSRPtr_.reset(nullptr);
if (upperPtr_)
{
lowerPtr_ = std::make_unique<scalarField>(*upperPtr_);
@ -338,6 +380,86 @@ Foam::scalarField& Foam::lduMatrix::lower(label nCoeffs)
}
const Foam::scalarField& Foam::lduMatrix::lowerCSR() const
{
if (!lowerCSRPtr_)
{
const label nLower = lduAddr().losortAddr().size();
lowerCSRPtr_ = std::make_unique<scalarField>(nLower);
if (lowerPtr_)
{
lduAddr().map(*lowerPtr_, *lowerCSRPtr_);
}
else if (upperPtr_)
{
lduAddr().map(*upperPtr_, *lowerCSRPtr_);
}
else
{
FatalErrorInFunction
<< "lowerPtr_ and upperPtr_ unallocated"
<< abort(FatalError);
}
}
return *lowerCSRPtr_;
}
Foam::scalarField& Foam::lduMatrix::lowerCSR()
{
if (!lowerCSRPtr_)
{
const label nLower = lduAddr().losortAddr().size();
lowerCSRPtr_ = std::make_unique<scalarField>(nLower);
if (lowerPtr_)
{
lduAddr().map(*lowerPtr_, *lowerCSRPtr_);
}
else if (upperPtr_)
{
lduAddr().map(*upperPtr_, *lowerCSRPtr_);
}
else
{
FatalErrorInFunction
<< "lowerPtr_ and upperPtr_ unallocated"
<< abort(FatalError);
}
}
return *lowerCSRPtr_;
}
Foam::solveScalarField& Foam::lduMatrix::work(label size) const
{
if (!workPtr_ || workPtr_->size() != size)
{
workPtr_ = std::make_unique<solveScalarField>(size, Foam::zero{});
}
return *workPtr_;
}
const Foam::solveScalarField& Foam::lduMatrix::work() const
{
if (!workPtr_)
{
FatalErrorInFunction
<< "workPtr_ unallocated"
<< abort(FatalError);
}
return *workPtr_;
}
void Foam::lduMatrix::setResidualField
(
const scalarField& residual,
@ -389,7 +511,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& mat)
{
os << mat.hasLower() << token::SPACE
<< mat.hasDiag() << token::SPACE
<< mat.hasUpper() << token::SPACE;
<< mat.hasUpper() << token::SPACE
<< mat.hasLowerCSR() << token::SPACE;
if (mat.hasLower())
{
@ -406,6 +529,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& mat)
os << mat.upper();
}
if (mat.hasLowerCSR() && !mat.hasLower())
{
// Only send over if can not be reconstructed locally
os << mat.lowerCSR();
}
os.check(FUNCTION_NAME);
return os;
@ -422,7 +551,9 @@ Foam::Ostream& Foam::operator<<
os << "Lower:" << Switch::name(mat.hasLower())
<< " Diag:" << Switch::name(mat.hasDiag())
<< " Upper:" << Switch::name(mat.hasUpper()) << endl;
<< " Upper:" << Switch::name(mat.hasUpper())
<< " lowerCSR:" << Switch::name(mat.hasLowerCSR())
<< endl;
if (mat.hasLower())
{

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

@ -99,6 +99,12 @@ class lduMatrix
//- Off-diagonal coefficients (not including interfaces)
std::unique_ptr<scalarField> upperPtr_;
//- Off-diagonal coefficients (not including interfaces) in CSR ordering
mutable std::unique_ptr<scalarField> lowerCSRPtr_;
//- Work space
mutable std::unique_ptr<solveScalarField> workPtr_;
public:
@ -671,6 +677,21 @@ public:
scalarField& lower(label nCoeffs);
// CSR
bool hasLowerCSR() const noexcept { return bool(lowerCSRPtr_); }
const scalarField& lowerCSR() const;
scalarField& lowerCSR();
//- Work array
const solveScalarField& work() const;
//- Work array
solveScalarField& work(label size) const;
// Characteristics
//- The matrix type (empty, diagonal, symmetric, ...)
@ -683,19 +704,19 @@ public:
//- Matrix has diagonal only
bool diagonal() const noexcept
{
return (diagPtr_ && !lowerPtr_ && !upperPtr_);
return (diagPtr_ && !lowerPtr_ && !lowerCSRPtr_ && !upperPtr_);
}
//- Matrix is symmetric
bool symmetric() const noexcept
{
return (diagPtr_ && !lowerPtr_ && upperPtr_);
return (diagPtr_ && !(lowerPtr_ || lowerCSRPtr_) && upperPtr_);
}
//- Matrix is asymmetric (ie, full)
bool asymmetric() const noexcept
{
return (diagPtr_ && lowerPtr_ && upperPtr_);
return (diagPtr_ && (lowerPtr_ || lowerCSRPtr_) && upperPtr_);
}

83
src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2017-2019 OpenCFD Ltd.
Copyright (C) 2017-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -43,6 +43,8 @@ void Foam::lduMatrix::Amul
const direction cmpt
) const
{
const auto& addr = lduAddr();
solveScalar* __restrict__ ApsiPtr = Apsi.begin();
const solveScalarField& psi = tpsi();
@ -50,8 +52,8 @@ void Foam::lduMatrix::Amul
const scalar* const __restrict__ diagPtr = diag().begin();
const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
const label* const __restrict__ uPtr = addr.upperAddr().begin();
const label* const __restrict__ lPtr = addr.lowerAddr().begin();
const scalar* const __restrict__ upperPtr = upper().begin();
const scalar* const __restrict__ lowerPtr = lower().begin();
@ -70,18 +72,75 @@ void Foam::lduMatrix::Amul
);
const label nCells = diag().size();
for (label cell=0; cell<nCells; cell++)
if (hasLowerCSR())
{
ApsiPtr[cell] = diagPtr[cell]*psiPtr[cell];
// Use cell-based looping
if (debug == 2) PoutInFunction<< "cell-based looping" << endl;
const label* const __restrict__ oStartPtr =
addr.ownerStartAddr().begin();
const label* const __restrict__ loStartPtr =
addr.losortStartAddr().begin();
const label* const __restrict__ lcsrPtr =
addr.lowerCSRAddr().begin();
// Note: lowerCSR constructed from lower if available, upper otherwise
// so is handling symmetric()
const scalar* const __restrict__ lowercsrPtr = lowerCSR().begin();
for (label cell=0; cell<nCells; cell++)
{
scalar& val = ApsiPtr[cell];
//Pout<< "cell:" << cell << endl;
val = diagPtr[cell]*psiPtr[cell];
// Add lower contributions
{
const label start = loStartPtr[cell];
const label end = loStartPtr[cell+1];
for (label i = start; i < end; i++)
{
const label nbrCell = lcsrPtr[i];
//Pout<< " adding from " << nbrCell
// << " coeff:" << lowercsrPtr[i] << endl;
val += lowercsrPtr[i]*psiPtr[nbrCell];
}
}
// Add upper contributions
{
const label start = oStartPtr[cell];
const label end = oStartPtr[cell+1];
for (label i = start; i < end; i++)
{
const label nbrCell = uPtr[i];
//Pout<< " adding from " << nbrCell
// << " coeff:" << upperPtr[i] << endl;
val += upperPtr[i]*psiPtr[nbrCell];
}
}
//Pout<< "cell:" << cell << " val after:" << val << endl;
}
}
const label nFaces = upper().size();
for (label face=0; face<nFaces; face++)
else
{
ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];
ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];
for (label cell=0; cell<nCells; cell++)
{
ApsiPtr[cell] = diagPtr[cell]*psiPtr[cell];
}
const label nFaces = upper().size();
for (label face=0; face<nFaces; face++)
{
ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];
ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];
}
}
// Update interface interfaces

4
src/OpenFOAM/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.C
View File

@ -83,7 +83,9 @@ void Foam::GaussSeidelSmoother::smooth
const label nCells = psi.size();
solveScalarField bPrime(nCells);
//solveScalarField bPrime(nCells);
solveScalarField& bPrime = matrix_.work(nCells);
solveScalar* __restrict__ bPrimePtr = bPrime.begin();
const scalar* const __restrict__ diagPtr = matrix_.diag().begin();

4
src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C
View File

@ -83,7 +83,9 @@ void Foam::symGaussSeidelSmoother::smooth
const label nCells = psi.size();
solveScalarField bPrime(nCells);
//solveScalarField bPrime(nCells);
solveScalarField& bPrime = matrix_.work(nCells);
solveScalar* __restrict__ bPrimePtr = bPrime.begin();
const scalar* const __restrict__ diagPtr = matrix_.diag().begin();