cellMDLimitedGrads.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 "cellMDLimitedGrad.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(cellMDLimitedGrad)
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  scalarField maxVsf(vsf.primitiveField());
66  scalarField minVsf(vsf.primitiveField());
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::Boundary& 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  forAll(owner, facei)
133  {
134  label own = owner[facei];
135  label nei = neighbour[facei];
136 
137  // owner side
138  limitFace
139  (
140  g[own],
141  maxVsf[own],
142  minVsf[own],
143  Cf[facei] - C[own]
144  );
145 
146  // neighbour side
147  limitFace
148  (
149  g[nei],
150  maxVsf[nei],
151  minVsf[nei],
152  Cf[facei] - C[nei]
153  );
154  }
155 
156 
157  forAll(bsf, patchi)
158  {
159  const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
160  const vectorField& pCf = Cf.boundaryField()[patchi];
161 
162  forAll(pOwner, pFacei)
163  {
164  label own = pOwner[pFacei];
165 
166  limitFace
167  (
168  g[own],
169  maxVsf[own],
170  minVsf[own],
171  pCf[pFacei] - C[own]
172  );
173  }
174  }
175 
176  g.correctBoundaryConditions();
178 
179  return tGrad;
180 }
181 
182 
183 template<>
186 (
187  const volVectorField& vsf,
188  const word& name
189 ) const
190 {
191  const fvMesh& mesh = vsf.mesh();
192 
193  tmp<volTensorField> tGrad = basicGradScheme_().calcGrad(vsf, name);
194 
195  if (k_ < small)
196  {
197  return tGrad;
198  }
199 
200  volTensorField& g = tGrad.ref();
201 
202  const labelUList& owner = mesh.owner();
203  const labelUList& neighbour = mesh.neighbour();
204 
205  const volVectorField& C = mesh.C();
206  const surfaceVectorField& Cf = mesh.Cf();
207 
208  vectorField maxVsf(vsf.primitiveField());
209  vectorField minVsf(vsf.primitiveField());
210 
211  forAll(owner, facei)
212  {
213  label own = owner[facei];
214  label nei = neighbour[facei];
215 
216  const vector& vsfOwn = vsf[own];
217  const vector& vsfNei = vsf[nei];
218 
219  maxVsf[own] = max(maxVsf[own], vsfNei);
220  minVsf[own] = min(minVsf[own], vsfNei);
221 
222  maxVsf[nei] = max(maxVsf[nei], vsfOwn);
223  minVsf[nei] = min(minVsf[nei], vsfOwn);
224  }
225 
226 
227  const volVectorField::Boundary& bsf = vsf.boundaryField();
228 
229  forAll(bsf, patchi)
230  {
231  const fvPatchVectorField& psf = bsf[patchi];
232  const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
233 
234  if (psf.coupled())
235  {
236  const vectorField psfNei(psf.patchNeighbourField());
237 
238  forAll(pOwner, pFacei)
239  {
240  label own = pOwner[pFacei];
241  const vector& vsfNei = psfNei[pFacei];
242 
243  maxVsf[own] = max(maxVsf[own], vsfNei);
244  minVsf[own] = min(minVsf[own], vsfNei);
245  }
246  }
247  else
248  {
249  forAll(pOwner, pFacei)
250  {
251  label own = pOwner[pFacei];
252  const vector& vsfNei = psf[pFacei];
253 
254  maxVsf[own] = max(maxVsf[own], vsfNei);
255  minVsf[own] = min(minVsf[own], vsfNei);
256  }
257  }
258  }
259 
260  maxVsf -= vsf;
261  minVsf -= vsf;
262 
263  if (k_ < 1.0)
264  {
265  const vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
266  maxVsf += maxMinVsf;
267  minVsf -= maxMinVsf;
268 
269  // maxVsf *= 1.0/k_;
270  // minVsf *= 1.0/k_;
271  }
272 
273 
274  forAll(owner, facei)
275  {
276  label own = owner[facei];
277  label nei = neighbour[facei];
278 
279  // owner side
280  limitFace
281  (
282  g[own],
283  maxVsf[own],
284  minVsf[own],
285  Cf[facei] - C[own]
286  );
287 
288  // neighbour side
289  limitFace
290  (
291  g[nei],
292  maxVsf[nei],
293  minVsf[nei],
294  Cf[facei] - C[nei]
295  );
296  }
297 
298 
299  forAll(bsf, patchi)
300  {
301  const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
302  const vectorField& pCf = Cf.boundaryField()[patchi];
303 
304  forAll(pOwner, pFacei)
305  {
306  label own = pOwner[pFacei];
307 
308  limitFace
309  (
310  g[own],
311  maxVsf[own],
312  minVsf[own],
313  pCf[pFacei] - C[own]
314  );
315  }
316  }
317 
320 
321  return tGrad;
322 }
323 
324 
325 // ************************************************************************* //
static const Foam::dimensionedScalar C("C", Foam::dimTemperature, 234.5)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
makeFvGradScheme(cellMDLimitedGrad)
Graphite solid properties.
Definition: C.H:51
const Mesh & mesh() const
Return mesh.
Generic GeometricBoundaryField class.
Generic GeometricField class.
const Internal::FieldType & primitiveField() const
Return a const-reference to the primitive field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
void correctBoundaryConditions()
Correct boundary field.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:99
const volVectorField & C() const
Return cell centres.
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:469
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:857
const surfaceVectorField & Cf() const
Return face centres.
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:475
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:88
virtual bool coupled() const
Return true if this patch field is coupled.
Definition: fvPatchField.H:328
virtual tmp< Field< Type > > patchNeighbourField(const Pstream::commsTypes commsType=Pstream::commsTypes::blocking) const
Return patchField on the opposite patch of a coupled patch.
Definition: fvPatchField.H:433
virtual tmp< VolField< typename outerProduct< vector, Type >::type > > calcGrad(const VolField< Type > &vsf, const word &name) const
Return the gradient of the given field to the gradScheme::grad.
A class for managing temporary objects.
Definition: tmp.H:55
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
A class for handling words, derived from string.
Definition: word.H:62
volScalarField scalarField(fieldObject, mesh)
volVectorField vectorField(fieldObject, mesh)
label patchi
U correctBoundaryConditions()
VolField< vector > volVectorField
Definition: volFieldsFwd.H:65
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
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:64
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
SurfaceField< vector > surfaceVectorField
UList< label > labelUList
Definition: UList.H:65
fvPatchField< scalar > fvPatchScalarField