FreeFOAM The Cross-Platform CFD Toolkit
cellLimitedGrads.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 "cellLimitedGrad.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<>
51 inline void cellLimitedGrad<scalar>::limitFace
52 (
53  scalar& limiter,
54  const scalar& maxDelta,
55  const scalar& minDelta,
56  const scalar& extrapolate
57 )
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 template<class Type>
71 (
72  Type& limiter,
73  const Type& maxDelta,
74  const Type& minDelta,
75  const Type& extrapolate
76 )
77 {
78  for(direction cmpt=0; cmpt<Type::nComponents; cmpt++)
79  {
81  (
82  limiter.component(cmpt),
83  maxDelta.component(cmpt),
84  minDelta.component(cmpt),
85  extrapolate.component(cmpt)
86  );
87  }
88 }
89 
90 
91 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
92 
93 template<>
95 (
96  const volScalarField& vsf
97 ) const
98 {
99  const fvMesh& mesh = vsf.mesh();
100 
101  tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
102 
103  if (k_ < SMALL)
104  {
105  return tGrad;
106  }
107 
108  volVectorField& g = tGrad();
109 
110  const unallocLabelList& owner = mesh.owner();
111  const unallocLabelList& neighbour = mesh.neighbour();
112 
113  const volVectorField& C = mesh.C();
114  const surfaceVectorField& Cf = mesh.Cf();
115 
116  scalarField maxVsf(vsf.internalField());
117  scalarField minVsf(vsf.internalField());
118 
119  forAll(owner, facei)
120  {
121  label own = owner[facei];
122  label nei = neighbour[facei];
123 
124  scalar vsfOwn = vsf[own];
125  scalar vsfNei = vsf[nei];
126 
127  maxVsf[own] = max(maxVsf[own], vsfNei);
128  minVsf[own] = min(minVsf[own], vsfNei);
129 
130  maxVsf[nei] = max(maxVsf[nei], vsfOwn);
131  minVsf[nei] = min(minVsf[nei], vsfOwn);
132  }
133 
134 
135  const volScalarField::GeometricBoundaryField& bsf = vsf.boundaryField();
136 
137  forAll(bsf, patchi)
138  {
139  const fvPatchScalarField& psf = bsf[patchi];
140 
141  const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
142 
143  if (psf.coupled())
144  {
145  scalarField psfNei = psf.patchNeighbourField();
146 
147  forAll(pOwner, pFacei)
148  {
149  label own = pOwner[pFacei];
150  scalar vsfNei = psfNei[pFacei];
151 
152  maxVsf[own] = max(maxVsf[own], vsfNei);
153  minVsf[own] = min(minVsf[own], vsfNei);
154  }
155  }
156  else
157  {
158  forAll(pOwner, pFacei)
159  {
160  label own = pOwner[pFacei];
161  scalar vsfNei = psf[pFacei];
162 
163  maxVsf[own] = max(maxVsf[own], vsfNei);
164  minVsf[own] = min(minVsf[own], vsfNei);
165  }
166  }
167  }
168 
169  maxVsf -= vsf;
170  minVsf -= vsf;
171 
172  if (k_ < 1.0)
173  {
174  scalarField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf);
175  maxVsf += maxMinVsf;
176  minVsf -= maxMinVsf;
177 
178  //maxVsf *= 1.0/k_;
179  //minVsf *= 1.0/k_;
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  limitFace
193  (
194  limiter[own],
195  maxVsf[own],
196  minVsf[own],
197  (Cf[facei] - C[own]) & g[own]
198  );
199 
200  // neighbour side
201  limitFace
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 unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
213  const vectorField& pCf = Cf.boundaryField()[patchi];
214 
215  forAll(pOwner, pFacei)
216  {
217  label own = pOwner[pFacei];
218 
219  limitFace
220  (
221  limiter[own],
222  maxVsf[own],
223  minVsf[own],
224  (pCf[pFacei] - C[own]) & g[own]
225  );
226  }
227  }
228 
229  if (fv::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.internalField() *= limiter;
240 
241  return tGrad;
242 }
243 
244 
245 template<>
247 (
248  const volVectorField& vsf
249 ) const
250 {
251  const fvMesh& mesh = vsf.mesh();
252 
253  tmp<volTensorField> tGrad = basicGradScheme_().grad(vsf);
254 
255  if (k_ < SMALL)
256  {
257  return tGrad;
258  }
259 
260  volTensorField& g = tGrad();
261 
262  const unallocLabelList& owner = mesh.owner();
263  const unallocLabelList& neighbour = mesh.neighbour();
264 
265  const volVectorField& C = mesh.C();
266  const surfaceVectorField& Cf = mesh.Cf();
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 volVectorField::GeometricBoundaryField& bsf = vsf.boundaryField();
288 
289  forAll(bsf, patchi)
290  {
291  const fvPatchVectorField& psf = bsf[patchi];
292  const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
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
336 
337  forAll(owner, facei)
338  {
339  label own = owner[facei];
340  label nei = neighbour[facei];
341 
342  // owner side
343  limitFace
344  (
345  limiter[own],
346  maxVsf[own],
347  minVsf[own],
348  (Cf[facei] - C[own]) & g[own]
349  );
350 
351  // neighbour side
352  limitFace
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 unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
364  const vectorField& pCf = Cf.boundaryField()[patchi];
365 
366  forAll(pOwner, pFacei)
367  {
368  label own = pOwner[pFacei];
369 
370  limitFace
371  (
372  limiter[own],
373  maxVsf[own],
374  minVsf[own],
375  ((pCf[pFacei] - C[own]) & g[own])
376  );
377  }
378  }
379 
380  if (fv::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.internalField();
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 
402 
403  return tGrad;
404 }
405 
406 
407 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
408 
409 } // End namespace fv
410 
411 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
412 
413 } // End namespace Foam
414 
415 // ************************ vim: set sw=4 sts=4 et: ************************ //