edgeLimitedFaGrads.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) 2016-2017 Wikki Ltd
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "edgeLimitedFaGrad.H"
29#include "gaussFaGrad.H"
30#include "faMesh.H"
31#include "areaFaMesh.H"
32#include "edgeFaMesh.H"
33#include "areaFields.H"
35
36// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37
38makeFaGradScheme(edgeLimitedGrad)
39
40// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41
42namespace Foam
43{
44
45// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
46
47namespace fa
48{
49
50// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
51
52template<class Type>
53inline void edgeLimitedGrad<Type>::limitEdge
54(
55 scalar& limiter,
56 const scalar maxDelta,
57 const scalar minDelta,
58 const scalar extrapolate
59) const
60{
61 if (extrapolate > maxDelta + VSMALL)
62 {
63 limiter = min(limiter, maxDelta/extrapolate);
64 }
65 else if (extrapolate < minDelta - VSMALL)
66 {
67 limiter = min(limiter, minDelta/extrapolate);
68 }
69}
70
71
72// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
73
74template<>
75tmp<areaVectorField> edgeLimitedGrad<scalar>::grad
76(
77 const areaScalarField& vsf
78) const
79{
80 const faMesh& mesh = vsf.mesh();
81
82 tmp<areaVectorField> tGrad = basicGradScheme_().grad(vsf);
83
84 if (k_ < SMALL)
85 {
86 return tGrad;
87 }
88
89 areaVectorField& g = tGrad.ref();
90
91 const labelUList& owner = mesh.owner();
92 const labelUList& neighbour = mesh.neighbour();
93
94 const areaVectorField& C = mesh.areaCentres();
95 const edgeVectorField& Cf = mesh.edgeCentres();
96
97 // create limiter
98 scalarField limiter(vsf.internalField().size(), 1.0);
99
100 scalar rk = (1.0/k_ - 1.0);
101
102 forAll(owner, edgei)
103 {
104 label own = owner[edgei];
105 label nei = neighbour[edgei];
106
107 scalar vsfOwn = vsf[own];
108 scalar vsfNei = vsf[nei];
109
110 scalar maxEdge = max(vsfOwn, vsfNei);
111 scalar minEdge = min(vsfOwn, vsfNei);
112 scalar maxMinEdge = rk*(maxEdge - minEdge);
113 maxEdge += maxMinEdge;
114 minEdge -= maxMinEdge;
115
116 // owner side
117 limitEdge
118 (
119 limiter[own],
120 maxEdge - vsfOwn, minEdge - vsfOwn,
121 (Cf[edgei] - C[own]) & g[own]
122 );
123
124 // neighbour side
125 limitEdge
126 (
127 limiter[nei],
128 maxEdge - vsfNei, minEdge - vsfNei,
129 (Cf[edgei] - C[nei]) & g[nei]
130 );
131 }
132
133 const areaScalarField::Boundary& bsf = vsf.boundaryField();
134
135 forAll(bsf, patchi)
136 {
137 const faPatchScalarField& psf = bsf[patchi];
138
139 const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
140 const vectorField& pCf = Cf.boundaryField()[patchi];
141
142 if (psf.coupled())
143 {
144 const scalarField psfNei(psf.patchNeighbourField());
145
146 forAll(pOwner, pEdgei)
147 {
148 label own = pOwner[pEdgei];
149
150 scalar vsfOwn = vsf[own];
151 scalar vsfNei = psfNei[pEdgei];
152
153 scalar maxEdge = max(vsfOwn, vsfNei);
154 scalar minEdge = min(vsfOwn, vsfNei);
155 scalar maxMinEdge = rk*(maxEdge - minEdge);
156 maxEdge += maxMinEdge;
157 minEdge -= maxMinEdge;
158
159 limitEdge
160 (
161 limiter[own],
162 maxEdge - vsfOwn, minEdge - vsfOwn,
163 (pCf[pEdgei] - C[own]) & g[own]
164 );
165 }
166 }
167 else if (psf.fixesValue())
168 {
169 forAll(pOwner, pEdgei)
170 {
171 label own = pOwner[pEdgei];
172
173 scalar vsfOwn = vsf[own];
174 scalar vsfNei = psf[pEdgei];
175
176 scalar maxEdge = max(vsfOwn, vsfNei);
177 scalar minEdge = min(vsfOwn, vsfNei);
178 scalar maxMinEdge = rk*(maxEdge - minEdge);
179 maxEdge += maxMinEdge;
180 minEdge -= maxMinEdge;
181
182 limitEdge
183 (
184 limiter[own],
185 maxEdge - vsfOwn, minEdge - vsfOwn,
186 (pCf[pEdgei] - C[own]) & g[own]
187 );
188 }
189 }
190 }
191
192 if (fa::debug)
193 {
194 Info<< "gradient limiter for: " << vsf.name()
195 << " max = " << gMax(limiter)
196 << " min = " << gMin(limiter)
197 << " average: " << gAverage(limiter) << endl;
198 }
199
200 g.primitiveFieldRef() *= limiter;
201 g.correctBoundaryConditions();
203
204 return tGrad;
205}
206
207
208template<>
209tmp<areaTensorField> edgeLimitedGrad<vector>::grad
210(
211 const areaVectorField& vvf
212) const
213{
214 const faMesh& mesh = vvf.mesh();
215
216 tmp<areaTensorField> tGrad = basicGradScheme_().grad(vvf);
217
218 if (k_ < SMALL)
219 {
220 return tGrad;
221 }
222
223 areaTensorField& g = tGrad.ref();
224
225 const labelUList& owner = mesh.owner();
226 const labelUList& neighbour = mesh.neighbour();
227
228 const areaVectorField& C = mesh.areaCentres();
229 const edgeVectorField& Cf = mesh.edgeCentres();
230
231 // create limiter
232 scalarField limiter(vvf.internalField().size(), 1.0);
233
234 scalar rk = (1.0/k_ - 1.0);
235
236 forAll(owner, edgei)
237 {
238 label own = owner[edgei];
239 label nei = neighbour[edgei];
240
241 vector vvfOwn = vvf[own];
242 vector vvfNei = vvf[nei];
243
244 // owner side
245 vector gradf = (Cf[edgei] - C[own]) & g[own];
246
247 scalar vsfOwn = gradf & vvfOwn;
248 scalar vsfNei = gradf & vvfNei;
249
250 scalar maxEdge = max(vsfOwn, vsfNei);
251 scalar minEdge = min(vsfOwn, vsfNei);
252 scalar maxMinEdge = rk*(maxEdge - minEdge);
253 maxEdge += maxMinEdge;
254 minEdge -= maxMinEdge;
255
256 limitEdge
257 (
258 limiter[own],
259 maxEdge - vsfOwn, minEdge - vsfOwn,
260 magSqr(gradf)
261 );
262
263
264 // neighbour side
265 gradf = (Cf[edgei] - C[nei]) & g[nei];
266
267 vsfOwn = gradf & vvfOwn;
268 vsfNei = gradf & vvfNei;
269
270 maxEdge = max(vsfOwn, vsfNei);
271 minEdge = min(vsfOwn, vsfNei);
272
273 limitEdge
274 (
275 limiter[nei],
276 maxEdge - vsfNei, minEdge - vsfNei,
277 magSqr(gradf)
278 );
279 }
280
281
282 const areaVectorField::Boundary& bvf = vvf.boundaryField();
283
284 forAll(bvf, patchi)
285 {
286 const faPatchVectorField& psf = bvf[patchi];
287
288 const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
289 const vectorField& pCf = Cf.boundaryField()[patchi];
290
291 if (psf.coupled())
292 {
293 const vectorField psfNei(psf.patchNeighbourField());
294
295 forAll(pOwner, pEdgei)
296 {
297 label own = pOwner[pEdgei];
298
299 vector vvfOwn = vvf[own];
300 vector vvfNei = psfNei[pEdgei];
301
302 vector gradf = (pCf[pEdgei] - C[own]) & g[own];
303
304 scalar vsfOwn = gradf & vvfOwn;
305 scalar vsfNei = gradf & vvfNei;
306
307 scalar maxEdge = max(vsfOwn, vsfNei);
308 scalar minEdge = min(vsfOwn, vsfNei);
309 scalar maxMinEdge = rk*(maxEdge - minEdge);
310 maxEdge += maxMinEdge;
311 minEdge -= maxMinEdge;
312
313 limitEdge
314 (
315 limiter[own],
316 maxEdge - vsfOwn, minEdge - vsfOwn,
317 magSqr(gradf)
318 );
319 }
320 }
321 else if (psf.fixesValue())
322 {
323 forAll(pOwner, pEdgei)
324 {
325 label own = pOwner[pEdgei];
326
327 vector vvfOwn = vvf[own];
328 vector vvfNei = psf[pEdgei];
329
330 vector gradf = (pCf[pEdgei] - C[own]) & g[own];
331
332 scalar vsfOwn = gradf & vvfOwn;
333 scalar vsfNei = gradf & vvfNei;
334
335 scalar maxEdge = max(vsfOwn, vsfNei);
336 scalar minEdge = min(vsfOwn, vsfNei);
337 scalar maxMinEdge = rk*(maxEdge - minEdge);
338 maxEdge += maxMinEdge;
339 minEdge -= maxMinEdge;
340
341 limitEdge
342 (
343 limiter[own],
344 maxEdge - vsfOwn, minEdge - vsfOwn,
345 magSqr(gradf)
346 );
347 }
348 }
349 }
350
351 if (fa::debug)
352 {
353 Info<< "gradient limiter for: " << vvf.name()
354 << " max = " << gMax(limiter)
355 << " min = " << gMin(limiter)
356 << " average: " << gAverage(limiter) << endl;
357 }
358
359 g.primitiveFieldRef() *= limiter;
360 g.correctBoundaryConditions();
362
363 return tGrad;
364}
365
366
367// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
368
369} // End namespace fa
370
371// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
372
373} // End namespace Foam
374
375// ************************************************************************* //
const uniformDimensionedVectorField & g
GeometricBoundaryField< scalar, faPatchField, areaMesh > Boundary
Type of boundary fields.
void correctBoundaryConditions()
Correct boundary field.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
dynamicFvMesh & mesh
#define makeFaGradScheme(SS)
Definition: faGradScheme.H:164
Namespace for OpenFOAM.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
GeometricField< vector, faePatchField, edgeMesh > edgeVectorField
Definition: edgeFieldsFwd.H:64
messageStream Info
Information stream (stdout output on master, null elsewhere)
GeometricField< tensor, faPatchField, areaMesh > areaTensorField
Definition: areaFieldsFwd.H:83
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
faPatchField< vector > faPatchVectorField
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
GeometricField< vector, faPatchField, areaMesh > areaVectorField
Definition: areaFieldsFwd.H:79
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
GeometricField< scalar, faPatchField, areaMesh > areaScalarField
Definition: areaFieldsFwd.H:78
Field< vector > vectorField
Specialisation of Field<T> for vector.
Type gAverage(const FieldField< Field, Type > &f)
Type gMin(const FieldField< Field, Type > &f)
UList< label > labelUList
A UList of labels.
Definition: UList.H:85
Type gMax(const FieldField< Field, Type > &f)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
faPatchField< scalar > faPatchScalarField
tmp< areaScalarField > limiter(const areaScalarField &phi)
Definition: faNVDscheme.C:38
volScalarField & C
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333