FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

cellLimitedGrads.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "cellLimitedGrad.H"
00027 #include <finiteVolume/gaussGrad.H>
00028 #include <finiteVolume/fvMesh.H>
00029 #include <finiteVolume/volMesh.H>
00030 #include <finiteVolume/surfaceMesh.H>
00031 #include <finiteVolume/volFields.H>
00032 #include <finiteVolume/fixedValueFvPatchFields.H>
00033 
00034 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00035 
00036 namespace Foam
00037 {
00038 
00039 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00040 
00041 namespace fv
00042 {
00043 
00044 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00045 
00046 makeFvGradScheme(cellLimitedGrad)
00047 
00048 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00049 
00050 template<>
00051 inline void cellLimitedGrad<scalar>::limitFace
00052 (
00053     scalar& limiter,
00054     const scalar& maxDelta,
00055     const scalar& minDelta,
00056     const scalar& extrapolate
00057 )
00058 {
00059     if (extrapolate > maxDelta + VSMALL)
00060     {
00061         limiter = min(limiter, maxDelta/extrapolate);
00062     }
00063     else if (extrapolate < minDelta - VSMALL)
00064     {
00065         limiter = min(limiter, minDelta/extrapolate);
00066     }
00067 }
00068 
00069 template<class Type>
00070 inline void cellLimitedGrad<Type>::limitFace
00071 (
00072     Type& limiter,
00073     const Type& maxDelta,
00074     const Type& minDelta,
00075     const Type& extrapolate
00076 )
00077 {
00078     for(direction cmpt=0; cmpt<Type::nComponents; cmpt++)
00079     {
00080         cellLimitedGrad<scalar>::limitFace
00081         (
00082             limiter.component(cmpt),
00083             maxDelta.component(cmpt),
00084             minDelta.component(cmpt),
00085             extrapolate.component(cmpt)
00086         );
00087     }
00088 }
00089 
00090 
00091 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00092 
00093 template<>
00094 tmp<volVectorField> cellLimitedGrad<scalar>::grad
00095 (
00096     const volScalarField& vsf
00097 ) const
00098 {
00099     const fvMesh& mesh = vsf.mesh();
00100 
00101     tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
00102 
00103     if (k_ < SMALL)
00104     {
00105         return tGrad;
00106     }
00107 
00108     volVectorField& g = tGrad();
00109 
00110     const unallocLabelList& owner = mesh.owner();
00111     const unallocLabelList& neighbour = mesh.neighbour();
00112 
00113     const volVectorField& C = mesh.C();
00114     const surfaceVectorField& Cf = mesh.Cf();
00115 
00116     scalarField maxVsf(vsf.internalField());
00117     scalarField minVsf(vsf.internalField());
00118 
00119     forAll(owner, facei)
00120     {
00121         label own = owner[facei];
00122         label nei = neighbour[facei];
00123 
00124         scalar vsfOwn = vsf[own];
00125         scalar vsfNei = vsf[nei];
00126 
00127         maxVsf[own] = max(maxVsf[own], vsfNei);
00128         minVsf[own] = min(minVsf[own], vsfNei);
00129 
00130         maxVsf[nei] = max(maxVsf[nei], vsfOwn);
00131         minVsf[nei] = min(minVsf[nei], vsfOwn);
00132     }
00133 
00134 
00135     const volScalarField::GeometricBoundaryField& bsf = vsf.boundaryField();
00136 
00137     forAll(bsf, patchi)
00138     {
00139         const fvPatchScalarField& psf = bsf[patchi];
00140 
00141         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
00142 
00143         if (psf.coupled())
00144         {
00145             scalarField psfNei = psf.patchNeighbourField();
00146 
00147             forAll(pOwner, pFacei)
00148             {
00149                 label own = pOwner[pFacei];
00150                 scalar vsfNei = psfNei[pFacei];
00151 
00152                 maxVsf[own] = max(maxVsf[own], vsfNei);
00153                 minVsf[own] = min(minVsf[own], vsfNei);
00154             }
00155         }
00156         else
00157         {
00158             forAll(pOwner, pFacei)
00159             {
00160                 label own = pOwner[pFacei];
00161                 scalar vsfNei = psf[pFacei];
00162 
00163                 maxVsf[own] = max(maxVsf[own], vsfNei);
00164                 minVsf[own] = min(minVsf[own], vsfNei);
00165             }
00166         }
00167     }
00168 
00169     maxVsf -= vsf;
00170     minVsf -= vsf;
00171 
00172     if (k_ < 1.0)
00173     {
00174         scalarField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf);
00175         maxVsf += maxMinVsf;
00176         minVsf -= maxMinVsf;
00177 
00178         //maxVsf *= 1.0/k_;
00179         //minVsf *= 1.0/k_;
00180     }
00181 
00182 
00183     // create limiter
00184     scalarField limiter(vsf.internalField().size(), 1.0);
00185 
00186     forAll(owner, facei)
00187     {
00188         label own = owner[facei];
00189         label nei = neighbour[facei];
00190 
00191         // owner side
00192         limitFace
00193         (
00194             limiter[own],
00195             maxVsf[own],
00196             minVsf[own],
00197             (Cf[facei] - C[own]) & g[own]
00198         );
00199 
00200         // neighbour side
00201         limitFace
00202         (
00203             limiter[nei],
00204             maxVsf[nei],
00205             minVsf[nei],
00206             (Cf[facei] - C[nei]) & g[nei]
00207         );
00208     }
00209 
00210     forAll(bsf, patchi)
00211     {
00212         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
00213         const vectorField& pCf = Cf.boundaryField()[patchi];
00214 
00215         forAll(pOwner, pFacei)
00216         {
00217             label own = pOwner[pFacei];
00218 
00219             limitFace
00220             (
00221                 limiter[own],
00222                 maxVsf[own],
00223                 minVsf[own],
00224                 (pCf[pFacei] - C[own]) & g[own]
00225             );
00226         }
00227     }
00228 
00229     if (fv::debug)
00230     {
00231         Info<< "gradient limiter for: " << vsf.name()
00232             << " max = " << gMax(limiter)
00233             << " min = " << gMin(limiter)
00234             << " average: " << gAverage(limiter) << endl;
00235     }
00236 
00237     g.internalField() *= limiter;
00238     g.correctBoundaryConditions();
00239     gaussGrad<scalar>::correctBoundaryConditions(vsf, g);
00240 
00241     return tGrad;
00242 }
00243 
00244 
00245 template<>
00246 tmp<volTensorField> cellLimitedGrad<vector>::grad
00247 (
00248     const volVectorField& vsf
00249 ) const
00250 {
00251     const fvMesh& mesh = vsf.mesh();
00252 
00253     tmp<volTensorField> tGrad = basicGradScheme_().grad(vsf);
00254 
00255     if (k_ < SMALL)
00256     {
00257         return tGrad;
00258     }
00259 
00260     volTensorField& g = tGrad();
00261 
00262     const unallocLabelList& owner = mesh.owner();
00263     const unallocLabelList& neighbour = mesh.neighbour();
00264 
00265     const volVectorField& C = mesh.C();
00266     const surfaceVectorField& Cf = mesh.Cf();
00267 
00268     vectorField maxVsf(vsf.internalField());
00269     vectorField minVsf(vsf.internalField());
00270 
00271     forAll(owner, facei)
00272     {
00273         label own = owner[facei];
00274         label nei = neighbour[facei];
00275 
00276         const vector& vsfOwn = vsf[own];
00277         const vector& vsfNei = vsf[nei];
00278 
00279         maxVsf[own] = max(maxVsf[own], vsfNei);
00280         minVsf[own] = min(minVsf[own], vsfNei);
00281 
00282         maxVsf[nei] = max(maxVsf[nei], vsfOwn);
00283         minVsf[nei] = min(minVsf[nei], vsfOwn);
00284     }
00285 
00286 
00287     const volVectorField::GeometricBoundaryField& bsf = vsf.boundaryField();
00288 
00289     forAll(bsf, patchi)
00290     {
00291         const fvPatchVectorField& psf = bsf[patchi];
00292         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
00293 
00294         if (psf.coupled())
00295         {
00296             vectorField psfNei = psf.patchNeighbourField();
00297 
00298             forAll(pOwner, pFacei)
00299             {
00300                 label own = pOwner[pFacei];
00301                 const vector& vsfNei = psfNei[pFacei];
00302 
00303                 maxVsf[own] = max(maxVsf[own], vsfNei);
00304                 minVsf[own] = min(minVsf[own], vsfNei);
00305             }
00306         }
00307         else
00308         {
00309             forAll(pOwner, pFacei)
00310             {
00311                 label own = pOwner[pFacei];
00312                 const vector& vsfNei = psf[pFacei];
00313 
00314                 maxVsf[own] = max(maxVsf[own], vsfNei);
00315                 minVsf[own] = min(minVsf[own], vsfNei);
00316             }
00317         }
00318     }
00319 
00320     maxVsf -= vsf;
00321     minVsf -= vsf;
00322 
00323     if (k_ < 1.0)
00324     {
00325         vectorField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf);
00326         maxVsf += maxMinVsf;
00327         minVsf -= maxMinVsf;
00328 
00329         //maxVsf *= 1.0/k_;
00330         //minVsf *= 1.0/k_;
00331     }
00332 
00333 
00334     // create limiter
00335     vectorField limiter(vsf.internalField().size(), vector::one);
00336 
00337     forAll(owner, facei)
00338     {
00339         label own = owner[facei];
00340         label nei = neighbour[facei];
00341 
00342         // owner side
00343         limitFace
00344         (
00345             limiter[own],
00346             maxVsf[own],
00347             minVsf[own],
00348             (Cf[facei] - C[own]) & g[own]
00349         );
00350 
00351         // neighbour side
00352         limitFace
00353         (
00354             limiter[nei],
00355             maxVsf[nei],
00356             minVsf[nei],
00357             (Cf[facei] - C[nei]) & g[nei]
00358         );
00359     }
00360 
00361     forAll(bsf, patchi)
00362     {
00363         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
00364         const vectorField& pCf = Cf.boundaryField()[patchi];
00365 
00366         forAll(pOwner, pFacei)
00367         {
00368             label own = pOwner[pFacei];
00369 
00370             limitFace
00371             (
00372                 limiter[own],
00373                 maxVsf[own],
00374                 minVsf[own],
00375                 ((pCf[pFacei] - C[own]) & g[own])
00376             );
00377         }
00378     }
00379 
00380     if (fv::debug)
00381     {
00382         Info<< "gradient limiter for: " << vsf.name()
00383             << " max = " << gMax(limiter)
00384             << " min = " << gMin(limiter)
00385             << " average: " << gAverage(limiter) << endl;
00386     }
00387 
00388     tensorField& gIf = g.internalField();
00389 
00390     forAll(gIf, celli)
00391     {
00392         gIf[celli] = tensor
00393         (
00394             cmptMultiply(limiter[celli], gIf[celli].x()),
00395             cmptMultiply(limiter[celli], gIf[celli].y()),
00396             cmptMultiply(limiter[celli], gIf[celli].z())
00397         );
00398     }
00399 
00400     g.correctBoundaryConditions();
00401     gaussGrad<vector>::correctBoundaryConditions(vsf, g);
00402 
00403     return tGrad;
00404 }
00405 
00406 
00407 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00408 
00409 } // End namespace fv
00410 
00411 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00412 
00413 } // End namespace Foam
00414 
00415 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines