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);
159 
160  switch (mode_)
161  {
162  case modeType::FUNCTION:
163  {
164  velocityPtr_ = Function1<vector>::New("velocity", dict);
165  break;
166  }
167  case modeType::ROTATION:
168  {
169  omegaPtr_ = Function1<scalar>::New("omega", dict);
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 
295  // Calculating incompressible flux based on stream function
296  // Note: R_ rotation not implemented in phi calculation
297  const scalarField xp(mesh_.points().component(0) - origin_[0]);
298  const scalarField yp(mesh_.points().component(1) - origin_[1]);
299  const scalarField zp(mesh_.points().component(2) - origin_[2]);
300  const scalarField psi((1.0/pi)*sqr(sin(pi*xp))*sqr(sin(pi*zp)));
301 
302  scalarField& phic = phi.primitiveFieldRef();
303  forAll(phic, fi)
304  {
305  phic[fi] = 0;
306  const face& f = mesh_.faces()[fi];
307  const label nPoints = f.size();
308 
309  forAll(f, fpi)
310  {
311  const label p1 = f[fpi];
312  const label p2 = f[(fpi + 1) % nPoints];
313  phic[fi] += 0.5*(psi[p1] + psi[p2])*(yp[p2] - yp[p1]);
314  }
315  }
316 
317  surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
318  forAll(phibf, patchi)
319  {
320  scalarField& phif = phibf[patchi];
321  const label start = mesh_.boundaryMesh()[patchi].start();
322 
323  forAll(phif, fi)
324  {
325  phif[fi] = 0;
326  const face& f = mesh_.faces()[start + fi];
327  const label nPoints = f.size();
328 
329  forAll(f, fpi)
330  {
331  const label p1 = f[fpi];
332  const label p2 = f[(fpi + 1) % nPoints];
333  phif[fi] += 0.5*(psi[p1] + psi[p2])*(yp[p2] - yp[p1]);
334  }
335  }
336  }
337 
338  break;
339  }
340  case modeType::VORTEX3D:
341  {
342  const scalar pi = Foam::constant::mathematical::pi;
343  const volVectorField& C = mesh_.C();
344 
345  const volVectorField d
346  (
347  typeName + ":d",
348  C - dimensionedVector("origin", dimLength, origin_)
349  );
350  const scalarField x(d.component(vector::X));
351  const scalarField y(d.component(vector::Y));
352  const scalarField z(d.component(vector::Z));
353 
354  vectorField& Uc = U.primitiveFieldRef();
355  Uc.replace(vector::X, 2*sqr(sin(pi*x))*sin(2*pi*y)*sin(2*pi*z));
356  Uc.replace(vector::Y, -sin(2*pi*x)*sqr(sin(pi*y))*sin(2*pi*z));
357  Uc.replace(vector::Z, -sin(2*pi*x)*sin(2*pi*y)*sqr(sin(pi*z)));
358 
359  U = U & R_;
360  U.correctBoundaryConditions();
361 
362  // Calculating phi
363  // Note: R_ rotation not implemented in phi calculation
364  const vectorField Cf(mesh_.Cf().primitiveField() - origin_);
365  const scalarField Xf(Cf.component(vector::X));
366  const scalarField Yf(Cf.component(vector::Y));
367  const scalarField Zf(Cf.component(vector::Z));
368  vectorField Uf(Xf.size());
369  Uf.replace(0, 2*sqr(sin(pi*Xf))*sin(2*pi*Yf)*sin(2*pi*Zf));
370  Uf.replace(1, -sin(2*pi*Xf)*sqr(sin(pi*Yf))*sin(2*pi*Zf));
371  Uf.replace(2, -sin(2*pi*Xf)*sin(2*pi*Yf)*sqr(sin(pi*Zf)));
372 
373  scalarField& phic = phi.primitiveFieldRef();
374  const vectorField& Sfc = mesh_.Sf().primitiveField();
375  phic = Uf & Sfc;
376 
377  surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
378  const surfaceVectorField::Boundary& Sfbf =
379  mesh_.Sf().boundaryField();
380  const surfaceVectorField::Boundary& Cfbf =
381  mesh_.Cf().boundaryField();
382 
383  forAll(phibf, patchi)
384  {
385  scalarField& phif = phibf[patchi];
386  const vectorField& Sff = Sfbf[patchi];
387  const vectorField& Cff = Cfbf[patchi];
388  const scalarField xf(Cff.component(vector::X));
389  const scalarField yf(Cff.component(vector::Y));
390  const scalarField zf(Cff.component(vector::Z));
391  vectorField Uf(xf.size());
392  Uf.replace(0, 2*sqr(sin(pi*xf))*sin(2*pi*yf)*sin(2*pi*zf));
393  Uf.replace(1, -sin(2*pi*xf)*sqr(sin(pi*yf))*sin(2*pi*zf));
394  Uf.replace(2, -sin(2*pi*xf)*sin(2*pi*yf)*sqr(sin(pi*zf)));
395  phif = Uf & Sff;
396  }
397 
398  break;
399  }
400  }
401 
402  if (t > reverseTime_)
403  {
404  Log << " flow direction: reverse" << nl;
405  U.negate();
406  phi.negate();
407  }
408 
409  // Apply scaling
410  const scalar s = scalePtr_->value(t);
411  U *= s;
412  phi *= s;
413 
414  U.correctBoundaryConditions();
415 
416  const scalarField sumPhi(fvc::surfaceIntegrate(phi));
417  Log << " Continuity error: max(mag(sum(phi))) = "
418  << gMax(mag(sumPhi)) << nl << endl;
419 
420  return true;
421 }
422 
423 
424 bool Foam::functionObjects::setFlow::write()
425 {
426  const auto* Uptr = mesh_.findObject<volVectorField>(UName_);
427  if (Uptr)
428  {
429  Uptr->write();
430  }
431 
432  const auto* phiptr = mesh_.findObject<surfaceScalarField>(phiName_);
433  if (phiptr)
434  {
435  phiptr->write();
436  }
437 
438  return true;
439 }
440 
441 
442 // ************************************************************************* //
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::Z
Definition: Vector.H:81
Foam::Vector::Y
Definition: Vector.H:81
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:62
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:53
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:350
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: Function1.H:86
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
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
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:424
Foam::Field< scalar >
Foam::Info
messageStream Info
Information stream (uses stdout - output is on the master only)
Foam::name
word name(const complex &c)
Return string representation of complex.
Definition: complex.C:76
Foam::Vector::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:563
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:121
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:551
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:381
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:385
setFlow.H
Foam::objectRegistry::names
wordList names() const
The names of all objects.
Definition: objectRegistry.C:164
f
labelList f(nPoints)
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::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:72
Foam::C
Graphite solid properties.
Definition: C.H:50
Foam::GeometricField< vector, fvPatchField, volMesh >
psi
const volScalarField & psi
Definition: createFieldRefs.H:1
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