src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2011 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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "adjustPhi.H"
  27 #include <finiteVolume/volFields.H>
  28 #include <finiteVolume/surfaceFields.H>
  29 #include <finiteVolume/processorFvsPatchFields.H>
  30 #include <finiteVolume/inletOutletFvPatchFields.H>
  31
  32 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
  33
  34 bool Foam::adjustPhi
  35 (
  36     surfaceScalarField& phi,
  37     const volVectorField& U,
  38     volScalarField& p
  39 )
  40 {
  41     if (p.needReference())
  42     {
  43         // p coefficients should not be updated here
  44         // p.boundaryField().updateCoeffs();
  45
  46         scalar massIn = 0.0;
  47         scalar fixedMassOut = 0.0;
  48         scalar adjustableMassOut = 0.0;
  49
  50         forAll (phi.boundaryField(), patchi)
  51         {
  52             const fvPatchVectorField& Up = U.boundaryField()[patchi];
  53             const fvsPatchScalarField& phip = phi.boundaryField()[patchi];
  54
  55             if (!isA<processorFvsPatchScalarField>(phip))
  56             {
  57                 if
  58                 (
  59                     Up.fixesValue()
  60                  && !isA<inletOutletFvPatchVectorField>(Up)
  61                 )
  62                 {
  63                     forAll(phip, i)
  64                     {
  65                         if (phip[i] < 0.0)
  66                         {
  67                             massIn -= phip[i];
  68                         }
  69                         else
  70                         {
  71                             fixedMassOut += phip[i];
  72                         }
  73                     }
  74                 }
  75                 else
  76                 {
  77                     forAll(phip, i)
  78                     {
  79                         if (phip[i] < 0.0)
  80                         {
  81                             massIn -= phip[i];
  82                         }
  83                         else
  84                         {
  85                             adjustableMassOut += phip[i];
  86                         }
  87                     }
  88                 }
  89             }
  90         }
  91
  92         // Calculate the total flux in the domain, used for normalisation
  93         scalar totalFlux = VSMALL + sum(mag(phi)).value();
  94
  95         reduce(massIn, sumOp<scalar>());
  96         reduce(fixedMassOut, sumOp<scalar>());
  97         reduce(adjustableMassOut, sumOp<scalar>());
  98
  99         scalar massCorr = 1.0;
 100         scalar magAdjustableMassOut = mag(adjustableMassOut);
 101
 102         if
 103         (
 104             magAdjustableMassOut > VSMALL
 105          && magAdjustableMassOut/totalFlux > SMALL
 106         )
 107         {
 108             massCorr = (massIn - fixedMassOut)/adjustableMassOut;
 109         }
 110         else if(mag(fixedMassOut - massIn)/totalFlux > 1e-10)
 111         {
 112             FatalErrorIn
 113             (
 114                 "adjustPhi(surfaceScalarField& phi, const volVectorField& U,"
 115                 "const volScalarField& p"
 116             )   << "Continuity error cannot be removed by adjusting the"
 117                    " outflow.\nPlease check the velocity boundary conditions"
 118                    " and/or run potentialFoam to initialise the outflow." << nl
 119                 << "Total flux              : " << totalFlux << nl
 120                 << "Specified mass inflow   : " << massIn << nl
 121                 << "Specified mass outflow  : " << fixedMassOut << nl
 122                 << "Adjustable mass outflow : " << adjustableMassOut << nl
 123                 << exit(FatalError);
 124         }
 125
 126         forAll (phi.boundaryField(), patchi)
 127         {
 128             const fvPatchVectorField& Up = U.boundaryField()[patchi];
 129             fvsPatchScalarField& phip = phi.boundaryField()[patchi];
 130
 131             if (!isA<processorFvsPatchScalarField>(phip))
 132             {
 133                 if
 134                 (
 135                     !Up.fixesValue()
 136                  || isA<inletOutletFvPatchVectorField>(Up)
 137                 )
 138                 {
 139                     forAll(phip, i)
 140                     {
 141                         if (phip[i] > 0.0)
 142                         {
 143                             phip[i] *= massCorr;
 144                         }
 145                     }
 146                 }
 147             }
 148         }
 149
 150         return mag(massIn)/totalFlux < SMALL
 151             && mag(fixedMassOut)/totalFlux < SMALL
 152             && mag(adjustableMassOut)/totalFlux < SMALL;
 153     }
 154     else
 155     {
 156         return false;
 157     }
 158 }
 159
 160
 161 // ************************ vim: set sw=4 sts=4 et: ************************ //