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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2017-2019 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "LduMatrix.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 template<class Type, class LUType>
39 :
41 {
42  const Field<LUType>& A_;
43 
44 public:
45 
47  :
48  A_(A)
49  {}
50 
51  virtual ~Amultiplier() = default;
52 
53  virtual void addAmul(Field<Type>& Apsi, const Field<Type>& psi) const
54  {
55  Apsi += A_*psi;
56  }
57 };
58 
59 }
60 
61 
62 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
63 
64 template<class Type, class DType, class LUType>
66 (
67  Field<Type>& Apsi,
68  const tmp<Field<Type>>& tpsi
69 ) const
70 {
71  Type* __restrict__ ApsiPtr = Apsi.begin();
72 
73  const Field<Type>& psi = tpsi();
74  const Type* const __restrict__ psiPtr = psi.begin();
75 
76  const DType* const __restrict__ diagPtr = diag().begin();
77 
78  const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
79  const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
80 
81  const LUType* const __restrict__ upperPtr = upper().begin();
82  const LUType* const __restrict__ lowerPtr = lower().begin();
83 
84  // Initialise the update of interfaced interfaces
85  initMatrixInterfaces
86  (
87  true,
88  interfacesUpper_,
89  psi,
90  Apsi
91  );
92 
93  const label nCells = diag().size();
94  for (label cell=0; cell<nCells; cell++)
95  {
96  ApsiPtr[cell] = dot(diagPtr[cell], psiPtr[cell]);
97  }
98 
99 
100  const label nFaces = upper().size();
101  for (label face=0; face<nFaces; face++)
102  {
103  ApsiPtr[uPtr[face]] += dot(lowerPtr[face], psiPtr[lPtr[face]]);
104  ApsiPtr[lPtr[face]] += dot(upperPtr[face], psiPtr[uPtr[face]]);
105  }
106 
107  // Update interface interfaces
108  updateMatrixInterfaces
109  (
110  true,
111  interfacesUpper_,
112  psi,
113  Apsi
114  );
115 
116  tpsi.clear();
117 }
118 
119 
120 template<class Type, class DType, class LUType>
122 (
123  Field<Type>& Tpsi,
124  const tmp<Field<Type>>& tpsi
125 ) const
126 {
127  Type* __restrict__ TpsiPtr = Tpsi.begin();
128 
129  const Field<Type>& psi = tpsi();
130  const Type* const __restrict__ psiPtr = psi.begin();
131 
132  const DType* const __restrict__ diagPtr = diag().begin();
133 
134  const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
135  const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
136 
137  const LUType* const __restrict__ lowerPtr = lower().begin();
138  const LUType* const __restrict__ upperPtr = upper().begin();
139 
140  // Initialise the update of interfaced interfaces
141  initMatrixInterfaces
142  (
143  true,
144  interfacesLower_,
145  psi,
146  Tpsi
147  );
148 
149  const label nCells = diag().size();
150  for (label cell=0; cell<nCells; cell++)
151  {
152  TpsiPtr[cell] = dot(diagPtr[cell], psiPtr[cell]);
153  }
154 
155  const label nFaces = upper().size();
156  for (label face=0; face<nFaces; face++)
157  {
158  TpsiPtr[uPtr[face]] += dot(upperPtr[face], psiPtr[lPtr[face]]);
159  TpsiPtr[lPtr[face]] += dot(lowerPtr[face], psiPtr[uPtr[face]]);
160  }
161 
162  // Update interface interfaces
163  updateMatrixInterfaces
164  (
165  true,
166  interfacesLower_,
167  psi,
168  Tpsi
169  );
170 
171  tpsi.clear();
172 }
173 
174 
175 template<class Type, class DType, class LUType>
177 (
178  Field<Type>& sumA
179 ) const
180 {
181  Type* __restrict__ sumAPtr = sumA.begin();
182 
183  const DType* __restrict__ diagPtr = diag().begin();
184 
185  const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
186  const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
187 
188  const LUType* __restrict__ lowerPtr = lower().begin();
189  const LUType* __restrict__ upperPtr = upper().begin();
190 
191  const label nCells = diag().size();
192  const label nFaces = upper().size();
193 
194  for (label cell=0; cell<nCells; cell++)
195  {
196  sumAPtr[cell] = dot(diagPtr[cell], pTraits<Type>::one);
197  }
198 
199  for (label face=0; face<nFaces; face++)
200  {
201  sumAPtr[uPtr[face]] += dot(lowerPtr[face], pTraits<Type>::one);
202  sumAPtr[lPtr[face]] += dot(upperPtr[face], pTraits<Type>::one);
203  }
204 
205  // Add the interface internal coefficients to diagonal
206  // and the interface boundary coefficients to the sum-off-diagonal
207  forAll(interfaces_, patchi)
208  {
209  if (interfaces_.set(patchi))
210  {
211  const labelUList& pa = lduAddr().patchAddr(patchi);
212  const Field<LUType>& pCoeffs = interfacesUpper_[patchi];
213 
214  forAll(pa, face)
215  {
216  sumAPtr[pa[face]] -= dot(pCoeffs[face], pTraits<Type>::one);
217  }
218  }
219  }
220 }
221 
222 
223 template<class Type, class DType, class LUType>
225 (
226  Field<Type>& rA,
227  const Field<Type>& psi
228 ) const
229 {
230  Type* __restrict__ rAPtr = rA.begin();
231 
232  const Type* const __restrict__ psiPtr = psi.begin();
233  const DType* const __restrict__ diagPtr = diag().begin();
234  const Type* const __restrict__ sourcePtr = source().begin();
235 
236  const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
237  const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
238 
239  const LUType* const __restrict__ upperPtr = upper().begin();
240  const LUType* const __restrict__ lowerPtr = lower().begin();
241 
242  // Parallel boundary initialisation.
243  // Note: there is a change of sign in the coupled
244  // interface update to add the contibution to the r.h.s.
245 
246  // Initialise the update of interfaced interfaces
247  initMatrixInterfaces
248  (
249  false, // negate interface contributions
250  interfacesUpper_,
251  psi,
252  rA
253  );
254 
255  const label nCells = diag().size();
256  for (label cell=0; cell<nCells; cell++)
257  {
258  rAPtr[cell] = sourcePtr[cell] - dot(diagPtr[cell], psiPtr[cell]);
259  }
260 
261 
262  const label nFaces = upper().size();
263  for (label face=0; face<nFaces; face++)
264  {
265  rAPtr[uPtr[face]] -= dot(lowerPtr[face], psiPtr[lPtr[face]]);
266  rAPtr[lPtr[face]] -= dot(upperPtr[face], psiPtr[uPtr[face]]);
267  }
268 
269  // Update interface interfaces
270  updateMatrixInterfaces
271  (
272  false,
273  interfacesUpper_,
274  psi,
275  rA
276  );
277 }
278 
279 
280 template<class Type, class DType, class LUType>
282 (
283  const Field<Type>& psi
284 ) const
285 {
286  tmp<Field<Type>> trA(new Field<Type>(psi.size()));
287  residual(trA.ref(), psi);
288  return trA;
289 }
290 
291 
292 // ************************************************************************* //
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::diag
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
Definition: pointPatchFieldFunctions.H:287
Foam::LduMatrix::sumA
void sumA(Field< Type > &) const
Sum the coefficients on each row of the matrix.
Definition: LduMatrixATmul.C:177
Foam::dot
void dot(FieldField< Field1, typename innerProduct< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
Definition: FieldFieldFunctions.C:944
A
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Foam::LduMatrix::Amul
void Amul(Field< Type > &, const tmp< Field< Type >> &) const
Matrix multiplication.
Definition: LduMatrixATmul.C:66
LduMatrix.H
Foam::Amultiplier::addAmul
virtual void addAmul(Field< Type > &Apsi, const Field< Type > &psi) const
Definition: LduMatrixATmul.C:53
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::LduMatrix::Tmul
void Tmul(Field< Type > &, const tmp< Field< Type >> &) const
Matrix transpose multiplication.
Definition: LduMatrixATmul.C:122
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
Foam::stringOps::lower
string lower(const std::string &s)
Return string copy transformed with std::tolower on each character.
Definition: stringOps.C:1184
Foam::Field< LUType >
Foam::LduMatrix::residual
void residual(Field< Type > &rA, const Field< Type > &psi) const
Definition: LduMatrixATmul.C:225
LduInterfaceFieldPtrsList.H
List of coupled interface fields to be used in coupling.
Foam::LduInterfaceField
An abstract base class for implicitly-coupled interface fields e.g. processor and cyclic patch fields...
Definition: LduInterfaceField.H:58
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::Amultiplier
Definition: LduMatrixATmul.C:38
Foam::pTraits
A traits class, which is primarily used for primitives.
Definition: pTraits.H:56
Foam::UList< label >
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:72
Foam::Amultiplier::~Amultiplier
virtual ~Amultiplier()=default
Foam::stringOps::upper
string upper(const std::string &s)
Return string copy transformed with std::toupper on each character.
Definition: stringOps.C:1200
psi
const volScalarField & psi
Definition: createFieldRefs.H:1
Foam::cell
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:54
Foam::Amultiplier::Amultiplier
Amultiplier(const Field< LUType > &A)
Definition: LduMatrixATmul.C:46