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