68 comm_(ldum.
mesh().comm())
133 nCells += lduMatrices[i].
size();
138 convert(lduMatrices);
146 convert(ldum, interfaceCoeffs, interfaces);
153 const label numRows =
m();
154 const label numCols =
n();
156 Pout<<
"LUscalarMatrix : size:" << numRows <<
endl;
157 for (label rowi = 0; rowi < numRows; ++rowi)
161 Pout<<
"cell:" << rowi <<
" diagCoeff:" << row[rowi] <<
endl;
163 Pout<<
" connects to upper cells :";
164 for (label coli = rowi+1; coli < numCols; ++coli)
166 if (
mag(row[coli]) > SMALL)
168 Pout<<
' ' << coli <<
" (coeff:" << row[coli] <<
')';
172 Pout<<
" connects to lower cells :";
173 for (label coli = 0; coli < rowi; ++coli)
175 if (
mag(row[coli]) > SMALL)
177 Pout<<
' ' << coli <<
" (coeff:" << row[coli] <<
')';
193void Foam::LUscalarMatrix::convert
203 const scalar* __restrict__ diagPtr = ldum.
diag().
begin();
204 const scalar* __restrict__ upperPtr = ldum.
upper().
begin();
205 const scalar* __restrict__ lowerPtr = ldum.
lower().
begin();
207 const label nCells = ldum.
diag().
size();
215 for (label face=0; face<
nFaces; face++)
217 label uCell = uPtr[face];
218 label lCell = lPtr[face];
220 operator[](uCell)[lCell] = lowerPtr[face];
221 operator[](lCell)[uCell] = upperPtr[face];
226 if (interfaces.
set(inti))
228 const lduInterface&
interface = interfaces[inti].
interface();
232 const label* __restrict__ lPtr =
interface.faceCells().begin();
234 const cyclicLduInterface& cycInterface =
235 refCast<const cyclicLduInterface>(
interface);
236 label nbrInt = cycInterface.neighbPatchID();
237 const label* __restrict__ uPtr =
238 interfaces[nbrInt].interface().faceCells().
begin();
240 const scalar* __restrict__ nbrUpperLowerPtr =
241 interfaceCoeffs[nbrInt].
begin();
243 label inFaces =
interface.faceCells().size();
245 for (label face=0; face<inFaces; face++)
247 label uCell = lPtr[face];
248 label lCell = uPtr[face];
250 operator[](uCell)[lCell] -= nbrUpperLowerPtr[face];
257void Foam::LUscalarMatrix::convert
259 const PtrList<procLduMatrix>& lduMatrices
262 procOffsets_.setSize(lduMatrices.size() + 1);
265 forAll(lduMatrices, ldumi)
267 procOffsets_[ldumi+1] = procOffsets_[ldumi] + lduMatrices[ldumi].size();
270 forAll(lduMatrices, ldumi)
272 const procLduMatrix& lduMatrixi = lduMatrices[ldumi];
273 label offset = procOffsets_[ldumi];
275 const label* __restrict__ uPtr = lduMatrixi.upperAddr_.begin();
276 const label* __restrict__ lPtr = lduMatrixi.lowerAddr_.begin();
278 const scalar* __restrict__ diagPtr = lduMatrixi.diag_.begin();
279 const scalar* __restrict__ upperPtr = lduMatrixi.upper_.begin();
280 const scalar* __restrict__ lowerPtr = lduMatrixi.lower_.begin();
282 const label nCells = lduMatrixi.
size();
283 const label
nFaces = lduMatrixi.upper_.
size();
285 for (label cell=0; cell<nCells; cell++)
287 label globalCell = cell + offset;
288 operator[](globalCell)[globalCell] = diagPtr[cell];
291 for (label face=0; face<
nFaces; face++)
293 label uCell = uPtr[face] + offset;
294 label lCell = lPtr[face] + offset;
296 operator[](uCell)[lCell] = lowerPtr[face];
297 operator[](lCell)[uCell] = upperPtr[face];
300 const PtrList<procLduInterface>& interfaces =
301 lduMatrixi.interfaces_;
305 const procLduInterface&
interface = interfaces[inti];
309 const label* __restrict__ ulPtr =
interface.faceCells_.begin();
311 const scalar* __restrict__ upperLowerPtr =
314 label inFaces =
interface.faceCells_.size()/2;
316 for (label face=0; face<inFaces; face++)
318 label uCell = ulPtr[face] + offset;
319 label lCell = ulPtr[face + inFaces] + offset;
321 operator[](uCell)[lCell] -= upperLowerPtr[face + inFaces];
322 operator[](lCell)[uCell] -= upperLowerPtr[face];
332 const PtrList<procLduInterface>& neiInterfaces =
333 lduMatrices[
interface.neighbProcNo_].interfaces_;
335 label neiInterfacei = -1;
337 forAll(neiInterfaces, ninti)
342 neiInterfaces[ninti].neighbProcNo_
345 && (neiInterfaces[ninti].tag_ ==
interface.tag_)
348 neiInterfacei = ninti;
353 if (neiInterfacei == -1)
358 const procLduInterface& neiInterface =
359 neiInterfaces[neiInterfacei];
361 const label* __restrict__ uPtr =
interface.faceCells_.begin();
362 const label* __restrict__ lPtr =
363 neiInterface.faceCells_.begin();
365 const scalar* __restrict__ upperPtr =
interface.coeffs_.begin();
366 const scalar* __restrict__ lowerPtr =
367 neiInterface.coeffs_.begin();
369 label inFaces =
interface.faceCells_.size();
370 label neiOffset = procOffsets_[
interface.neighbProcNo_];
372 for (label face=0; face<inFaces; face++)
374 label uCell = uPtr[face] + offset;
375 label lCell = lPtr[face] + neiOffset;
377 operator[](uCell)[lCell] -= lowerPtr[face];
378 operator[](lCell)[uCell] -= upperPtr[face];
386void Foam::LUscalarMatrix::printDiagonalDominance()
const
388 for (label i=0; i<m(); i++)
391 for (label j=0; j<m(); j++)
395 sum += operator[](i)[j];
406 pivotIndices_.setSize(m());
415 for (label j=0; j<m(); j++)
420 for (label i=0; i<m(); i++)
A field of fields is a PtrList of fields with reference counting.
Input inter-processor communications stream.
Class to perform the LU decomposition on a symmetric matrix.
LUscalarMatrix()
Construct null.
void setSize(const label n)
Alias for resize()
label n() const noexcept
The number of columns.
void transfer(Matrix< Form, Type > &mat)
const Type * operator[](const label irow) const
Return const pointer to data in the specified row - rowData().
label m() const noexcept
The number of rows.
Output inter-processor communications stream.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
UPstream::rangeType subProcs() const noexcept
Range of sub-processes indices associated with PstreamBuffers.
Inter-processor communications stream.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
const T * set(const label i) const
SquareMatrix & operator=(const SquareMatrix &)=default
Copy assignment.
iterator begin() noexcept
Return an iterator to begin traversing the UList.
void size(const label n)
Older name for setAddressableSize.
static int & msgType() noexcept
Message tag of standard messages.
static constexpr int masterNo() noexcept
Process index of the master (always 0)
static bool & parRun() noexcept
Test if this a parallel run.
const T * set(const label i) const
label size() const noexcept
The number of elements in the list.
iterator begin() noexcept
Return an iterator to begin of UPtrList traversal.
A cell is defined as a list of faces with extra functionality.
gradingDescriptor inv() const
Return the inverse gradingDescriptor with 1/expansionRatio.
virtual const labelUList & upperAddr() const =0
Return upper addressing.
label size() const
Return number of equations.
virtual const labelUList & lowerAddr() const =0
Return lower addressing.
lduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
const lduAddressing & lduAddr() const
Return the LDU addressing.
I/O for lduMatrix and interface values.
splitCell * master() const
bool decompose() const noexcept
Query the decompose flag (normally off)
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
static constexpr const zero Zero
Global zero (0)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a)
interfaceProperties interface(alpha1, U, thermo->transportPropertiesDict())
#define forAll(list, i)
Loop across all elements in list.