src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/PhiScheme/PhiScheme.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 #include "surfaceInterpolate.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 template<class Type, class PhiLimiter>
  41 tmp<surfaceScalarField> PhiScheme<Type, PhiLimiter>::limiter
  42 (
  43     const GeometricField<Type, fvPatchField, volMesh>& phi
  44 ) const
  45 {
  46     const fvMesh& mesh = this->mesh();
  47
  48     tmp<surfaceScalarField> tLimiter
  49     (
  50         new surfaceScalarField
  51         (
  52             IOobject
  53             (
  54                 "PhiLimiter",
  55                 mesh.time().timeName(),
  56                 mesh
  57             ),
  58             mesh,
  59             dimless
  60         )
  61     );
  62     surfaceScalarField& Limiter = tLimiter();
  63
  64     const surfaceScalarField& CDweights = mesh.surfaceInterpolation::weights();
  65
  66     const surfaceVectorField& Sf = mesh.Sf();
  67     const surfaceScalarField& magSf = mesh.magSf();
  68
  69     const unallocLabelList& owner = mesh.owner();
  70     const unallocLabelList& neighbour = mesh.neighbour();
  71
  72     tmp<surfaceScalarField> tUflux = this->faceFlux_;
  73
  74     if (this->faceFlux_.dimensions() == dimDensity*dimVelocity*dimArea)
  75     {
  76         const volScalarField& rho =
  77             phi.db().objectRegistry::lookupObject<volScalarField>("rho");
  78         tUflux = this->faceFlux_/fvc::interpolate(rho);
  79     }
  80     else if (this->faceFlux_.dimensions() != dimVelocity*dimArea)
  81     {
  82         FatalErrorIn
  83         (
  84             "PhiScheme<PhiLimiter>::limiter"
  85             "(const GeometricField<Type, fvPatchField, volMesh>& phi)"
  86         )   << "dimensions of faceFlux are not correct"
  87             << exit(FatalError);
  88     }
  89
  90     const surfaceScalarField& Uflux = tUflux();
  91
  92     scalarField& pLimiter = Limiter.internalField();
  93
  94     forAll(pLimiter, face)
  95     {
  96         pLimiter[face] = PhiLimiter::limiter
  97         (
  98             CDweights[face],
  99             Uflux[face],
 100             phi[owner[face]],
 101             phi[neighbour[face]],
 102             Sf[face],
 103             magSf[face]
 104         );
 105     }
 106
 107
 108     surfaceScalarField::GeometricBoundaryField& bLimiter =
 109         Limiter.boundaryField();
 110
 111     forAll(bLimiter, patchI)
 112     {
 113         scalarField& pLimiter = bLimiter[patchI];
 114
 115         if (bLimiter[patchI].coupled())
 116         {
 117             const scalarField& pCDweights = CDweights.boundaryField()[patchI];
 118             const vectorField& pSf = Sf.boundaryField()[patchI];
 119             const scalarField& pmagSf = magSf.boundaryField()[patchI];
 120             const scalarField& pFaceFlux = Uflux.boundaryField()[patchI];
 121             Field<Type> pphiP =
 122                 phi.boundaryField()[patchI].patchInternalField();
 123             Field<Type> pphiN =
 124                 phi.boundaryField()[patchI].patchNeighbourField();
 125
 126             forAll(pLimiter, face)
 127             {
 128                 pLimiter[face] = PhiLimiter::limiter
 129                 (
 130                     pCDweights[face],
 131                     pFaceFlux[face],
 132                     pphiP[face],
 133                     pphiN[face],
 134                     pSf[face],
 135                     pmagSf[face]
 136                 );
 137             }
 138         }
 139         else
 140         {
 141             pLimiter = 1.0;
 142         }
 143     }
 144
 145     return tLimiter;
 146 }
 147
 148
 149 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 150
 151 } // End namespace Foam
 152
 153 // ************************************************************************* //