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 | Copyright (C) 2011-2015 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 
26 #include "TGaussSeidelSmoother.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class Type, class DType, class LUType>
32 (
33  const word& fieldName,
34  const LduMatrix<Type, DType, LUType>& matrix
35 )
36 :
38  (
39  fieldName,
40  matrix
41  ),
42  rD_(matrix.diag().size())
43 {
44  const label nCells = matrix.diag().size();
45  const DType* const __restrict__ diagPtr = matrix.diag().begin();
46  DType* __restrict__ rDPtr = rD_.begin();
47 
48  for (label cellI=0; cellI<nCells; cellI++)
49  {
50  rDPtr[cellI] = inv(diagPtr[cellI]);
51  }
52 }
53 
54 
55 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
56 
57 template<class Type, class DType, class LUType>
59 (
60  const word& fieldName_,
61  Field<Type>& psi,
62  const LduMatrix<Type, DType, LUType>& matrix_,
63  const Field<DType>& rD_,
64  const label nSweeps
65 )
66 {
67  Type* __restrict__ psiPtr = psi.begin();
68 
69  const label nCells = psi.size();
70 
71  Field<Type> bPrime(nCells);
72  Type* __restrict__ bPrimePtr = bPrime.begin();
73 
74  const DType* const __restrict__ rDPtr = rD_.begin();
75 
76  const LUType* const __restrict__ upperPtr =
77  matrix_.upper().begin();
78 
79  const LUType* const __restrict__ lowerPtr =
80  matrix_.lower().begin();
81 
82  const label* const __restrict__ uPtr =
83  matrix_.lduAddr().upperAddr().begin();
84 
85  const label* const __restrict__ ownStartPtr =
86  matrix_.lduAddr().ownerStartAddr().begin();
87 
88 
89  // Parallel boundary initialisation. The parallel boundary is treated
90  // as an effective jacobi interface in the boundary.
91  // Note: there is a change of sign in the coupled
92  // interface update to add the contibution to the r.h.s.
93 
94  FieldField<Field, LUType> mBouCoeffs(matrix_.interfacesUpper().size());
95 
96  forAll(mBouCoeffs, patchi)
97  {
98  if (matrix_.interfaces().set(patchi))
99  {
100  mBouCoeffs.set(patchi, -matrix_.interfacesUpper()[patchi]);
101  }
102  }
103 
104  for (label sweep=0; sweep<nSweeps; sweep++)
105  {
106  bPrime = matrix_.source();
107 
108  matrix_.initMatrixInterfaces
109  (
110  mBouCoeffs,
111  psi,
112  bPrime
113  );
114 
115  matrix_.updateMatrixInterfaces
116  (
117  mBouCoeffs,
118  psi,
119  bPrime
120  );
121 
122  Type curPsi;
123  label fStart;
124  label fEnd = ownStartPtr[0];
125 
126  for (label cellI=0; cellI<nCells; cellI++)
127  {
128  // Start and end of this row
129  fStart = fEnd;
130  fEnd = ownStartPtr[cellI + 1];
131 
132  // Get the accumulated neighbour side
133  curPsi = bPrimePtr[cellI];
134 
135  // Accumulate the owner product side
136  for (label curFace=fStart; curFace<fEnd; curFace++)
137  {
138  curPsi -= dot(upperPtr[curFace], psiPtr[uPtr[curFace]]);
139  }
140 
141  // Finish current psi
142  curPsi = dot(rDPtr[cellI], curPsi);
143 
144  // Distribute the neighbour side using current psi
145  for (label curFace=fStart; curFace<fEnd; curFace++)
146  {
147  bPrimePtr[uPtr[curFace]] -= dot(lowerPtr[curFace], curPsi);
148  }
149 
150  psiPtr[cellI] = curPsi;
151  }
152  }
153 }
154 
155 
156 template<class Type, class DType, class LUType>
158 (
159  Field<Type>& psi,
160  const label nSweeps
161 ) const
162 {
163  smooth(this->fieldName_, psi, this->matrix_, rD_, nSweeps);
164 }
165 
166 
167 // ************************************************************************* //
const LduInterfaceFieldPtrsList< Type > & interfaces() const
Return interfaces.
Definition: LduMatrix.H:511
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.
Generic field type.
Definition: FieldField.H:51
void updateMatrixInterfaces(const FieldField< Field, LUType > &interfaceCoeffs, const Field< Type > &psiif, Field< Type > &result) const
Update interfaced interfaces for matrix operations.
void initMatrixInterfaces(const FieldField< Field, LUType > &interfaceCoeffs, const Field< Type > &psiif, Field< Type > &result) const
Initialise the update of interfaced interfaces.
Field< DType > & diag()
Definition: LduMatrix.C:190
TGaussSeidelSmoother(const word &fieldName, const LduMatrix< Type, DType, LUType > &matrix)
Construct from components.
A class for handling words, derived from string.
Definition: word.H:59
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Abstract base-class for LduMatrix smoothers.
Definition: LduMatrix.H:258
const lduAddressing & lduAddr() const
Return the LDU addressing.
Definition: LduMatrix.H:499
void sweep(volScalarField &field, const volScalarField &alpha, const label nLayers, const scalar alphaDiff=0.2)
Definition: fvcSmooth.C:222
Field< LUType > & lower()
Definition: LduMatrix.C:225
Field< LUType > & upper()
Definition: LduMatrix.C:202
Field< Type > & source()
Definition: LduMatrix.C:248
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
void dot(FieldField< Field1, typename innerProduct< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
#define forAll(list, i)
Definition: UList.H:421
label patchi
Pre-declare SubField and related Field type.
Definition: Field.H:57
LduMatrix is a general matrix class in which the coefficients are stored as three arrays...
Definition: LduMatrix.H:69
FieldField< Field, LUType > & interfacesUpper()
Definition: LduMatrix.H:530
void smooth(volScalarField &field, const scalar coeff)
Definition: fvcSmooth.C:35