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 -------------------------------------------------------------------------------
10 License
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 
38 makeFaGradScheme(edgeLimitedGrad)
39 
40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41 
42 namespace Foam
43 {
44 
45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
46 
47 namespace fa
48 {
49 
50 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
51 
52 template<class Type>
53 inline 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 
74 template<>
75 tmp<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 
208 template<>
209 tmp<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 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
gaussFaGrad.H
Foam::edgeVectorField
GeometricField< vector, faePatchField, edgeMesh > edgeVectorField
Definition: edgeFieldsFwd.H:57
Foam::gAverage
Type gAverage(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:604
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
faMesh.H
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
C
volScalarField & C
Definition: readThermalProperties.H:102
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:54
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::faPatchVectorField
faPatchField< vector > faPatchVectorField
Definition: faPatchFieldsFwd.H:43
edgeLimitedFaGrad.H
Foam::areaVectorField
GeometricField< vector, faPatchField, areaMesh > areaVectorField
Definition: areaFieldsFwd.H:58
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
fixedValueFaPatchFields.H
edgeFaMesh.H
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
areaFields.H
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
g
const uniformDimensionedVectorField & g
Definition: createFluidFields.H:26
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::areaScalarField
GeometricField< scalar, faPatchField, areaMesh > areaScalarField
Definition: areaFieldsFwd.H:53
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::fa::gaussGrad::correctBoundaryConditions
static void correctBoundaryConditions(const GeometricField< Type, faPatchField, areaMesh > &, GeometricField< typename outerProduct< vector, Type >::type, faPatchField, areaMesh > &)
Correct the boundary values of the gradient using the patchField.
Definition: gaussFaGrad.C:82
areaFaMesh.H
Foam::faPatchScalarField
faPatchField< scalar > faPatchScalarField
Definition: faPatchFieldsFwd.H:40
makeFaGradScheme
makeFaGradScheme(edgeLimitedGrad)
Definition: edgeLimitedFaGrads.C:38
Foam::gMin
Type gMin(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:593
Foam::labelUList
UList< label > labelUList
A UList of labels.
Definition: UList.H:85
Foam::limiter
tmp< areaScalarField > limiter(const areaScalarField &phi)
Definition: faNVDscheme.C:37
Foam::fa::edgeLimitedGrad::grad
tmp< GeometricField< typename outerProduct< vector, Type >::type, faPatchField, areaMesh > > grad(const GeometricField< Type, faPatchField, areaMesh > &) const
Calculate and return the grad of the given field.
Foam::gMax
Type gMax(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:592
Foam::areaTensorField
GeometricField< tensor, faPatchField, areaMesh > areaTensorField
Definition: areaFieldsFwd.H:62