src/dbns/numericFlux/fineNumericFlux.C

   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 \*---------------------------------------------------------------------------*/
  25
  26 #include "fineNumericFlux.H"
  27
  28 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  29
  30 // Construct from components
  31 template<class Flux, class Limiter>
  32 Foam::fineNumericFlux<Flux, Limiter>::fineNumericFlux
  33 (
  34     const mgMeshLevel& meshLevel,
  35     const mgFieldLevel& fieldLevel
  36     basicThermo& thermo,
  37 )
  38 :
  39     meshLevel_(meshLevel),
  40     fieldLevel_(fieldLevel),
  41     p_(fieldLevel.p()),
  42     U_(fieldLevel.U()),
  43     T_(fieldLevel.T()),
  44     thermo_(thermo),
  45     rhoFlux_(fieldLevel.rhoFlux()),
  46     rhoUFlux_(fieldLevel.rhoUFlux()),
  47     rhoEFlux_(fieldLevel.rhoEFlux()),
  48     gradP_(fvc::grad(p_)),
  49     gradU_(fvc::grad(U_)),
  50     gradT_(fvc::grad(T_))
  51 {}
  52
  53
  54 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  55
  56 template<class Flux, class Limiter>
  57 void Foam::fineNumericFlux<Flux, Limiter>::computeInteriorFlux()
  58 {
  59     // Get face-to-cell addressing: face area point from owner to neighbour
  60     const unallocLabelList& owner = meshLevel_.owner();
  61     const unallocLabelList& neighbour = meshLevel_.neighbour();
  62
  63     // Get the face area vector
  64     const surfaceVectorField& Sf = mesh().Sf();
  65     const surfaceScalarField& magSf = mesh().magSf();
  66
  67     const scalarField& Cv = fieldLevel_.CvVar();
  68     const scalarField& R  = fieldLevel_.R();
  69
  70     gradP_ = fvc::grad(p_);
  71     gradP_.correctBoundaryConditions();
  72
  73     gradU_ = fvc::grad(U_);
  74     gradU_.correctBoundaryConditions();
  75
  76     gradT_ = fvc::grad(T_);
  77     gradT_.correctBoundaryConditions();
  78
  79     MDLimiter<scalar, Limiter> scalarPLimiter
  80     (
  81         this->p_,
  82         this->gradP_
  83     );
  84
  85     MDLimiter<vector, Limiter> vectorULimiter
  86     (
  87         this->U_,
  88         this->gradU_
  89     );
  90
  91     MDLimiter<scalar, Limiter> scalarTLimiter
  92     (
  93         this->T_,
  94         this->gradT_
  95     );
  96
  97     // Get limiters
  98     const volScalarField& pLimiter = scalarPLimiter.phiLimiter();
  99     const volVectorField& ULimiter = vectorULimiter.phiLimiter();
 100     const volScalarField& TLimiter = scalarTLimiter.phiLimiter();
 101
 102     // Calculate fluxes at internal faces
 103     forAll (owner, faceI)
 104     {
 105         const label own = owner[faceI];
 106         const label nei = neighbour[faceI];
 107
 108         const vector deltaRLeft = faceCentre[faceI] - cellCentre[own];
 109         const vector deltaRRight = faceCentre[faceI] - cellCentre[nei];
 110
 111         // calculate fluxes with reconstructed primitive variables at faces
 112         Flux::evaluateFlux
 113         (
 114             rhoFlux_[faceI],
 115             rhoUFlux_[faceI],
 116             rhoEFlux_[faceI],
 117             p_[own] + pLimiter[own]*(deltaRLeft & gradP_[own]),
 118             p_[nei] + pLimiter[nei]*(deltaRRight & gradP_[nei]),
 119             U_[own] + cmptMultiply(ULimiter[own], (deltaRLeft & gradU_[own])),
 120             U_[nei] + cmptMultiply(ULimiter[nei], (deltaRRight & gradU_[nei])),
 121             T_[own] + TLimiter[own]*(deltaRLeft & gradT_[own]),
 122             T_[nei] + TLimiter[nei]*(deltaRRight & gradT_[nei]),
 123             R[own],
 124             R[nei],
 125             Cv[own],
 126             Cv[nei],
 127             Sf[faceI],
 128             magSf[faceI]
 129         );
 130     }
 131
 132     // Update boundary field and values
 133     forAll (rhoFlux_.boundaryField(), patchi)
 134     {
 135         const fvPatch& curPatch = p_.boundaryField()[patchi].patch();
 136
 137         // Fluxes
 138         fvsPatchScalarField& pRhoFlux  = rhoFlux_.boundaryField()[patchi];
 139         fvsPatchVectorField& pRhoUFlux = rhoUFlux_.boundaryField()[patchi];
 140         fvsPatchScalarField& pRhoEFlux = rhoEFlux_.boundaryField()[patchi];
 141
 142         const scalarField& pp = fieldLevel_.patchP(patchi);
 143         const vectorField& pU = fieldLevel_.patchU(patchi);
 144         const scalarField& pT = fieldLevel_.patchT(patchi);
 145
 146         const scalarField& pCv = fieldLevel_.patchCv(patchi);
 147         const scalarField& pR = fieldLevel_.patchR(patchi);
 148
 149         // Gradients
 150         const fvPatchVectorField& pGradP = gradP_.boundaryField()[patchi];
 151         const fvPatchTensorField& pGradU = gradU_.boundaryField()[patchi];
 152         const fvPatchVectorField& pGradT = gradT_.boundaryField()[patchi];
 153
 154         // Limiters
 155         const fvPatchScalarField& pPatchLim = pLimiter.boundaryField()[patchi];
 156         const fvPatchVectorField& UPatchLim = ULimiter.boundaryField()[patchi];
 157         const fvPatchScalarField& TPatchLim = TLimiter.boundaryField()[patchi];
 158
 159         // Face areas
 160         const vectorField& pSf  = meshLevel_.patchFaceAreas(patchi);
 161         const scalarField& pMagSf = meshLevel_.magPatchFaceAreas(patchi);
 162
 163         if (p_.boundaryField()[patchi].coupled())
 164         {
 165             // Coupled patch
 166             const scalarField ppLeft  =
 167                 p_.boundaryField()[patchi].patchInternalField();
 168
 169             const scalarField ppRight =
 170                 p_.boundaryField()[patchi].patchNeighbourField();
 171
 172             const vectorField pULeft  =
 173                 U_.boundaryField()[patchi].patchInternalField();
 174
 175             const vectorField pURight =
 176                 U_.boundaryField()[patchi].patchNeighbourField();
 177
 178             const scalarField pTLeft  =
 179                 T_.boundaryField()[patchi].patchInternalField();
 180
 181             const scalarField pTRight =
 182                 T_.boundaryField()[patchi].patchNeighbourField();
 183
 184             // Gradients
 185             const vectorField pgradPLeft = pGradP.patchInternalField();
 186             const vectorField pgradPRight = pGradP.patchNeighbourField();
 187
 188             const tensorField pgradULeft = pGradU.patchInternalField();
 189             const tensorField pgradURight = pGradU.patchNeighbourField();
 190
 191             const vectorField pgradTLeft = pGradT.patchInternalField();
 192             const vectorField pgradTRight = pGradT.patchNeighbourField();
 193
 194             // Geometry
 195             vectorField pDeltaRLeft = meshLevel_.patchDeltaR(patchi);
 196             vectorField pDdeltaRRight =
 197                 meshLevel_.patchDeltaRNeighbour(patchi);
 198
 199             // Limiters
 200
 201             const scalarField ppLimiterLeft = pPatchLim.patchInternalField();
 202             const scalarField ppLimiterRight = pPatchLim.patchNeighbourField();
 203
 204             const vectorField pULimiterLeft = UPatchLim.patchInternalField();
 205             const vectorField pULimiterRight = UPatchLim.patchNeighbourField();
 206
 207             const scalarField pTLimiterLeft = TPatchLim.patchInternalField();
 208             const scalarField pTLimiterRight = TPatchLim.patchNeighbourField();
 209
 210             forAll (pp, facei)
 211             {
 212                 Flux::evaluateFlux
 213                 (
 214                     pRhoFlux[facei],
 215                     pRhoUFlux[facei],
 216                     pRhoEFlux[facei],
 217
 218                     ppLeft[facei]
 219                   + ppLimiterLeft[facei]*
 220                     (pDeltaRLeft[facei] & pgradPLeft[facei]),
 221
 222                     ppRight[facei]
 223                   + ppLimiterRight[facei]*
 224                     (pDdeltaRRight[facei] & pgradPRight[facei]),
 225
 226                     pULeft[facei]
 227                   + cmptMultiply
 228                     (
 229                         pULimiterLeft[facei],
 230                         pDeltaRLeft[facei] & pgradULeft[facei]
 231                     ),
 232
 233                     pURight[facei]
 234                   + cmptMultiply
 235                     (
 236                         pULimiterRight[facei],
 237                         pDdeltaRRight[facei] & pgradURight[facei]
 238                     ),
 239
 240                     pTLeft[facei]
 241                   + pTLimiterLeft[facei]*
 242                     (pDeltaRLeft[facei] & pgradTLeft[facei]),
 243
 244                     pTRight[facei]
 245                   + pTLimiterRight[facei]*
 246                     (pDdeltaRRight[facei] & pgradTRight[facei]),
 247
 248                     pR[facei],
 249                     pR[facei],
 250                     pCv[facei],
 251                     pCv[facei],
 252                     pSf[facei],
 253                     pMagSf[facei]
 254                 );
 255             }
 256         }
 257         else
 258         {
 259             forAll (pp, facei)
 260             {
 261                 // Calculate fluxes
 262                 Flux::evaluateFlux
 263                 (
 264                     pRhoFlux[facei],
 265                     pRhoUFlux[facei],
 266                     pRhoEFlux[facei],
 267                     pp[facei],
 268                     pp[facei],
 269                     pU[facei],
 270                     pU[facei],
 271                     pT[facei],
 272                     pT[facei],
 273                     pR[facei],
 274                     pR[facei],
 275                     pCv[facei],
 276                     pCv[facei],
 277                     pSf[facei],
 278                     pMagSf[facei]
 279                 );
 280             }
 281         }
 282     }
 283 }
 284
 285
 286 // ************************************************************************* //