src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellMDLimitedGrad/cellMDLimitedGrads.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 "cellMDLimitedGrad.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(cellMDLimitedGrad)
  48
  49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  50
  51 template<>
  52 tmp<volVectorField> cellMDLimitedGrad<scalar>::grad
  53 (
  54     const volScalarField& vsf
  55 ) const
  56 {
  57     const fvMesh& mesh = vsf.mesh();
  58
  59     tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
  60
  61     if (k_ < SMALL)
  62     {
  63         return tGrad;
  64     }
  65
  66     volVectorField& g = tGrad();
  67
  68     const unallocLabelList& owner = mesh.owner();
  69     const unallocLabelList& neighbour = mesh.neighbour();
  70
  71     const volVectorField& C = mesh.C();
  72     const surfaceVectorField& Cf = mesh.Cf();
  73
  74     scalarField maxVsf(vsf.internalField());
  75     scalarField minVsf(vsf.internalField());
  76
  77     forAll(owner, facei)
  78     {
  79         label own = owner[facei];
  80         label nei = neighbour[facei];
  81
  82         scalar vsfOwn = vsf[own];
  83         scalar vsfNei = vsf[nei];
  84
  85         maxVsf[own] = max(maxVsf[own], vsfNei);
  86         minVsf[own] = min(minVsf[own], vsfNei);
  87
  88         maxVsf[nei] = max(maxVsf[nei], vsfOwn);
  89         minVsf[nei] = min(minVsf[nei], vsfOwn);
  90     }
  91
  92
  93     const volScalarField::GeometricBoundaryField& bsf = vsf.boundaryField();
  94
  95     forAll(bsf, patchi)
  96     {
  97         const fvPatchScalarField& psf = bsf[patchi];
  98
  99         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
 100
 101         if (psf.coupled())
 102         {
 103             scalarField psfNei = psf.patchNeighbourField();
 104
 105             forAll(pOwner, pFacei)
 106             {
 107                 label own = pOwner[pFacei];
 108                 scalar vsfNei = psfNei[pFacei];
 109
 110                 maxVsf[own] = max(maxVsf[own], vsfNei);
 111                 minVsf[own] = min(minVsf[own], vsfNei);
 112             }
 113         }
 114         else
 115         {
 116             forAll(pOwner, pFacei)
 117             {
 118                 label own = pOwner[pFacei];
 119                 scalar vsfNei = psf[pFacei];
 120
 121                 maxVsf[own] = max(maxVsf[own], vsfNei);
 122                 minVsf[own] = min(minVsf[own], vsfNei);
 123             }
 124         }
 125     }
 126
 127     maxVsf -= vsf;
 128     minVsf -= vsf;
 129
 130     if (k_ < 1.0)
 131     {
 132         scalarField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf);
 133         maxVsf += maxMinVsf;
 134         minVsf -= maxMinVsf;
 135
 136         //maxVsf *= 1.0/k_;
 137         //minVsf *= 1.0/k_;
 138     }
 139
 140
 141     forAll(owner, facei)
 142     {
 143         label own = owner[facei];
 144         label nei = neighbour[facei];
 145
 146         // owner side
 147         limitFace
 148         (
 149             g[own],
 150             maxVsf[own],
 151             minVsf[own],
 152             Cf[facei] - C[own]
 153         );
 154
 155         // neighbour side
 156         limitFace
 157         (
 158             g[nei],
 159             maxVsf[nei],
 160             minVsf[nei],
 161             Cf[facei] - C[nei]
 162         );
 163     }
 164
 165
 166     forAll(bsf, patchi)
 167     {
 168         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
 169         const vectorField& pCf = Cf.boundaryField()[patchi];
 170
 171         forAll(pOwner, pFacei)
 172         {
 173             label own = pOwner[pFacei];
 174
 175             limitFace
 176             (
 177                 g[own],
 178                 maxVsf[own],
 179                 minVsf[own],
 180                 pCf[pFacei] - C[own]
 181             );
 182         }
 183     }
 184
 185     g.correctBoundaryConditions();
 186     gaussGrad<scalar>::correctBoundaryConditions(vsf, g);
 187
 188     return tGrad;
 189 }
 190
 191
 192 template<>
 193 tmp<volTensorField> cellMDLimitedGrad<vector>::grad
 194 (
 195     const volVectorField& vsf
 196 ) const
 197 {
 198     const fvMesh& mesh = vsf.mesh();
 199
 200     tmp<volTensorField> tGrad = basicGradScheme_().grad(vsf);
 201
 202     if (k_ < SMALL)
 203     {
 204         return tGrad;
 205     }
 206
 207     volTensorField& g = tGrad();
 208
 209     const unallocLabelList& owner = mesh.owner();
 210     const unallocLabelList& neighbour = mesh.neighbour();
 211
 212     const volVectorField& C = mesh.C();
 213     const surfaceVectorField& Cf = mesh.Cf();
 214
 215     vectorField maxVsf(vsf.internalField());
 216     vectorField minVsf(vsf.internalField());
 217
 218     forAll(owner, facei)
 219     {
 220         label own = owner[facei];
 221         label nei = neighbour[facei];
 222
 223         const vector& vsfOwn = vsf[own];
 224         const vector& vsfNei = vsf[nei];
 225
 226         maxVsf[own] = max(maxVsf[own], vsfNei);
 227         minVsf[own] = min(minVsf[own], vsfNei);
 228
 229         maxVsf[nei] = max(maxVsf[nei], vsfOwn);
 230         minVsf[nei] = min(minVsf[nei], vsfOwn);
 231     }
 232
 233
 234     const volVectorField::GeometricBoundaryField& bsf = vsf.boundaryField();
 235
 236     forAll(bsf, patchi)
 237     {
 238         const fvPatchVectorField& psf = bsf[patchi];
 239         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
 240
 241         if (psf.coupled())
 242         {
 243             vectorField psfNei = psf.patchNeighbourField();
 244
 245             forAll(pOwner, pFacei)
 246             {
 247                 label own = pOwner[pFacei];
 248                 const vector& vsfNei = psfNei[pFacei];
 249
 250                 maxVsf[own] = max(maxVsf[own], vsfNei);
 251                 minVsf[own] = min(minVsf[own], vsfNei);
 252             }
 253         }
 254         else
 255         {
 256             forAll(pOwner, pFacei)
 257             {
 258                 label own = pOwner[pFacei];
 259                 const vector& vsfNei = psf[pFacei];
 260
 261                 maxVsf[own] = max(maxVsf[own], vsfNei);
 262                 minVsf[own] = min(minVsf[own], vsfNei);
 263             }
 264         }
 265     }
 266
 267     maxVsf -= vsf;
 268     minVsf -= vsf;
 269
 270     if (k_ < 1.0)
 271     {
 272         vectorField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf);
 273         maxVsf += maxMinVsf;
 274         minVsf -= maxMinVsf;
 275
 276         //maxVsf *= 1.0/k_;
 277         //minVsf *= 1.0/k_;
 278     }
 279
 280
 281     forAll(owner, facei)
 282     {
 283         label own = owner[facei];
 284         label nei = neighbour[facei];
 285
 286         // owner side
 287         limitFace
 288         (
 289             g[own],
 290             maxVsf[own],
 291             minVsf[own],
 292             Cf[facei] - C[own]
 293         );
 294
 295         // neighbour side
 296         limitFace
 297         (
 298             g[nei],
 299             maxVsf[nei],
 300             minVsf[nei],
 301             Cf[facei] - C[nei]
 302         );
 303     }
 304
 305
 306     forAll(bsf, patchi)
 307     {
 308         const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
 309         const vectorField& pCf = Cf.boundaryField()[patchi];
 310
 311         forAll(pOwner, pFacei)
 312         {
 313             label own = pOwner[pFacei];
 314
 315             limitFace
 316             (
 317                 g[own],
 318                 maxVsf[own],
 319                 minVsf[own],
 320                 pCf[pFacei] - C[own]
 321             );
 322         }
 323     }
 324
 325     g.correctBoundaryConditions();
 326     gaussGrad<vector>::correctBoundaryConditions(vsf, g);
 327
 328     return tGrad;
 329 }
 330
 331
 332 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 333
 334 } // End namespace fv
 335
 336 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 337
 338 } // End namespace Foam
 339
 340 // ************************************************************************* //