faNVDscheme.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2016-2017 Wikki Ltd
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "areaFields.H"
29 #include "edgeFields.H"
30 #include "facGrad.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
38 {
39  return phi;
40 }
41 
43 {
44  return magSqr(phi);
45 }
46 
47 inline tmp<areaScalarField> limiter(const areaTensorField& phi)
48 {
49  return magSqr(phi);
50 }
51 
52 }
53 
54 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
55 
56 template<class Type, class NVDweight>
58 (
60 ) const
61 {
62  const faMesh& mesh = this->mesh();
63 
64  tmp<edgeScalarField> tWeightingFactors
65  (
66  new edgeScalarField(mesh.edgeInterpolation::weights())
67  );
68  edgeScalarField& weightingFactors = tWeightingFactors.ref();
69 
70  scalarField& weights = weightingFactors.primitiveFieldRef();
71 
73  const areaScalarField& vf = tvf();
74 
75  const areaVectorField gradc(fac::grad(vf));
76 
77 // edgeVectorField d
78 // (
79 // mesh.Le()
80 // /(mesh.magLe()*mesh.edgeInterpolation::deltaCoeffs())
81 // );
82 
83 // if (!mesh.orthogonal())
84 // {
85 // d -=
86 // mesh.edgeInterpolation::correctionVectors()
87 // /mesh.edgeInterpolation::deltaCoeffs();
88 // }
89 
90  const labelUList& owner = mesh.owner();
91  const labelUList& neighbour = mesh.neighbour();
92  const vectorField& n = mesh.faceAreaNormals().internalField();
93  const vectorField& c = mesh.areaCentres().internalField();
94 
95  forAll(weights, edge)
96  {
97  vector d(Zero);
98 
99  if (edgeFlux_[edge] > 0)
100  {
101  d = c[neighbour[edge]] - c[owner[edge]];
102  d -= n[owner[edge]]*(n[owner[edge]]&d);
103  d /= mag(d)/mesh.edgeInterpolation::lPN().internalField()[edge];
104  }
105  else
106  {
107  d = c[neighbour[edge]] - c[owner[edge]];
108  d -= n[neighbour[edge]]*(n[neighbour[edge]]&d);
109  d /= mag(d)/mesh.edgeInterpolation::lPN().internalField()[edge];
110  }
111 
112  weights[edge] =
113  this->weight
114  (
115  weights[edge],
116  edgeFlux_[edge],
117  vf[owner[edge]],
118  vf[neighbour[edge]],
119  gradc[owner[edge]],
120  gradc[neighbour[edge]],
121  d
122  );
123  }
124 
125 
127  bWeights = weightingFactors.boundaryFieldRef();
128 
129  forAll(bWeights, patchI)
130  {
131  if (bWeights[patchI].coupled())
132  {
133  scalarField& pWeights = bWeights[patchI];
134 
135  const scalarField& pEdgeFlux = edgeFlux_.boundaryField()[patchI];
136 
137  scalarField pVfP(vf.boundaryField()[patchI].patchInternalField());
138 
139  scalarField pVfN(vf.boundaryField()[patchI].patchNeighbourField());
140 
141  vectorField pGradcP
142  (
143  gradc.boundaryField()[patchI].patchInternalField()
144  );
145 
146  vectorField pGradcN
147  (
148  gradc.boundaryField()[patchI].patchNeighbourField()
149  );
150 
151  vectorField CP
152  (
153  mesh.areaCentres().boundaryField()[patchI].patchInternalField()
154  );
155 
156  vectorField CN
157  (
158  mesh.areaCentres().boundaryField()[patchI]
159  .patchNeighbourField()
160  );
161 
162  vectorField nP
163  (
164  mesh.faceAreaNormals().boundaryField()[patchI]
165  .patchInternalField()
166  );
167 
168  vectorField nN
169  (
170  mesh.faceAreaNormals().boundaryField()[patchI]
171  .patchNeighbourField()
172  );
173 
174  scalarField pLPN
175  (
176  mesh.edgeInterpolation::lPN().boundaryField()[patchI]
177  );
178 
179  forAll(pWeights, edgeI)
180  {
181  vector d(Zero);
182 
183  if (pEdgeFlux[edgeI] > 0)
184  {
185  d = CN[edgeI] - CP[edgeI];
186  d -= nP[edgeI]*(nP[edgeI]&d);
187  d /= mag(d)/pLPN[edgeI];
188  }
189  else
190  {
191  d = CN[edgeI] - CP[edgeI];
192  d -= nN[edgeI]*(nN[edgeI]&d);
193  d /= mag(d)/pLPN[edgeI];
194  }
195 
196  pWeights[edgeI] =
197  this->weight
198  (
199  pWeights[edgeI],
200  pEdgeFlux[edgeI],
201  pVfP[edgeI],
202  pVfN[edgeI],
203  pGradcP[edgeI],
204  pGradcN[edgeI],
205  d
206  );
207  }
208  }
209  }
210 
211  return tWeightingFactors;
212 }
213 
214 
215 // ************************************************************************* //
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::edge
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:63
Foam::fac::grad
tmp< GeometricField< typename outerProduct< vector, Type >::type, faPatchField, areaMesh >> grad(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
Definition: facGrad.C:56
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::faNVDscheme::weights
virtual tmp< edgeScalarField > weights(const GeometricField< Type, faPatchField, areaMesh > &) const
Return the interpolation weighting factors.
Definition: faNVDscheme.C:58
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
Foam::Field< scalar >
edgeFields.H
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::edgeScalarField
GeometricField< scalar, faePatchField, edgeMesh > edgeScalarField
Definition: edgeFieldsFwd.H:52
areaFields.H
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::GeometricField::primitiveFieldRef
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
Definition: GeometricField.C:766
Foam::limiter
tmp< areaScalarField > limiter(const areaVectorField &phi)
Definition: faNVDscheme.C:42
facGrad.H
Calculate the gradient of the given field.
Foam::GeometricField::boundaryFieldRef
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
Definition: GeometricField.C:783
Foam::Vector< scalar >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::UList< label >
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::faMesh
Finite area mesh. Used for 2-D non-Euclidian finite area method.
Definition: faMesh.H:82
coupled
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
Foam::GeometricField< scalar, faPatchField, areaMesh >
Foam::limiter
tmp< areaScalarField > limiter(const areaScalarField &phi)
Definition: faNVDscheme.C:37
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:62
Foam::areaTensorField
GeometricField< tensor, faPatchField, areaMesh > areaTensorField
Definition: areaFieldsFwd.H:62