FreeFOAM The Cross-Platform CFD Toolkit
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 | Copyright (C) 1991-2010 OpenCFD Ltd.
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 "faceLimitedGrad.H"
27 #include <finiteVolume/gaussGrad.H>
28 #include <finiteVolume/fvMesh.H>
29 #include <finiteVolume/volMesh.H>
31 #include <finiteVolume/volFields.H>
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 
39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
40 
41 namespace fv
42 {
43 
44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
45 
47 
48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
49 
50 template<class Type>
51 inline void faceLimitedGrad<Type>::limitFace
52 (
53  scalar& limiter,
54  const scalar maxDelta,
55  const scalar minDelta,
56  const scalar extrapolate
57 ) const
58 {
59  if (extrapolate > maxDelta + VSMALL)
60  {
61  limiter = min(limiter, maxDelta/extrapolate);
62  }
63  else if (extrapolate < minDelta - VSMALL)
64  {
65  limiter = min(limiter, minDelta/extrapolate);
66  }
67 }
68 
69 
70 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
71 
72 template<>
74 (
75  const volScalarField& vsf
76 ) const
77 {
78  const fvMesh& mesh = vsf.mesh();
79 
80  tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
81 
82  if (k_ < SMALL)
83  {
84  return tGrad;
85  }
86 
87  volVectorField& g = tGrad();
88 
89  const unallocLabelList& owner = mesh.owner();
90  const unallocLabelList& neighbour = mesh.neighbour();
91 
92  const volVectorField& C = mesh.C();
93  const surfaceVectorField& Cf = mesh.Cf();
94 
95  // create limiter
96  scalarField limiter(vsf.internalField().size(), 1.0);
97 
98  scalar rk = (1.0/k_ - 1.0);
99 
100  forAll(owner, facei)
101  {
102  label own = owner[facei];
103  label nei = neighbour[facei];
104 
105  scalar vsfOwn = vsf[own];
106  scalar vsfNei = vsf[nei];
107 
108  scalar maxFace = max(vsfOwn, vsfNei);
109  scalar minFace = min(vsfOwn, vsfNei);
110  scalar maxMinFace = rk*(maxFace - minFace);
111  maxFace += maxMinFace;
112  minFace -= maxMinFace;
113 
114  // owner side
115  limitFace
116  (
117  limiter[own],
118  maxFace - vsfOwn, minFace - vsfOwn,
119  (Cf[facei] - C[own]) & g[own]
120  );
121 
122  // neighbour side
123  limitFace
124  (
125  limiter[nei],
126  maxFace - vsfNei, minFace - vsfNei,
127  (Cf[facei] - C[nei]) & g[nei]
128  );
129  }
130 
131  const volScalarField::GeometricBoundaryField& bsf = vsf.boundaryField();
132 
133  forAll(bsf, patchi)
134  {
135  const fvPatchScalarField& psf = bsf[patchi];
136 
137  const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
138  const vectorField& pCf = Cf.boundaryField()[patchi];
139 
140  if (psf.coupled())
141  {
142  scalarField psfNei = psf.patchNeighbourField();
143 
144  forAll(pOwner, pFacei)
145  {
146  label own = pOwner[pFacei];
147 
148  scalar vsfOwn = vsf[own];
149  scalar vsfNei = psfNei[pFacei];
150 
151  scalar maxFace = max(vsfOwn, vsfNei);
152  scalar minFace = min(vsfOwn, vsfNei);
153  scalar maxMinFace = rk*(maxFace - minFace);
154  maxFace += maxMinFace;
155  minFace -= maxMinFace;
156 
157  limitFace
158  (
159  limiter[own],
160  maxFace - vsfOwn, minFace - vsfOwn,
161  (pCf[pFacei] - C[own]) & g[own]
162  );
163  }
164  }
165  else if (psf.fixesValue())
166  {
167  forAll(pOwner, pFacei)
168  {
169  label own = pOwner[pFacei];
170 
171  scalar vsfOwn = vsf[own];
172  scalar vsfNei = psf[pFacei];
173 
174  scalar maxFace = max(vsfOwn, vsfNei);
175  scalar minFace = min(vsfOwn, vsfNei);
176  scalar maxMinFace = rk*(maxFace - minFace);
177  maxFace += maxMinFace;
178  minFace -= maxMinFace;
179 
180  limitFace
181  (
182  limiter[own],
183  maxFace - vsfOwn, minFace - vsfOwn,
184  (pCf[pFacei] - C[own]) & g[own]
185  );
186  }
187  }
188  }
189 
190  if (fv::debug)
191  {
192  Info<< "gradient limiter for: " << vsf.name()
193  << " max = " << gMax(limiter)
194  << " min = " << gMin(limiter)
195  << " average: " << gAverage(limiter) << endl;
196  }
197 
198  g.internalField() *= limiter;
201 
202  return tGrad;
203 }
204 
205 
206 template<>
208 (
209  const volVectorField& vvf
210 ) const
211 {
212  const fvMesh& mesh = vvf.mesh();
213 
214  tmp<volTensorField> tGrad = basicGradScheme_().grad(vvf);
215 
216  if (k_ < SMALL)
217  {
218  return tGrad;
219  }
220 
221  volTensorField& g = tGrad();
222 
223  const unallocLabelList& owner = mesh.owner();
224  const unallocLabelList& neighbour = mesh.neighbour();
225 
226  const volVectorField& C = mesh.C();
227  const surfaceVectorField& Cf = mesh.Cf();
228 
229  // create limiter
230  scalarField limiter(vvf.internalField().size(), 1.0);
231 
232  scalar rk = (1.0/k_ - 1.0);
233 
234  forAll(owner, facei)
235  {
236  label own = owner[facei];
237  label nei = neighbour[facei];
238 
239  vector vvfOwn = vvf[own];
240  vector vvfNei = vvf[nei];
241 
242  // owner side
243  vector gradf = (Cf[facei] - C[own]) & g[own];
244 
245  scalar vsfOwn = gradf & vvfOwn;
246  scalar vsfNei = gradf & vvfNei;
247 
248  scalar maxFace = max(vsfOwn, vsfNei);
249  scalar minFace = min(vsfOwn, vsfNei);
250  scalar maxMinFace = rk*(maxFace - minFace);
251  maxFace += maxMinFace;
252  minFace -= maxMinFace;
253 
254  limitFace
255  (
256  limiter[own],
257  maxFace - vsfOwn, minFace - vsfOwn,
258  magSqr(gradf)
259  );
260 
261 
262  // neighbour side
263  gradf = (Cf[facei] - C[nei]) & g[nei];
264 
265  vsfOwn = gradf & vvfOwn;
266  vsfNei = gradf & vvfNei;
267 
268  maxFace = max(vsfOwn, vsfNei);
269  minFace = min(vsfOwn, vsfNei);
270 
271  limitFace
272  (
273  limiter[nei],
274  maxFace - vsfNei, minFace - vsfNei,
275  magSqr(gradf)
276  );
277  }
278 
279 
280  const volVectorField::GeometricBoundaryField& bvf = vvf.boundaryField();
281 
282  forAll(bvf, patchi)
283  {
284  const fvPatchVectorField& psf = bvf[patchi];
285 
286  const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
287  const vectorField& pCf = Cf.boundaryField()[patchi];
288 
289  if (psf.coupled())
290  {
291  vectorField psfNei = psf.patchNeighbourField();
292 
293  forAll(pOwner, pFacei)
294  {
295  label own = pOwner[pFacei];
296 
297  vector vvfOwn = vvf[own];
298  vector vvfNei = psfNei[pFacei];
299 
300  vector gradf = (pCf[pFacei] - C[own]) & g[own];
301 
302  scalar vsfOwn = gradf & vvfOwn;
303  scalar vsfNei = gradf & vvfNei;
304 
305  scalar maxFace = max(vsfOwn, vsfNei);
306  scalar minFace = min(vsfOwn, vsfNei);
307  scalar maxMinFace = rk*(maxFace - minFace);
308  maxFace += maxMinFace;
309  minFace -= maxMinFace;
310 
311  limitFace
312  (
313  limiter[own],
314  maxFace - vsfOwn, minFace - vsfOwn,
315  magSqr(gradf)
316  );
317  }
318  }
319  else if (psf.fixesValue())
320  {
321  forAll(pOwner, pFacei)
322  {
323  label own = pOwner[pFacei];
324 
325  vector vvfOwn = vvf[own];
326  vector vvfNei = psf[pFacei];
327 
328  vector gradf = (pCf[pFacei] - C[own]) & g[own];
329 
330  scalar vsfOwn = gradf & vvfOwn;
331  scalar vsfNei = gradf & vvfNei;
332 
333  scalar maxFace = max(vsfOwn, vsfNei);
334  scalar minFace = min(vsfOwn, vsfNei);
335  scalar maxMinFace = rk*(maxFace - minFace);
336  maxFace += maxMinFace;
337  minFace -= maxMinFace;
338 
339  limitFace
340  (
341  limiter[own],
342  maxFace - vsfOwn, minFace - vsfOwn,
343  magSqr(gradf)
344  );
345  }
346  }
347  }
348 
349  if (fv::debug)
350  {
351  Info<< "gradient limiter for: " << vvf.name()
352  << " max = " << gMax(limiter)
353  << " min = " << gMin(limiter)
354  << " average: " << gAverage(limiter) << endl;
355  }
356 
357  g.internalField() *= limiter;
360 
361  return tGrad;
362 }
363 
364 
365 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
366 
367 } // End namespace fv
368 
369 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
370 
371 } // End namespace Foam
372 
373 // ************************ vim: set sw=4 sts=4 et: ************************ //