EulerDdtScheme.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 | Copyright (C) 2011-2016 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 "EulerDdtScheme.H"
27 #include "surfaceInterpolate.H"
28 #include "fvcDiv.H"
29 #include "fvMatrices.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 namespace fv
39 {
40 
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 
43 template<class Type>
44 tmp<GeometricField<Type, fvPatchField, volMesh>>
45 EulerDdtScheme<Type>::fvcDdt
46 (
47  const dimensioned<Type>& dt
48 )
49 {
50  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
51 
52  IOobject ddtIOobject
53  (
54  "ddt("+dt.name()+')',
55  mesh().time().timeName(),
56  mesh()
57  );
58 
59  if (mesh().moving())
60  {
61  tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
62  (
63  new GeometricField<Type, fvPatchField, volMesh>
64  (
65  ddtIOobject,
66  mesh(),
67  dimensioned<Type>
68  (
69  "0",
70  dt.dimensions()/dimTime,
71  Zero
72  )
73  )
74  );
75 
76  tdtdt.ref().primitiveFieldRef() =
77  rDeltaT.value()*dt.value()*(1.0 - mesh().Vsc0()/mesh().Vsc());
78 
79  return tdtdt;
80  }
81  else
82  {
83  return tmp<GeometricField<Type, fvPatchField, volMesh>>
84  (
85  new GeometricField<Type, fvPatchField, volMesh>
86  (
87  ddtIOobject,
88  mesh(),
89  dimensioned<Type>
90  (
91  "0",
92  dt.dimensions()/dimTime,
93  Zero
94  ),
95  calculatedFvPatchField<Type>::typeName
96  )
97  );
98  }
99 }
100 
101 
102 template<class Type>
103 tmp<GeometricField<Type, fvPatchField, volMesh>>
104 EulerDdtScheme<Type>::fvcDdt
105 (
106  const GeometricField<Type, fvPatchField, volMesh>& vf
107 )
108 {
109  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
110 
111  IOobject ddtIOobject
112  (
113  "ddt("+vf.name()+')',
114  mesh().time().timeName(),
115  mesh()
116  );
117 
118  if (mesh().moving())
119  {
120  return tmp<GeometricField<Type, fvPatchField, volMesh>>
121  (
122  new GeometricField<Type, fvPatchField, volMesh>
123  (
124  ddtIOobject,
125  rDeltaT*
126  (
127  vf()
128  - vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
129  ),
130  rDeltaT.value()*
131  (
132  vf.boundaryField() - vf.oldTime().boundaryField()
133  )
134  )
135  );
136  }
137  else
138  {
139  return tmp<GeometricField<Type, fvPatchField, volMesh>>
140  (
141  new GeometricField<Type, fvPatchField, volMesh>
142  (
143  ddtIOobject,
144  rDeltaT*(vf - vf.oldTime())
145  )
146  );
147  }
148 }
149 
150 
151 template<class Type>
152 tmp<GeometricField<Type, fvPatchField, volMesh>>
153 EulerDdtScheme<Type>::fvcDdt
154 (
155  const dimensionedScalar& rho,
156  const GeometricField<Type, fvPatchField, volMesh>& vf
157 )
158 {
159  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
160 
161  IOobject ddtIOobject
162  (
163  "ddt("+rho.name()+','+vf.name()+')',
164  mesh().time().timeName(),
165  mesh()
166  );
167 
168  if (mesh().moving())
169  {
170  return tmp<GeometricField<Type, fvPatchField, volMesh>>
171  (
172  new GeometricField<Type, fvPatchField, volMesh>
173  (
174  ddtIOobject,
175  rDeltaT*rho*
176  (
177  vf()
178  - vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
179  ),
180  rDeltaT.value()*rho.value()*
181  (
182  vf.boundaryField() - vf.oldTime().boundaryField()
183  )
184  )
185  );
186  }
187  else
188  {
189  return tmp<GeometricField<Type, fvPatchField, volMesh>>
190  (
191  new GeometricField<Type, fvPatchField, volMesh>
192  (
193  ddtIOobject,
194  rDeltaT*rho*(vf - vf.oldTime())
195  )
196  );
197  }
198 }
199 
200 
201 template<class Type>
202 tmp<GeometricField<Type, fvPatchField, volMesh>>
203 EulerDdtScheme<Type>::fvcDdt
204 (
205  const volScalarField& rho,
206  const GeometricField<Type, fvPatchField, volMesh>& vf
207 )
208 {
209  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
210 
211  IOobject ddtIOobject
212  (
213  "ddt("+rho.name()+','+vf.name()+')',
214  mesh().time().timeName(),
215  mesh()
216  );
217 
218  if (mesh().moving())
219  {
220  return tmp<GeometricField<Type, fvPatchField, volMesh>>
221  (
222  new GeometricField<Type, fvPatchField, volMesh>
223  (
224  ddtIOobject,
225  rDeltaT*
226  (
227  rho()*vf()
228  - rho.oldTime()()
229  *vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
230  ),
231  rDeltaT.value()*
232  (
233  rho.boundaryField()*vf.boundaryField()
234  - rho.oldTime().boundaryField()
235  *vf.oldTime().boundaryField()
236  )
237  )
238  );
239  }
240  else
241  {
242  return tmp<GeometricField<Type, fvPatchField, volMesh>>
243  (
244  new GeometricField<Type, fvPatchField, volMesh>
245  (
246  ddtIOobject,
247  rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
248  )
249  );
250  }
251 }
252 
253 
254 template<class Type>
255 tmp<GeometricField<Type, fvPatchField, volMesh>>
256 EulerDdtScheme<Type>::fvcDdt
257 (
258  const volScalarField& alpha,
259  const volScalarField& rho,
260  const GeometricField<Type, fvPatchField, volMesh>& vf
261 )
262 {
263  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
264 
265  IOobject ddtIOobject
266  (
267  "ddt("+alpha.name()+','+rho.name()+','+vf.name()+')',
268  mesh().time().timeName(),
269  mesh()
270  );
271 
272  if (mesh().moving())
273  {
274  return tmp<GeometricField<Type, fvPatchField, volMesh>>
275  (
276  new GeometricField<Type, fvPatchField, volMesh>
277  (
278  ddtIOobject,
279  rDeltaT.value()*
280  (
281  alpha()
282  *rho()
283  *vf()
284 
285  - alpha.oldTime()()
286  *rho.oldTime()()
287  *vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
288  ),
289  rDeltaT.value()*
290  (
291  alpha.boundaryField()
292  *rho.boundaryField()
293  *vf.boundaryField()
294 
295  - alpha.oldTime().boundaryField()
296  *rho.oldTime().boundaryField()
297  *vf.oldTime().boundaryField()
298  )
299  )
300  );
301  }
302  else
303  {
304  return tmp<GeometricField<Type, fvPatchField, volMesh>>
305  (
306  new GeometricField<Type, fvPatchField, volMesh>
307  (
308  ddtIOobject,
309  rDeltaT
310  *(
311  alpha*rho*vf
312  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
313  )
314  )
315  );
316  }
317 }
318 
319 
320 template<class Type>
321 tmp<fvMatrix<Type>>
322 EulerDdtScheme<Type>::fvmDdt
323 (
324  const GeometricField<Type, fvPatchField, volMesh>& vf
325 )
326 {
327  tmp<fvMatrix<Type>> tfvm
328  (
329  new fvMatrix<Type>
330  (
331  vf,
332  vf.dimensions()*dimVol/dimTime
333  )
334  );
335 
336  fvMatrix<Type>& fvm = tfvm.ref();
337 
338  scalar rDeltaT = 1.0/mesh().time().deltaTValue();
339 
340  fvm.diag() = rDeltaT*mesh().Vsc();
341 
342  if (mesh().moving())
343  {
344  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc0();
345  }
346  else
347  {
348  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc();
349  }
350 
351  return tfvm;
352 }
353 
354 
355 template<class Type>
356 tmp<fvMatrix<Type>>
357 EulerDdtScheme<Type>::fvmDdt
358 (
359  const dimensionedScalar& rho,
360  const GeometricField<Type, fvPatchField, volMesh>& vf
361 )
362 {
363  tmp<fvMatrix<Type>> tfvm
364  (
365  new fvMatrix<Type>
366  (
367  vf,
368  rho.dimensions()*vf.dimensions()*dimVol/dimTime
369  )
370  );
371  fvMatrix<Type>& fvm = tfvm.ref();
372 
373  scalar rDeltaT = 1.0/mesh().time().deltaTValue();
374 
375  fvm.diag() = rDeltaT*rho.value()*mesh().Vsc();
376 
377  if (mesh().moving())
378  {
379  fvm.source() = rDeltaT
380  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc0();
381  }
382  else
383  {
384  fvm.source() = rDeltaT
385  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc();
386  }
387 
388  return tfvm;
389 }
390 
391 
392 template<class Type>
393 tmp<fvMatrix<Type>>
394 EulerDdtScheme<Type>::fvmDdt
395 (
396  const volScalarField& rho,
397  const GeometricField<Type, fvPatchField, volMesh>& vf
398 )
399 {
400  tmp<fvMatrix<Type>> tfvm
401  (
402  new fvMatrix<Type>
403  (
404  vf,
405  rho.dimensions()*vf.dimensions()*dimVol/dimTime
406  )
407  );
408  fvMatrix<Type>& fvm = tfvm.ref();
409 
410  scalar rDeltaT = 1.0/mesh().time().deltaTValue();
411 
412  fvm.diag() = rDeltaT*rho.primitiveField()*mesh().Vsc();
413 
414  if (mesh().moving())
415  {
416  fvm.source() = rDeltaT
417  *rho.oldTime().primitiveField()
418  *vf.oldTime().primitiveField()*mesh().Vsc0();
419  }
420  else
421  {
422  fvm.source() = rDeltaT
423  *rho.oldTime().primitiveField()
424  *vf.oldTime().primitiveField()*mesh().Vsc();
425  }
426 
427  return tfvm;
428 }
429 
430 
431 template<class Type>
432 tmp<fvMatrix<Type>>
433 EulerDdtScheme<Type>::fvmDdt
434 (
435  const volScalarField& alpha,
436  const volScalarField& rho,
437  const GeometricField<Type, fvPatchField, volMesh>& vf
438 )
439 {
440  tmp<fvMatrix<Type>> tfvm
441  (
442  new fvMatrix<Type>
443  (
444  vf,
445  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
446  )
447  );
448  fvMatrix<Type>& fvm = tfvm.ref();
449 
450  scalar rDeltaT = 1.0/mesh().time().deltaTValue();
451 
452  fvm.diag() =
453  rDeltaT*alpha.primitiveField()*rho.primitiveField()*mesh().Vsc();
454 
455  if (mesh().moving())
456  {
457  fvm.source() = rDeltaT
458  *alpha.oldTime().primitiveField()
459  *rho.oldTime().primitiveField()
460  *vf.oldTime().primitiveField()*mesh().Vsc0();
461  }
462  else
463  {
464  fvm.source() = rDeltaT
465  *alpha.oldTime().primitiveField()
466  *rho.oldTime().primitiveField()
467  *vf.oldTime().primitiveField()*mesh().Vsc();
468  }
469 
470  return tfvm;
471 }
472 
473 
474 template<class Type>
475 tmp<typename EulerDdtScheme<Type>::fluxFieldType>
476 EulerDdtScheme<Type>::fvcDdtUfCorr
477 (
478  const GeometricField<Type, fvPatchField, volMesh>& U,
479  const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
480 )
481 {
482  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
483 
484  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
485  fluxFieldType phiCorr
486  (
487  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
488  );
489 
490  return tmp<fluxFieldType>
491  (
492  new fluxFieldType
493  (
494  IOobject
495  (
496  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
497  mesh().time().timeName(),
498  mesh()
499  ),
500  this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
501  *rDeltaT*phiCorr
502  )
503  );
504 }
505 
506 
507 template<class Type>
508 tmp<typename EulerDdtScheme<Type>::fluxFieldType>
509 EulerDdtScheme<Type>::fvcDdtPhiCorr
510 (
511  const GeometricField<Type, fvPatchField, volMesh>& U,
512  const fluxFieldType& phi
513 )
514 {
515  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
516 
517  fluxFieldType phiCorr
518  (
519  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
520  );
521 
522  return tmp<fluxFieldType>
523  (
524  new fluxFieldType
525  (
526  IOobject
527  (
528  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
529  mesh().time().timeName(),
530  mesh()
531  ),
532  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
533  *rDeltaT*phiCorr
534  )
535  );
536 }
537 
538 
539 template<class Type>
540 tmp<typename EulerDdtScheme<Type>::fluxFieldType>
541 EulerDdtScheme<Type>::fvcDdtUfCorr
542 (
543  const volScalarField& rho,
544  const GeometricField<Type, fvPatchField, volMesh>& U,
545  const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
546 )
547 {
548  if
549  (
550  U.dimensions() == dimVelocity
551  && Uf.dimensions() == rho.dimensions()*dimVelocity
552  )
553  {
554  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
555 
556  GeometricField<Type, fvPatchField, volMesh> rhoU0
557  (
558  rho.oldTime()*U.oldTime()
559  );
560 
561  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
562  fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
563 
564  return tmp<fluxFieldType>
565  (
566  new fluxFieldType
567  (
568  IOobject
569  (
570  "ddtCorr("
571  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
572  mesh().time().timeName(),
573  mesh()
574  ),
575  this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr)*rDeltaT*phiCorr
576  )
577  );
578  }
579  else if
580  (
581  U.dimensions() == rho.dimensions()*dimVelocity
582  && Uf.dimensions() == rho.dimensions()*dimVelocity
583  )
584  {
585  return fvcDdtUfCorr(U, Uf);
586  }
587  else
588  {
589  FatalErrorInFunction
590  << "dimensions of Uf are not correct"
591  << abort(FatalError);
592 
593  return fluxFieldType::null();
594  }
595 }
596 
597 
598 template<class Type>
599 tmp<typename EulerDdtScheme<Type>::fluxFieldType>
600 EulerDdtScheme<Type>::fvcDdtPhiCorr
601 (
602  const volScalarField& rho,
603  const GeometricField<Type, fvPatchField, volMesh>& U,
604  const fluxFieldType& phi
605 )
606 {
607  if
608  (
609  U.dimensions() == dimVelocity
610  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
611  )
612  {
613  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
614 
615  GeometricField<Type, fvPatchField, volMesh> rhoU0
616  (
617  rho.oldTime()*U.oldTime()
618  );
619 
620  fluxFieldType phiCorr
621  (
622  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
623  );
624 
625  return tmp<fluxFieldType>
626  (
627  new fluxFieldType
628  (
629  IOobject
630  (
631  "ddtCorr("
632  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
633  mesh().time().timeName(),
634  mesh()
635  ),
636  this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr)
637  *rDeltaT*phiCorr
638  )
639  );
640  }
641  else if
642  (
643  U.dimensions() == rho.dimensions()*dimVelocity
644  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
645  )
646  {
647  return fvcDdtPhiCorr(U, phi);
648  }
649  else
650  {
651  FatalErrorInFunction
652  << "dimensions of phi are not correct"
653  << abort(FatalError);
654 
655  return fluxFieldType::null();
656  }
657 }
658 
659 
660 template<class Type>
661 tmp<surfaceScalarField> EulerDdtScheme<Type>::meshPhi
662 (
663  const GeometricField<Type, fvPatchField, volMesh>&
664 )
665 {
666  return mesh().phi();
667 }
668 
669 
670 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
671 
672 } // End namespace fv
673 
674 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
675 
676 } // End namespace Foam
677 
678 // ************************************************************************* //
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const
Return const reference to dimensions.
Definition: dimensionedType.C:248
Foam::FatalError
error FatalError
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Foam::fvc::dotInterpolate
static tmp< GeometricField< typename innerProduct< vector, Type >::type, fvsPatchField, surfaceMesh > > dotInterpolate(const surfaceVectorField &Sf, const GeometricField< Type, fvPatchField, volMesh > &tvf)
Interpolate field onto faces.
Foam::dimVol
const dimensionSet dimVol(dimVolume)
Definition: dimensionSets.H:59
Foam::fv::EulerDdtScheme::fvcDdt
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
Definition: EulerDdtScheme.C:46
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Definition: GeometricFieldI.H:51
Foam::dimensioned::value
const Type & value() const
Return const reference to value.
Definition: dimensionedType.C:261
Foam::GeometricField
Generic GeometricField class.
Definition: surfaceFieldsFwd.H:52
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedScalarFwd.H:41
Foam::dimensioned::name
const word & name() const
Return const reference to name.
Definition: dimensionedType.C:235
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:809
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:18
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:60
Foam::fv::EulerDdtScheme::fvcDdtPhiCorr
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
Definition: EulerDdtScheme.C:510
Foam::fv::EulerDdtScheme::meshPhi
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: EulerDdtScheme.C:662
fv
labelList fv(nPoints)
Foam::calculatedFvPatchField
This boundary condition is not designed to be evaluated; it is assmued that the value is assigned via...
Definition: calculatedFvPatchField.H:64
Foam::fvMatrix
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:71
Foam::Zero
static const zero Zero
Definition: zero.H:91
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::fv::EulerDdtScheme::fvcDdtUfCorr
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition: EulerDdtScheme.C:477
Foam::DimensionedField::dimensions
const dimensionSet & dimensions() const
Return dimensions.
Definition: DimensionedFieldI.H:46
fvcDiv.H
Calculate the divergence of the given field.
Foam::fvMatrix::source
Field< Type > & source()
Definition: fvMatrix.H:291
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
fvMatrices.H
A special matrix type and solver, designed for finite volume solutions of scalar equations.
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:54
Foam::lduMatrix::diag
scalarField & diag()
Definition: lduMatrix.C:183
Foam::tmp::ref
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
Foam::constant::atomic::alpha
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Definition: readThermalProperties.H:213
Foam::IOobject
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
Foam::IOobject::name
const word & name() const
Return name.
Definition: IOobject.H:260
Foam::fv::EulerDdtScheme::fvmDdt
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: EulerDdtScheme.C:323
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Foam
Namespace for OpenFOAM.
Definition: combustionModel.C:30
Foam::dimVelocity
const dimensionSet dimVelocity
Generated by   doxygen 1.8.11