faceLimitedFaGrads.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 "faceLimitedFaGrad.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(faceLimitedGrad)
39
40// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41
42namespace Foam
43{
44
45// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
46
47namespace fa
48{
49
50// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
51
52template<>
54(
55 scalar& limiter,
56 const scalar& maxDelta,
57 const scalar& minDelta,
58 const scalar& extrapolate
59)
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
72template<class Type>
74(
75 Type& limiter,
76 const Type& maxDelta,
77 const Type& minDelta,
78 const Type& extrapolate
79)
80{
81 for (direction cmpt = 0; cmpt < Type::nComponents; ++cmpt)
82 {
84 (
85 limiter.component(cmpt),
86 maxDelta.component(cmpt),
87 minDelta.component(cmpt),
88 extrapolate.component(cmpt)
89 );
90 }
91}
92
93
94// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
95
96template<>
97tmp<areaVectorField> faceLimitedGrad<scalar>::grad
98(
99 const areaScalarField& vsf
100) const
101{
102 const faMesh& mesh = vsf.mesh();
103
104 tmp<areaVectorField> tGrad = basicGradScheme_().grad(vsf);
105
106 if (k_ < SMALL)
107 {
108 return tGrad;
109 }
110
111 areaVectorField& g = tGrad.ref();
112
113 const labelUList& owner = mesh.owner();
114 const labelUList& neighbour = mesh.neighbour();
115
116 const areaVectorField& C = mesh.areaCentres();
117 const edgeVectorField& Cf = mesh.edgeCentres();
118
119 scalarField maxVsf(vsf.internalField());
120 scalarField minVsf(vsf.internalField());
121
122 forAll(owner, facei)
123 {
124 label own = owner[facei];
125 label nei = neighbour[facei];
126
127 scalar vsfOwn = vsf[own];
128 scalar vsfNei = vsf[nei];
129
130 maxVsf[own] = max(maxVsf[own], vsfNei);
131 minVsf[own] = min(minVsf[own], vsfNei);
132
133 maxVsf[nei] = max(maxVsf[nei], vsfOwn);
134 minVsf[nei] = min(minVsf[nei], vsfOwn);
135 }
136
137
138 const areaScalarField::Boundary& bsf = vsf.boundaryField();
139
140 forAll(bsf, patchi)
141 {
142 const faPatchScalarField& psf = bsf[patchi];
143
144 const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
145
146 if (psf.coupled())
147 {
148 scalarField psfNei(psf.patchNeighbourField());
149
150 forAll(pOwner, pFacei)
151 {
152 label own = pOwner[pFacei];
153 scalar vsfNei = psfNei[pFacei];
154
155 maxVsf[own] = max(maxVsf[own], vsfNei);
156 minVsf[own] = min(minVsf[own], vsfNei);
157 }
158 }
159 else
160 {
161 forAll(pOwner, pFacei)
162 {
163 label own = pOwner[pFacei];
164 scalar vsfNei = psf[pFacei];
165
166 maxVsf[own] = max(maxVsf[own], vsfNei);
167 minVsf[own] = min(minVsf[own], vsfNei);
168 }
169 }
170 }
171
172 maxVsf -= vsf;
173 minVsf -= vsf;
174
175 if (k_ < 1.0)
176 {
177 scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
178 maxVsf += maxMinVsf;
179 minVsf -= maxMinVsf;
180 }
181
182
183 // create limiter
184 scalarField limiter(vsf.internalField().size(), 1.0);
185
186 forAll(owner, facei)
187 {
188 label own = owner[facei];
189 label nei = neighbour[facei];
190
191 // owner side
192 limitEdge
193 (
194 limiter[own],
195 maxVsf[own],
196 minVsf[own],
197 (Cf[facei] - C[own]) & g[own]
198 );
199
200 // neighbour side
201 limitEdge
202 (
203 limiter[nei],
204 maxVsf[nei],
205 minVsf[nei],
206 (Cf[facei] - C[nei]) & g[nei]
207 );
208 }
209
210 forAll(bsf, patchi)
211 {
212 const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
213 const vectorField& pCf = Cf.boundaryField()[patchi];
214
215 forAll(pOwner, pFacei)
216 {
217 label own = pOwner[pFacei];
218
219 limitEdge
220 (
221 limiter[own],
222 maxVsf[own],
223 minVsf[own],
224 (pCf[pFacei] - C[own]) & g[own]
225 );
226 }
227 }
228
229 if (fa::debug)
230 {
231 Info<< "gradient limiter for: " << vsf.name()
232 << " max = " << gMax(limiter)
233 << " min = " << gMin(limiter)
234 << " average: " << gAverage(limiter) << endl;
235 }
236
237 g.primitiveFieldRef() *= limiter;
238 g.correctBoundaryConditions();
240
241 return tGrad;
242}
243
244
245template<>
246tmp<areaTensorField> faceLimitedGrad<vector>::grad
247(
248 const areaVectorField& vsf
249) const
250{
251 const faMesh& mesh = vsf.mesh();
252
253 tmp<areaTensorField> tGrad = basicGradScheme_().grad(vsf);
254
255 if (k_ < SMALL)
256 {
257 return tGrad;
258 }
259
260 areaTensorField& g = tGrad.ref();
261
262 const labelUList& owner = mesh.owner();
263 const labelUList& neighbour = mesh.neighbour();
264
265 const areaVectorField& C = mesh.areaCentres();
266 const edgeVectorField& Cf = mesh.edgeCentres();
267
268 vectorField maxVsf(vsf.internalField());
269 vectorField minVsf(vsf.internalField());
270
271 forAll(owner, facei)
272 {
273 label own = owner[facei];
274 label nei = neighbour[facei];
275
276 const vector& vsfOwn = vsf[own];
277 const vector& vsfNei = vsf[nei];
278
279 maxVsf[own] = max(maxVsf[own], vsfNei);
280 minVsf[own] = min(minVsf[own], vsfNei);
281
282 maxVsf[nei] = max(maxVsf[nei], vsfOwn);
283 minVsf[nei] = min(minVsf[nei], vsfOwn);
284 }
285
286
287 const areaVectorField::Boundary& bsf = vsf.boundaryField();
288
289 forAll(bsf, patchi)
290 {
291 const faPatchVectorField& psf = bsf[patchi];
292 const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
293
294 if (psf.coupled())
295 {
296 vectorField psfNei(psf.patchNeighbourField());
297
298 forAll(pOwner, pFacei)
299 {
300 label own = pOwner[pFacei];
301 const vector& vsfNei = psfNei[pFacei];
302
303 maxVsf[own] = max(maxVsf[own], vsfNei);
304 minVsf[own] = min(minVsf[own], vsfNei);
305 }
306 }
307 else
308 {
309 forAll(pOwner, pFacei)
310 {
311 label own = pOwner[pFacei];
312 const vector& vsfNei = psf[pFacei];
313
314 maxVsf[own] = max(maxVsf[own], vsfNei);
315 minVsf[own] = min(minVsf[own], vsfNei);
316 }
317 }
318 }
319
320 maxVsf -= vsf;
321 minVsf -= vsf;
322
323 if (k_ < 1.0)
324 {
325 vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
326 maxVsf += maxMinVsf;
327 minVsf -= maxMinVsf;
328
329 //maxVsf *= 1.0/k_;
330 //minVsf *= 1.0/k_;
331 }
332
333
334 // create limiter
335 vectorField limiter(vsf.internalField().size(), vector::one);
336
337 forAll(owner, facei)
338 {
339 label own = owner[facei];
340 label nei = neighbour[facei];
341
342 // owner side
343 limitEdge
344 (
345 limiter[own],
346 maxVsf[own],
347 minVsf[own],
348 (Cf[facei] - C[own]) & g[own]
349 );
350
351 // neighbour side
352 limitEdge
353 (
354 limiter[nei],
355 maxVsf[nei],
356 minVsf[nei],
357 (Cf[facei] - C[nei]) & g[nei]
358 );
359 }
360
361 forAll(bsf, patchi)
362 {
363 const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
364 const vectorField& pCf = Cf.boundaryField()[patchi];
365
366 forAll(pOwner, pFacei)
367 {
368 label own = pOwner[pFacei];
369
370 limitEdge
371 (
372 limiter[own],
373 maxVsf[own],
374 minVsf[own],
375 ((pCf[pFacei] - C[own]) & g[own])
376 );
377 }
378 }
379
380 if (fa::debug)
381 {
382 Info<< "gradient limiter for: " << vsf.name()
383 << " max = " << gMax(limiter)
384 << " min = " << gMin(limiter)
385 << " average: " << gAverage(limiter) << endl;
386 }
387
388 tensorField& gIf = g.primitiveFieldRef();
389
390 forAll(gIf, celli)
391 {
392 gIf[celli] = tensor
393 (
394 cmptMultiply(limiter[celli], gIf[celli].x()),
395 cmptMultiply(limiter[celli], gIf[celli].y()),
396 cmptMultiply(limiter[celli], gIf[celli].z())
397 );
398 }
399
400 g.correctBoundaryConditions();
402
403 return tGrad;
404}
405
406
407// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
408
409} // End namespace fa
410
411// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
412
413} // End namespace Foam
414
415// ************************************************************************* //
scalar y
const uniformDimensionedVectorField & g
GeometricBoundaryField< scalar, faPatchField, areaMesh > Boundary
Type of boundary fields.
void correctBoundaryConditions()
Correct boundary field.
static void limitEdge(Type &limiter, const Type &maxDelta, const Type &minDelta, const Type &extrapolate)
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
Tensor< scalar > tensor
Definition: symmTensor.H:61
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.
uint8_t direction
Definition: direction.H:56
Type gAverage(const FieldField< Field, Type > &f)
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
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)
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