porosityModel.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) 2012-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 "porosityModel.H"
27 #include "volFields.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  defineTypeNameAndDebug(porosityModel, 0);
34  defineRunTimeSelectionTable(porosityModel, mesh);
35 }
36 
37 
38 // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
39 
41 {
42  scalar maxCmpt = max(0, cmptMax(resist.value()));
43 
44  if (maxCmpt < 0)
45  {
47  << "Negative resistances are invalid, resistance = " << resist
48  << exit(FatalError);
49  }
50  else
51  {
52  vector& val = resist.value();
53  for (label cmpt = 0; cmpt < vector::nComponents; cmpt++)
54  {
55  if (val[cmpt] < 0)
56  {
57  val[cmpt] *= -maxCmpt;
58  }
59  }
60  }
61 }
62 
63 
65 {
66  label index = 0;
67  if (!coordSys_.R().uniform())
68  {
69  index = i;
70  }
71  return index;
72 }
73 
74 
75 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
76 
77 Foam::porosityModel::porosityModel
78 (
79  const word& name,
80  const word& modelType,
81  const fvMesh& mesh,
82  const dictionary& dict,
83  const word& cellZoneName
84 )
85 :
87  (
88  IOobject
89  (
90  name,
91  mesh.time().timeName(),
92  mesh,
95  )
96  ),
97  name_(name),
98  mesh_(mesh),
99  dict_(dict),
100  coeffs_(dict.optionalSubDict(modelType + "Coeffs")),
101  active_(true),
102  zoneName_(cellZoneName),
103  cellZoneIDs_(),
104  coordSys_(coordinateSystem::New(mesh, coeffs_))
105 {
106  if (zoneName_ == word::null)
107  {
108  dict.readIfPresent("active", active_);
109  dict_.lookup("cellZone") >> zoneName_;
110  }
111 
112  cellZoneIDs_ = mesh_.cellZones().findIndices(zoneName_);
113 
114  Info<< " creating porous zone: " << zoneName_ << endl;
115 
116  bool foundZone = !cellZoneIDs_.empty();
117  reduce(foundZone, orOp<bool>());
118 
119  if (!foundZone && Pstream::master())
120  {
122  << "cannot find porous cellZone " << zoneName_
123  << exit(FatalError);
124  }
125 }
126 
127 
128 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
129 
131 {}
132 
133 
134 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
135 
137 {
138  if (!mesh_.upToDatePoints(*this))
139  {
140  calcTransformModelData();
141 
142  // set model up-to-date wrt points
143  mesh_.setUpToDatePoints(*this);
144  }
145 }
146 
147 
148 Foam::tmp<Foam::vectorField> Foam::porosityModel::porosityModel::force
149 (
150  const volVectorField& U,
151  const volScalarField& rho,
152  const volScalarField& mu
153 )
154 {
155  transformModelData();
156 
157  tmp<vectorField> tforce(new vectorField(U.size(), Zero));
158 
159  if (!cellZoneIDs_.empty())
160  {
161  this->calcForce(U, rho, mu, tforce.ref());
162  }
163 
164  return tforce;
165 }
166 
167 
169 {
170  if (cellZoneIDs_.empty())
171  {
172  return;
173  }
174 
175  transformModelData();
176  this->correct(UEqn);
177 }
178 
179 
181 (
182  fvVectorMatrix& UEqn,
183  const volScalarField& rho,
184  const volScalarField& mu
185 )
186 {
187  if (cellZoneIDs_.empty())
188  {
189  return;
190  }
191 
192  transformModelData();
193  this->correct(UEqn, rho, mu);
194 }
195 
196 
198 (
199  const fvVectorMatrix& UEqn,
200  volTensorField& AU,
201  bool correctAUprocBC
202 )
203 {
204  if (cellZoneIDs_.empty())
205  {
206  return;
207  }
208 
209  transformModelData();
210  this->correct(UEqn, AU);
211 
212  if (correctAUprocBC)
213  {
214  // Correct the boundary conditions of the tensorial diagonal to ensure
215  // processor boundaries are correctly handled when AU^-1 is interpolated
216  // for the pressure equation.
218  }
219 }
220 
221 
223 {
224  return true;
225 }
226 
227 
229 {
230  dict.readIfPresent("active", active_);
231 
232  coeffs_ = dict.optionalSubDict(type() + "Coeffs");
233 
234  dict.lookup("cellZone") >> zoneName_;
235  cellZoneIDs_ = mesh_.cellZones().findIndices(zoneName_);
236 
237  return true;
238 }
239 
240 
241 // ************************************************************************* //
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
virtual bool writeData(Ostream &os) const
Write.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
virtual bool read()
Read object.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:421
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
label fieldIndex(const label index) const
Return label index.
Definition: porosityModel.C:64
void adjustNegativeResistance(dimensionedVector &resist)
Adjust negative resistance values to be multiplier of max value.
Definition: porosityModel.C:40
static const direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:96
dynamicFvMesh & mesh
const dictionary & optionalSubDict(const word &) const
Find and return a sub-dictionary if found.
Definition: dictionary.C:759
A class for handling words, derived from string.
Definition: word.H:59
const Type & value() const
Return const reference to value.
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
static const word null
An empty word.
Definition: word.H:77
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:72
static const zero Zero
Definition: zero.H:97
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
virtual void transformModelData()
Transform the model data wrt mesh changes.
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil) - psi) *pSat, rhoMin);# 1 "/home/ubuntu/OpenFOAM-6/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
Definition: pEqn.H:72
defineTypeNameAndDebug(combustionModel, 0)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
fileName::Type type(const fileName &, const bool followLink=true)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:481
virtual ~porosityModel()
Destructor.
static autoPtr< coordinateSystem > New(const objectRegistry &obr, const dictionary &dict)
Select constructed from dictionary and objectRegistry.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
void correctBoundaryConditions()
Correct boundary field.
regIOobject is an abstract class derived from IOobject to handle automatic object registration with t...
Definition: regIOobject.H:65
messageStream Info
Field< vector > vectorField
Specialisation of Field<T> for vector.
A class for managing temporary objects.
Definition: PtrList.H:53
virtual void addResistance(fvVectorMatrix &UEqn)
Add resistance.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:576