FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/lduMatrix.H>
28 #include "procLduMatrix.H"
29 #include "procLduInterface.H"
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
34 :
35  scalarSquareMatrix(matrix),
36  pivotIndices_(n())
37 {
38  LUDecompose(*this, pivotIndices_);
39 }
40 
41 
43 (
44  const lduMatrix& ldum,
45  const FieldField<Field, scalar>& interfaceCoeffs,
46  const lduInterfaceFieldPtrsList& interfaces
47 )
48 {
49  if (Pstream::parRun())
50  {
52 
53  label lduMatrixi = 0;
54 
55  lduMatrices.set
56  (
57  lduMatrixi++,
58  new procLduMatrix
59  (
60  ldum,
61  interfaceCoeffs,
62  interfaces
63  )
64  );
65 
66  if (Pstream::master())
67  {
68  for
69  (
70  int slave=Pstream::firstSlave();
71  slave<=Pstream::lastSlave();
72  slave++
73  )
74  {
75  lduMatrices.set
76  (
77  lduMatrixi++,
79  );
80  }
81  }
82  else
83  {
85  procLduMatrix cldum
86  (
87  ldum,
88  interfaceCoeffs,
89  interfaces
90  );
91  toMaster<< cldum;
92 
93  }
94 
95  if (Pstream::master())
96  {
97  label nCells = 0;
98  forAll(lduMatrices, i)
99  {
100  nCells += lduMatrices[i].size();
101  }
102 
103  scalarSquareMatrix m(nCells, nCells, 0.0);
104  transfer(m);
105  convert(lduMatrices);
106  }
107  }
108  else
109  {
110  label nCells = ldum.lduAddr().size();
111  scalarSquareMatrix m(nCells, nCells, 0.0);
112  transfer(m);
113  convert(ldum, interfaceCoeffs, interfaces);
114  }
115 
116  if (Pstream::master())
117  {
118  pivotIndices_.setSize(n());
119  LUDecompose(*this, pivotIndices_);
120  }
121 }
122 
123 
124 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
125 
126 void Foam::LUscalarMatrix::convert
127 (
128  const lduMatrix& ldum,
129  const FieldField<Field, scalar>& interfaceCoeffs,
130  const lduInterfaceFieldPtrsList& interfaces
131 )
132 {
133  const label* __restrict__ uPtr = ldum.lduAddr().upperAddr().begin();
134  const label* __restrict__ lPtr = ldum.lduAddr().lowerAddr().begin();
135 
136  const scalar* __restrict__ diagPtr = ldum.diag().begin();
137  const scalar* __restrict__ upperPtr = ldum.upper().begin();
138  const scalar* __restrict__ lowerPtr = ldum.lower().begin();
139 
140  register const label nCells = ldum.diag().size();
141  register const label nFaces = ldum.upper().size();
142 
143  for (register label cell=0; cell<nCells; cell++)
144  {
145  operator[](cell)[cell] = diagPtr[cell];
146  }
147 
148  for (register label face=0; face<nFaces; face++)
149  {
150  label uCell = uPtr[face];
151  label lCell = lPtr[face];
152 
153  operator[](uCell)[lCell] = lowerPtr[face];
154  operator[](lCell)[uCell] = upperPtr[face];
155  }
156 
157  forAll(interfaces, inti)
158  {
159  if (interfaces.set(inti))
160  {
161  const lduInterface& interface = interfaces[inti].interface();
162 
163  const label* __restrict__ ulPtr = interface.faceCells().begin();
164  const scalar* __restrict__ upperLowerPtr =
165  interfaceCoeffs[inti].begin();
166 
167  register label inFaces = interface.faceCells().size()/2;
168 
169  for (register label face=0; face<inFaces; face++)
170  {
171  label uCell = ulPtr[face];
172  label lCell = ulPtr[face + inFaces];
173 
174  operator[](uCell)[lCell] -= upperLowerPtr[face + inFaces];
175  operator[](lCell)[uCell] -= upperLowerPtr[face];
176  }
177  }
178  }
179 
180  //printDiagonalDominance();
181 }
182 
183 
184 void Foam::LUscalarMatrix::convert
185 (
186  const PtrList<procLduMatrix>& lduMatrices
187 )
188 {
189  procOffsets_.setSize(lduMatrices.size() + 1);
190  procOffsets_[0] = 0;
191 
192  forAll(lduMatrices, ldumi)
193  {
194  procOffsets_[ldumi+1] = procOffsets_[ldumi] + lduMatrices[ldumi].size();
195  }
196 
197  forAll(lduMatrices, ldumi)
198  {
199  const procLduMatrix& lduMatrixi = lduMatrices[ldumi];
200  label offset = procOffsets_[ldumi];
201 
202  const label* __restrict__ uPtr = lduMatrixi.upperAddr_.begin();
203  const label* __restrict__ lPtr = lduMatrixi.lowerAddr_.begin();
204 
205  const scalar* __restrict__ diagPtr = lduMatrixi.diag_.begin();
206  const scalar* __restrict__ upperPtr = lduMatrixi.upper_.begin();
207  const scalar* __restrict__ lowerPtr = lduMatrixi.lower_.begin();
208 
209  register const label nCells = lduMatrixi.size();
210  register const label nFaces = lduMatrixi.upper_.size();
211 
212  for (register label cell=0; cell<nCells; cell++)
213  {
214  label globalCell = cell + offset;
215  operator[](globalCell)[globalCell] = diagPtr[cell];
216  }
217 
218  for (register label face=0; face<nFaces; face++)
219  {
220  label uCell = uPtr[face] + offset;
221  label lCell = lPtr[face] + offset;
222 
223  operator[](uCell)[lCell] = lowerPtr[face];
224  operator[](lCell)[uCell] = upperPtr[face];
225  }
226 
227  const PtrList<procLduInterface>& interfaces =
228  lduMatrixi.interfaces_;
229 
230  forAll(interfaces, inti)
231  {
232  const procLduInterface& interface = interfaces[inti];
233 
234  if (interface.myProcNo_ == interface.neighbProcNo_)
235  {
236  const label* __restrict__ ulPtr = interface.faceCells_.begin();
237 
238  const scalar* __restrict__ upperLowerPtr =
239  interface.coeffs_.begin();
240 
241  register label inFaces = interface.faceCells_.size()/2;
242 
243  for (register label face=0; face<inFaces; face++)
244  {
245  label uCell = ulPtr[face] + offset;
246  label lCell = ulPtr[face + inFaces] + offset;
247 
248  operator[](uCell)[lCell] -= upperLowerPtr[face + inFaces];
249  operator[](lCell)[uCell] -= upperLowerPtr[face];
250  }
251  }
252  else if (interface.myProcNo_ < interface.neighbProcNo_)
253  {
254  const PtrList<procLduInterface>& neiInterfaces =
255  lduMatrices[interface.neighbProcNo_].interfaces_;
256 
257  label neiInterfacei = -1;
258 
259  forAll(neiInterfaces, ninti)
260  {
261  if
262  (
263  neiInterfaces[ninti].neighbProcNo_
264  == interface.myProcNo_
265  )
266  {
267  neiInterfacei = ninti;
268  break;
269  }
270  }
271 
272  if (neiInterfacei == -1)
273  {
274  FatalErrorIn("LUscalarMatrix::convert") << exit(FatalError);
275  }
276 
277  const procLduInterface& neiInterface =
278  neiInterfaces[neiInterfacei];
279 
280  const label* __restrict__ uPtr = interface.faceCells_.begin();
281  const label* __restrict__ lPtr = neiInterface.faceCells_.begin();
282 
283  const scalar* __restrict__ upperPtr = interface.coeffs_.begin();
284  const scalar* __restrict__ lowerPtr = neiInterface.coeffs_.begin();
285 
286  register label inFaces = interface.faceCells_.size();
287  label neiOffset = procOffsets_[interface.neighbProcNo_];
288 
289  for (register label face=0; face<inFaces; face++)
290  {
291  label uCell = uPtr[face] + offset;
292  label lCell = lPtr[face] + neiOffset;
293 
294  operator[](uCell)[lCell] -= lowerPtr[face];
295  operator[](lCell)[uCell] -= upperPtr[face];
296  }
297  }
298  }
299  }
300 
301  //printDiagonalDominance();
302 }
303 
304 
305 void Foam::LUscalarMatrix::printDiagonalDominance() const
306 {
307  for (label i=0; i<n(); i++)
308  {
309  scalar sum = 0.0;
310  for (label j=0; j<n(); j++)
311  {
312  if (i != j)
313  {
314  sum += operator[](i)[j];
315  }
316  }
317  Info<< mag(sum)/mag(operator[](i)[i]) << endl;
318  }
319 }
320 
321 
322 // ************************ vim: set sw=4 sts=4 et: ************************ //