src/dbns/limiter/MDLimiter.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Class
  25     MDLimiter
  26
  27 Description
  28     Multidimensional limiter for density based solver
  29
  30 Author
  31     Aleksandar Jemcov
  32     Updated by Hrvoje Jasak
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #ifndef MDLimiter_H
  37 #define MDLimiter_H
  38
  39 #include "volFieldsFwd.H"
  40 #include "surfaceFieldsFwd.H"
  41 #include "BarthJespersenLimiter.H"
  42 #include "VenkatakrishnanLimiter.H"
  43
  44 namespace Foam
  45 {
  46
  47 /*---------------------------------------------------------------------------*\
  48                           Class MDLimiter Declaration
  49 \*---------------------------------------------------------------------------*/
  50
  51 template
  52 <
  53     class Type,
  54     typename LimiterFunctionType
  55 >
  56 class MDLimiter
  57 {
  58 public:
  59
  60     // Public typedefs
  61
  62     typedef Field<Type> FieldType;
  63     typedef GeometricField<Type, fvPatchField, volMesh> GeoFieldType;
  64
  65     typedef Field<typename outerProduct<vector, Type>::type> GradFieldType;
  66     typedef GeometricField
  67     <
  68         typename outerProduct<vector, Type>::type,
  69         fvPatchField,
  70         volMesh
  71     > GeoGradFieldType;
  72
  73
  74 private:
  75
  76     // Private data
  77
  78         //- Reference to the mesh
  79         const fvMesh& mesh_;
  80
  81         //- Field
  82         const GeoFieldType& phi_;
  83
  84         //- Grad field
  85         const GeoGradFieldType& gradPhi_;
  86
  87         //- Limiter
  88         GeoFieldType phiLimiter_;
  89
  90
  91 public:
  92
  93     // Constructors
  94
  95         //- Construct from field and gradient field
  96         MDLimiter
  97         (
  98             const GeoFieldType& phi,
  99             const GeoGradFieldType& gradPhi
 100         )
 101         :
 102             mesh_(phi.mesh()),
 103             phi_(phi),
 104             gradPhi_(gradPhi),
 105             phiLimiter_
 106             (
 107                 IOobject
 108                 (
 109                     "phiLimiter",
 110                     mesh_.time().timeName(),
 111                     mesh_,
 112                     IOobject::NO_READ,
 113                     IOobject::NO_WRITE
 114                 ),
 115                 mesh_,
 116                 dimensioned<Type>("pLimiter", dimless, pTraits<Type>::one)
 117             )
 118         {
 119             // Collect min and max values from neighbourhood
 120             GeoFieldType phiMinValue
 121             (
 122                 "phiMinValue",
 123                 phi
 124             );
 125
 126             GeoFieldType phiMaxValue
 127             (
 128                 "phiMaxValue",
 129                 phi
 130             );
 131
 132             const unallocLabelList& owner = mesh_.owner();
 133             const unallocLabelList& neighbour = mesh_.neighbour();
 134
 135             Field<Type>& phiMinIn = phiMinValue.internalField();
 136             Field<Type>& phiMaxIn = phiMaxValue.internalField();
 137
 138             forAll (owner, faceI)
 139             {
 140                 const label& own = owner[faceI];
 141                 const label& nei = neighbour[faceI];
 142
 143                 // min values
 144                 phiMinIn[own] = min(phiMinIn[own], phi_[nei]);
 145                 phiMinIn[nei] = min(phiMinIn[nei], phi_[own]);
 146
 147                 // max values
 148                 phiMaxIn[own] = max(phiMaxIn[own], phi_[nei]);
 149                 phiMaxIn[nei] = max(phiMaxIn[nei], phi_[own]);
 150             }
 151
 152             // Coupled boundaries
 153             forAll (phi.boundaryField(), patchI)
 154             {
 155                 if (phi.boundaryField()[patchI].coupled())
 156                 {
 157                     const Field<Type> pNei =
 158                         phi.boundaryField()[patchI].patchNeighbourField();
 159
 160                     const labelList& fc =
 161                         phi.boundaryField()[patchI].patch().faceCells();
 162
 163                     forAll (fc, faceI)
 164                     {
 165                         const label& curFC = fc[faceI];
 166
 167                         // min value
 168                         phiMinIn[curFC] =
 169                             min(phiMinIn[curFC], pNei[faceI]);
 170
 171                         // max value
 172                         phiMaxIn[curFC] =
 173                             max(phiMaxIn[curFC], pNei[faceI]);
 174                     }
 175                 }
 176             }
 177
 178             // Calculate limiter
 179
 180             // Get geometrical information
 181
 182             const DimensionedField<scalar, volMesh>& cellVolume = mesh_.V();
 183             const volVectorField& cellCentre = mesh_.C();
 184             const surfaceVectorField& faceCentre = mesh_.Cf();
 185
 186             // Create limiter function
 187             LimiterFunctionType limitFunction;
 188
 189             FieldType& phiLimiterIn = phiLimiter_.internalField();
 190
 191             const GradFieldType& gradPhiIn = gradPhi.internalField();
 192
 193             // Compute limiter values, internal faces
 194
 195             // Dummy flux
 196             scalar upwindFlux = 1;
 197
 198             forAll (owner, faceI)
 199             {
 200                 const label& own = owner[faceI];
 201                 const label& nei = neighbour[faceI];
 202
 203                 vector deltaRLeft = faceCentre[faceI] - cellCentre[own];
 204                 vector deltaRRight = faceCentre[faceI] - cellCentre[nei];
 205
 206                 Type phiOwnerLimiter = limitFunction.limiter
 207                 (
 208                     cellVolume[own],
 209                     upwindFlux,
 210                     phiMaxIn[own] - phi_[own],
 211                     phiMinIn[own] - phi_[own],
 212                     gradPhiIn[own],
 213                     gradPhiIn[nei],
 214                     deltaRLeft
 215                 );
 216
 217                 Type phiNeighbourLimiter = limitFunction.limiter
 218                 (
 219                     cellVolume[nei],
 220                     upwindFlux,
 221                     phiMaxIn[nei] - phi_[nei],
 222                     phiMinIn[nei] - phi_[nei],
 223                     gradPhiIn[nei],
 224                     gradPhiIn[own],
 225                     deltaRRight
 226                 );
 227
 228                 // find minimal limiter value in each cell
 229                 phiLimiterIn[own] =
 230                     min(phiLimiterIn[own], phiOwnerLimiter);
 231
 232                 phiLimiterIn[nei] =
 233                     min(phiLimiterIn[nei], phiNeighbourLimiter);
 234             }
 235
 236             // Coupled boundaries
 237             forAll (phi.boundaryField(), patchI)
 238             {
 239                 if (phi.boundaryField()[patchI].coupled())
 240                 {
 241                     // Get patch
 242                     const fvPatch& p = phi.boundaryField()[patchI].patch();
 243
 244                     const FieldType pNei =
 245                         phi.boundaryField()[patchI].patchNeighbourField();
 246
 247                     const labelList& fc = p.faceCells();
 248
 249                     const vectorField deltaR = p.Cf() - p.Cn();
 250
 251                     // Get gradients
 252                     const GradFieldType gradPhiLeft =
 253                         gradPhi.boundaryField()[patchI].patchInternalField();
 254
 255                     const GradFieldType gradPhiRight =
 256                         gradPhi.boundaryField()[patchI].patchNeighbourField();
 257
 258                     forAll (fc, faceI)
 259                     {
 260                         const label& curFC = fc[faceI];
 261
 262                         Type phiFCLimiter = limitFunction.limiter
 263                         (
 264                             cellVolume[curFC],
 265                             upwindFlux,
 266                             phiMaxIn[curFC] - phi_[curFC],
 267                             phiMinIn[curFC] - phi_[curFC],
 268                             gradPhiLeft[faceI],
 269                             gradPhiRight[faceI],
 270                             deltaR[faceI]
 271                         );
 272
 273                         // Find minimal limiter value in each cell
 274                         phiLimiterIn[curFC] =
 275                             min(phiLimiterIn[curFC], phiFCLimiter);
 276                     }
 277                 }
 278             }
 279
 280             // Do parallel communication to correct limiter on
 281             // coupled boundaries
 282             phiLimiter_.correctBoundaryConditions();
 283         }
 284
 285
 286     // Destructor - default
 287
 288
 289     // Member functions
 290
 291         //- Return limiter
 292         const GeoFieldType& phiLimiter() const
 293         {
 294             return phiLimiter_;
 295         }
 296 };
 297
 298
 299 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 300
 301 } // End namespace Foam
 302
 303 #endif
 304
 305 // ************************************************************************* //