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