src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "volFields.H"
  28 #include "surfaceFields.H"
  29 #include "fvcGrad.H"
  30 #include "coupledFvPatchFields.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 template<class Type, class Limiter, template<class> class LimitFunc>
  40 tmp<surfaceScalarField> LimitedScheme<Type, Limiter, LimitFunc>::limiter
  41 (
  42     const GeometricField<Type, fvPatchField, volMesh>& phi
  43 ) const
  44 {
  45     const fvMesh& mesh = this->mesh();
  46
  47     tmp<surfaceScalarField> tLimiter
  48     (
  49         new surfaceScalarField
  50         (
  51             IOobject
  52             (
  53                 type() + "Limiter(" + phi.name() + ')',
  54                 mesh.time().timeName(),
  55                 mesh
  56             ),
  57             mesh,
  58             dimless
  59         )
  60     );
  61     surfaceScalarField& lim = tLimiter();
  62
  63     tmp<GeometricField<typename Limiter::phiType, fvPatchField, volMesh> >
  64         tlPhi = LimitFunc<Type>()(phi);
  65
  66     const GeometricField<typename Limiter::phiType, fvPatchField, volMesh>&
  67         lPhi = tlPhi();
  68
  69     GeometricField<typename Limiter::gradPhiType, fvPatchField, volMesh>
  70         gradc(fvc::grad(lPhi));
  71
  72     const surfaceScalarField& CDweights = mesh.surfaceInterpolation::weights();
  73
  74     const unallocLabelList& owner = mesh.owner();
  75     const unallocLabelList& neighbour = mesh.neighbour();
  76
  77     const vectorField& C = mesh.C();
  78
  79     scalarField& pLim = lim.internalField();
  80
  81     forAll(pLim, face)
  82     {
  83         label own = owner[face];
  84         label nei = neighbour[face];
  85
  86         pLim[face] = Limiter::limiter
  87         (
  88             CDweights[face],
  89             this->faceFlux_[face],
  90             lPhi[own],
  91             lPhi[nei],
  92             gradc[own],
  93             gradc[nei],
  94             C[nei] - C[own]
  95         );
  96     }
  97
  98     surfaceScalarField::GeometricBoundaryField& bLim = lim.boundaryField();
  99
 100     forAll(bLim, patchi)
 101     {
 102         scalarField& pLim = bLim[patchi];
 103
 104         if (bLim[patchi].coupled())
 105         {
 106             const scalarField& pCDweights = CDweights.boundaryField()[patchi];
 107             const scalarField& pFaceFlux =
 108                 this->faceFlux_.boundaryField()[patchi];
 109             Field<typename Limiter::phiType> plPhiP =
 110                 lPhi.boundaryField()[patchi].patchInternalField();
 111             Field<typename Limiter::phiType> plPhiN =
 112                 lPhi.boundaryField()[patchi].patchNeighbourField();
 113             Field<typename Limiter::gradPhiType> pGradcP =
 114                 gradc.boundaryField()[patchi].patchInternalField();
 115             Field<typename Limiter::gradPhiType> pGradcN =
 116                 gradc.boundaryField()[patchi].patchNeighbourField();
 117
 118             // Build the d-vectors
 119             vectorField pd =
 120                 mesh.Sf().boundaryField()[patchi]
 121                /(
 122                    mesh.magSf().boundaryField()[patchi]
 123                   *mesh.deltaCoeffs().boundaryField()[patchi]
 124                 );
 125
 126             if (!mesh.orthogonal())
 127             {
 128                 pd -= mesh.correctionVectors().boundaryField()[patchi]
 129                     /mesh.deltaCoeffs().boundaryField()[patchi];
 130             }
 131
 132             forAll(pLim, face)
 133             {
 134                 pLim[face] = Limiter::limiter
 135                 (
 136                     pCDweights[face],
 137                     pFaceFlux[face],
 138                     plPhiP[face],
 139                     plPhiN[face],
 140                     pGradcP[face],
 141                     pGradcN[face],
 142                     pd[face]
 143                 );
 144             }
 145         }
 146         else
 147         {
 148             pLim = 1.0;
 149         }
 150     }
 151
 152     return tLimiter;
 153 }
 154
 155
 156 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 157
 158 } // End namespace Foam
 159
 160 // ************************************************************************* //