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