setFlow.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2017-2020 OpenCFD Ltd.
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "setFlow.H"
29 #include "volFields.H"
30 #include "surfaceFields.H"
31 #include "fvcFlux.H"
32 #include "addToRunTimeSelectionTable.H"
33 #include "fvcSurfaceIntegrate.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 namespace functionObjects
40 {
41  defineTypeNameAndDebug(setFlow, 0);
42  addToRunTimeSelectionTable(functionObject, setFlow, dictionary);
43 }
44 }
45 
46 
47 const Foam::Enum
48 <
49  Foam::functionObjects::setFlow::modeType
50 >
51 Foam::functionObjects::setFlow::modeTypeNames
52 ({
53  { functionObjects::setFlow::modeType::FUNCTION, "function" },
54  { functionObjects::setFlow::modeType::ROTATION, "rotation" },
55  { functionObjects::setFlow::modeType::VORTEX2D, "vortex2D" },
56  { functionObjects::setFlow::modeType::VORTEX3D, "vortex3D" },
57 });
58 
59 
60 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
61 
62 void Foam::functionObjects::setFlow::setPhi(const volVectorField& U)
63 {
64  surfaceScalarField* phiptr =
65  mesh_.getObjectPtr<surfaceScalarField>(phiName_);
66 
67  if (!phiptr)
68  {
69  return;
70  }
71 
72  if (rhoName_ != "none")
73  {
74  const volScalarField* rhoptr =
75  mesh_.findObject<volScalarField>(rhoName_);
76 
77  if (rhoptr)
78  {
79  const volScalarField& rho = *rhoptr;
80  *phiptr = fvc::flux(rho*U);
81  }
82  else
83  {
84  FatalErrorInFunction
85  << "Unable to find rho field'" << rhoName_
86  << "' in the mesh database. Available fields are:"
87  << mesh_.names<volScalarField>()
88  << exit(FatalError);
89  }
90  }
91  else
92  {
93  *phiptr = fvc::flux(U);
94  }
95 }
96 
97 
98 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
99 
100 Foam::functionObjects::setFlow::setFlow
101 (
102  const word& name,
103  const Time& runTime,
104  const dictionary& dict
105 )
106 :
107  fvMeshFunctionObject(name, runTime, dict),
108  mode_(modeType::FUNCTION),
109  UName_("U"),
110  rhoName_("none"),
111  phiName_("phi"),
112  reverseTime_(VGREAT),
113  scalePtr_(nullptr),
114  origin_(Zero),
115  R_(tensor::I),
116  omegaPtr_(nullptr),
117  velocityPtr_(nullptr)
118 {
119  read(dict);
120 }
121 
122 
123 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
124 
125 bool Foam::functionObjects::setFlow::read(const dictionary& dict)
126 {
127  if (fvMeshFunctionObject::read(dict))
128  {
129  Info<< name() << ":" << endl;
130 
131  modeTypeNames.readEntry("mode", dict, mode_);
132 
133  Info<< " operating mode: " << modeTypeNames[mode_] << endl;
134 
135  if (dict.readIfPresent("U", UName_))
136  {
137  Info<< " U field name: " << UName_ << endl;
138  }
139 
140  if (dict.readIfPresent("rho", rhoName_))
141  {
142  Info<< " rho field name: " << rhoName_ << endl;
143  }
144 
145  if (dict.readIfPresent("phi", phiName_))
146  {
147  Info<< " phi field name: " << phiName_ << endl;
148  }
149 
150  if (dict.readIfPresent("reverseTime", reverseTime_))
151  {
152  Info<< " reverse flow direction at time: " << reverseTime_
153  << endl;
154  reverseTime_ = mesh_.time().userTimeToTime(reverseTime_);
155  }
156 
157  // Scaling is applied across all modes
158  scalePtr_ = Function1<scalar>::New("scale", dict, &mesh_);
159 
160  switch (mode_)
161  {
162  case modeType::FUNCTION:
163  {
164  velocityPtr_ = Function1<vector>::New("velocity", dict, &mesh_);
165  break;
166  }
167  case modeType::ROTATION:
168  {
169  omegaPtr_ = Function1<scalar>::New("omega", dict, &mesh_);
170 
171  dict.readEntry("origin", origin_);
172  const vector refDir(dict.get<vector>("refDir").normalise());
173  const vector axis(dict.get<vector>("axis").normalise());
174 
175  R_ = tensor(refDir, axis, refDir^axis);
176  break;
177  }
178  case modeType::VORTEX2D:
179  case modeType::VORTEX3D:
180  {
181  dict.readEntry("origin", origin_);
182  const vector refDir(dict.get<vector>("refDir").normalise());
183  const vector axis(dict.get<vector>("axis").normalise());
184 
185  R_ = tensor(refDir, axis, refDir^axis);
186  break;
187  }
188  }
189 
190  Info<< endl;
191 
192  return true;
193  }
194 
195  return false;
196 }
197 
198 
199 bool Foam::functionObjects::setFlow::execute()
200 {
201  volVectorField* Uptr =
202  mesh_.getObjectPtr<volVectorField>(UName_);
203 
204  surfaceScalarField* phiptr =
205  mesh_.getObjectPtr<surfaceScalarField>(phiName_);
206 
207  Log << nl << name() << ":" << nl;
208 
209  if (!Uptr || !phiptr)
210  {
211  Info<< " Either field " << UName_ << " or " << phiName_
212  << " not found in the mesh database" << nl;
213 
214  return true;
215  }
216 
217  const scalar t = mesh_.time().timeOutputValue();
218 
219  Log << " setting " << UName_ << " and " << phiName_ << nl;
220 
221  volVectorField& U = *Uptr;
222  surfaceScalarField& phi = *phiptr;
223 
224  switch (mode_)
225  {
226  case modeType::FUNCTION:
227  {
228  const vector Uc = velocityPtr_->value(t);
229  U == dimensionedVector("Uc", dimVelocity, Uc);
230  U.correctBoundaryConditions();
231  setPhi(U);
232 
233  break;
234  }
235  case modeType::ROTATION:
236  {
237  const volVectorField& C = mesh_.C();
238  const volVectorField d
239  (
240  typeName + ":d",
241  C - dimensionedVector("origin", dimLength, origin_)
242  );
243  const scalarField x(d.component(vector::X));
244  const scalarField z(d.component(vector::Z));
245 
246  const scalar omega = omegaPtr_->value(t);
247  vectorField& Uc = U.primitiveFieldRef();
248  Uc.replace(vector::X, -omega*z);
249  Uc.replace(vector::Y, scalar(0));
250  Uc.replace(vector::Z, omega*x);
251 
252  volVectorField::Boundary& Ubf = U.boundaryFieldRef();
253  forAll(Ubf, patchi)
254  {
255  vectorField& Uf = Ubf[patchi];
256  if (Uf.size())
257  {
258  const vectorField& Cp = C.boundaryField()[patchi];
259  const vectorField dp(Cp - origin_);
260  const scalarField xp(dp.component(vector::X));
261  const scalarField zp(dp.component(vector::Z));
262  Uf.replace(vector::X, -omega*zp);
263  Uf.replace(vector::Y, scalar(0));
264  Uf.replace(vector::Z, omega*xp);
265  }
266  }
267 
268  U = U & R_;
269  U.correctBoundaryConditions();
270  setPhi(U);
271 
272  break;
273  }
274  case modeType::VORTEX2D:
275  {
276  const scalar pi = Foam::constant::mathematical::pi;
277 
278  const volVectorField& C = mesh_.C();
279 
280  const volVectorField d
281  (
282  typeName + ":d",
283  C - dimensionedVector("origin", dimLength, origin_)
284  );
285  const scalarField x(d.component(vector::X));
286  const scalarField z(d.component(vector::Z));
287 
288  vectorField& Uc = U.primitiveFieldRef();
289  Uc.replace(vector::X, -sin(2*pi*z)*sqr(sin(pi*x)));
290  Uc.replace(vector::Y, scalar(0));
291  Uc.replace(vector::Z, sin(2*pi*x)*sqr(sin(pi*z)));
292 
293  U = U & R_;
294  U.correctBoundaryConditions();
295 
296  // Calculating phi
297  // Note: R_ rotation not implemented in phi calculation
298  const vectorField Cf(mesh_.Cf().primitiveField() - origin_);
299  const scalarField Xf(Cf.component(vector::X));
300  const scalarField Yf(Cf.component(vector::Y));
301  const scalarField Zf(Cf.component(vector::Z));
302  vectorField Uf(Xf.size());
303  Uf.replace(vector::X, -sin(2*pi*Zf)*sqr(sin(pi*Xf)));
304  Uf.replace(vector::Y, scalar(0));
305  Uf.replace(vector::Z, sin(2*pi*Xf)*sqr(sin(pi*Zf)));
306 
307  scalarField& phic = phi.primitiveFieldRef();
308  const vectorField& Sfc = mesh_.Sf().primitiveField();
309  phic = Uf & Sfc;
310 
311  surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
312  const surfaceVectorField::Boundary& Sfbf =
313  mesh_.Sf().boundaryField();
314  const surfaceVectorField::Boundary& Cfbf =
315  mesh_.Cf().boundaryField();
316 
317  forAll(phibf, patchi)
318  {
319  scalarField& phif = phibf[patchi];
320  const vectorField& Sff = Sfbf[patchi];
321  const vectorField& Cff = Cfbf[patchi];
322  const scalarField xf(Cff.component(vector::X));
323  const scalarField yf(Cff.component(vector::Y));
324  const scalarField zf(Cff.component(vector::Z));
325  vectorField Ufb(xf.size());
326  Ufb.replace(vector::X, -sin(2*pi*zf)*sqr(sin(pi*xf)));
327  Ufb.replace(vector::Y, scalar(0));
328  Ufb.replace(vector::Z, sin(2*pi*xf)*sqr(sin(pi*zf)));
329  phif = Ufb & Sff;
330  }
331 
332  break;
333  }
334  case modeType::VORTEX3D:
335  {
336  const scalar pi = Foam::constant::mathematical::pi;
337  const volVectorField& C = mesh_.C();
338 
339  const volVectorField d
340  (
341  typeName + ":d",
342  C - dimensionedVector("origin", dimLength, origin_)
343  );
344  const scalarField x(d.component(vector::X));
345  const scalarField y(d.component(vector::Y));
346  const scalarField z(d.component(vector::Z));
347 
348  vectorField& Uc = U.primitiveFieldRef();
349  Uc.replace(vector::X, 2*sqr(sin(pi*x))*sin(2*pi*y)*sin(2*pi*z));
350  Uc.replace(vector::Y, -sin(2*pi*x)*sqr(sin(pi*y))*sin(2*pi*z));
351  Uc.replace(vector::Z, -sin(2*pi*x)*sin(2*pi*y)*sqr(sin(pi*z)));
352 
353  U = U & R_;
354  U.correctBoundaryConditions();
355 
356  // Calculating phi
357  // Note: R_ rotation not implemented in phi calculation
358  const vectorField Cf(mesh_.Cf().primitiveField() - origin_);
359  const scalarField Xf(Cf.component(vector::X));
360  const scalarField Yf(Cf.component(vector::Y));
361  const scalarField Zf(Cf.component(vector::Z));
362  vectorField Uf(Xf.size());
363  Uf.replace(0, 2*sqr(sin(pi*Xf))*sin(2*pi*Yf)*sin(2*pi*Zf));
364  Uf.replace(1, -sin(2*pi*Xf)*sqr(sin(pi*Yf))*sin(2*pi*Zf));
365  Uf.replace(2, -sin(2*pi*Xf)*sin(2*pi*Yf)*sqr(sin(pi*Zf)));
366 
367  scalarField& phic = phi.primitiveFieldRef();
368  const vectorField& Sfc = mesh_.Sf().primitiveField();
369  phic = Uf & Sfc;
370 
371  surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
372  const surfaceVectorField::Boundary& Sfbf =
373  mesh_.Sf().boundaryField();
374  const surfaceVectorField::Boundary& Cfbf =
375  mesh_.Cf().boundaryField();
376 
377  forAll(phibf, patchi)
378  {
379  scalarField& phif = phibf[patchi];
380  const vectorField& Sff = Sfbf[patchi];
381  const vectorField& Cff = Cfbf[patchi];
382  const scalarField xf(Cff.component(vector::X));
383  const scalarField yf(Cff.component(vector::Y));
384  const scalarField zf(Cff.component(vector::Z));
385  vectorField Uf(xf.size());
386  Uf.replace(0, 2*sqr(sin(pi*xf))*sin(2*pi*yf)*sin(2*pi*zf));
387  Uf.replace(1, -sin(2*pi*xf)*sqr(sin(pi*yf))*sin(2*pi*zf));
388  Uf.replace(2, -sin(2*pi*xf)*sin(2*pi*yf)*sqr(sin(pi*zf)));
389  phif = Uf & Sff;
390  }
391 
392  break;
393  }
394  }
395 
396  if (t > reverseTime_)
397  {
398  Log << " flow direction: reverse" << nl;
399  U.negate();
400  phi.negate();
401  }
402 
403  // Apply scaling
404  const scalar s = scalePtr_->value(t);
405  U *= s;
406  phi *= s;
407 
408  U.correctBoundaryConditions();
409 
410  const scalarField sumPhi(fvc::surfaceIntegrate(phi));
411  Log << " Continuity error: max(mag(sum(phi))) = "
412  << gMax(mag(sumPhi)) << nl << endl;
413 
414  return true;
415 }
416 
417 
418 bool Foam::functionObjects::setFlow::write()
419 {
420  const auto* Uptr = mesh_.findObject<volVectorField>(UName_);
421  if (Uptr)
422  {
423  Uptr->write();
424  }
425 
426  const auto* phiptr = mesh_.findObject<surfaceScalarField>(phiName_);
427  if (phiptr)
428  {
429  phiptr->write();
430  }
431 
432  return true;
433 }
434 
435 
436 // ************************************************************************* //
Foam::functionObjects::setFlow::read
virtual bool read(const dictionary &dict)
Read the setFlow data.
Definition: setFlow.C:125
volFields.H
Foam::objectRegistry::getObjectPtr
Type * getObjectPtr(const word &name, const bool recursive=false) const
Definition: objectRegistryTemplates.C:423
runTime
engineTime & runTime
Definition: createEngineTime.H:13
Foam::GeometricField::component
tmp< GeometricField< cmptType, PatchField, GeoMesh > > component(const direction) const
Return a component of the field.
Foam::Enum
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: IOstreamOption.H:57
Log
#define Log
Definition: PDRblock.C:35
Foam::fvc::flux
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:73
Foam::Vector< scalar >::Z
Definition: Vector.H:81
Foam::Vector< scalar >::Y
Definition: Vector.H:81
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::Tensor::I
static const Tensor I
Definition: Tensor.H:85
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Definition: gmvOutputSpray.H:25
Foam::dimLength
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::dimVelocity
const dimensionSet dimVelocity
Foam::sin
dimensionedScalar sin(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:264
Foam::read
bool read(const char *buf, int32_t &val)
Same as readInt32.
Definition: int32.H:108
Foam::functionObjects::setFlow::setFlow
setFlow(const word &name, const Time &runTime, const dictionary &dict)
Construct from Time and dictionary.
Definition: setFlow.C:101
Uf
autoPtr< surfaceVectorField > Uf
Definition: createUfIfPresent.H:33
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
surfaceFields.H
Foam::surfaceFields.
Foam::Vector::normalise
Vector< Cmpt > & normalise()
Normalise the vector by its magnitude.
Definition: VectorI.H:123
Foam::functionObjects::fvMeshFunctionObject
Specialization of Foam::functionObject for an Foam::fvMesh, providing a reference to the Foam::fvMesh...
Definition: fvMeshFunctionObject.H:64
rho
rho
Definition: readInitialConditions.H:88
Foam::Function1
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
Definition: propellerInfo.H:291
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
phic
surfaceScalarField phic(mixture.cAlpha() *mag(alphaPhic/mesh.magSf()))
Foam::dimensionedVector
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Definition: dimensionedVector.H:50
Foam::functionObjects::setFlow::execute
virtual bool execute()
Do nothing.
Definition: setFlow.C:199
Foam::functionObjects::setFlow::write
virtual bool write()
Calculate the setFlow and write.
Definition: setFlow.C:418
Foam::Field< scalar >
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::Vector< scalar >::X
Definition: Vector.H:81
Foam::Field::replace
void replace(const direction, const UList< cmptType > &)
Replace a component field of the field.
Definition: Field.C:551
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::volScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
fvcSurfaceIntegrate.H
Surface integrate surfaceField creating a volField. Surface sum a surfaceField creating a volField.
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
Foam::functionObjects::regionFunctionObject::read
virtual bool read(const dictionary &dict)
Read optional controls.
Definition: regionFunctionObject.C:173
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
Foam::C::C
C()
Construct null.
Definition: C.C:43
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::Field::component
tmp< Field< cmptType > > component(const direction) const
Return a component field of the field.
Definition: Field.C:539
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::volVectorField
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:62
U
U
Definition: pEqn.H:72
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::objectRegistry::findObject
const Type * findObject(const word &name, const bool recursive=false) const
Return const pointer to the object of the given Type.
Definition: objectRegistryTemplates.C:401
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::nl
constexpr char nl
Definition: Ostream.H:404
setFlow.H
Foam::objectRegistry::names
wordList names() const
The names of all objects.
Definition: objectRegistry.C:147
Foam::surfaceScalarField
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Definition: surfaceFieldsFwd.H:54
Foam::functionObjects::addToRunTimeSelectionTable
addToRunTimeSelectionTable(functionObject, ObukhovLength, dictionary)
Foam::Vector< scalar >
Foam::fvc::surfaceIntegrate
void surfaceIntegrate(Field< Type > &ivf, const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcSurfaceIntegrate.C:46
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
fvcFlux.H
Calculate the face-flux of the given field.
Foam::functionObjects::fvMeshFunctionObject::mesh_
const fvMesh & mesh_
Reference to the fvMesh.
Definition: fvMeshFunctionObject.H:73
Cp
const volScalarField & Cp
Definition: EEqn.H:7
Foam::functionObjects::defineTypeNameAndDebug
defineTypeNameAndDebug(ObukhovLength, 0)
x
x
Definition: LISASMDCalcMethod2.H:52
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Foam::C
Graphite solid properties.
Definition: C.H:50
Foam::GeometricField< vector, fvPatchField, volMesh >
Foam::tensor
Tensor< scalar > tensor
Tensor of scalars, i.e. Tensor<scalar>.
Definition: symmTensor.H:61
Foam::gMax
Type gMax(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:592
y
scalar y
Definition: LISASMDCalcMethod1.H:14