nonBlockingGaussSeidelSmoother.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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 
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32  defineTypeNameAndDebug(nonBlockingGaussSeidelSmoother, 0);
33 
34  lduMatrix::smoother::
35  addsymMatrixConstructorToTable<nonBlockingGaussSeidelSmoother>
37 
38  lduMatrix::smoother::
39  addasymMatrixConstructorToTable<nonBlockingGaussSeidelSmoother>
41 }
42 
43 
44 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
45 
47 (
48  const word& fieldName,
49  const lduMatrix& matrix,
50  const FieldField<Field, scalar>& interfaceBouCoeffs,
51  const FieldField<Field, scalar>& interfaceIntCoeffs,
52  const lduInterfaceFieldPtrsList& interfaces
53 )
54 :
56  (
57  fieldName,
58  matrix,
59  interfaceBouCoeffs,
60  interfaceIntCoeffs,
61  interfaces
62  )
63 {
64  // Check that all interface addressing is sorted to be after the
65  // non-interface addressing.
66 
67  const label nCells = matrix.diag().size();
68 
69  blockStart_ = nCells;
70 
71  labelList startCelli(interfaceBouCoeffs.size(), -1);
72  forAll(interfaces, patchi)
73  {
74  if (interfaces.set(patchi))
75  {
76  const labelUList& faceCells = matrix_.lduAddr().patchAddr(patchi);
77 
78  blockStart_ = min(blockStart_, min(faceCells));
79  }
80  }
81 
82  if (debug)
83  {
84  Pout<< "nonBlockingGaussSeidelSmoother :"
85  << " Starting block on cell " << blockStart_
86  << " out of " << nCells << endl;
87  }
88 }
89 
90 
91 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
92 
94 (
95  const word& fieldName_,
96  scalarField& psi,
97  const lduMatrix& matrix_,
98  const label blockStart,
99  const scalarField& source,
100  const FieldField<Field, scalar>& interfaceBouCoeffs_,
101  const lduInterfaceFieldPtrsList& interfaces_,
102  const direction cmpt,
103  const label nSweeps
104 )
105 {
106  scalar* __restrict__ psiPtr = psi.begin();
107 
108  const label nCells = psi.size();
109 
110  scalarField bPrime(nCells);
111  scalar* __restrict__ bPrimePtr = bPrime.begin();
112 
113  const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
114  const scalar* const __restrict__ upperPtr =
115  matrix_.upper().begin();
116  const scalar* const __restrict__ lowerPtr =
117  matrix_.lower().begin();
118 
119  const label* const __restrict__ uPtr =
120  matrix_.lduAddr().upperAddr().begin();
121 
122  const label* const __restrict__ ownStartPtr =
123  matrix_.lduAddr().ownerStartAddr().begin();
124 
125  // Parallel boundary initialisation. The parallel boundary is treated
126  // as an effective jacobi interface in the boundary.
127  // Note: there is a change of sign in the coupled
128  // interface update. The reason for this is that the
129  // internal coefficients are all located at the l.h.s. of
130  // the matrix whereas the "implicit" coefficients on the
131  // coupled boundaries are all created as if the
132  // coefficient contribution is of a source-kind (i.e. they
133  // have a sign as if they are on the r.h.s. of the matrix.
134  // To compensate for this, it is necessary to turn the
135  // sign of the contribution.
136 
137  FieldField<Field, scalar>& mBouCoeffs =
138  const_cast<FieldField<Field, scalar>&>
139  (
140  interfaceBouCoeffs_
141  );
142  forAll(mBouCoeffs, patchi)
143  {
144  if (interfaces_.set(patchi))
145  {
146  mBouCoeffs[patchi].negate();
147  }
148  }
149 
150  for (label sweep=0; sweep<nSweeps; sweep++)
151  {
152  bPrime = source;
153 
154  matrix_.initMatrixInterfaces
155  (
156  mBouCoeffs,
157  interfaces_,
158  psi,
159  bPrime,
160  cmpt
161  );
162 
163  scalar curPsi;
164  label fStart;
165  label fEnd = ownStartPtr[0];
166 
167  for (label celli=0; celli<blockStart; celli++)
168  {
169  // Start and end of this row
170  fStart = fEnd;
171  fEnd = ownStartPtr[celli + 1];
172 
173  // Get the accumulated neighbour side
174  curPsi = bPrimePtr[celli];
175 
176  // Accumulate the owner product side
177  for (label curFace=fStart; curFace<fEnd; curFace++)
178  {
179  curPsi -= upperPtr[curFace]*psiPtr[uPtr[curFace]];
180  }
181 
182  // Finish current psi
183  curPsi /= diagPtr[celli];
184 
185  // Distribute the neighbour side using current psi
186  for (label curFace=fStart; curFace<fEnd; curFace++)
187  {
188  bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curPsi;
189  }
190 
191  psiPtr[celli] = curPsi;
192  }
193 
194  matrix_.updateMatrixInterfaces
195  (
196  mBouCoeffs,
197  interfaces_,
198  psi,
199  bPrime,
200  cmpt
201  );
202 
203  // Update rest of the cells
204  for (label celli=blockStart; celli < nCells; celli++)
205  {
206  // Start and end of this row
207  fStart = fEnd;
208  fEnd = ownStartPtr[celli + 1];
209 
210  // Get the accumulated neighbour side
211  curPsi = bPrimePtr[celli];
212 
213  // Accumulate the owner product side
214  for (label curFace=fStart; curFace<fEnd; curFace++)
215  {
216  curPsi -= upperPtr[curFace]*psiPtr[uPtr[curFace]];
217  }
218 
219  // Finish current psi
220  curPsi /= diagPtr[celli];
221 
222  // Distribute the neighbour side using current psi
223  for (label curFace=fStart; curFace<fEnd; curFace++)
224  {
225  bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curPsi;
226  }
227 
228  psiPtr[celli] = curPsi;
229  }
230  }
231 
232  // Restore interfaceBouCoeffs_
233  forAll(mBouCoeffs, patchi)
234  {
235  if (interfaces_.set(patchi))
236  {
237  mBouCoeffs[patchi].negate();
238  }
239  }
240 }
241 
242 
244 (
245  scalarField& psi,
246  const scalarField& source,
247  const direction cmpt,
248  const label nSweeps
249 ) const
250 {
251  smooth
252  (
253  fieldName_,
254  psi,
255  matrix_,
256  blockStart_,
257  source,
258  interfaceBouCoeffs_,
259  interfaces_,
260  cmpt,
261  nSweeps
262  );
263 }
264 
265 
266 // ************************************************************************* //
nonBlockingGaussSeidelSmoother(const word &fieldName, const lduMatrix &matrix, const FieldField< Field, scalar > &interfaceBouCoeffs, const FieldField< Field, scalar > &interfaceIntCoeffs, const lduInterfaceFieldPtrsList &interfaces)
Construct from components.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
lduMatrix::smoother::addsymMatrixConstructorToTable< nonBlockingGaussSeidelSmoother > addnonBlockingGaussSeidelSmootherSymMatrixConstructorToTable_
uint8_t direction
Definition: direction.H:45
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void initMatrixInterfaces(const FieldField< Field, scalar > &interfaceCoeffs, const lduInterfaceFieldPtrsList &interfaces, const scalarField &psiif, scalarField &result, const direction cmpt) const
Initialise the update of interfaced interfaces.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
scalarField & upper()
Definition: lduMatrix.C:197
Generic field type.
Definition: FieldField.H:51
void smooth(volScalarField &field, const scalar coeff)
Definition: fvcSmooth.C:35
void sweep(volScalarField &field, const volScalarField &alpha, const label nLayers, const scalar alphaDiff=0.2)
Definition: fvcSmooth.C:240
Abstract base-class for lduMatrix smoothers.
Definition: lduMatrix.H:270
void negate()
Negate this field.
Definition: FieldField.C:218
A class for handling words, derived from string.
Definition: word.H:59
virtual const labelUList & upperAddr() const =0
Return upper addressing.
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
virtual const labelUList & patchAddr(const label patchNo) const =0
Return patch to internal addressing given patch number.
void updateMatrixInterfaces(const FieldField< Field, scalar > &interfaceCoeffs, const lduInterfaceFieldPtrsList &interfaces, const scalarField &psiif, scalarField &result, const direction cmpt) const
Update interfaced interfaces for matrix operations.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
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
bool set(const label) const
Is element set.
Definition: UPtrListI.H:78
defineTypeNameAndDebug(combustionModel, 0)
lduMatrix::smoother::addasymMatrixConstructorToTable< nonBlockingGaussSeidelSmoother > addnonBlockingGaussSeidelSmootherAsymMatrixConstructorToTable_
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
lduMatrix is a general matrix class in which the coefficients are stored as three arrays...
Definition: lduMatrix.H:79
label patchi
const lduAddressing & lduAddr() const
Return the LDU addressing.
Definition: lduMatrix.H:550
scalarField & lower()
Definition: lduMatrix.C:168
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
scalarField & diag()
Definition: lduMatrix.C:186
static void smooth(const word &fieldName, scalarField &psi, const lduMatrix &matrix, const label blockStart, const scalarField &source, const FieldField< Field, scalar > &interfaceBouCoeffs, const lduInterfaceFieldPtrsList &interfaces, const direction cmpt, const label nSweeps)
Smooth for the given number of sweeps.
const labelUList & ownerStartAddr() const
Return owner start addressing.
Namespace for OpenFOAM.