67 comm_(ldum.
mesh().comm())
91 auto& mat = lduMatrices.emplace_set(proci);
111 convert(ldum, interfaceCoeffs, interfaces);
117 const label numRows =
nRows();
118 const label numCols =
nCols();
120 Pout<<
"LUscalarMatrix : size:" << numRows <<
endl;
121 for (label rowi = 0; rowi < numRows; ++rowi)
125 Pout<<
"cell:" << rowi <<
" diagCoeff:" << row[rowi] <<
nl;
127 Pout<<
" connects to upper cells :";
128 for (label coli = rowi+1; coli < numCols; ++coli)
130 if (
mag(row[coli]) > SMALL)
132 Pout<<
' ' << coli <<
" (coeff:" << row[coli] <<
')';
136 Pout<<
" connects to lower cells :";
137 for (label coli = 0; coli < rowi; ++coli)
139 if (
mag(row[coli]) > SMALL)
141 Pout<<
' ' << coli <<
" (coeff:" << row[coli] <<
')';
158 void Foam::LUscalarMatrix::convert
160 const lduMatrix& ldum,
161 const FieldField<Field, scalar>& interfaceCoeffs,
169 const label* __restrict__ uPtr = ldum.lduAddr().upperAddr().
begin();
170 const label* __restrict__ lPtr = ldum.lduAddr().lowerAddr().begin();
172 const scalar* __restrict__ diagPtr = ldum.diag().begin();
173 const scalar* __restrict__ upperPtr = ldum.upper().begin();
174 const scalar* __restrict__ lowerPtr = ldum.lower().begin();
176 const label nCells = ldum.diag().size();
177 const label nFaces = ldum.upper().size();
179 for (label cell=0; cell<nCells; cell++)
181 operator[](cell)[cell] = diagPtr[cell];
184 for (label face=0; face<nFaces; face++)
186 label uCell = uPtr[face];
187 label lCell = lPtr[face];
189 operator[](uCell)[lCell] = lowerPtr[face];
190 operator[](lCell)[uCell] = upperPtr[face];
195 if (interfaces.set(inti))
197 const lduInterface&
interface = interfaces[inti].
interface();
201 const label* __restrict__ lPtr =
interface.faceCells().begin();
203 const cyclicLduInterface& cycInterface =
204 refCast<const cyclicLduInterface>(
interface);
205 label nbrInt = cycInterface.neighbPatchID();
206 const label* __restrict__ uPtr =
207 interfaces[nbrInt].interface().faceCells().begin();
209 const scalar* __restrict__ nbrUpperLowerPtr =
210 interfaceCoeffs[nbrInt].begin();
212 label inFaces =
interface.faceCells().size();
214 for (label face=0; face<inFaces; face++)
216 label uCell = lPtr[face];
217 label lCell = uPtr[face];
219 operator[](uCell)[lCell] -= nbrUpperLowerPtr[face];
226 void Foam::LUscalarMatrix::convert
228 const PtrList<procLduMatrix>& lduMatrices
231 procOffsets_.resize_nocopy(lduMatrices.size() + 1);
234 auto iter = procOffsets_.begin();
236 label nCellsTotal = 0;
237 *iter++ = nCellsTotal;
239 for (
const auto& mat : lduMatrices)
241 nCellsTotal += mat.size();
242 *iter++ = nCellsTotal;
251 forAll(lduMatrices, ldumi)
253 const procLduMatrix& lduMatrixi = lduMatrices[ldumi];
254 label offset = procOffsets_[ldumi];
256 const label* __restrict__ uPtr = lduMatrixi.upperAddr_.begin();
257 const label* __restrict__ lPtr = lduMatrixi.lowerAddr_.begin();
259 const scalar* __restrict__ diagPtr = lduMatrixi.diag_.begin();
260 const scalar* __restrict__ upperPtr = lduMatrixi.upper_.begin();
261 const scalar* __restrict__ lowerPtr = lduMatrixi.lower_.begin();
263 const label nCells = lduMatrixi.size();
264 const label nFaces = lduMatrixi.upper_.size();
266 for (label cell=0; cell<nCells; cell++)
268 label globalCell = cell + offset;
269 operator[](globalCell)[globalCell] = diagPtr[cell];
272 for (label face=0; face<nFaces; face++)
274 label uCell = uPtr[face] + offset;
275 label lCell = lPtr[face] + offset;
277 operator[](uCell)[lCell] = lowerPtr[face];
278 operator[](lCell)[uCell] = upperPtr[face];
281 const PtrList<procLduInterface>& interfaces =
282 lduMatrixi.interfaces_;
286 const procLduInterface&
interface = interfaces[inti];
290 const label* __restrict__ ulPtr =
interface.faceCells_.begin();
292 const scalar* __restrict__ upperLowerPtr =
295 label inFaces =
interface.faceCells_.size()/2;
297 for (label face=0; face<inFaces; face++)
299 label uCell = ulPtr[face] + offset;
300 label lCell = ulPtr[face + inFaces] + offset;
302 operator[](uCell)[lCell] -= upperLowerPtr[face + inFaces];
303 operator[](lCell)[uCell] -= upperLowerPtr[face];
313 const PtrList<procLduInterface>& neiInterfaces =
314 lduMatrices[
interface.neighbProcNo_].interfaces_;
316 label neiInterfacei = -1;
318 forAll(neiInterfaces, ninti)
323 neiInterfaces[ninti].neighbProcNo_
326 && (neiInterfaces[ninti].tag_ ==
interface.tag_)
329 neiInterfacei = ninti;
334 if (neiInterfacei == -1)
339 const procLduInterface& neiInterface =
340 neiInterfaces[neiInterfacei];
342 const label* __restrict__ uPtr =
interface.faceCells_.begin();
343 const label* __restrict__ lPtr =
344 neiInterface.faceCells_.begin();
346 const scalar* __restrict__ upperPtr =
interface.coeffs_.begin();
347 const scalar* __restrict__ lowerPtr =
348 neiInterface.coeffs_.begin();
350 label inFaces =
interface.faceCells_.size();
351 label neiOffset = procOffsets_[
interface.neighbProcNo_];
353 for (label face=0; face<inFaces; face++)
355 label uCell = uPtr[face] + offset;
356 label lCell = lPtr[face] + neiOffset;
358 operator[](uCell)[lCell] -= lowerPtr[face];
359 operator[](lCell)[uCell] -= upperPtr[face];
367 void Foam::LUscalarMatrix::printDiagonalDominance()
const 369 for (label i=0; i<m(); i++)
372 for (label j=0; j<m(); j++)
376 sum += operator[](i)[j];
395 for (label j=0; j<m(); j++)
400 for (label i=0; i<m(); i++)
UPtrList< const lduInterfaceField > lduInterfaceFieldPtrsList
List of coupled interface fields to be used in coupling.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
LUscalarMatrix() noexcept
Default construct.
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)
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().
iterator begin() noexcept
Return an iterator to begin traversing a Matrix.
A field of fields is a PtrList of fields with reference counting.
#define forAll(list, i)
Loop across all elements in list.
label nCols() const noexcept
The number of columns.
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run. ...
static void recv(Type &value, const int fromProcNo, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm, IOstreamOption::streamFormat fmt=IOstreamOption::BINARY)
Receive and deserialize a value. Uses operator>> for de-serialization.
label nRows() const noexcept
The number of rows.
void resize_nocopy(const label n)
Resize the matrix without preserving existing content.
static constexpr int masterNo() noexcept
Relative rank for the master process - is always 0.
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
lduMatrix is a general matrix class in which the coefficients are stored as three arrays...
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers...
void decompose(const scalarSquareMatrix &mat)
Perform the LU decomposition of the matrix.
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
messageStream Info
Information stream (stdout output on master, null elsewhere)
A class representing the concept of 0 (zero) that can be used to avoid manipulating objects known to ...
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.
bool send()
Send buffer contents now and not in destructor [advanced usage]. Returns true on success.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Inter-processor communications stream.
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution in the source.
void inv(scalarSquareMatrix &M) const
Set M to the inverse of this square matrix.
static constexpr const zero Zero
Global zero (0)