nonBlockingGaussSeidelSmoother.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 
30 #include "PrecisionAdaptor.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(nonBlockingGaussSeidelSmoother, 0);
37 
38  lduMatrix::smoother::
39  addsymMatrixConstructorToTable<nonBlockingGaussSeidelSmoother>
41 
42  lduMatrix::smoother::
43  addasymMatrixConstructorToTable<nonBlockingGaussSeidelSmoother>
45 }
46 
47 
48 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
49 
51 (
52  const word& fieldName,
53  const lduMatrix& matrix,
54  const FieldField<Field, scalar>& interfaceBouCoeffs,
55  const FieldField<Field, scalar>& interfaceIntCoeffs,
56  const lduInterfaceFieldPtrsList& interfaces
57 )
58 :
60  (
61  fieldName,
62  matrix,
63  interfaceBouCoeffs,
64  interfaceIntCoeffs,
65  interfaces
66  )
67 {
68  // Check that all interface addressing is sorted to be after the
69  // non-interface addressing.
70 
71  const label nCells = matrix.diag().size();
72 
73  blockStart_ = nCells;
74 
75  labelList startCelli(interfaceBouCoeffs.size(), -1);
76  forAll(interfaces, patchi)
77  {
78  if (interfaces.set(patchi))
79  {
80  const labelUList& faceCells = matrix_.lduAddr().patchAddr(patchi);
81 
82  blockStart_ = min(blockStart_, min(faceCells));
83  }
84  }
85 
86  if (debug)
87  {
88  Pout<< "nonBlockingGaussSeidelSmoother :"
89  << " Starting block on cell " << blockStart_
90  << " out of " << nCells << endl;
91  }
92 }
93 
94 
95 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
96 
98 (
99  const word& fieldName_,
101  const lduMatrix& matrix_,
102  const label blockStart,
103  const solveScalarField& source,
104  const FieldField<Field, scalar>& interfaceBouCoeffs_,
105  const lduInterfaceFieldPtrsList& interfaces_,
106  const direction cmpt,
107  const label nSweeps
108 )
109 {
110  solveScalar* __restrict__ psiPtr = psi.begin();
111 
112  const label nCells = psi.size();
113 
114  solveScalarField bPrime(nCells);
115  solveScalar* __restrict__ bPrimePtr = bPrime.begin();
116 
117  const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
118  const scalar* const __restrict__ upperPtr =
119  matrix_.upper().begin();
120  const scalar* const __restrict__ lowerPtr =
121  matrix_.lower().begin();
122 
123  const label* const __restrict__ uPtr =
124  matrix_.lduAddr().upperAddr().begin();
125 
126  const label* const __restrict__ ownStartPtr =
127  matrix_.lduAddr().ownerStartAddr().begin();
128 
129  // Parallel boundary initialisation. The parallel boundary is treated
130  // as an effective jacobi interface in the boundary.
131  // Note: there is a change of sign in the coupled
132  // interface update. The reason for this is that the
133  // internal coefficients are all located at the l.h.s. of
134  // the matrix whereas the "implicit" coefficients on the
135  // coupled boundaries are all created as if the
136  // coefficient contribution is of a source-kind (i.e. they
137  // have a sign as if they are on the r.h.s. of the matrix.
138  // To compensate for this, it is necessary to turn the
139  // sign of the contribution.
140 
141  for (label sweep=0; sweep<nSweeps; sweep++)
142  {
143  bPrime = source;
144 
145  const label startRequest = Pstream::nRequests();
146 
147  matrix_.initMatrixInterfaces
148  (
149  false,
150  interfaceBouCoeffs_,
151  interfaces_,
152  psi,
153  bPrime,
154  cmpt
155  );
156 
157  solveScalar curPsi;
158  label fStart;
159  label fEnd = ownStartPtr[0];
160 
161  for (label celli=0; celli<blockStart; celli++)
162  {
163  // Start and end of this row
164  fStart = fEnd;
165  fEnd = ownStartPtr[celli + 1];
166 
167  // Get the accumulated neighbour side
168  curPsi = bPrimePtr[celli];
169 
170  // Accumulate the owner product side
171  for (label curFace=fStart; curFace<fEnd; curFace++)
172  {
173  curPsi -= upperPtr[curFace]*psiPtr[uPtr[curFace]];
174  }
175 
176  // Finish current psi
177  curPsi /= diagPtr[celli];
178 
179  // Distribute the neighbour side using current psi
180  for (label curFace=fStart; curFace<fEnd; curFace++)
181  {
182  bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curPsi;
183  }
184 
185  psiPtr[celli] = curPsi;
186  }
187 
188  matrix_.updateMatrixInterfaces
189  (
190  false,
191  interfaceBouCoeffs_,
192  interfaces_,
193  psi,
194  bPrime,
195  cmpt,
196  startRequest
197  );
198 
199  // Update rest of the cells
200  for (label celli=blockStart; celli < nCells; celli++)
201  {
202  // Start and end of this row
203  fStart = fEnd;
204  fEnd = ownStartPtr[celli + 1];
205 
206  // Get the accumulated neighbour side
207  curPsi = bPrimePtr[celli];
208 
209  // Accumulate the owner product side
210  for (label curFace=fStart; curFace<fEnd; curFace++)
211  {
212  curPsi -= upperPtr[curFace]*psiPtr[uPtr[curFace]];
213  }
214 
215  // Finish current psi
216  curPsi /= diagPtr[celli];
217 
218  // Distribute the neighbour side using current psi
219  for (label curFace=fStart; curFace<fEnd; curFace++)
220  {
221  bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curPsi;
222  }
223 
224  psiPtr[celli] = curPsi;
225  }
226  }
227 }
228 
229 
231 (
233  const solveScalarField& source,
234  const direction cmpt,
235  const label nSweeps
236 ) const
237 {
238  smooth
239  (
240  fieldName_,
241  psi,
242  matrix_,
243  blockStart_,
244  source,
245  interfaceBouCoeffs_,
246  interfaces_,
247  cmpt,
248  nSweeps
249  );
250 }
251 
252 
254 (
256  const scalarField& source,
257  const direction cmpt,
258  const label nSweeps
259 ) const
260 {
261  scalarSmooth
262  (
263  psi,
265  cmpt,
266  nSweeps
267  );
268 }
269 
270 
271 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::lduMatrix::lduAddr
const lduAddressing & lduAddr() const
Return the LDU addressing.
Definition: lduMatrix.H:578
Foam::addnonBlockingGaussSeidelSmootherSymMatrixConstructorToTable_
lduMatrix::smoother::addsymMatrixConstructorToTable< nonBlockingGaussSeidelSmoother > addnonBlockingGaussSeidelSmootherSymMatrixConstructorToTable_
Definition: nonBlockingGaussSeidelSmoother.C:40
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::FieldField
A field of fields is a PtrList of fields with reference counting.
Definition: FieldField.H:53
Foam::nonBlockingGaussSeidelSmoother::scalarSmooth
virtual void scalarSmooth(solveScalarField &psi, const solveScalarField &source, const direction cmpt, const label nSweeps) const
Smooth the solution for a given number of sweeps.
Definition: nonBlockingGaussSeidelSmoother.C:231
Foam::lduMatrix::upper
scalarField & upper()
Definition: lduMatrix.C:203
Foam::lduMatrix
lduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: lduMatrix.H:83
Foam::nonBlockingGaussSeidelSmoother::nonBlockingGaussSeidelSmoother
nonBlockingGaussSeidelSmoother(const word &fieldName, const lduMatrix &matrix, const FieldField< Field, scalar > &interfaceBouCoeffs, const FieldField< Field, scalar > &interfaceIntCoeffs, const lduInterfaceFieldPtrsList &interfaces)
Construct from components.
Definition: nonBlockingGaussSeidelSmoother.C:51
PrecisionAdaptor.H
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::lduMatrix::updateMatrixInterfaces
void updateMatrixInterfaces(const bool add, const FieldField< Field, scalar > &interfaceCoeffs, const lduInterfaceFieldPtrsList &interfaces, const solveScalarField &psiif, solveScalarField &result, const direction cmpt, const label startRequest) const
Update interfaced interfaces for matrix operations.
Definition: lduMatrixUpdateMatrixInterfaces.C:107
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::lduAddressing::ownerStartAddr
const labelUList & ownerStartAddr() const
Return owner start addressing.
Definition: lduAddressing.C:199
Foam::lduAddressing::upperAddr
virtual const labelUList & upperAddr() const =0
Return upper addressing.
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::UList::begin
iterator begin() noexcept
Return an iterator to begin traversing the UList.
Definition: UListI.H:329
Foam::Field< solveScalar >
Foam::UPtrList< const lduInterfaceField >
Foam::lduMatrix::smoother
Abstract base-class for lduMatrix smoothers.
Definition: lduMatrix.H:287
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::lduMatrix::diag
scalarField & diag()
Definition: lduMatrix.C:192
Foam::lduMatrix::initMatrixInterfaces
void initMatrixInterfaces(const bool add, const FieldField< Field, scalar > &interfaceCoeffs, const lduInterfaceFieldPtrsList &interfaces, const solveScalarField &psiif, solveScalarField &result, const direction cmpt) const
Definition: lduMatrixUpdateMatrixInterfaces.C:34
Foam::UPtrList::set
const T * set(const label i) const
Definition: UPtrList.H:176
Foam::List< label >
Foam::UList< label >
Foam::ConstPrecisionAdaptor
A const Field/List wrapper with possible data conversion.
Definition: PrecisionAdaptor.H:57
Foam::direction
uint8_t direction
Definition: direction.H:52
Foam::lduAddressing::patchAddr
virtual const labelUList & patchAddr(const label patchNo) const =0
Return patch to internal addressing given patch number.
Foam::addnonBlockingGaussSeidelSmootherAsymMatrixConstructorToTable_
lduMatrix::smoother::addasymMatrixConstructorToTable< nonBlockingGaussSeidelSmoother > addnonBlockingGaussSeidelSmootherAsymMatrixConstructorToTable_
Definition: nonBlockingGaussSeidelSmoother.C:44
psi
const volScalarField & psi
Definition: createFieldRefs.H:1
Foam::nonBlockingGaussSeidelSmoother::smooth
static void smooth(const word &fieldName, solveScalarField &psi, const lduMatrix &matrix, const label blockStart, const solveScalarField &source, const FieldField< Field, scalar > &interfaceBouCoeffs, const lduInterfaceFieldPtrsList &interfaces, const direction cmpt, const label nSweeps)
Smooth for the given number of sweeps.
Definition: nonBlockingGaussSeidelSmoother.C:98
Foam::lduMatrix::lower
scalarField & lower()
Definition: lduMatrix.C:174
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
Foam::faceCells
Smooth ATC in cells next to a set of patches supplied by type.
Definition: faceCells.H:56
Foam::fvc::sweep
void sweep(volScalarField &field, const volScalarField &alpha, const label nLayers, const scalar alphaDiff=0.2)
Definition: fvcSmooth.C:231
nonBlockingGaussSeidelSmoother.H
Foam::fvc::smooth
void smooth(volScalarField &field, const scalar coeff)
Definition: fvcSmooth.C:44