TGaussSeidelSmoother.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 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 "TGaussSeidelSmoother.H"
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
33 template<class Type, class DType, class LUType>
35 (
36  const word& fieldName,
37  const LduMatrix<Type, DType, LUType>& matrix
38 )
39 :
41  (
42  fieldName,
43  matrix
44  ),
45  rD_(matrix.diag().size())
46 {
47  const label nCells = matrix.diag().size();
48  const DType* const __restrict__ diagPtr = matrix.diag().begin();
49  DType* __restrict__ rDPtr = rD_.begin();
50 
51  for (label celli=0; celli<nCells; celli++)
52  {
53  rDPtr[celli] = inv(diagPtr[celli]);
54  }
55 }
56 
57 
58 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
59 
60 template<class Type, class DType, class LUType>
62 (
63  const word& fieldName_,
65  const LduMatrix<Type, DType, LUType>& matrix_,
66  const Field<DType>& rD_,
67  const label nSweeps
68 )
69 {
70  Type* __restrict__ psiPtr = psi.begin();
71 
72  const label nCells = psi.size();
73 
74  Field<Type> bPrime(nCells);
75  Type* __restrict__ bPrimePtr = bPrime.begin();
76 
77  const DType* const __restrict__ rDPtr = rD_.begin();
78 
79  const LUType* const __restrict__ upperPtr =
80  matrix_.upper().begin();
81 
82  const LUType* const __restrict__ lowerPtr =
83  matrix_.lower().begin();
84 
85  const label* const __restrict__ uPtr =
86  matrix_.lduAddr().upperAddr().begin();
87 
88  const label* const __restrict__ ownStartPtr =
89  matrix_.lduAddr().ownerStartAddr().begin();
90 
91 
92  // Parallel boundary initialisation. The parallel boundary is treated
93  // as an effective jacobi interface in the boundary.
94  // Note: there is a change of sign in the coupled
95  // interface update to add the contibution to the r.h.s.
96 
97  for (label sweep=0; sweep<nSweeps; sweep++)
98  {
99  bPrime = matrix_.source();
100 
101  matrix_.initMatrixInterfaces
102  (
103  false,
104  matrix_.interfacesUpper(),
105  psi,
106  bPrime
107  );
108 
109  matrix_.updateMatrixInterfaces
110  (
111  false,
112  matrix_.interfacesUpper(),
113  psi,
114  bPrime
115  );
116 
117  Type curPsi;
118  label fStart;
119  label fEnd = ownStartPtr[0];
120 
121  for (label celli=0; celli<nCells; celli++)
122  {
123  // Start and end of this row
124  fStart = fEnd;
125  fEnd = ownStartPtr[celli + 1];
126 
127  // Get the accumulated neighbour side
128  curPsi = bPrimePtr[celli];
129 
130  // Accumulate the owner product side
131  for (label curFace=fStart; curFace<fEnd; curFace++)
132  {
133  curPsi -= dot(upperPtr[curFace], psiPtr[uPtr[curFace]]);
134  }
135 
136  // Finish current psi
137  curPsi = dot(rDPtr[celli], curPsi);
138 
139  // Distribute the neighbour side using current psi
140  for (label curFace=fStart; curFace<fEnd; curFace++)
141  {
142  bPrimePtr[uPtr[curFace]] -= dot(lowerPtr[curFace], curPsi);
143  }
144 
145  psiPtr[celli] = curPsi;
146  }
147  }
148 }
149 
150 
151 template<class Type, class DType, class LUType>
153 (
154  Field<Type>& psi,
155  const label nSweeps
156 ) const
157 {
158  smooth(this->fieldName_, psi, this->matrix_, rD_, nSweeps);
159 }
160 
161 
162 // ************************************************************************* //
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::LduMatrix::initMatrixInterfaces
void initMatrixInterfaces(const bool add, const FieldField< Field, LUType > &interfaceCoeffs, const Field< Type > &psiif, Field< Type > &result) const
Initialise the update of interfaced interfaces.
Definition: LduMatrixUpdateMatrixInterfaces.C:36
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
Foam::LduMatrix
LduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: LduMatrix.H:72
Foam::TGaussSeidelSmoother::TGaussSeidelSmoother
TGaussSeidelSmoother(const word &fieldName, const LduMatrix< Type, DType, LUType > &matrix)
Construct from components.
Definition: TGaussSeidelSmoother.C:35
Foam::LduMatrix::updateMatrixInterfaces
void updateMatrixInterfaces(const bool add, const FieldField< Field, LUType > &interfaceCoeffs, const Field< Type > &psiif, Field< Type > &result) const
Update interfaced interfaces for matrix operations.
Definition: LduMatrixUpdateMatrixInterfaces.C:108
Foam::LduMatrix::source
Field< Type > & source()
Definition: LduMatrix.C:250
Foam::LduMatrix::upper
Field< LUType > & upper()
Definition: LduMatrix.C:204
Foam::Field
Generic templated field type.
Definition: Field.H:63
Foam::inv
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:73
Foam::LduMatrix::diag
Field< DType > & diag()
Definition: LduMatrix.C:192
Foam::LduMatrix::smoother
Abstract base-class for LduMatrix smoothers.
Definition: LduMatrix.H:265
TGaussSeidelSmoother.H
Foam::TGaussSeidelSmoother::smooth
static void smooth(const word &fieldName, Field< Type > &psi, const LduMatrix< Type, DType, LUType > &matrix, const Field< DType > &rD, const label nSweeps)
Smooth for the given number of sweeps.
Definition: TGaussSeidelSmoother.C:62
psi
const volScalarField & psi
Definition: createFieldRefs.H:1
Foam::LduMatrix::interfacesUpper
FieldField< Field, LUType > & interfacesUpper()
Definition: LduMatrix.H:529
Foam::LduMatrix::lduAddr
const lduAddressing & lduAddr() const
Return the LDU addressing.
Definition: LduMatrix.H:498
Foam::fvc::sweep
void sweep(volScalarField &field, const volScalarField &alpha, const label nLayers, const scalar alphaDiff=0.2)
Definition: fvcSmooth.C:231
Foam::LduMatrix::lower
Field< LUType > & lower()
Definition: LduMatrix.C:227
Foam::fvc::smooth
void smooth(volScalarField &field, const scalar coeff)
Definition: fvcSmooth.C:44