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] <<
')';
193 void Foam::LUscalarMatrix::convert
195 const lduMatrix& ldum,
196 const FieldField<Field, scalar>& interfaceCoeffs,
200 const label* __restrict__ uPtr = ldum.lduAddr().upperAddr().begin();
201 const label* __restrict__ lPtr = ldum.lduAddr().lowerAddr().begin();
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();
208 const label nFaces = ldum.upper().size();
210 for (label cell=0; cell<nCells; cell++)
212 operator[](cell)[cell] = diagPtr[cell];
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];
257 void 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];
386 void 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++)
UPtrList< const lduInterfaceField > lduInterfaceFieldPtrsList
List of coupled interface fields to be used in coupling.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Class to perform the LU decomposition on a symmetric matrix.
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
constexpr char nl
The newline '\n' character (0x0a)
label n() const noexcept
The number of columns.
Ostream & endl(Ostream &os)
Add newline and flush stream.
static bool & parRun() noexcept
Test if this a parallel run.
static int & msgType() noexcept
Message tag of standard messages.
interfaceProperties interface(alpha1, U, thermo->transportPropertiesDict())
const scalar * operator[](const label irow) const
Return const pointer to data in the specified row - rowData().
void transfer(Matrix< SquareMatrix< scalar >, scalar > &mat)
Transfer the contents of the argument Matrix into this Matrix and annul the argument Matrix...
A field of fields is a PtrList of fields with reference counting.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
#define forAll(list, i)
Loop across all elements in list.
Input inter-processor communications stream.
void decompose(const scalarSquareMatrix &M)
Perform the LU decomposition of the matrix M.
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator) is 1 for serial run.
void setSize(const label n)
Alias for resize()
"scheduled" : (MPI_Send, MPI_Recv)
Inter-processor communications stream.
static constexpr int masterNo() noexcept
Process index of the master (always 0)
LUscalarMatrix()
Construct null.
label m() const noexcept
The number of rows.
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
lduMatrix is a general matrix class in which the coefficients are stored as three arrays...
const lduAddressing & lduAddr() const
Return the LDU addressing.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers...
static bool master(const label communicator=worldComm)
Am I the master rank.
messageStream Info
Information stream (stdout output on master, null elsewhere)
SquareMatrix & operator=(const SquareMatrix &)=default
Copy assignment.
static rangeType subProcs(const label communicator=worldComm)
Range of process indices for sub-processes.
I/O for lduMatrix and interface values.
SquareMatrix< scalar > scalarSquareMatrix
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution in the source.
label size() const
Return number of equations.
void inv(scalarSquareMatrix &M) const
Set M to the inverse of this square matrix.
static constexpr const zero Zero
Global zero (0)