faceLimitedGrads.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) 2021 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
29#include "faceLimitedGrad.H"
30#include "gaussGrad.H"
31#include "fvMesh.H"
32#include "volMesh.H"
33#include "surfaceMesh.H"
34#include "volFields.H"
36
37// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38
39makeFvGradScheme(faceLimitedGrad)
40
41// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42
43template<>
46(
47 const volScalarField& vsf,
48 const word& name
49) const
50{
51 const fvMesh& mesh = vsf.mesh();
52
53 tmp<volVectorField> tGrad = basicGradScheme_().calcGrad(vsf, name);
54
55 if (k_ < SMALL)
56 {
57 return tGrad;
58 }
59
60 volVectorField& g = tGrad.ref();
61
62 const labelUList& owner = mesh.owner();
63 const labelUList& neighbour = mesh.neighbour();
64
65 const volVectorField& C = mesh.C();
66 const surfaceVectorField& Cf = mesh.Cf();
67
68 // create limiter
69 scalarField limiter(vsf.primitiveField().size(), 1.0);
70
71 const scalar rk = (1.0/k_ - 1.0);
72
73 forAll(owner, facei)
74 {
75 const label own = owner[facei];
76 const label nei = neighbour[facei];
77
78 const scalar vsfOwn = vsf[own];
79 const scalar vsfNei = vsf[nei];
80
81 scalar maxFace = max(vsfOwn, vsfNei);
82 scalar minFace = min(vsfOwn, vsfNei);
83 const scalar maxMinFace = rk*(maxFace - minFace);
84 maxFace += maxMinFace;
85 minFace -= maxMinFace;
86
87 // owner side
88 limitFace
89 (
90 limiter[own],
91 maxFace - vsfOwn, minFace - vsfOwn,
92 (Cf[facei] - C[own]) & g[own]
93 );
94
95 // neighbour side
96 limitFace
97 (
98 limiter[nei],
99 maxFace - vsfNei, minFace - vsfNei,
100 (Cf[facei] - C[nei]) & g[nei]
101 );
102 }
103
104 const volScalarField::Boundary& bsf = vsf.boundaryField();
105
106 forAll(bsf, patchi)
107 {
108 const fvPatchScalarField& psf = bsf[patchi];
109
110 const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
111 const vectorField& pCf = Cf.boundaryField()[patchi];
112
113 if (psf.coupled())
114 {
115 const scalarField psfNei(psf.patchNeighbourField());
116
117 forAll(pOwner, pFacei)
118 {
119 const label own = pOwner[pFacei];
120
121 const scalar vsfOwn = vsf[own];
122 const scalar vsfNei = psfNei[pFacei];
123
124 scalar maxFace = max(vsfOwn, vsfNei);
125 scalar minFace = min(vsfOwn, vsfNei);
126 const scalar maxMinFace = rk*(maxFace - minFace);
127 maxFace += maxMinFace;
128 minFace -= maxMinFace;
129
130 limitFace
131 (
132 limiter[own],
133 maxFace - vsfOwn, minFace - vsfOwn,
134 (pCf[pFacei] - C[own]) & g[own]
135 );
136 }
137 }
138 else if (psf.fixesValue())
139 {
140 forAll(pOwner, pFacei)
141 {
142 const label own = pOwner[pFacei];
143
144 const scalar vsfOwn = vsf[own];
145 const scalar vsfNei = psf[pFacei];
146
147 scalar maxFace = max(vsfOwn, vsfNei);
148 scalar minFace = min(vsfOwn, vsfNei);
149 const scalar maxMinFace = rk*(maxFace - minFace);
150 maxFace += maxMinFace;
151 minFace -= maxMinFace;
152
153 limitFace
154 (
155 limiter[own],
156 maxFace - vsfOwn, minFace - vsfOwn,
157 (pCf[pFacei] - C[own]) & g[own]
158 );
159 }
160 }
161 }
162
163 if (fv::debug)
164 {
165 Info<< "gradient limiter for: " << vsf.name()
166 << " max = " << gMax(limiter)
167 << " min = " << gMin(limiter)
168 << " average: " << gAverage(limiter) << endl;
169 }
170
171 g.primitiveFieldRef() *= limiter;
172 g.correctBoundaryConditions();
173 gaussGrad<scalar>::correctBoundaryConditions(vsf, g);
174
175 return tGrad;
176}
177
178
179template<>
182(
183 const volVectorField& vvf,
184 const word& name
185) const
186{
187 const fvMesh& mesh = vvf.mesh();
188
189 tmp<volTensorField> tGrad = basicGradScheme_().calcGrad(vvf, name);
190
191 if (k_ < SMALL)
192 {
193 return tGrad;
194 }
195
196 volTensorField& g = tGrad.ref();
197
198 const labelUList& owner = mesh.owner();
199 const labelUList& neighbour = mesh.neighbour();
200
201 const volVectorField& C = mesh.C();
202 const surfaceVectorField& Cf = mesh.Cf();
203
204 // create limiter
206
207 const scalar rk = (1.0/k_ - 1.0);
208
209 forAll(owner, facei)
210 {
211 const label own = owner[facei];
212 const label nei = neighbour[facei];
213
214 const vector& vvfOwn = vvf[own];
215 const vector& vvfNei = vvf[nei];
216
217 // owner side
218 vector gradf((Cf[facei] - C[own]) & g[own]);
219
220 scalar vsfOwn = gradf & vvfOwn;
221 scalar vsfNei = gradf & vvfNei;
222
223 scalar maxFace = max(vsfOwn, vsfNei);
224 scalar minFace = min(vsfOwn, vsfNei);
225 const scalar maxMinFace = rk*(maxFace - minFace);
226 maxFace += maxMinFace;
227 minFace -= maxMinFace;
228
229 limitFace
230 (
231 limiter[own],
232 maxFace - vsfOwn, minFace - vsfOwn,
233 magSqr(gradf)
234 );
235
236
237 // neighbour side
238 gradf = (Cf[facei] - C[nei]) & g[nei];
239
240 vsfOwn = gradf & vvfOwn;
241 vsfNei = gradf & vvfNei;
242
243 maxFace = max(vsfOwn, vsfNei);
244 minFace = min(vsfOwn, vsfNei);
245
246 limitFace
247 (
248 limiter[nei],
249 maxFace - vsfNei, minFace - vsfNei,
250 magSqr(gradf)
251 );
252 }
253
254
255 const volVectorField::Boundary& bvf = vvf.boundaryField();
256
257 forAll(bvf, patchi)
258 {
259 const fvPatchVectorField& psf = bvf[patchi];
260
261 const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
262 const vectorField& pCf = Cf.boundaryField()[patchi];
263
264 if (psf.coupled())
265 {
266 const vectorField psfNei(psf.patchNeighbourField());
267
268 forAll(pOwner, pFacei)
269 {
270 const label own = pOwner[pFacei];
271
272 const vector& vvfOwn = vvf[own];
273 const vector& vvfNei = psfNei[pFacei];
274
275 const vector gradf((pCf[pFacei] - C[own]) & g[own]);
276
277 const scalar vsfOwn = gradf & vvfOwn;
278 const scalar vsfNei = gradf & vvfNei;
279
280 scalar maxFace = max(vsfOwn, vsfNei);
281 scalar minFace = min(vsfOwn, vsfNei);
282 const scalar maxMinFace = rk*(maxFace - minFace);
283 maxFace += maxMinFace;
284 minFace -= maxMinFace;
285
286 limitFace
287 (
288 limiter[own],
289 maxFace - vsfOwn, minFace - vsfOwn,
290 magSqr(gradf)
291 );
292 }
293 }
294 else if (psf.fixesValue())
295 {
296 forAll(pOwner, pFacei)
297 {
298 const label own = pOwner[pFacei];
299
300 const vector& vvfOwn = vvf[own];
301 const vector& vvfNei = psf[pFacei];
302
303 const vector gradf((pCf[pFacei] - C[own]) & g[own]);
304
305 const scalar vsfOwn = gradf & vvfOwn;
306 const scalar vsfNei = gradf & vvfNei;
307
308 scalar maxFace = max(vsfOwn, vsfNei);
309 scalar minFace = min(vsfOwn, vsfNei);
310 const scalar maxMinFace = rk*(maxFace - minFace);
311 maxFace += maxMinFace;
312 minFace -= maxMinFace;
313
314 limitFace
315 (
316 limiter[own],
317 maxFace - vsfOwn, minFace - vsfOwn,
318 magSqr(gradf)
319 );
320 }
321 }
322 }
323
324 if (fv::debug)
325 {
326 Info<< "gradient limiter for: " << vvf.name()
327 << " max = " << gMax(limiter)
328 << " min = " << gMin(limiter)
329 << " average: " << gAverage(limiter) << endl;
330 }
331
332 g.primitiveFieldRef() *= limiter;
333 g.correctBoundaryConditions();
335
336 return tGrad;
337}
338
339
340// ************************************************************************* //
static const Foam::dimensionedScalar C("", Foam::dimTemperature, 234.5)
Y[inertIndex] max(0.0)
const uniformDimensionedVectorField & g
Graphite solid properties.
Definition: C.H:53
const Mesh & mesh() const
Return mesh.
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
virtual bool fixesValue() const
Return true if this patch field fixes a value.
Definition: fvPatchField.H:337
virtual tmp< Field< Type > > patchNeighbourField() const
Return patchField on the opposite patch of a coupled patch.
Definition: fvPatchField.H:453
virtual bool coupled() const
Return true if this patch field is coupled.
Definition: fvPatchField.H:350
faceLimitedGrad gradient scheme applied to a runTime selected base gradient scheme.
Basic second-order gradient scheme using face-interpolation and Gauss' theorem.
Definition: gaussGrad.H:66
A class for managing temporary objects.
Definition: tmp.H:65
T & ref() const
Definition: tmpI.H:227
A class for handling words, derived from Foam::string.
Definition: word.H:68
dynamicFvMesh & mesh
#define makeFvGradScheme(SS)
Definition: gradScheme.H:207
messageStream Info
Information stream (stdout output on master, null elsewhere)
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
Type gAverage(const FieldField< Field, Type > &f)
Type gMin(const FieldField< Field, Type > &f)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Type gMax(const FieldField< Field, Type > &f)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
tmp< areaScalarField > limiter(const areaScalarField &phi)
Definition: faNVDscheme.C:38
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333