lduMatrixATmul.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-2016 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 Description
25  Multiply a given vector (second argument) by the matrix or its transpose
26  and return the result in the first argument.
27 
28 \*---------------------------------------------------------------------------*/
29 
30 #include "lduMatrix.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
35 (
36  scalarField& Apsi,
37  const tmp<scalarField>& tpsi,
38  const FieldField<Field, scalar>& interfaceBouCoeffs,
39  const lduInterfaceFieldPtrsList& interfaces,
40  const direction cmpt
41 ) const
42 {
43  scalar* __restrict__ ApsiPtr = Apsi.begin();
44 
45  const scalarField& psi = tpsi();
46  const scalar* const __restrict__ psiPtr = psi.begin();
47 
48  const scalar* const __restrict__ diagPtr = diag().begin();
49 
50  const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
51  const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
52 
53  const scalar* const __restrict__ upperPtr = upper().begin();
54  const scalar* const __restrict__ lowerPtr = lower().begin();
55 
56  // Initialise the update of interfaced interfaces
57  initMatrixInterfaces
58  (
59  interfaceBouCoeffs,
60  interfaces,
61  psi,
62  Apsi,
63  cmpt
64  );
65 
66  const label nCells = diag().size();
67  for (label cell=0; cell<nCells; cell++)
68  {
69  ApsiPtr[cell] = diagPtr[cell]*psiPtr[cell];
70  }
71 
72 
73  const label nFaces = upper().size();
74 
75  for (label face=0; face<nFaces; face++)
76  {
77  ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];
78  ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];
79  }
80 
81  // Update interface interfaces
82  updateMatrixInterfaces
83  (
84  interfaceBouCoeffs,
85  interfaces,
86  psi,
87  Apsi,
88  cmpt
89  );
90 
91  tpsi.clear();
92 }
93 
94 
96 (
97  scalarField& Tpsi,
98  const tmp<scalarField>& tpsi,
99  const FieldField<Field, scalar>& interfaceIntCoeffs,
100  const lduInterfaceFieldPtrsList& interfaces,
101  const direction cmpt
102 ) const
103 {
104  scalar* __restrict__ TpsiPtr = Tpsi.begin();
105 
106  const scalarField& psi = tpsi();
107  const scalar* const __restrict__ psiPtr = psi.begin();
108 
109  const scalar* const __restrict__ diagPtr = diag().begin();
110 
111  const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
112  const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
113 
114  const scalar* const __restrict__ lowerPtr = lower().begin();
115  const scalar* const __restrict__ upperPtr = upper().begin();
116 
117  // Initialise the update of interfaced interfaces
118  initMatrixInterfaces
119  (
120  interfaceIntCoeffs,
121  interfaces,
122  psi,
123  Tpsi,
124  cmpt
125  );
126 
127  const label nCells = diag().size();
128  for (label cell=0; cell<nCells; cell++)
129  {
130  TpsiPtr[cell] = diagPtr[cell]*psiPtr[cell];
131  }
132 
133  const label nFaces = upper().size();
134  for (label face=0; face<nFaces; face++)
135  {
136  TpsiPtr[uPtr[face]] += upperPtr[face]*psiPtr[lPtr[face]];
137  TpsiPtr[lPtr[face]] += lowerPtr[face]*psiPtr[uPtr[face]];
138  }
139 
140  // Update interface interfaces
141  updateMatrixInterfaces
142  (
143  interfaceIntCoeffs,
144  interfaces,
145  psi,
146  Tpsi,
147  cmpt
148  );
149 
150  tpsi.clear();
151 }
152 
153 
155 (
156  scalarField& sumA,
157  const FieldField<Field, scalar>& interfaceBouCoeffs,
158  const lduInterfaceFieldPtrsList& interfaces
159 ) const
160 {
161  scalar* __restrict__ sumAPtr = sumA.begin();
162 
163  const scalar* __restrict__ diagPtr = diag().begin();
164 
165  const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
166  const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
167 
168  const scalar* __restrict__ lowerPtr = lower().begin();
169  const scalar* __restrict__ upperPtr = upper().begin();
170 
171  const label nCells = diag().size();
172  const label nFaces = upper().size();
173 
174  for (label cell=0; cell<nCells; cell++)
175  {
176  sumAPtr[cell] = diagPtr[cell];
177  }
178 
179  for (label face=0; face<nFaces; face++)
180  {
181  sumAPtr[uPtr[face]] += lowerPtr[face];
182  sumAPtr[lPtr[face]] += upperPtr[face];
183  }
184 
185  // Add the interface internal coefficients to diagonal
186  // and the interface boundary coefficients to the sum-off-diagonal
187  forAll(interfaces, patchi)
188  {
189  if (interfaces.set(patchi))
190  {
191  const labelUList& pa = lduAddr().patchAddr(patchi);
192  const scalarField& pCoeffs = interfaceBouCoeffs[patchi];
193 
194  forAll(pa, face)
195  {
196  sumAPtr[pa[face]] -= pCoeffs[face];
197  }
198  }
199  }
200 }
201 
202 
204 (
205  scalarField& rA,
206  const scalarField& psi,
207  const scalarField& source,
208  const FieldField<Field, scalar>& interfaceBouCoeffs,
209  const lduInterfaceFieldPtrsList& interfaces,
210  const direction cmpt
211 ) const
212 {
213  scalar* __restrict__ rAPtr = rA.begin();
214 
215  const scalar* const __restrict__ psiPtr = psi.begin();
216  const scalar* const __restrict__ diagPtr = diag().begin();
217  const scalar* const __restrict__ sourcePtr = source.begin();
218 
219  const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
220  const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
221 
222  const scalar* const __restrict__ upperPtr = upper().begin();
223  const scalar* const __restrict__ lowerPtr = lower().begin();
224 
225  // Parallel boundary initialisation.
226  // Note: there is a change of sign in the coupled
227  // interface update. The reason for this is that the
228  // internal coefficients are all located at the l.h.s. of
229  // the matrix whereas the "implicit" coefficients on the
230  // coupled boundaries are all created as if the
231  // coefficient contribution is of a source-kind (i.e. they
232  // have a sign as if they are on the r.h.s. of the matrix.
233  // To compensate for this, it is necessary to turn the
234  // sign of the contribution.
235 
236  FieldField<Field, scalar> mBouCoeffs(interfaceBouCoeffs.size());
237 
238  forAll(mBouCoeffs, patchi)
239  {
240  if (interfaces.set(patchi))
241  {
242  mBouCoeffs.set(patchi, -interfaceBouCoeffs[patchi]);
243  }
244  }
245 
246  // Initialise the update of interfaced interfaces
247  initMatrixInterfaces
248  (
249  mBouCoeffs,
250  interfaces,
251  psi,
252  rA,
253  cmpt
254  );
255 
256  const label nCells = diag().size();
257  for (label cell=0; cell<nCells; cell++)
258  {
259  rAPtr[cell] = sourcePtr[cell] - diagPtr[cell]*psiPtr[cell];
260  }
261 
262 
263  const label nFaces = upper().size();
264 
265  for (label face=0; face<nFaces; face++)
266  {
267  rAPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
268  rAPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
269  }
270 
271  // Update interface interfaces
272  updateMatrixInterfaces
273  (
274  mBouCoeffs,
275  interfaces,
276  psi,
277  rA,
278  cmpt
279  );
280 }
281 
282 
284 (
285  const scalarField& psi,
286  const scalarField& source,
287  const FieldField<Field, scalar>& interfaceBouCoeffs,
288  const lduInterfaceFieldPtrsList& interfaces,
289  const direction cmpt
290 ) const
291 {
292  tmp<scalarField> trA(new scalarField(psi.size()));
293  residual(trA.ref(), psi, source, interfaceBouCoeffs, interfaces, cmpt);
294  return trA;
295 }
296 
297 
299 {
301  (
302  new scalarField(lduAddr().size(), 0.0)
303  );
304 
305  if (lowerPtr_ || upperPtr_)
306  {
307  scalarField& H1_ = tH1.ref();
308 
309  scalar* __restrict__ H1Ptr = H1_.begin();
310 
311  const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
312  const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
313 
314  const scalar* __restrict__ lowerPtr = lower().begin();
315  const scalar* __restrict__ upperPtr = upper().begin();
316 
317  const label nFaces = upper().size();
318 
319  for (label face=0; face<nFaces; face++)
320  {
321  H1Ptr[uPtr[face]] -= lowerPtr[face];
322  H1Ptr[lPtr[face]] -= upperPtr[face];
323  }
324  }
325 
326  return tH1;
327 }
328 
329 
330 // ************************************************************************* //
#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
uint8_t direction
Definition: direction.H:46
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
tmp< scalarField > H1() const
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
void Tmul(scalarField &, const tmp< scalarField > &, const FieldField< Field, scalar > &, const lduInterfaceFieldPtrsList &, const direction cmpt) const
Matrix transpose multiplication with updated interfaces.
scalarField & upper()
Definition: lduMatrix.C:194
virtual const labelUList & lowerAddr() const =0
Return lower addressing.
Generic field type.
Definition: FieldField.H:51
bool set(const label) const
Is element set.
Definition: UPtrListI.H:78
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:230
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
void residual(scalarField &rA, const scalarField &psi, const scalarField &source, const FieldField< Field, scalar > &interfaceBouCoeffs, const lduInterfaceFieldPtrsList &interfaces, const direction cmpt) const
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual const labelUList & upperAddr() const =0
Return upper addressing.
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
void Amul(scalarField &, const tmp< scalarField > &, const FieldField< Field, scalar > &, const lduInterfaceFieldPtrsList &, const direction cmpt) const
Matrix multiplication with updated interfaces.
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
volScalarField scalarField(fieldObject, mesh)
const lduAddressing & lduAddr() const
Return the LDU addressing.
Definition: lduMatrix.H:547
label patchi
scalarField & lower()
Definition: lduMatrix.C:165
void sumA(scalarField &, const FieldField< Field, scalar > &, const lduInterfaceFieldPtrsList &) const
Sum the coefficients on each row of the matrix.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
const volScalarField & psi
A class for managing temporary objects.
Definition: PtrList.H:54
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29