FreeFOAM The Cross-Platform CFD Toolkit
LimitedScheme.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 <finiteVolume/volFields.H>
28 #include <finiteVolume/fvcGrad.H>
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 template<class Type, class Limiter, template<class> class LimitFunc>
40 (
42 ) const
43 {
44  const fvMesh& mesh = this->mesh();
45 
47  (
49  (
50  IOobject
51  (
52  type() + "Limiter(" + phi.name() + ')',
53  mesh.time().timeName(),
54  mesh
55  ),
56  mesh,
57  dimless
58  )
59  );
60  surfaceScalarField& lim = tLimiter();
61 
63  tlPhi = LimitFunc<Type>()(phi);
64 
66  lPhi = tlPhi();
67 
69  gradc(fvc::grad(lPhi));
70 
71  const surfaceScalarField& CDweights = mesh.surfaceInterpolation::weights();
72 
73  const unallocLabelList& owner = mesh.owner();
74  const unallocLabelList& neighbour = mesh.neighbour();
75 
76  const vectorField& C = mesh.C();
77 
78  scalarField& pLim = lim.internalField();
79 
80  forAll(pLim, face)
81  {
82  label own = owner[face];
83  label nei = neighbour[face];
84 
85  pLim[face] = Limiter::limiter
86  (
87  CDweights[face],
88  this->faceFlux_[face],
89  lPhi[own],
90  lPhi[nei],
91  gradc[own],
92  gradc[nei],
93  C[nei] - C[own]
94  );
95  }
96 
97  surfaceScalarField::GeometricBoundaryField& bLim = lim.boundaryField();
98 
99  forAll(bLim, patchi)
100  {
101  scalarField& pLim = bLim[patchi];
102 
103  if (bLim[patchi].coupled())
104  {
105  const scalarField& pCDweights = CDweights.boundaryField()[patchi];
106  const scalarField& pFaceFlux =
107  this->faceFlux_.boundaryField()[patchi];
109  lPhi.boundaryField()[patchi].patchInternalField();
111  lPhi.boundaryField()[patchi].patchNeighbourField();
113  gradc.boundaryField()[patchi].patchInternalField();
115  gradc.boundaryField()[patchi].patchNeighbourField();
116 
117  // Build the d-vectors
118  vectorField pd =
119  mesh.Sf().boundaryField()[patchi]
120  /(
121  mesh.magSf().boundaryField()[patchi]
122  *mesh.deltaCoeffs().boundaryField()[patchi]
123  );
124 
125  if (!mesh.orthogonal())
126  {
127  pd -= mesh.correctionVectors().boundaryField()[patchi]
128  /mesh.deltaCoeffs().boundaryField()[patchi];
129  }
130 
131  forAll(pLim, face)
132  {
133  pLim[face] = Limiter::limiter
134  (
135  pCDweights[face],
136  pFaceFlux[face],
137  plPhiP[face],
138  plPhiN[face],
139  pGradcP[face],
140  pGradcN[face],
141  pd[face]
142  );
143  }
144  }
145  else
146  {
147  pLim = 1.0;
148  }
149  }
150 
151  return tLimiter;
152 }
153 
154 
155 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
156 
157 } // End namespace Foam
158 
159 // ************************ vim: set sw=4 sts=4 et: ************************ //