src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.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 "faceLimitedGrad.H"
  28 #include "gaussGrad.H"
  29 #include "fvMesh.H"
  30 #include "volMesh.H"
  31 #include "surfaceMesh.H"
  32 #include "volFields.H"
  33 #include "fixedValueFvPatchFields.H"
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 namespace fv
  43 {
  44
  45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  46
  47 makeFvGradScheme(faceLimitedGrad)
  48
  49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  50
  51 template<class Type>
  52 inline void faceLimitedGrad<Type>::limitFace
  53 (
  54     scalar& limiter,
  55     const scalar maxDelta,
  56     const scalar minDelta,
  57     const scalar extrapolate
  58 ) const
  59 {
  60     if (extrapolate > maxDelta + VSMALL)
  61     {
  62         limiter = min(limiter, maxDelta/extrapolate);
  63     }
  64     else if (extrapolate < minDelta - VSMALL)
  65     {
  66         limiter = min(limiter, minDelta/extrapolate);
  67     }
  68 }
  69
  70
  71 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  72
  73 template<>
  74 tmp<volVectorField> faceLimitedGrad<scalar>::grad
  75 (
  76     const volScalarField& vsf
  77 ) const
  78 {
  79     const fvMesh& mesh = vsf.mesh();
  80
  81     tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
  82
  83     if (k_ < SMALL)
  84     {
  85         return tGrad;
  86     }
  87
  88     volVectorField& g = tGrad();
  89
  90     const unallocLabelList& owner = mesh.owner();
  91     const unallocLabelList& neighbour = mesh.neighbour();
  92
  93     const volVectorField& C = mesh.C();
  94     const surfaceVectorField& Cf = mesh.Cf();
  95
  96     // create limiter
  97     scalarField limiter(vsf.internalField().size(), 1.0);
  98
  99     scalar rk = (1.0/k_ - 1.0);
 100
 101     forAll(owner, facei)
 102     {
 103         label own = owner[facei];
 104         label nei = neighbour[facei];
 105
 106         scalar vsfOwn = vsf[own];
 107         scalar vsfNei = vsf[nei];
 108
 109         scalar maxFace = max(vsfOwn, vsfNei);
 110         scalar minFace = min(vsfOwn, vsfNei);
 111         scalar maxMinFace = rk*(maxFace - minFace);
 112         maxFace += maxMinFace;
 113         minFace -= maxMinFace;
 114
 115         // owner side
 116         limitFace
 117         (
 118             limiter[own],
 119             maxFace - vsfOwn, minFace - vsfOwn,
 120             (Cf[facei] - C[own]) & g[own]
 121         );
 122
 123         // neighbour side
 124         limitFace
 125         (
 126             limiter[nei],
 127             maxFace - vsfNei, minFace - vsfNei,
 128             (Cf[facei] - C[nei]) & g[nei]
 129         );
 130     }
 131
 132     const volScalarField::GeometricBoundaryField& bsf = vsf.boundaryField();
 133
 134     forAll(bsf, patchi)
 135     {
 136         const fvPatchScalarField& psf = bsf[patchi];
 137
 138         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
 139         const vectorField& pCf = Cf.boundaryField()[patchi];
 140
 141         if (psf.coupled())
 142         {
 143             scalarField psfNei = psf.patchNeighbourField();
 144
 145             forAll(pOwner, pFacei)
 146             {
 147                 label own = pOwner[pFacei];
 148
 149                 scalar vsfOwn = vsf[own];
 150                 scalar vsfNei = psfNei[pFacei];
 151
 152                 scalar maxFace = max(vsfOwn, vsfNei);
 153                 scalar minFace = min(vsfOwn, vsfNei);
 154                 scalar maxMinFace = rk*(maxFace - minFace);
 155                 maxFace += maxMinFace;
 156                 minFace -= maxMinFace;
 157
 158                 limitFace
 159                 (
 160                     limiter[own],
 161                     maxFace - vsfOwn, minFace - vsfOwn,
 162                     (pCf[pFacei] - C[own]) & g[own]
 163                 );
 164             }
 165         }
 166         else if (psf.fixesValue())
 167         {
 168             forAll(pOwner, pFacei)
 169             {
 170                 label own = pOwner[pFacei];
 171
 172                 scalar vsfOwn = vsf[own];
 173                 scalar vsfNei = psf[pFacei];
 174
 175                 scalar maxFace = max(vsfOwn, vsfNei);
 176                 scalar minFace = min(vsfOwn, vsfNei);
 177                 scalar maxMinFace = rk*(maxFace - minFace);
 178                 maxFace += maxMinFace;
 179                 minFace -= maxMinFace;
 180
 181                 limitFace
 182                 (
 183                     limiter[own],
 184                     maxFace - vsfOwn, minFace - vsfOwn,
 185                     (pCf[pFacei] - C[own]) & g[own]
 186                 );
 187             }
 188         }
 189     }
 190
 191     if (fv::debug)
 192     {
 193         Info<< "gradient limiter for: " << vsf.name()
 194             << " max = " << gMax(limiter)
 195             << " min = " << gMin(limiter)
 196             << " average: " << gAverage(limiter) << endl;
 197     }
 198
 199     g.internalField() *= limiter;
 200     g.correctBoundaryConditions();
 201     gaussGrad<scalar>::correctBoundaryConditions(vsf, g);
 202
 203     return tGrad;
 204 }
 205
 206
 207 template<>
 208 tmp<volTensorField> faceLimitedGrad<vector>::grad
 209 (
 210     const volVectorField& vvf
 211 ) const
 212 {
 213     const fvMesh& mesh = vvf.mesh();
 214
 215     tmp<volTensorField> tGrad = basicGradScheme_().grad(vvf);
 216
 217     if (k_ < SMALL)
 218     {
 219         return tGrad;
 220     }
 221
 222     volTensorField& g = tGrad();
 223
 224     const unallocLabelList& owner = mesh.owner();
 225     const unallocLabelList& neighbour = mesh.neighbour();
 226
 227     const volVectorField& C = mesh.C();
 228     const surfaceVectorField& Cf = mesh.Cf();
 229
 230     // create limiter
 231     scalarField limiter(vvf.internalField().size(), 1.0);
 232
 233     scalar rk = (1.0/k_ - 1.0);
 234
 235     forAll(owner, facei)
 236     {
 237         label own = owner[facei];
 238         label nei = neighbour[facei];
 239
 240         vector vvfOwn = vvf[own];
 241         vector vvfNei = vvf[nei];
 242
 243         // owner side
 244         vector gradf = (Cf[facei] - C[own]) & g[own];
 245
 246         scalar vsfOwn = gradf & vvfOwn;
 247         scalar vsfNei = gradf & vvfNei;
 248
 249         scalar maxFace = max(vsfOwn, vsfNei);
 250         scalar minFace = min(vsfOwn, vsfNei);
 251         scalar maxMinFace = rk*(maxFace - minFace);
 252         maxFace += maxMinFace;
 253         minFace -= maxMinFace;
 254
 255         limitFace
 256         (
 257             limiter[own],
 258             maxFace - vsfOwn, minFace - vsfOwn,
 259             magSqr(gradf)
 260         );
 261
 262
 263         // neighbour side
 264         gradf = (Cf[facei] - C[nei]) & g[nei];
 265
 266         vsfOwn = gradf & vvfOwn;
 267         vsfNei = gradf & vvfNei;
 268
 269         maxFace = max(vsfOwn, vsfNei);
 270         minFace = min(vsfOwn, vsfNei);
 271
 272         limitFace
 273         (
 274             limiter[nei],
 275             maxFace - vsfNei, minFace - vsfNei,
 276             magSqr(gradf)
 277         );
 278     }
 279
 280
 281     const volVectorField::GeometricBoundaryField& bvf = vvf.boundaryField();
 282
 283     forAll(bvf, patchi)
 284     {
 285         const fvPatchVectorField& psf = bvf[patchi];
 286
 287         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
 288         const vectorField& pCf = Cf.boundaryField()[patchi];
 289
 290         if (psf.coupled())
 291         {
 292             vectorField psfNei = psf.patchNeighbourField();
 293
 294             forAll(pOwner, pFacei)
 295             {
 296                 label own = pOwner[pFacei];
 297
 298                 vector vvfOwn = vvf[own];
 299                 vector vvfNei = psfNei[pFacei];
 300
 301                 vector gradf = (pCf[pFacei] - C[own]) & g[own];
 302
 303                 scalar vsfOwn = gradf & vvfOwn;
 304                 scalar vsfNei = gradf & vvfNei;
 305
 306                 scalar maxFace = max(vsfOwn, vsfNei);
 307                 scalar minFace = min(vsfOwn, vsfNei);
 308                 scalar maxMinFace = rk*(maxFace - minFace);
 309                 maxFace += maxMinFace;
 310                 minFace -= maxMinFace;
 311
 312                 limitFace
 313                 (
 314                     limiter[own],
 315                     maxFace - vsfOwn, minFace - vsfOwn,
 316                     magSqr(gradf)
 317                 );
 318             }
 319         }
 320         else if (psf.fixesValue())
 321         {
 322             forAll(pOwner, pFacei)
 323             {
 324                 label own = pOwner[pFacei];
 325
 326                 vector vvfOwn = vvf[own];
 327                 vector vvfNei = psf[pFacei];
 328
 329                 vector gradf = (pCf[pFacei] - C[own]) & g[own];
 330
 331                 scalar vsfOwn = gradf & vvfOwn;
 332                 scalar vsfNei = gradf & vvfNei;
 333
 334                 scalar maxFace = max(vsfOwn, vsfNei);
 335                 scalar minFace = min(vsfOwn, vsfNei);
 336                 scalar maxMinFace = rk*(maxFace - minFace);
 337                 maxFace += maxMinFace;
 338                 minFace -= maxMinFace;
 339
 340                 limitFace
 341                 (
 342                     limiter[own],
 343                     maxFace - vsfOwn, minFace - vsfOwn,
 344                     magSqr(gradf)
 345                 );
 346             }
 347         }
 348     }
 349
 350     if (fv::debug)
 351     {
 352         Info<< "gradient limiter for: " << vvf.name()
 353             << " max = " << gMax(limiter)
 354             << " min = " << gMin(limiter)
 355             << " average: " << gAverage(limiter) << endl;
 356     }
 357
 358     g.internalField() *= limiter;
 359     g.correctBoundaryConditions();
 360     gaussGrad<vector>::correctBoundaryConditions(vvf, g);
 361
 362     return tGrad;
 363 }
 364
 365
 366 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 367
 368 } // End namespace fv
 369
 370 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 371
 372 } // End namespace Foam
 373
 374 // ************************************************************************* //