GaussSeidelSmoother.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-2015 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 "GaussSeidelSmoother.H"
30 #include "PrecisionAdaptor.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(GaussSeidelSmoother, 0);
37 
38  lduMatrix::smoother::addsymMatrixConstructorToTable<GaussSeidelSmoother>
40 
41  lduMatrix::smoother::addasymMatrixConstructorToTable<GaussSeidelSmoother>
43 }
44 
45 
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47 
49 (
50  const word& fieldName,
51  const lduMatrix& matrix,
52  const FieldField<Field, scalar>& interfaceBouCoeffs,
53  const FieldField<Field, scalar>& interfaceIntCoeffs,
54  const lduInterfaceFieldPtrsList& interfaces
55 )
56 :
58  (
59  fieldName,
60  matrix,
61  interfaceBouCoeffs,
62  interfaceIntCoeffs,
63  interfaces
64  )
65 {}
66 
67 
68 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
69 
71 (
72  const word& fieldName_,
74  const lduMatrix& matrix_,
75  const solveScalarField& source,
76  const FieldField<Field, scalar>& interfaceBouCoeffs_,
77  const lduInterfaceFieldPtrsList& interfaces_,
78  const direction cmpt,
79  const label nSweeps
80 )
81 {
82  solveScalar* __restrict__ psiPtr = psi.begin();
83 
84  const label nCells = psi.size();
85 
86  solveScalarField bPrime(nCells);
87  solveScalar* __restrict__ bPrimePtr = bPrime.begin();
88 
89  const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
90  const scalar* const __restrict__ upperPtr =
91  matrix_.upper().begin();
92  const scalar* const __restrict__ lowerPtr =
93  matrix_.lower().begin();
94 
95  const label* const __restrict__ uPtr =
96  matrix_.lduAddr().upperAddr().begin();
97 
98  const label* const __restrict__ ownStartPtr =
99  matrix_.lduAddr().ownerStartAddr().begin();
100 
101 
102  // Parallel boundary initialisation. The parallel boundary is treated
103  // as an effective jacobi interface in the boundary.
104  // Note: there is a change of sign in the coupled
105  // interface update. The reason for this is that the
106  // internal coefficients are all located at the l.h.s. of
107  // the matrix whereas the "implicit" coefficients on the
108  // coupled boundaries are all created as if the
109  // coefficient contribution is of a source-kind (i.e. they
110  // have a sign as if they are on the r.h.s. of the matrix.
111  // To compensate for this, it is necessary to turn the
112  // sign of the contribution.
113 
114  for (label sweep=0; sweep<nSweeps; sweep++)
115  {
116  bPrime = source;
117 
118  const label startRequest = Pstream::nRequests();
119 
120  matrix_.initMatrixInterfaces
121  (
122  false,
123  interfaceBouCoeffs_,
124  interfaces_,
125  psi,
126  bPrime,
127  cmpt
128  );
129 
130  matrix_.updateMatrixInterfaces
131  (
132  false,
133  interfaceBouCoeffs_,
134  interfaces_,
135  psi,
136  bPrime,
137  cmpt,
138  startRequest
139  );
140 
141  solveScalar psii;
142  label fStart;
143  label fEnd = ownStartPtr[0];
144 
145  for (label celli=0; celli<nCells; celli++)
146  {
147  // Start and end of this row
148  fStart = fEnd;
149  fEnd = ownStartPtr[celli + 1];
150 
151  // Get the accumulated neighbour side
152  psii = bPrimePtr[celli];
153 
154  // Accumulate the owner product side
155  for (label facei=fStart; facei<fEnd; facei++)
156  {
157  psii -= upperPtr[facei]*psiPtr[uPtr[facei]];
158  }
159 
160  // Finish psi for this cell
161  psii /= diagPtr[celli];
162 
163  // Distribute the neighbour side using psi for this cell
164  for (label facei=fStart; facei<fEnd; facei++)
165  {
166  bPrimePtr[uPtr[facei]] -= lowerPtr[facei]*psii;
167  }
168 
169  psiPtr[celli] = psii;
170  }
171  }
172 }
173 
174 
176 (
178  const scalarField& source,
179  const direction cmpt,
180  const label nSweeps
181 ) const
182 {
183  smooth
184  (
185  fieldName_,
186  psi,
187  matrix_,
189  interfaceBouCoeffs_,
190  interfaces_,
191  cmpt,
192  nSweeps
193  );
194 }
195 
196 
198 (
200  const solveScalarField& source,
201  const direction cmpt,
202  const label nSweeps
203 ) const
204 {
205  smooth
206  (
207  fieldName_,
208  psi,
209  matrix_,
210  source,
211  interfaceBouCoeffs_,
212  interfaces_,
213  cmpt,
214  nSweeps
215  );
216 }
217 
218 
219 // ************************************************************************* //
Foam::lduMatrix::lduAddr
const lduAddressing & lduAddr() const
Return the LDU addressing.
Definition: lduMatrix.H:578
Foam::GaussSeidelSmoother::smooth
static void smooth(const word &fieldName, solveScalarField &psi, const lduMatrix &matrix, 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: GaussSeidelSmoother.C:71
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::addGaussSeidelSmootherAsymMatrixConstructorToTable_
lduMatrix::smoother::addasymMatrixConstructorToTable< GaussSeidelSmoother > addGaussSeidelSmootherAsymMatrixConstructorToTable_
Definition: GaussSeidelSmoother.C:42
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::addGaussSeidelSmootherSymMatrixConstructorToTable_
lduMatrix::smoother::addsymMatrixConstructorToTable< GaussSeidelSmoother > addGaussSeidelSmootherSymMatrixConstructorToTable_
Definition: GaussSeidelSmoother.C:39
PrecisionAdaptor.H
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::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::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::GaussSeidelSmoother::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: GaussSeidelSmoother.C:198
GaussSeidelSmoother.H
Foam::GaussSeidelSmoother::GaussSeidelSmoother
GaussSeidelSmoother(const word &fieldName, const lduMatrix &matrix, const FieldField< Field, scalar > &interfaceBouCoeffs, const FieldField< Field, scalar > &interfaceIntCoeffs, const lduInterfaceFieldPtrsList &interfaces)
Construct from components.
Definition: GaussSeidelSmoother.C:49
Foam::ConstPrecisionAdaptor
A const Field/List wrapper with possible data conversion.
Definition: PrecisionAdaptor.H:57
Foam::direction
uint8_t direction
Definition: direction.H:52
psi
const volScalarField & psi
Definition: createFieldRefs.H:1
Foam::lduMatrix::lower
scalarField & lower()
Definition: lduMatrix.C:174
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
Foam::fvc::sweep
void sweep(volScalarField &field, const volScalarField &alpha, const label nLayers, const scalar alphaDiff=0.2)
Definition: fvcSmooth.C:231
Foam::fvc::smooth
void smooth(volScalarField &field, const scalar coeff)
Definition: fvcSmooth.C:44