src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.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 "cellLimitedGrad.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(cellLimitedGrad)
  48
  49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  50
  51 template<>
  52 inline void cellLimitedGrad<scalar>::limitFace
  53 (
  54     scalar& limiter,
  55     const scalar& maxDelta,
  56     const scalar& minDelta,
  57     const scalar& extrapolate
  58 )
  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 template<class Type>
  71 inline void cellLimitedGrad<Type>::limitFace
  72 (
  73     Type& limiter,
  74     const Type& maxDelta,
  75     const Type& minDelta,
  76     const Type& extrapolate
  77 )
  78 {
  79     for(direction cmpt=0; cmpt<Type::nComponents; cmpt++)
  80     {
  81         cellLimitedGrad<scalar>::limitFace
  82         (
  83             limiter.component(cmpt),
  84             maxDelta.component(cmpt),
  85             minDelta.component(cmpt),
  86             extrapolate.component(cmpt)
  87         );
  88     }
  89 }
  90
  91
  92 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  93
  94 template<>
  95 tmp<volVectorField> cellLimitedGrad<scalar>::grad
  96 (
  97     const volScalarField& vsf
  98 ) const
  99 {
 100     const fvMesh& mesh = vsf.mesh();
 101
 102     tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
 103
 104     if (k_ < SMALL)
 105     {
 106         return tGrad;
 107     }
 108
 109     volVectorField& g = tGrad();
 110
 111     const unallocLabelList& owner = mesh.owner();
 112     const unallocLabelList& neighbour = mesh.neighbour();
 113
 114     const volVectorField& C = mesh.C();
 115     const surfaceVectorField& Cf = mesh.Cf();
 116
 117     scalarField maxVsf(vsf.internalField());
 118     scalarField minVsf(vsf.internalField());
 119
 120     forAll(owner, facei)
 121     {
 122         label own = owner[facei];
 123         label nei = neighbour[facei];
 124
 125         scalar vsfOwn = vsf[own];
 126         scalar vsfNei = vsf[nei];
 127
 128         maxVsf[own] = max(maxVsf[own], vsfNei);
 129         minVsf[own] = min(minVsf[own], vsfNei);
 130
 131         maxVsf[nei] = max(maxVsf[nei], vsfOwn);
 132         minVsf[nei] = min(minVsf[nei], vsfOwn);
 133     }
 134
 135
 136     const volScalarField::GeometricBoundaryField& bsf = vsf.boundaryField();
 137
 138     forAll(bsf, patchi)
 139     {
 140         const fvPatchScalarField& psf = bsf[patchi];
 141
 142         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
 143
 144         if (psf.coupled())
 145         {
 146             scalarField psfNei = psf.patchNeighbourField();
 147
 148             forAll(pOwner, pFacei)
 149             {
 150                 label own = pOwner[pFacei];
 151                 scalar vsfNei = psfNei[pFacei];
 152
 153                 maxVsf[own] = max(maxVsf[own], vsfNei);
 154                 minVsf[own] = min(minVsf[own], vsfNei);
 155             }
 156         }
 157         else
 158         {
 159             forAll(pOwner, pFacei)
 160             {
 161                 label own = pOwner[pFacei];
 162                 scalar vsfNei = psf[pFacei];
 163
 164                 maxVsf[own] = max(maxVsf[own], vsfNei);
 165                 minVsf[own] = min(minVsf[own], vsfNei);
 166             }
 167         }
 168     }
 169
 170     maxVsf -= vsf;
 171     minVsf -= vsf;
 172
 173     if (k_ < 1.0)
 174     {
 175         scalarField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf);
 176         maxVsf += maxMinVsf;
 177         minVsf -= maxMinVsf;
 178
 179         //maxVsf *= 1.0/k_;
 180         //minVsf *= 1.0/k_;
 181     }
 182
 183
 184     // create limiter
 185     scalarField limiter(vsf.internalField().size(), 1.0);
 186
 187     forAll(owner, facei)
 188     {
 189         label own = owner[facei];
 190         label nei = neighbour[facei];
 191
 192         // owner side
 193         limitFace
 194         (
 195             limiter[own],
 196             maxVsf[own],
 197             minVsf[own],
 198             (Cf[facei] - C[own]) & g[own]
 199         );
 200
 201         // neighbour side
 202         limitFace
 203         (
 204             limiter[nei],
 205             maxVsf[nei],
 206             minVsf[nei],
 207             (Cf[facei] - C[nei]) & g[nei]
 208         );
 209     }
 210
 211     forAll(bsf, patchi)
 212     {
 213         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
 214         const vectorField& pCf = Cf.boundaryField()[patchi];
 215
 216         forAll(pOwner, pFacei)
 217         {
 218             label own = pOwner[pFacei];
 219
 220             limitFace
 221             (
 222                 limiter[own],
 223                 maxVsf[own],
 224                 minVsf[own],
 225                 (pCf[pFacei] - C[own]) & g[own]
 226             );
 227         }
 228     }
 229
 230     if (fv::debug)
 231     {
 232         Info<< "gradient limiter for: " << vsf.name()
 233             << " max = " << gMax(limiter)
 234             << " min = " << gMin(limiter)
 235             << " average: " << gAverage(limiter) << endl;
 236     }
 237
 238     g.internalField() *= limiter;
 239     g.correctBoundaryConditions();
 240     gaussGrad<scalar>::correctBoundaryConditions(vsf, g);
 241
 242     return tGrad;
 243 }
 244
 245
 246 template<>
 247 tmp<volTensorField> cellLimitedGrad<vector>::grad
 248 (
 249     const volVectorField& vsf
 250 ) const
 251 {
 252     const fvMesh& mesh = vsf.mesh();
 253
 254     tmp<volTensorField> tGrad = basicGradScheme_().grad(vsf);
 255
 256     if (k_ < SMALL)
 257     {
 258         return tGrad;
 259     }
 260
 261     volTensorField& g = tGrad();
 262
 263     const unallocLabelList& owner = mesh.owner();
 264     const unallocLabelList& neighbour = mesh.neighbour();
 265
 266     const volVectorField& C = mesh.C();
 267     const surfaceVectorField& Cf = mesh.Cf();
 268
 269     vectorField maxVsf(vsf.internalField());
 270     vectorField minVsf(vsf.internalField());
 271
 272     forAll(owner, facei)
 273     {
 274         label own = owner[facei];
 275         label nei = neighbour[facei];
 276
 277         const vector& vsfOwn = vsf[own];
 278         const vector& vsfNei = vsf[nei];
 279
 280         maxVsf[own] = max(maxVsf[own], vsfNei);
 281         minVsf[own] = min(minVsf[own], vsfNei);
 282
 283         maxVsf[nei] = max(maxVsf[nei], vsfOwn);
 284         minVsf[nei] = min(minVsf[nei], vsfOwn);
 285     }
 286
 287
 288     const volVectorField::GeometricBoundaryField& bsf = vsf.boundaryField();
 289
 290     forAll(bsf, patchi)
 291     {
 292         const fvPatchVectorField& psf = bsf[patchi];
 293         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
 294
 295         if (psf.coupled())
 296         {
 297             vectorField psfNei = psf.patchNeighbourField();
 298
 299             forAll(pOwner, pFacei)
 300             {
 301                 label own = pOwner[pFacei];
 302                 const vector& vsfNei = psfNei[pFacei];
 303
 304                 maxVsf[own] = max(maxVsf[own], vsfNei);
 305                 minVsf[own] = min(minVsf[own], vsfNei);
 306             }
 307         }
 308         else
 309         {
 310             forAll(pOwner, pFacei)
 311             {
 312                 label own = pOwner[pFacei];
 313                 const vector& vsfNei = psf[pFacei];
 314
 315                 maxVsf[own] = max(maxVsf[own], vsfNei);
 316                 minVsf[own] = min(minVsf[own], vsfNei);
 317             }
 318         }
 319     }
 320
 321     maxVsf -= vsf;
 322     minVsf -= vsf;
 323
 324     if (k_ < 1.0)
 325     {
 326         vectorField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf);
 327         maxVsf += maxMinVsf;
 328         minVsf -= maxMinVsf;
 329
 330         //maxVsf *= 1.0/k_;
 331         //minVsf *= 1.0/k_;
 332     }
 333
 334
 335     // create limiter
 336     vectorField limiter(vsf.internalField().size(), vector::one);
 337
 338     forAll(owner, facei)
 339     {
 340         label own = owner[facei];
 341         label nei = neighbour[facei];
 342
 343         // owner side
 344         limitFace
 345         (
 346             limiter[own],
 347             maxVsf[own],
 348             minVsf[own],
 349             (Cf[facei] - C[own]) & g[own]
 350         );
 351
 352         // neighbour side
 353         limitFace
 354         (
 355             limiter[nei],
 356             maxVsf[nei],
 357             minVsf[nei],
 358             (Cf[facei] - C[nei]) & g[nei]
 359         );
 360     }
 361
 362     forAll(bsf, patchi)
 363     {
 364         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
 365         const vectorField& pCf = Cf.boundaryField()[patchi];
 366
 367         forAll(pOwner, pFacei)
 368         {
 369             label own = pOwner[pFacei];
 370
 371             limitFace
 372             (
 373                 limiter[own],
 374                 maxVsf[own],
 375                 minVsf[own],
 376                 ((pCf[pFacei] - C[own]) & g[own])
 377             );
 378         }
 379     }
 380
 381     if (fv::debug)
 382     {
 383         Info<< "gradient limiter for: " << vsf.name()
 384             << " max = " << gMax(limiter)
 385             << " min = " << gMin(limiter)
 386             << " average: " << gAverage(limiter) << endl;
 387     }
 388
 389     tensorField& gIf = g.internalField();
 390
 391     forAll(gIf, celli)
 392     {
 393         gIf[celli].x() = cmptMultiply(limiter[celli], gIf[celli].x());
 394         gIf[celli].y() = cmptMultiply(limiter[celli], gIf[celli].y());
 395         gIf[celli].z() = cmptMultiply(limiter[celli], gIf[celli].z());
 396     }
 397
 398     g.correctBoundaryConditions();
 399     gaussGrad<vector>::correctBoundaryConditions(vsf, g);
 400
 401     return tGrad;
 402 }
 403
 404
 405 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 406
 407 } // End namespace fv
 408
 409 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 410
 411 } // End namespace Foam
 412
 413 // ************************************************************************* //