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-2022 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"
29 #include "processorFvPatchFields.H"
31 
32 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
33 
34 bool Foam::adjustPhi
35 (
36  surfaceScalarField& phi,
37  const volVectorField& U,
39 )
40 {
41  if (p.needReference())
42  {
43  scalar massIn = 0.0;
44  scalar fixedMassOut = 0.0;
45  scalar adjustableMassOut = 0.0;
46 
48  phi.boundaryFieldRef();
49 
50  forAll(bphi, patchi)
51  {
52  const fvPatchVectorField& Up = U.boundaryField()[patchi];
53  const fvsPatchScalarField& phip = bphi[patchi];
54 
55  if (!isA<processorFvPatchVectorField>(Up))
56  {
57  if (Up.fixesValue() && !isA<inletOutletFvPatchVectorField>(Up))
58  {
59  forAll(phip, i)
60  {
61  if (phip[i] < 0.0)
62  {
63  massIn -= phip[i];
64  }
65  else
66  {
67  fixedMassOut += phip[i];
68  }
69  }
70  }
71  else
72  {
73  forAll(phip, i)
74  {
75  if (phip[i] < 0.0)
76  {
77  massIn -= phip[i];
78  }
79  else
80  {
81  adjustableMassOut += phip[i];
82  }
83  }
84  }
85  }
86  }
87 
88  // Calculate the total flux in the domain, used for normalisation
89  scalar totalFlux = vSmall + sum(mag(phi)).value();
90 
91  reduce(massIn, sumOp<scalar>());
92  reduce(fixedMassOut, sumOp<scalar>());
93  reduce(adjustableMassOut, sumOp<scalar>());
94 
95  scalar massCorr = 1.0;
96  scalar magAdjustableMassOut = mag(adjustableMassOut);
97 
98  if
99  (
100  magAdjustableMassOut > vSmall
101  && magAdjustableMassOut/totalFlux > small
102  )
103  {
104  massCorr = (massIn - fixedMassOut)/adjustableMassOut;
105  }
106  else if (mag(fixedMassOut - massIn)/totalFlux > 1e-8)
107  {
109  << "Continuity error cannot be removed by adjusting the"
110  " outflow.\nPlease check the velocity boundary conditions"
111  " and/or run potentialFoam to initialise the outflow." << nl
112  << "Total flux : " << totalFlux << nl
113  << "Specified mass inflow : " << massIn << nl
114  << "Specified mass outflow : " << fixedMassOut << nl
115  << "Adjustable mass outflow : " << adjustableMassOut << nl
116  << exit(FatalError);
117  }
118 
119  forAll(bphi, patchi)
120  {
121  const fvPatchVectorField& Up = U.boundaryField()[patchi];
122  fvsPatchScalarField& phip = bphi[patchi];
123 
124  if (!phip.coupled())
125  {
126  if
127  (
128  !Up.fixesValue()
129  || isA<inletOutletFvPatchVectorField>(Up)
130  )
131  {
132  forAll(phip, i)
133  {
134  if (phip[i] > 0.0)
135  {
136  phip[i] *= massCorr;
137  }
138  }
139  }
140  }
141  }
142 
143  return mag(massIn)/totalFlux < small
144  && mag(fixedMassOut)/totalFlux < small
145  && mag(adjustableMassOut)/totalFlux < small;
146  }
147  else
148  {
149  return false;
150  }
151 }
152 
153 
154 // ************************************************************************* //
Foam::surfaceFields.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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:310
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
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:35
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
bool needReference() const
Does the field need a reference level for solution.
virtual bool coupled() const
Return true if this patch field is coupled.
static const char nl
Definition: Ostream.H:260
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
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105