LUscalarMatrix.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "LUscalarMatrix.H"
27 #include "lduMatrix.H"
28 #include "procLduMatrix.H"
29 #include "procLduInterface.H"
30 #include "cyclicLduInterface.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
37 }
38 
39 
40 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
41 
43 :
44  comm_(Pstream::worldComm)
45 {}
46 
47 
49 :
51  comm_(Pstream::worldComm),
52  pivotIndices_(n)
53 {}
54 
55 
57 :
58  scalarSquareMatrix(matrix),
59  comm_(Pstream::worldComm),
60  pivotIndices_(m())
61 {
62  LUDecompose(*this, pivotIndices_);
63 }
64 
65 
67 (
68  const lduMatrix& ldum,
69  const FieldField<Field, scalar>& interfaceCoeffs,
70  const lduInterfaceFieldPtrsList& interfaces
71 )
72 :
73  comm_(ldum.mesh().comm())
74 {
75  if (Pstream::parRun())
76  {
77  PtrList<procLduMatrix> lduMatrices(Pstream::nProcs(comm_));
78 
79  label lduMatrixi = 0;
80 
81  lduMatrices.set
82  (
83  lduMatrixi++,
84  new procLduMatrix
85  (
86  ldum,
87  interfaceCoeffs,
88  interfaces
89  )
90  );
91 
92  if (Pstream::master(comm_))
93  {
94  for
95  (
96  int slave=Pstream::firstSlave();
97  slave<=Pstream::lastSlave(comm_);
98  slave++
99  )
100  {
101  lduMatrices.set
102  (
103  lduMatrixi++,
104  new procLduMatrix
105  (
106  IPstream
107  (
109  slave,
110  0, // bufSize
112  comm_
113  )()
114  )
115  );
116  }
117  }
118  else
119  {
120  OPstream toMaster
121  (
124  0, // bufSize
126  comm_
127  );
128  procLduMatrix cldum
129  (
130  ldum,
131  interfaceCoeffs,
132  interfaces
133  );
134  toMaster<< cldum;
135 
136  }
137 
138  if (Pstream::master(comm_))
139  {
140  label nCells = 0;
141  forAll(lduMatrices, i)
142  {
143  nCells += lduMatrices[i].size();
144  }
145 
146  scalarSquareMatrix m(nCells, 0.0);
147  transfer(m);
148  convert(lduMatrices);
149  }
150  }
151  else
152  {
153  label nCells = ldum.lduAddr().size();
154  scalarSquareMatrix m(nCells, 0.0);
155  transfer(m);
156  convert(ldum, interfaceCoeffs, interfaces);
157  }
158 
159  if (Pstream::master(comm_))
160  {
161  label mRows = m();
162  label nColumns = n();
163 
164  if (debug)
165  {
166  Pout<< "LUscalarMatrix : size:" << mRows << endl;
167  for (label rowI = 0; rowI < mRows; rowI++)
168  {
169  const scalar* row = operator[](rowI);
170 
171  Pout<< "cell:" << rowI << " diagCoeff:" << row[rowI] << endl;
172 
173  Pout<< " connects to upper cells :";
174  for (label columnI = rowI+1; columnI < nColumns; columnI++)
175  {
176  if (mag(row[columnI]) > small)
177  {
178  Pout<< ' ' << columnI << " (coeff:" << row[columnI]
179  << ")";
180  }
181  }
182  Pout<< endl;
183  Pout<< " connects to lower cells :";
184  for (label columnI = 0; columnI < rowI; columnI++)
185  {
186  if (mag(row[columnI]) > small)
187  {
188  Pout<< ' ' << columnI << " (coeff:" << row[columnI]
189  << ")";
190  }
191  }
192  Pout<< endl;
193  }
194  Pout<< endl;
195  }
196 
197  pivotIndices_.setSize(m());
198  LUDecompose(*this, pivotIndices_);
199  }
200 }
201 
202 
203 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
204 
205 void Foam::LUscalarMatrix::convert
206 (
207  const lduMatrix& ldum,
208  const FieldField<Field, scalar>& interfaceCoeffs,
209  const lduInterfaceFieldPtrsList& interfaces
210 )
211 {
212  const label* __restrict__ uPtr = ldum.lduAddr().upperAddr().begin();
213  const label* __restrict__ lPtr = ldum.lduAddr().lowerAddr().begin();
214 
215  const scalar* __restrict__ diagPtr = ldum.diag().begin();
216  const scalar* __restrict__ upperPtr = ldum.upper().begin();
217  const scalar* __restrict__ lowerPtr = ldum.lower().begin();
218 
219  const label nCells = ldum.diag().size();
220  const label nFaces = ldum.upper().size();
221 
222  for (label cell=0; cell<nCells; cell++)
223  {
224  operator[](cell)[cell] = diagPtr[cell];
225  }
226 
227  for (label face=0; face<nFaces; face++)
228  {
229  label uCell = uPtr[face];
230  label lCell = lPtr[face];
231 
232  operator[](uCell)[lCell] = lowerPtr[face];
233  operator[](lCell)[uCell] = upperPtr[face];
234  }
235 
236  forAll(interfaces, inti)
237  {
238  if (interfaces.set(inti))
239  {
240  const lduInterface& interface = interfaces[inti].interface();
241 
242  // Assume any interfaces are cyclic ones
243 
244  const label* __restrict__ lPtr = interface.faceCells().begin();
245 
246  const cyclicLduInterface& cycInterface =
247  refCast<const cyclicLduInterface>(interface);
248  label nbrInt = cycInterface.nbrPatchIndex();
249  const label* __restrict__ uPtr =
250  interfaces[nbrInt].interface().faceCells().begin();
251 
252  const scalar* __restrict__ nbrUpperLowerPtr =
253  interfaceCoeffs[nbrInt].begin();
254 
255  label inFaces = interface.faceCells().size();
256 
257  for (label face=0; face<inFaces; face++)
258  {
259  label uCell = lPtr[face];
260  label lCell = uPtr[face];
261 
262  operator[](uCell)[lCell] -= nbrUpperLowerPtr[face];
263  }
264  }
265  }
266 }
267 
268 
269 void Foam::LUscalarMatrix::convert
270 (
271  const PtrList<procLduMatrix>& lduMatrices
272 )
273 {
274  procOffsets_.setSize(lduMatrices.size() + 1);
275  procOffsets_[0] = 0;
276 
277  forAll(lduMatrices, ldumi)
278  {
279  procOffsets_[ldumi+1] = procOffsets_[ldumi] + lduMatrices[ldumi].size();
280  }
281 
282  forAll(lduMatrices, ldumi)
283  {
284  const procLduMatrix& lduMatrixi = lduMatrices[ldumi];
285  label offset = procOffsets_[ldumi];
286 
287  const label* __restrict__ uPtr = lduMatrixi.upperAddr_.begin();
288  const label* __restrict__ lPtr = lduMatrixi.lowerAddr_.begin();
289 
290  const scalar* __restrict__ diagPtr = lduMatrixi.diag_.begin();
291  const scalar* __restrict__ upperPtr = lduMatrixi.upper_.begin();
292  const scalar* __restrict__ lowerPtr = lduMatrixi.lower_.begin();
293 
294  const label nCells = lduMatrixi.size();
295  const label nFaces = lduMatrixi.upper_.size();
296 
297  for (label cell=0; cell<nCells; cell++)
298  {
299  label globalCell = cell + offset;
300  operator[](globalCell)[globalCell] = diagPtr[cell];
301  }
302 
303  for (label face=0; face<nFaces; face++)
304  {
305  label uCell = uPtr[face] + offset;
306  label lCell = lPtr[face] + offset;
307 
308  operator[](uCell)[lCell] = lowerPtr[face];
309  operator[](lCell)[uCell] = upperPtr[face];
310  }
311 
312  const PtrList<procLduInterface>& interfaces =
313  lduMatrixi.interfaces_;
314 
315  forAll(interfaces, inti)
316  {
317  const procLduInterface& interface = interfaces[inti];
318 
319  if (interface.myProcNo_ == interface.neighbProcNo_)
320  {
321  const label* __restrict__ ulPtr = interface.faceCells_.begin();
322 
323  const scalar* __restrict__ upperLowerPtr =
324  interface.coeffs_.begin();
325 
326  label inFaces = interface.faceCells_.size()/2;
327 
328  for (label face=0; face<inFaces; face++)
329  {
330  label uCell = ulPtr[face] + offset;
331  label lCell = ulPtr[face + inFaces] + offset;
332 
333  operator[](uCell)[lCell] -= upperLowerPtr[face + inFaces];
334  operator[](lCell)[uCell] -= upperLowerPtr[face];
335  }
336  }
337  else if (interface.myProcNo_ < interface.neighbProcNo_)
338  {
339  // Interface to neighbour proc. Find on neighbour proc the
340  // corresponding interface. The problem is that there can
341  // be multiple interfaces between two processors (from
342  // processorCyclics) so also compare the communication tag
343 
344  const PtrList<procLduInterface>& neiInterfaces =
345  lduMatrices[interface.neighbProcNo_].interfaces_;
346 
347  label neiInterfacei = -1;
348 
349  forAll(neiInterfaces, ninti)
350  {
351  if
352  (
353  (
354  neiInterfaces[ninti].neighbProcNo_
355  == interface.myProcNo_
356  )
357  && (neiInterfaces[ninti].tag_ == interface.tag_)
358  )
359  {
360  neiInterfacei = ninti;
361  break;
362  }
363  }
364 
365  if (neiInterfacei == -1)
366  {
368  }
369 
370  const procLduInterface& neiInterface =
371  neiInterfaces[neiInterfacei];
372 
373  const label* __restrict__ uPtr = interface.faceCells_.begin();
374  const label* __restrict__ lPtr =
375  neiInterface.faceCells_.begin();
376 
377  const scalar* __restrict__ upperPtr = interface.coeffs_.begin();
378  const scalar* __restrict__ lowerPtr =
379  neiInterface.coeffs_.begin();
380 
381  label inFaces = interface.faceCells_.size();
382  label neiOffset = procOffsets_[interface.neighbProcNo_];
383 
384  for (label face=0; face<inFaces; face++)
385  {
386  label uCell = uPtr[face] + offset;
387  label lCell = lPtr[face] + neiOffset;
388 
389  operator[](uCell)[lCell] -= lowerPtr[face];
390  operator[](lCell)[uCell] -= upperPtr[face];
391  }
392  }
393  }
394  }
395 }
396 
397 
398 void Foam::LUscalarMatrix::printDiagonalDominance() const
399 {
400  for (label i=0; i<m(); i++)
401  {
402  scalar sum = 0.0;
403  for (label j=0; j<m(); j++)
404  {
405  if (i != j)
406  {
407  sum += operator[](i)[j];
408  }
409  }
410  Info<< mag(sum)/mag(operator[](i)[i]) << endl;
411  }
412 }
413 
414 
416 {
417  pivotIndices_.setSize(m());
418  LUDecompose(*this, pivotIndices_);
419 }
420 
421 
423 {
425  pivotIndices_.setSize(m());
426  LUDecompose(*this, pivotIndices_);
427 }
428 
429 
431 {
432  scalarField source(m());
433 
434  for (label j=0; j<m(); j++)
435  {
436  source = Zero;
437  source[j] = 1;
438  LUBacksubstitute(*this, pivotIndices_, source);
439  for (label i=0; i<m(); i++)
440  {
441  M(i, j) = source[i];
442  }
443  }
444 }
445 
446 
447 // ************************************************************************* //
#define M(I)
label n
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Generic field type.
Definition: FieldField.H:77
Input inter-processor communications stream.
Definition: IPstream.H:54
Class to perform the LU decomposition on a symmetric matrix.
LUscalarMatrix()
Construct null.
void inv(scalarSquareMatrix &M) const
Set M to the inverse of this square matrix.
void decompose()
Perform the LU decomposition of the matrix.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void setSize(const label)
Reset size of List.
Definition: List.C:281
label n() const
Return the number of columns.
Definition: MatrixI.H:64
label m() const
Return the number of rows.
Definition: MatrixI.H:57
void transfer(mType &)
Transfer the contents of the argument Matrix into this Matrix.
Definition: Matrix.C:228
Type * operator[](const label)
Return subscript-checked row of Matrix.
Definition: MatrixI.H:328
Output inter-processor communications stream.
Definition: OPstream.H:54
Inter-processor communications stream.
Definition: Pstream.H:56
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: PtrList.H:75
bool set(const label) const
Is element set.
Definition: PtrListI.H:62
void operator=(const zero)
Assignment of all elements to zero.
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
static int masterNo()
Process index of the master.
Definition: UPstream.H:417
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
static int lastSlave(const label communicator=0)
Process index of last slave.
Definition: UPstream.H:452
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
static int firstSlave()
Process index of first slave.
Definition: UPstream.H:446
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
static int & msgType()
Message tag of standard messages.
Definition: UPstream.H:476
iterator begin()
Return an iterator to begin traversing the UPtrList.
Definition: UPtrListI.H:308
bool set(const label) const
Is element set.
Definition: UPtrListI.H:87
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:60
label size() const
Return number of equations.
virtual const labelUList & upperAddr() const =0
Return upper addressing.
virtual const labelUList & lowerAddr() const =0
Return lower addressing.
lduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: lduMatrix.H:80
scalarField & upper()
Definition: lduMatrix.C:197
scalarField & lower()
Definition: lduMatrix.C:168
scalarField & diag()
Definition: lduMatrix.C:186
const lduAddressing & lduAddr() const
Return the LDU addressing.
Definition: lduMatrix.H:550
I/O for lduMatrix and interface values.
Definition: procLduMatrix.H:62
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static const zero Zero
Definition: zero.H:97
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
error FatalError
void offset(label &lst, const label o)