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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
  41 Foam::DeferredCorrectionLimitedScheme<Type, Limiter, LimitFunc>::correction
  42 (
  43     const GeometricField<Type, fvPatchField, volMesh>& vf
  44 ) const
  45 {
  46     // Correction = full interpolation minus upwinded part
  47     return
  48         surfaceInterpolationScheme<Type>::interpolate
  49         (
  50             vf,
  51             limitedSurfaceInterpolationScheme<Type>::weights(vf)
  52         )
  53       - upwindScheme_.interpolate(vf);
  54 }
  55
  56
  57 template<class Type, class Limiter, template<class> class LimitFunc>
  58 tmp<surfaceScalarField>
  59 DeferredCorrectionLimitedScheme<Type, Limiter, LimitFunc>::limiter
  60 (
  61     const GeometricField<Type, fvPatchField, volMesh>& phi
  62 ) const
  63 {
  64     const fvMesh& mesh = this->mesh();
  65
  66     tmp<surfaceScalarField> tLimiter
  67     (
  68         new surfaceScalarField
  69         (
  70             IOobject
  71             (
  72                 type() + "Limiter(" + phi.name() + ')',
  73                 mesh.time().timeName(),
  74                 mesh
  75             ),
  76             mesh,
  77             dimless
  78         )
  79     );
  80     surfaceScalarField& lim = tLimiter();
  81
  82     tmp<GeometricField<typename Limiter::phiType, fvPatchField, volMesh> >
  83         tlPhi = LimitFunc<Type>()(phi);
  84
  85     const GeometricField<typename Limiter::phiType, fvPatchField, volMesh>&
  86         lPhi = tlPhi();
  87
  88     GeometricField<typename Limiter::gradPhiType, fvPatchField, volMesh>
  89         gradc(fvc::grad(lPhi));
  90
  91     // Note: in order for the patchNeighbourField to be correct on coupled
  92     // boundaries, correctBoundaryConditions needs to be called.
  93     // The call shall be moved into the library fvc operators
  94     gradc.correctBoundaryConditions();
  95
  96     const surfaceScalarField& CDweights = mesh.surfaceInterpolation::weights();
  97
  98     const unallocLabelList& owner = mesh.owner();
  99     const unallocLabelList& neighbour = mesh.neighbour();
 100
 101     const vectorField& C = mesh.C();
 102
 103     scalarField& pLim = lim.internalField();
 104
 105     forAll(pLim, face)
 106     {
 107         label own = owner[face];
 108         label nei = neighbour[face];
 109
 110         pLim[face] = Limiter::limiter
 111         (
 112             CDweights[face],
 113             this->faceFlux_[face],
 114             lPhi[own],
 115             lPhi[nei],
 116             gradc[own],
 117             gradc[nei],
 118             C[nei] - C[own]
 119         );
 120     }
 121
 122     surfaceScalarField::GeometricBoundaryField& bLim = lim.boundaryField();
 123
 124     forAll (bLim, patchi)
 125     {
 126         scalarField& pLim = bLim[patchi];
 127
 128         if (bLim[patchi].coupled())
 129         {
 130             const scalarField& pCDweights = CDweights.boundaryField()[patchi];
 131             const scalarField& pFaceFlux =
 132                 this->faceFlux_.boundaryField()[patchi];
 133             Field<typename Limiter::phiType> plPhiP =
 134                 lPhi.boundaryField()[patchi].patchInternalField();
 135             Field<typename Limiter::phiType> plPhiN =
 136                 lPhi.boundaryField()[patchi].patchNeighbourField();
 137             Field<typename Limiter::gradPhiType> pGradcP =
 138                 gradc.boundaryField()[patchi].patchInternalField();
 139             Field<typename Limiter::gradPhiType> pGradcN =
 140                 gradc.boundaryField()[patchi].patchNeighbourField();
 141
 142             // Build the d-vectors
 143             // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010
 144             vectorField pd = bLim[patchi].patch().delta();
 145
 146             forAll(pLim, face)
 147             {
 148                 pLim[face] = Limiter::limiter
 149                 (
 150                     pCDweights[face],
 151                     pFaceFlux[face],
 152                     plPhiP[face],
 153                     plPhiN[face],
 154                     pGradcP[face],
 155                     pGradcN[face],
 156                     pd[face]
 157                 );
 158             }
 159         }
 160         else
 161         {
 162             pLim = 1.0;
 163         }
 164     }
 165
 166     return tLimiter;
 167 }
 168
 169
 170 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 171
 172 } // End namespace Foam
 173
 174 // ************************************************************************* //