src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.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 "gaussLaplacianScheme.H"
  28 #include "surfaceInterpolate.H"
  29 #include "fvcDiv.H"
  30 #include "fvcGrad.H"
  31 #include "fvMatrices.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 namespace fv
  41 {
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 template<class Type, class GType>
  46 tmp<fvMatrix<Type> >
  47 gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
  48 (
  49     const surfaceScalarField& gammaMagSf,
  50     GeometricField<Type, fvPatchField, volMesh>& vf
  51 )
  52 {
  53     tmp<surfaceScalarField> tdeltaCoeffs =
  54         this->tsnGradScheme_().deltaCoeffs(vf);
  55     const surfaceScalarField& deltaCoeffs = tdeltaCoeffs();
  56
  57     tmp<fvMatrix<Type> > tfvm
  58     (
  59         new fvMatrix<Type>
  60         (
  61             vf,
  62             deltaCoeffs.dimensions()*gammaMagSf.dimensions()*vf.dimensions()
  63         )
  64     );
  65     fvMatrix<Type>& fvm = tfvm();
  66
  67     fvm.upper() = deltaCoeffs.internalField()*gammaMagSf.internalField();
  68     fvm.negSumDiag();
  69
  70     forAll(fvm.psi().boundaryField(), patchI)
  71     {
  72         const fvPatchField<Type>& psf = fvm.psi().boundaryField()[patchI];
  73         const fvsPatchScalarField& patchGamma =
  74             gammaMagSf.boundaryField()[patchI];
  75
  76         fvm.internalCoeffs()[patchI] = patchGamma*psf.gradientInternalCoeffs();
  77         fvm.boundaryCoeffs()[patchI] = -patchGamma*psf.gradientBoundaryCoeffs();
  78     }
  79
  80     return tfvm;
  81 }
  82
  83
  84 template<class Type, class GType>
  85 tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
  86 gaussLaplacianScheme<Type, GType>::gammaSnGradCorr
  87 (
  88     const surfaceVectorField& SfGammaCorr,
  89     const GeometricField<Type, fvPatchField, volMesh>& vf
  90 )
  91 {
  92     const fvMesh& mesh = this->mesh();
  93
  94     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tgammaSnGradCorr
  95     (
  96         new GeometricField<Type, fvsPatchField, surfaceMesh>
  97         (
  98             IOobject
  99             (
 100                 "gammaSnGradCorr("+vf.name()+')',
 101                 vf.instance(),
 102                 mesh,
 103                 IOobject::NO_READ,
 104                 IOobject::NO_WRITE
 105             ),
 106             mesh,
 107             SfGammaCorr.dimensions()
 108            *vf.dimensions()*mesh.deltaCoeffs().dimensions()
 109         )
 110     );
 111
 112     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
 113     {
 114         tgammaSnGradCorr().replace
 115         (
 116             cmpt,
 117             SfGammaCorr & fvc::interpolate(fvc::grad(vf.component(cmpt)))
 118         );
 119     }
 120
 121     return tgammaSnGradCorr;
 122 }
 123
 124
 125 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 126
 127 template<class Type, class GType>
 128 tmp<GeometricField<Type, fvPatchField, volMesh> >
 129 gaussLaplacianScheme<Type, GType>::fvcLaplacian
 130 (
 131     const GeometricField<Type, fvPatchField, volMesh>& vf
 132 )
 133 {
 134     const fvMesh& mesh = this->mesh();
 135
 136     tmp<GeometricField<Type, fvPatchField, volMesh> > tLaplacian
 137     (
 138         fvc::div(this->tsnGradScheme_().snGrad(vf)*mesh.magSf())
 139     );
 140
 141     tLaplacian().rename("laplacian(" + vf.name() + ')');
 142
 143     return tLaplacian;
 144 }
 145
 146
 147 template<class Type, class GType>
 148 tmp<fvMatrix<Type> >
 149 gaussLaplacianScheme<Type, GType>::fvmLaplacian
 150 (
 151     const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
 152     GeometricField<Type, fvPatchField, volMesh>& vf
 153 )
 154 {
 155     const fvMesh& mesh = this->mesh();
 156
 157     surfaceVectorField Sn = mesh.Sf()/mesh.magSf();
 158
 159     surfaceVectorField SfGamma = mesh.Sf() & gamma;
 160     GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn = SfGamma & Sn;
 161     surfaceVectorField SfGammaCorr = SfGamma - SfGammaSn*Sn;
 162
 163     tmp<fvMatrix<Type> > tfvm = fvmLaplacianUncorrected(SfGammaSn, vf);
 164     fvMatrix<Type>& fvm = tfvm();
 165
 166     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tfaceFluxCorrection
 167         = gammaSnGradCorr(SfGammaCorr, vf);
 168
 169     if (this->tsnGradScheme_().corrected())
 170     {
 171         tfaceFluxCorrection() +=
 172             SfGammaSn*this->tsnGradScheme_().correction(vf);
 173     }
 174
 175     fvm.source() -= mesh.V()*fvc::div(tfaceFluxCorrection())().internalField();
 176
 177     if (mesh.fluxRequired(vf.name()))
 178     {
 179         fvm.faceFluxCorrectionPtr() = tfaceFluxCorrection.ptr();
 180     }
 181
 182     return tfvm;
 183 }
 184
 185
 186 template<class Type, class GType>
 187 tmp<GeometricField<Type, fvPatchField, volMesh> >
 188 gaussLaplacianScheme<Type, GType>::fvcLaplacian
 189 (
 190     const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
 191     const GeometricField<Type, fvPatchField, volMesh>& vf
 192 )
 193 {
 194     const fvMesh& mesh = this->mesh();
 195
 196     surfaceVectorField Sn = mesh.Sf()/mesh.magSf();
 197
 198     surfaceVectorField SfGamma = mesh.Sf() & gamma;
 199     GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn = SfGamma & Sn;
 200     surfaceVectorField SfGammaCorr = SfGamma - SfGammaSn*Sn;
 201
 202     tmp<GeometricField<Type, fvPatchField, volMesh> > tLaplacian
 203     (
 204         fvc::div
 205         (
 206             SfGammaSn*this->tsnGradScheme_().snGrad(vf)
 207           + gammaSnGradCorr(SfGammaCorr, vf)
 208         )
 209     );
 210
 211     tLaplacian().rename("laplacian(" + gamma.name() + ',' + vf.name() + ')');
 212
 213     return tLaplacian;
 214 }
 215
 216
 217 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 218
 219 } // End namespace fv
 220
 221 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 222
 223 } // End namespace Foam
 224
 225 // ************************************************************************* //