faceLimitedGrads.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
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 "gaussGrad.H"
28 #include "fvMesh.H"
29 #include "volMesh.H"
30 #include "surfaceMesh.H"
31 #include "volFields.H"
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 makeFvGradScheme(faceLimitedGrad)
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 template<>
43 (
44  const volScalarField& vsf,
45  const word& name
46 ) const
47 {
48  const fvMesh& mesh = vsf.mesh();
49 
50  tmp<volVectorField> tGrad = basicGradScheme_().calcGrad(vsf, name);
51 
52  if (k_ < small)
53  {
54  return tGrad;
55  }
56 
57  volVectorField& g = tGrad.ref();
58 
59  const labelUList& owner = mesh.owner();
60  const labelUList& neighbour = mesh.neighbour();
61 
62  const volVectorField& C = mesh.C();
63  const surfaceVectorField& Cf = mesh.Cf();
64 
65  // create limiter
66  scalarField limiter(vsf.primitiveField().size(), 1.0);
67 
68  scalar rk = (1.0/k_ - 1.0);
69 
70  forAll(owner, facei)
71  {
72  label own = owner[facei];
73  label nei = neighbour[facei];
74 
75  scalar vsfOwn = vsf[own];
76  scalar vsfNei = vsf[nei];
77 
78  scalar maxFace = max(vsfOwn, vsfNei);
79  scalar minFace = min(vsfOwn, vsfNei);
80  scalar maxMinFace = rk*(maxFace - minFace);
81  maxFace += maxMinFace;
82  minFace -= maxMinFace;
83 
84  // owner side
85  limitFace
86  (
87  limiter[own],
88  maxFace - vsfOwn, minFace - vsfOwn,
89  (Cf[facei] - C[own]) & g[own]
90  );
91 
92  // neighbour side
93  limitFace
94  (
95  limiter[nei],
96  maxFace - vsfNei, minFace - vsfNei,
97  (Cf[facei] - C[nei]) & g[nei]
98  );
99  }
100 
101  const volScalarField::Boundary& bsf = vsf.boundaryField();
102 
103  forAll(bsf, patchi)
104  {
105  const fvPatchScalarField& psf = bsf[patchi];
106 
107  const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
108  const vectorField& pCf = Cf.boundaryField()[patchi];
109 
110  if (psf.coupled())
111  {
112  const scalarField psfNei(psf.patchNeighbourField());
113 
114  forAll(pOwner, pFacei)
115  {
116  label own = pOwner[pFacei];
117 
118  scalar vsfOwn = vsf[own];
119  scalar vsfNei = psfNei[pFacei];
120 
121  scalar maxFace = max(vsfOwn, vsfNei);
122  scalar minFace = min(vsfOwn, vsfNei);
123  scalar maxMinFace = rk*(maxFace - minFace);
124  maxFace += maxMinFace;
125  minFace -= maxMinFace;
126 
127  limitFace
128  (
129  limiter[own],
130  maxFace - vsfOwn, minFace - vsfOwn,
131  (pCf[pFacei] - C[own]) & g[own]
132  );
133  }
134  }
135  else if (psf.fixesValue())
136  {
137  forAll(pOwner, pFacei)
138  {
139  label own = pOwner[pFacei];
140 
141  scalar vsfOwn = vsf[own];
142  scalar vsfNei = psf[pFacei];
143 
144  scalar maxFace = max(vsfOwn, vsfNei);
145  scalar minFace = min(vsfOwn, vsfNei);
146  scalar maxMinFace = rk*(maxFace - minFace);
147  maxFace += maxMinFace;
148  minFace -= maxMinFace;
149 
150  limitFace
151  (
152  limiter[own],
153  maxFace - vsfOwn, minFace - vsfOwn,
154  (pCf[pFacei] - C[own]) & g[own]
155  );
156  }
157  }
158  }
159 
160  if (fv::debug)
161  {
162  Info<< "gradient limiter for: " << vsf.name()
163  << " max = " << gMax(limiter)
164  << " min = " << gMin(limiter)
165  << " average: " << gAverage(limiter) << endl;
166  }
167 
168  g.primitiveFieldRef() *= limiter;
169  g.correctBoundaryConditions();
171 
172  return tGrad;
173 }
174 
175 
176 template<>
179 (
180  const volVectorField& vvf,
181  const word& name
182 ) const
183 {
184  const fvMesh& mesh = vvf.mesh();
185 
186  tmp<volTensorField> tGrad = basicGradScheme_().calcGrad(vvf, name);
187 
188  if (k_ < small)
189  {
190  return tGrad;
191  }
192 
193  volTensorField& g = tGrad.ref();
194 
195  const labelUList& owner = mesh.owner();
196  const labelUList& neighbour = mesh.neighbour();
197 
198  const volVectorField& C = mesh.C();
199  const surfaceVectorField& Cf = mesh.Cf();
200 
201  // create limiter
202  scalarField limiter(vvf.primitiveField().size(), 1.0);
203 
204  scalar rk = (1.0/k_ - 1.0);
205 
206  forAll(owner, facei)
207  {
208  label own = owner[facei];
209  label nei = neighbour[facei];
210 
211  vector vvfOwn = vvf[own];
212  vector vvfNei = vvf[nei];
213 
214  // owner side
215  vector gradf = (Cf[facei] - C[own]) & g[own];
216 
217  scalar vsfOwn = gradf & vvfOwn;
218  scalar vsfNei = gradf & vvfNei;
219 
220  scalar maxFace = max(vsfOwn, vsfNei);
221  scalar minFace = min(vsfOwn, vsfNei);
222  scalar maxMinFace = rk*(maxFace - minFace);
223  maxFace += maxMinFace;
224  minFace -= maxMinFace;
225 
226  limitFace
227  (
228  limiter[own],
229  maxFace - vsfOwn, minFace - vsfOwn,
230  magSqr(gradf)
231  );
232 
233 
234  // neighbour side
235  gradf = (Cf[facei] - C[nei]) & g[nei];
236 
237  vsfOwn = gradf & vvfOwn;
238  vsfNei = gradf & vvfNei;
239 
240  maxFace = max(vsfOwn, vsfNei);
241  minFace = min(vsfOwn, vsfNei);
242 
243  limitFace
244  (
245  limiter[nei],
246  maxFace - vsfNei, minFace - vsfNei,
247  magSqr(gradf)
248  );
249  }
250 
251 
252  const volVectorField::Boundary& bvf = vvf.boundaryField();
253 
254  forAll(bvf, patchi)
255  {
256  const fvPatchVectorField& psf = bvf[patchi];
257 
258  const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
259  const vectorField& pCf = Cf.boundaryField()[patchi];
260 
261  if (psf.coupled())
262  {
263  const vectorField psfNei(psf.patchNeighbourField());
264 
265  forAll(pOwner, pFacei)
266  {
267  label own = pOwner[pFacei];
268 
269  vector vvfOwn = vvf[own];
270  vector vvfNei = psfNei[pFacei];
271 
272  vector gradf = (pCf[pFacei] - C[own]) & g[own];
273 
274  scalar vsfOwn = gradf & vvfOwn;
275  scalar vsfNei = gradf & vvfNei;
276 
277  scalar maxFace = max(vsfOwn, vsfNei);
278  scalar minFace = min(vsfOwn, vsfNei);
279  scalar maxMinFace = rk*(maxFace - minFace);
280  maxFace += maxMinFace;
281  minFace -= maxMinFace;
282 
283  limitFace
284  (
285  limiter[own],
286  maxFace - vsfOwn, minFace - vsfOwn,
287  magSqr(gradf)
288  );
289  }
290  }
291  else if (psf.fixesValue())
292  {
293  forAll(pOwner, pFacei)
294  {
295  label own = pOwner[pFacei];
296 
297  vector vvfOwn = vvf[own];
298  vector vvfNei = psf[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  }
320 
321  if (fv::debug)
322  {
323  Info<< "gradient limiter for: " << vvf.name()
324  << " max = " << gMax(limiter)
325  << " min = " << gMin(limiter)
326  << " average: " << gAverage(limiter) << endl;
327  }
328 
329  g.primitiveFieldRef() *= limiter;
332 
333  return tGrad;
334 }
335 
336 
337 // ************************************************************************* //
Graphite solid properties.
Definition: C.H:48
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const word & name() const
Return name.
Definition: IOobject.H:303
rDeltaT correctBoundaryConditions()
virtual bool fixesValue() const
Return true if this patch field fixes a value.
Definition: fvPatchField.H:296
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const surfaceVectorField & Cf() const
Return face centres as surfaceVectorField.
Type gMin(const FieldField< Field, Type > &f)
const Boundary & boundaryField() const
Return const-reference to the boundary field.
makeFvGradScheme(faceLimitedGrad)
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:66
volVectorField vectorField(fieldObject, mesh)
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:58
UList< label > labelUList
Definition: UList.H:65
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:284
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
dynamicFvMesh & mesh
fvPatchField< scalar > fvPatchScalarField
void limiter(scalarField &allLambda, const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phiBD, const surfaceScalarField &phiCorr, const SpType &Sp, const SuType &Su, const PsiMaxType &psiMax, const PsiMinType &psiMin)
A class for handling words, derived from string.
Definition: word.H:59
virtual tmp< Field< Type > > patchNeighbourField() const
Return patchField on the opposite patch of a coupled patch.
Definition: fvPatchField.H:411
virtual bool coupled() const
Return true if this patch field is coupled.
Definition: fvPatchField.H:309
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
dimensioned< scalar > magSqr(const dimensioned< Type > &)
volScalarField scalarField(fieldObject, mesh)
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
const Mesh & mesh() const
Return mesh.
Type gMax(const FieldField< Field, Type > &f)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:278
label patchi
Type gAverage(const FieldField< Field, Type > &f)
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
void correctBoundaryConditions()
Correct boundary field.
messageStream Info
virtual tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > calcGrad(const GeometricField< Type, fvPatchField, volMesh > &vsf, const word &name) const
Return the gradient of the given field to the gradScheme::grad.
const volVectorField & C() const
Return cell centres as volVectorField.
A class for managing temporary objects.
Definition: PtrList.H:53
const dimensionedVector & g
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:540