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-------------------------------------------------------------------------------
11License
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
34namespace Foam
35{
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:
59 lduMatrix::smoother
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
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
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// ************************************************************************* //
A const Field/List wrapper with possible data conversion.
A field of fields is a PtrList of fields with reference counting.
Definition: FieldField.H:80
iterator begin() noexcept
Return an iterator to begin traversing the UList.
Definition: UListI.H:329
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
static label nRequests()
Get number of outstanding requests.
Definition: UPstream.C:90
const T * set(const label i) const
Definition: UPtrList.H:248
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
Smooth ATC in cells next to a set of patches supplied by type.
Definition: faceCells.H:59
const labelUList & ownerStartAddr() const
Return owner start addressing.
virtual const labelUList & upperAddr() const =0
Return upper addressing.
virtual const labelUList & patchAddr(const label patchNo) const =0
Return patch to internal addressing given patch number.
const lduMatrix & matrix_
Definition: lduMatrix.H:294
const lduInterfaceFieldPtrsList & interfaces() const noexcept
Definition: lduMatrix.H:406
const lduMatrix & matrix() const noexcept
Definition: lduMatrix.H:391
const FieldField< Field, scalar > & interfaceBouCoeffs() const noexcept
Definition: lduMatrix.H:396
lduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: lduMatrix.H:84
scalarField & upper()
Definition: lduMatrix.C:203
const lduAddressing & lduAddr() const
Return the LDU addressing.
Definition: lduMatrix.H:578
void initMatrixInterfaces(const bool add, const FieldField< Field, scalar > &interfaceCoeffs, const lduInterfaceFieldPtrsList &interfaces, const solveScalarField &psiif, solveScalarField &result, const direction cmpt) const
scalarField & lower()
Definition: lduMatrix.C:174
scalarField & diag()
Definition: lduMatrix.C:192
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.
Variant of gaussSeidelSmoother that expects processor boundary cells to be sorted last and so can blo...
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.
virtual void scalarSmooth(solveScalarField &psi, const solveScalarField &source, const direction cmpt, const label nSweeps) const
Smooth the solution for a given number of sweeps.
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const volScalarField & psi
Namespace for OpenFOAM.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
lduMatrix::smoother::addasymMatrixConstructorToTable< nonBlockingGaussSeidelSmoother > addnonBlockingGaussSeidelSmootherAsymMatrixConstructorToTable_
uint8_t direction
Definition: direction.H:56
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
lduMatrix::smoother::addsymMatrixConstructorToTable< nonBlockingGaussSeidelSmoother > addnonBlockingGaussSeidelSmootherSymMatrixConstructorToTable_
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333