src/dbns/numericFlux/FASNumericFlux.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 Class
  25     numericFlux
  26
  27 Description
  28     Generic Godunov flux class for coarse level multigrid
  29
  30 Author
  31     Aleksandar Jemcov
  32     Rewrite by Hrvoje Jasak
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #include "FASNumericFlux.H"
  37
  38 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  39
  40 // Construct from components
  41 template<class Flux>
  42 Foam::FASNumericFlux<Flux>::FASNumericFlux
  43 (
  44     mgMeshLevel const& meshLevel,
  45     mgFieldLevel const& fieldLevel
  46 )
  47 :
  48     meshLevel_(meshLevel),
  49     fieldLevel_(fieldLevel),
  50     p_(fieldLevel.pVar()),
  51     U_(fieldLevel.UVar()),
  52     T_(fieldLevel.TVar()),
  53     fineRhoFlux_(fieldLevel.rhoFlux()),
  54     fineRhoUFlux_(fieldLevel.rhoUFlux()),
  55     fineRhoEFlux_(fieldLevel.rhoEFlux()),
  56     rhoFlux_(scalarField(meshLevel.nInternalFaces())),
  57     rhoUFlux_(vectorField(meshLevel.nInternalFaces())),
  58     rhoEFlux_(scalarField(meshLevel.nInternalFaces())),
  59     rhoResidual_(scalarField(meshLevel.nCells())),
  60     rhoUResidual_(vectorField(meshLevel.nCells())),
  61     rhoEResidual_(scalarField(meshLevel.nCells()))
  62 {
  63     rhoFlux_ = 0;
  64     rhoUFlux_ = vector::zero;
  65     rhoEFlux_ = 0;
  66     rhoResidual_ = 0;
  67     rhoUResidual_ = vector::zero;
  68     rhoEResidual_ = 0;
  69     Info<< "mgLevel = " << fieldLevel.level() << endl;
  70     Info<< "p_ size = " << p_.size() << endl;
  71     Info<< "U_ size = " << U_.size() << endl;
  72     Info<< "T_ size = " << T_.size() << endl;
  73     Info<< "fineRhoFlux_ size = " << fineRhoFlux_.size() << endl;
  74     Info<< "fineRhoUFlux_ size = " << fineRhoUFlux_.size() << endl;
  75     Info<< "fineRhoEFlux_ size = " << fineRhoEFlux_.size() << endl;
  76     Info<< "rhoFlux_ size = " << rhoFlux_.size() << endl;
  77     Info<< "rhoUFlux_ size = " << rhoUFlux_.size() << endl;
  78     Info<< "rhoEFlux_ size = " << rhoEFlux_.size() << endl;
  79     Info<< "rhoResidual_ size = " << rhoResidual_.size() << endl;
  80     Info<< "rhoUResidual_ size = " << rhoUResidual_.size() << endl;
  81     Info<< "rhoEResidual_ size = " << rhoEResidual_.size() << endl;
  82 }
  83
  84
  85 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  86
  87 template<class Flux>
  88 inline
  89 void Foam::FASNumericFlux<Flux>::computeFlux()
  90 {
  91     // Get face-to-cell addressing: face area point from owner to neighbour
  92     const unallocLabelList& owner     = meshLevel_.owner();
  93     const unallocLabelList& neighbour = meshLevel_.neighbour();
  94
  95     // Get the face area vector
  96     const vectorField& Sf    = meshLevel_.faceAreas();
  97     const scalarField& magSf = meshLevel_.magFaceAreas();
  98
  99     const scalarField& Cv = fieldLevel_.CvVar();
 100     const scalarField& R  = fieldLevel_.R();
 101
 102     const vectorField& cellCenter = meshLevel_.cellCentres();
 103     const vectorField& faceCenter = meshLevel_.faceCentres();
 104
 105     // Calculate fluxes at internal faces
 106     forAll(owner, faceI)
 107     {
 108         label own = owner[faceI];
 109         label nei = neighbour[faceI];
 110
 111         vector deltaRLeft  = faceCenter[faceI] - cellCenter[own];
 112         vector deltaRRight = faceCenter[faceI] - cellCenter[nei];
 113
 114         // calculate fluxes with reconstructed primitive variables at faces
 115         Flux::evaluateFlux
 116         (
 117             rhoFlux_[faceI],
 118             rhoUFlux_[faceI],
 119             rhoEFlux_[faceI],
 120             p_[own],
 121             p_[nei],
 122             U_[own],
 123             U_[nei],
 124             T_[own],
 125             T_[nei],
 126             R [own],
 127             R [nei],
 128             Cv[own],
 129             Cv[nei],
 130             Sf[faceI],
 131             magSf[faceI]
 132         );
 133         rhoResidual_[own] += rhoFlux_[faceI];
 134         rhoResidual_[nei] -= rhoFlux_[faceI];
 135         rhoUResidual_[own] += rhoUFlux_[faceI];
 136         rhoUResidual_[nei] -= rhoUFlux_[faceI];
 137         rhoEResidual_[own] += rhoEFlux_[faceI];
 138         rhoEResidual_[nei] -= rhoEFlux_[faceI];
 139     }
 140     Info<< "Done with interior on level " << fieldLevel_.level() << endl;
 141
 142     // Update boundary field and values
 143     forAll(fineRhoFlux_.boundaryField(), patchi)
 144     {
 145         unallocLabelList const& owner = meshLevel_.faceCells(patchi);
 146
 147         const scalarField& pp = fieldLevel_.patchP(patchi);
 148         const vectorField& pU = fieldLevel_.patchU(patchi);
 149         const scalarField& pT = fieldLevel_.patchT(patchi);
 150
 151         const scalarField& pCv = fieldLevel_.patchCv(patchi);
 152         const scalarField& pR  = fieldLevel_.patchR(patchi);
 153
 154         const vectorField& pSf    = meshLevel_.patchFaceAreas(patchi);
 155         const scalarField& pMagSf = meshLevel_.magPatchFaceAreas(patchi);
 156
 157         if (fieldLevel_.p().boundaryField()[patchi].coupled())
 158         {
 159             const scalarField ppLeft  =
 160                 fieldLevel_.p().boundaryField()[patchi].patchInternalField();
 161             const scalarField ppRight =
 162                 fieldLevel_.p().boundaryField()[patchi].patchNeighbourField();
 163             const vectorField pULeft  =
 164                 fieldLevel_.U().boundaryField()[patchi].patchInternalField();
 165             const vectorField pURight =
 166                 fieldLevel_.U().boundaryField()[patchi].patchNeighbourField();
 167             const scalarField pTLeft  =
 168                 fieldLevel_.T().boundaryField()[patchi].patchInternalField();
 169             const scalarField pTRight =
 170                 fieldLevel_.T().boundaryField()[patchi].patchNeighbourField();
 171
 172             forAll(owner, facei)
 173             {
 174                 label own = owner[facei];
 175                 Flux::evaluateFlux
 176                 (
 177                     rhoFlux_[facei],
 178                     rhoUFlux_[facei],
 179                     rhoEFlux_[facei],
 180                     ppLeft[facei],
 181                     ppRight[facei],
 182                     pULeft[facei],
 183                     pURight[facei],
 184                     pTLeft[facei],
 185                     pTRight[facei],
 186                     pR[facei],
 187                     pR[facei],
 188                     pCv[facei],
 189                     pCv[facei],
 190                     pSf[facei],
 191                     pMagSf[facei]
 192                 );
 193                 rhoResidual_[own]  += rhoFlux_[facei];
 194                 rhoUResidual_[own] += rhoUFlux_[facei];
 195                 rhoEResidual_[own] += rhoEFlux_[facei];
 196             }
 197         }
 198         else
 199         {
 200             forAll(owner, facei)
 201             {
 202                 label own = owner[facei];
 203                 Flux::evaluateFlux
 204                 (
 205                     rhoFlux_[facei],
 206                     rhoUFlux_[facei],
 207                     rhoEFlux_[facei],
 208                     pp[facei],
 209                     pp[facei],
 210                     pU[facei],
 211                     pU[facei],
 212                     pT[facei],
 213                     pT[facei],
 214                     pR[facei],
 215                     pR[facei],
 216                     pCv[facei],
 217                     pCv[facei],
 218                     pSf[facei],
 219                     pMagSf[facei]
 220                 );
 221                 rhoResidual_[own]  += rhoFlux_[facei];
 222                 rhoUResidual_[own] += rhoUFlux_[facei];
 223                 rhoEResidual_[own] += rhoEFlux_[facei];
 224             }
 225         }
 226     }
 227     Info<< "Done with the boundary on level " << fieldLevel_.level() << endl;
 228 }
 229
 230
 231 // ************************************************************************* //
 232