FreeFOAM The Cross-Platform CFD Toolkit
adjustPhi.C
Go to the documentation of this file.
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>
31 
32 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
33 
34 bool Foam::adjustPhi
35 (
37  const volVectorField& U,
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  {
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];
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: ************************ //