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) 2011-2015 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 "cellLimitedGrad.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(cellLimitedGrad)
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 template<>
42 Foam::fv::cellLimitedGrad<Foam::scalar>::calcGrad
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();
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  scalarField maxVsf(vsf.internalField());
66  scalarField minVsf(vsf.internalField());
67 
68  forAll(owner, facei)
69  {
70  label own = owner[facei];
71  label nei = neighbour[facei];
72 
73  scalar vsfOwn = vsf[own];
74  scalar vsfNei = vsf[nei];
75 
76  maxVsf[own] = max(maxVsf[own], vsfNei);
77  minVsf[own] = min(minVsf[own], vsfNei);
78 
79  maxVsf[nei] = max(maxVsf[nei], vsfOwn);
80  minVsf[nei] = min(minVsf[nei], vsfOwn);
81  }
82 
83 
84  const volScalarField::GeometricBoundaryField& bsf = vsf.boundaryField();
85 
86  forAll(bsf, patchi)
87  {
88  const fvPatchScalarField& psf = bsf[patchi];
89 
90  const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
91 
92  if (psf.coupled())
93  {
94  const scalarField psfNei(psf.patchNeighbourField());
95 
96  forAll(pOwner, pFacei)
97  {
98  label own = pOwner[pFacei];
99  scalar vsfNei = psfNei[pFacei];
100 
101  maxVsf[own] = max(maxVsf[own], vsfNei);
102  minVsf[own] = min(minVsf[own], vsfNei);
103  }
104  }
105  else
106  {
107  forAll(pOwner, pFacei)
108  {
109  label own = pOwner[pFacei];
110  scalar vsfNei = psf[pFacei];
111 
112  maxVsf[own] = max(maxVsf[own], vsfNei);
113  minVsf[own] = min(minVsf[own], vsfNei);
114  }
115  }
116  }
117 
118  maxVsf -= vsf;
119  minVsf -= vsf;
120 
121  if (k_ < 1.0)
122  {
123  const scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
124  maxVsf += maxMinVsf;
125  minVsf -= maxMinVsf;
126 
127  //maxVsf *= 1.0/k_;
128  //minVsf *= 1.0/k_;
129  }
130 
131 
132  // create limiter
133  scalarField limiter(vsf.internalField().size(), 1.0);
134 
135  forAll(owner, facei)
136  {
137  label own = owner[facei];
138  label nei = neighbour[facei];
139 
140  // owner side
141  limitFace
142  (
143  limiter[own],
144  maxVsf[own],
145  minVsf[own],
146  (Cf[facei] - C[own]) & g[own]
147  );
148 
149  // neighbour side
150  limitFace
151  (
152  limiter[nei],
153  maxVsf[nei],
154  minVsf[nei],
155  (Cf[facei] - C[nei]) & g[nei]
156  );
157  }
158 
159  forAll(bsf, patchi)
160  {
161  const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
162  const vectorField& pCf = Cf.boundaryField()[patchi];
163 
164  forAll(pOwner, pFacei)
165  {
166  label own = pOwner[pFacei];
167 
168  limitFace
169  (
170  limiter[own],
171  maxVsf[own],
172  minVsf[own],
173  (pCf[pFacei] - C[own]) & g[own]
174  );
175  }
176  }
177 
178  if (fv::debug)
179  {
180  Info<< "gradient limiter for: " << vsf.name()
181  << " max = " << gMax(limiter)
182  << " min = " << gMin(limiter)
183  << " average: " << gAverage(limiter) << endl;
184  }
185 
186  g.internalField() *= limiter;
187  g.correctBoundaryConditions();
189 
190  return tGrad;
191 }
192 
193 
194 template<>
197 (
198  const volVectorField& vsf,
199  const word& name
200 ) const
201 {
202  const fvMesh& mesh = vsf.mesh();
203 
204  tmp<volTensorField> tGrad = basicGradScheme_().calcGrad(vsf, name);
205 
206  if (k_ < SMALL)
207  {
208  return tGrad;
209  }
210 
211  volTensorField& g = tGrad();
212 
213  const labelUList& owner = mesh.owner();
214  const labelUList& neighbour = mesh.neighbour();
215 
216  const volVectorField& C = mesh.C();
217  const surfaceVectorField& Cf = mesh.Cf();
218 
219  vectorField maxVsf(vsf.internalField());
220  vectorField minVsf(vsf.internalField());
221 
222  forAll(owner, facei)
223  {
224  label own = owner[facei];
225  label nei = neighbour[facei];
226 
227  const vector& vsfOwn = vsf[own];
228  const vector& vsfNei = vsf[nei];
229 
230  maxVsf[own] = max(maxVsf[own], vsfNei);
231  minVsf[own] = min(minVsf[own], vsfNei);
232 
233  maxVsf[nei] = max(maxVsf[nei], vsfOwn);
234  minVsf[nei] = min(minVsf[nei], vsfOwn);
235  }
236 
237 
238  const volVectorField::GeometricBoundaryField& bsf = vsf.boundaryField();
239 
240  forAll(bsf, patchi)
241  {
242  const fvPatchVectorField& psf = bsf[patchi];
243  const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
244 
245  if (psf.coupled())
246  {
247  const vectorField psfNei(psf.patchNeighbourField());
248 
249  forAll(pOwner, pFacei)
250  {
251  label own = pOwner[pFacei];
252  const vector& vsfNei = psfNei[pFacei];
253 
254  maxVsf[own] = max(maxVsf[own], vsfNei);
255  minVsf[own] = min(minVsf[own], vsfNei);
256  }
257  }
258  else
259  {
260  forAll(pOwner, pFacei)
261  {
262  label own = pOwner[pFacei];
263  const vector& vsfNei = psf[pFacei];
264 
265  maxVsf[own] = max(maxVsf[own], vsfNei);
266  minVsf[own] = min(minVsf[own], vsfNei);
267  }
268  }
269  }
270 
271  maxVsf -= vsf;
272  minVsf -= vsf;
273 
274  if (k_ < 1.0)
275  {
276  const vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
277  maxVsf += maxMinVsf;
278  minVsf -= maxMinVsf;
279 
280  //maxVsf *= 1.0/k_;
281  //minVsf *= 1.0/k_;
282  }
283 
284 
285  // create limiter
286  vectorField limiter(vsf.internalField().size(), vector::one);
287 
288  forAll(owner, facei)
289  {
290  label own = owner[facei];
291  label nei = neighbour[facei];
292 
293  // owner side
294  limitFace
295  (
296  limiter[own],
297  maxVsf[own],
298  minVsf[own],
299  (Cf[facei] - C[own]) & g[own]
300  );
301 
302  // neighbour side
303  limitFace
304  (
305  limiter[nei],
306  maxVsf[nei],
307  minVsf[nei],
308  (Cf[facei] - C[nei]) & g[nei]
309  );
310  }
311 
312  forAll(bsf, patchi)
313  {
314  const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
315  const vectorField& pCf = Cf.boundaryField()[patchi];
316 
317  forAll(pOwner, pFacei)
318  {
319  label own = pOwner[pFacei];
320 
321  limitFace
322  (
323  limiter[own],
324  maxVsf[own],
325  minVsf[own],
326  ((pCf[pFacei] - C[own]) & g[own])
327  );
328  }
329  }
330 
331  if (fv::debug)
332  {
333  Info<< "gradient limiter for: " << vsf.name()
334  << " max = " << gMax(limiter)
335  << " min = " << gMin(limiter)
336  << " average: " << gAverage(limiter) << endl;
337  }
338 
339  tensorField& gIf = g.internalField();
340 
341  forAll(gIf, celli)
342  {
343  gIf[celli] = tensor
344  (
345  cmptMultiply(limiter[celli], gIf[celli].x()),
346  cmptMultiply(limiter[celli], gIf[celli].y()),
347  cmptMultiply(limiter[celli], gIf[celli].z())
348  );
349  }
350 
353 
354  return tGrad;
355 }
356 
357 
358 // ************************************************************************* //
UList< label > labelUList
Definition: UList.H:63
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
A class for handling words, derived from string.
Definition: word.H:59
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
InternalField & internalField()
Return internal field.
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:282
Graphite solid properties.
Definition: C.H:57
messageStream Info
const Mesh & mesh() const
Return mesh.
dynamicFvMesh & mesh
const surfaceVectorField & Cf() const
Return face centres as surfaceVectorField.
virtual bool coupled() const
Return true if this patch field is coupled.
Definition: fvPatchField.H:342
Namespace for OpenFOAM.
virtual tmp< Field< Type > > patchNeighbourField() const
Return patchField on the opposite patch of a coupled patch.
Definition: fvPatchField.H:411
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const volVectorField & C() const
Return cell centres as volVectorField.
volVectorField vectorField(fieldObject, mesh)
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
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.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
label patchi
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:65
fvPatchField< scalar > fvPatchScalarField
const word & name() const
Return name.
Definition: IOobject.H:260
Type gMin(const FieldField< Field, Type > &f)
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
const dimensionedVector & g
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
scalar y
void correctBoundaryConditions()
Correct boundary field.
makeFvGradScheme(cellLimitedGrad)
Type gAverage(const FieldField< Field, Type > &f)
labelList fv(nPoints)
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
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:288
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
U correctBoundaryConditions()
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 scalar psiMax, const scalar psiMin)
volScalarField scalarField(fieldObject, mesh)
volScalarField & C
Type gMax(const FieldField< Field, Type > &f)
A class for managing temporary objects.
Definition: PtrList.H:118
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:552