adjustPhi.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
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 "volFields.H"
28 #include "surfaceFields.H"
30 
31 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
32 
33 bool Foam::adjustPhi
34 (
35  surfaceScalarField& phi,
36  const volVectorField& U,
38 )
39 {
40  if (p.needReference())
41  {
42  scalar massIn = 0.0;
43  scalar fixedMassOut = 0.0;
44  scalar adjustableMassOut = 0.0;
45 
46  surfaceScalarField::Boundary& bphi =
47  phi.boundaryFieldRef();
48 
49  forAll(bphi, patchi)
50  {
51  const fvPatchVectorField& Up = U.boundaryField()[patchi];
52  const fvsPatchScalarField& phip = bphi[patchi];
53 
54  if (!phip.coupled())
55  {
56  if (Up.fixesValue() && !isA<inletOutletFvPatchVectorField>(Up))
57  {
58  forAll(phip, i)
59  {
60  if (phip[i] < 0.0)
61  {
62  massIn -= phip[i];
63  }
64  else
65  {
66  fixedMassOut += phip[i];
67  }
68  }
69  }
70  else
71  {
72  forAll(phip, i)
73  {
74  if (phip[i] < 0.0)
75  {
76  massIn -= phip[i];
77  }
78  else
79  {
80  adjustableMassOut += phip[i];
81  }
82  }
83  }
84  }
85  }
86 
87  // Calculate the total flux in the domain, used for normalisation
88  scalar totalFlux = vSmall + sum(mag(phi)).value();
89 
90  reduce(massIn, sumOp<scalar>());
91  reduce(fixedMassOut, sumOp<scalar>());
92  reduce(adjustableMassOut, sumOp<scalar>());
93 
94  scalar massCorr = 1.0;
95  scalar magAdjustableMassOut = mag(adjustableMassOut);
96 
97  if
98  (
99  magAdjustableMassOut > vSmall
100  && magAdjustableMassOut/totalFlux > small
101  )
102  {
103  massCorr = (massIn - fixedMassOut)/adjustableMassOut;
104  }
105  else if (mag(fixedMassOut - massIn)/totalFlux > 1e-8)
106  {
108  << "Continuity error cannot be removed by adjusting the"
109  " outflow.\nPlease check the velocity boundary conditions"
110  " and/or run potentialFoam to initialise the outflow." << nl
111  << "Total flux : " << totalFlux << nl
112  << "Specified mass inflow : " << massIn << nl
113  << "Specified mass outflow : " << fixedMassOut << nl
114  << "Adjustable mass outflow : " << adjustableMassOut << nl
115  << exit(FatalError);
116  }
117 
118  forAll(bphi, patchi)
119  {
120  const fvPatchVectorField& Up = U.boundaryField()[patchi];
121  fvsPatchScalarField& phip = bphi[patchi];
122 
123  if (!phip.coupled())
124  {
125  if
126  (
127  !Up.fixesValue()
128  || isA<inletOutletFvPatchVectorField>(Up)
129  )
130  {
131  forAll(phip, i)
132  {
133  if (phip[i] > 0.0)
134  {
135  phip[i] *= massCorr;
136  }
137  }
138  }
139  }
140  }
141 
142  return mag(massIn)/totalFlux < small
143  && mag(fixedMassOut)/totalFlux < small
144  && mag(adjustableMassOut)/totalFlux < small;
145  }
146  else
147  {
148  return false;
149  }
150 }
151 
152 
153 // ************************************************************************* //
Foam::surfaceFields.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
virtual bool fixesValue() const
Return true if this patch field fixes a value.
Definition: fvPatchField.H:315
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:66
bool adjustPhi(surfaceScalarField &phi, const volVectorField &U, volScalarField &p)
Adjust the balance of fluxes to obey continuity.
Definition: adjustPhi.C:34
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
bool needReference() const
Does the field need a reference level for solution.
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:98
virtual bool coupled() const
Return true if this patch field is coupled.
static const char nl
Definition: Ostream.H:265
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
label patchi
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
For cases which do no have a pressure boundary adjust the balance of fluxes to obey continuity...
dimensioned< scalar > mag(const dimensioned< Type > &)
An abstract base class with a fat-interface to all derived classes covering all possible ways in whic...
Definition: fvsPatchField.H:65