cellLimitedGrad.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2018-2021 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 "cellLimitedGrad.H"
27 #include "gaussGrad.H"
28 
29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
30 
31 template<class Type, class Limiter>
33 (
34  const Field<scalar>& limiter,
35  Field<vector>& gIf
36 ) const
37 {
38  gIf *= limiter;
39 }
40 
41 
42 template<class Type, class Limiter>
44 (
45  const Field<vector>& limiter,
46  Field<tensor>& gIf
47 ) const
48 {
49  forAll(gIf, celli)
50  {
51  gIf[celli] = tensor
52  (
53  cmptMultiply(limiter[celli], gIf[celli].x()),
54  cmptMultiply(limiter[celli], gIf[celli].y()),
55  cmptMultiply(limiter[celli], gIf[celli].z())
56  );
57  }
58 }
59 
60 
61 template<class Type, class Limiter>
63 <
65  <
66  typename Foam::outerProduct<Foam::vector, Type>::type,
67  Foam::fvPatchField,
69  >
70 >
72 (
73  const GeometricField<Type, fvPatchField, volMesh>& vsf,
74  const word& name
75 ) const
76 {
77  const fvMesh& mesh = vsf.mesh();
78 
79  tmp
80  <
81  GeometricField
82  <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
83  > tGrad = basicGradScheme_().calcGrad(vsf, name);
84 
85  if (k_ < small)
86  {
87  return tGrad;
88  }
89 
90  GeometricField
91  <
92  typename outerProduct<vector, Type>::type,
93  fvPatchField,
94  volMesh
95  >& g = tGrad.ref();
96 
97  const labelUList& owner = mesh.owner();
98  const labelUList& neighbour = mesh.neighbour();
99 
100  const volVectorField& C = mesh.C();
101  const surfaceVectorField& Cf = mesh.Cf();
102 
103  Field<Type> maxVsf(vsf.primitiveField());
104  Field<Type> minVsf(vsf.primitiveField());
105 
106  forAll(owner, facei)
107  {
108  label own = owner[facei];
109  label nei = neighbour[facei];
110 
111  const Type& vsfOwn = vsf[own];
112  const Type& vsfNei = vsf[nei];
113 
114  maxVsf[own] = max(maxVsf[own], vsfNei);
115  minVsf[own] = min(minVsf[own], vsfNei);
116 
117  maxVsf[nei] = max(maxVsf[nei], vsfOwn);
118  minVsf[nei] = min(minVsf[nei], vsfOwn);
119  }
120 
121 
122  const typename GeometricField<Type, fvPatchField, volMesh>::Boundary& bsf =
123  vsf.boundaryField();
124 
125  forAll(bsf, patchi)
126  {
127  const fvPatchField<Type>& psf = bsf[patchi];
128  const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
129 
130  if (psf.coupled())
131  {
132  const Field<Type> psfNei(psf.patchNeighbourField());
133 
134  forAll(pOwner, pFacei)
135  {
136  label own = pOwner[pFacei];
137  const Type& vsfNei = psfNei[pFacei];
138 
139  maxVsf[own] = max(maxVsf[own], vsfNei);
140  minVsf[own] = min(minVsf[own], vsfNei);
141  }
142  }
143  else
144  {
145  forAll(pOwner, pFacei)
146  {
147  label own = pOwner[pFacei];
148  const Type& vsfNei = psf[pFacei];
149 
150  maxVsf[own] = max(maxVsf[own], vsfNei);
151  minVsf[own] = min(minVsf[own], vsfNei);
152  }
153  }
154  }
155 
156  maxVsf -= vsf;
157  minVsf -= vsf;
158 
159  if (k_ < 1.0)
160  {
161  const Field<Type> maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
162  maxVsf += maxMinVsf;
163  minVsf -= maxMinVsf;
164  }
165 
166 
167  // Create limiter initialised to 1
168  // Note: the limiter is not permitted to be > 1
169  Field<Type> limiter(vsf.primitiveField().size(), pTraits<Type>::one);
170 
171  forAll(owner, facei)
172  {
173  label own = owner[facei];
174  label nei = neighbour[facei];
175 
176  // owner side
177  limitFace
178  (
179  limiter[own],
180  maxVsf[own],
181  minVsf[own],
182  (Cf[facei] - C[own]) & g[own]
183  );
184 
185  // neighbour side
186  limitFace
187  (
188  limiter[nei],
189  maxVsf[nei],
190  minVsf[nei],
191  (Cf[facei] - C[nei]) & g[nei]
192  );
193  }
194 
195  forAll(bsf, patchi)
196  {
197  const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
198  const vectorField& pCf = Cf.boundaryField()[patchi];
199 
200  forAll(pOwner, pFacei)
201  {
202  label own = pOwner[pFacei];
203 
204  limitFace
205  (
206  limiter[own],
207  maxVsf[own],
208  minVsf[own],
209  ((pCf[pFacei] - C[own]) & g[own])
210  );
211  }
212  }
213 
214  if (fv::debug)
215  {
216  Info<< "gradient limiter for: " << vsf.name()
217  << " max = " << gMax(limiter)
218  << " min = " << gMin(limiter)
219  << " average: " << gAverage(limiter) << endl;
220  }
221 
222  limitGradient(limiter, g);
223  g.correctBoundaryConditions();
225 
226  return tGrad;
227 }
228 
229 
230 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
virtual tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > calcGrad(const GeometricField< Type, fvPatchField, volMesh > &vsf, const word &name) const
Return the gradient of the given field to the gradScheme::grad.
rDeltaT correctBoundaryConditions()
Type gMin(const FieldField< Field, Type > &f)
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
volVectorField vectorField(fieldObject, mesh)
Generic GeometricField class.
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:59
UList< label > labelUList
Definition: UList.H:65
scalar y
Mesh data needed to do the Finite Volume discretisation.
Definition: volMesh.H:53
cellLimitedGrad gradient scheme applied to a runTime selected base gradient scheme.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
void limiter(surfaceScalarField &lambda, const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phiBD, const surfaceScalarField &phiCorr, const SpType &Sp, const SuType &Su, const PsiMaxType &psiMax, const PsiMinType &psiMin)
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
Type gMax(const FieldField< Field, Type > &f)
label patchi
Type gAverage(const FieldField< Field, Type > &f)
messageStream Info
A class for managing temporary objects.
Definition: PtrList.H:53
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51