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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2018 OpenFOAM Foundation
9  Copyright (C) 2023 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "EulerDdtScheme.H"
30 #include "surfaceInterpolate.H"
31 #include "fvcDiv.H"
32 #include "fvMatrices.H"
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 
39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
40 
41 namespace fv
42 {
43 
44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
45 
46 template<class Type>
49 (
50  const dimensioned<Type>& dt
51 )
52 {
53  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
54 
55  IOobject ddtIOobject
56  (
57  "ddt("+dt.name()+')',
58  mesh().time().timeName(),
59  mesh().thisDb()
60  );
61 
62  if (mesh().moving())
63  {
65  (
67  (
68  ddtIOobject,
69  mesh(),
71  )
72  );
73 
74  tdtdt.ref().primitiveFieldRef() =
75  rDeltaT.value()*dt.value()*(1.0 - mesh().Vsc0()/mesh().Vsc());
76 
77  return tdtdt;
78  }
79 
81  (
82  ddtIOobject,
83  mesh(),
86  );
87 }
88 
89 
90 template<class Type>
93 (
95 )
96 {
97  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
98 
99  IOobject ddtIOobject
100  (
101  "ddt("+vf.name()+')',
102  mesh().time().timeName(),
103  mesh().thisDb()
104  );
105 
106  if (mesh().moving())
107  {
109  (
111  (
112  ddtIOobject,
113  rDeltaT*
114  (
115  vf()
116  - vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
117  ),
118  rDeltaT.value()*
119  (
120  vf.boundaryField() - vf.oldTime().boundaryField()
121  )
122  )
123  );
124 
125  // Different operation on boundary v.s. internal so re-evaluate
126  // coupled boundaries
127  tdtdt.ref().boundaryFieldRef().
128  template evaluateCoupled<coupledFvPatch>();
129 
130  return tdtdt;
131  }
132  else
133  {
135  (
137  (
138  ddtIOobject,
139  rDeltaT*(vf - vf.oldTime())
140  )
141  );
142  }
143 }
144 
145 
146 template<class Type>
149 (
150  const dimensionedScalar& rho,
152 )
153 {
154  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
155 
156  IOobject ddtIOobject
157  (
158  "ddt("+rho.name()+','+vf.name()+')',
159  mesh().time().timeName(),
160  mesh().thisDb()
161  );
162 
163  if (mesh().moving())
164  {
166  (
168  (
169  ddtIOobject,
170  rDeltaT*rho*
171  (
172  vf()
173  - vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
174  ),
175  rDeltaT.value()*rho.value()*
176  (
177  vf.boundaryField() - vf.oldTime().boundaryField()
178  )
179  )
180  );
181 
182  // Different operation on boundary v.s. internal so re-evaluate
183  // coupled boundaries
184  tdtdt.ref().boundaryFieldRef().
185  template evaluateCoupled<coupledFvPatch>();
186 
187  return tdtdt;
188  }
189  else
190  {
191  return tmp<GeometricField<Type, fvPatchField, volMesh>>
192  (
193  new GeometricField<Type, fvPatchField, volMesh>
194  (
195  ddtIOobject,
196  rDeltaT*rho*(vf - vf.oldTime())
197  )
198  );
199  }
200 }
201 
202 
203 template<class Type>
206 (
207  const volScalarField& rho,
209 )
210 {
211  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
212 
213  IOobject ddtIOobject
214  (
215  "ddt("+rho.name()+','+vf.name()+')',
216  mesh().time().timeName(),
217  mesh().thisDb()
218  );
219 
220  if (mesh().moving())
221  {
223  (
225  (
226  ddtIOobject,
227  rDeltaT*
228  (
229  rho()*vf()
230  - rho.oldTime()()
231  *vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
232  ),
233  rDeltaT.value()*
234  (
235  rho.boundaryField()*vf.boundaryField()
236  - rho.oldTime().boundaryField()
237  *vf.oldTime().boundaryField()
238  )
239  )
240  );
241 
242  // Different operation on boundary v.s. internal so re-evaluate
243  // coupled boundaries
244  tdtdt.ref().boundaryFieldRef().
245  template evaluateCoupled<coupledFvPatch>();
246 
247  return tdtdt;
248  }
249  else
250  {
251  return tmp<GeometricField<Type, fvPatchField, volMesh>>
252  (
253  new GeometricField<Type, fvPatchField, volMesh>
254  (
255  ddtIOobject,
256  rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
257  )
258  );
259  }
260 }
261 
262 
263 template<class Type>
266 (
267  const volScalarField& alpha,
268  const volScalarField& rho,
270 )
271 {
272  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
273 
274  IOobject ddtIOobject
275  (
276  "ddt("+alpha.name()+','+rho.name()+','+vf.name()+')',
277  mesh().time().timeName(),
278  mesh().thisDb()
279  );
280 
281  if (mesh().moving())
282  {
284  (
286  (
287  ddtIOobject,
288  rDeltaT*
289  (
290  alpha()
291  *rho()
292  *vf()
293 
294  - alpha.oldTime()()
295  *rho.oldTime()()
296  *vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
297  ),
298  rDeltaT.value()*
299  (
300  alpha.boundaryField()
301  *rho.boundaryField()
302  *vf.boundaryField()
303 
304  - alpha.oldTime().boundaryField()
305  *rho.oldTime().boundaryField()
306  *vf.oldTime().boundaryField()
307  )
308  )
309  );
310 
311  // Different operation on boundary v.s. internal so re-evaluate
312  // coupled boundaries
313  tdtdt.ref().boundaryFieldRef().
314  template evaluateCoupled<coupledFvPatch>();
315 
316  return tdtdt;
317  }
318  else
319  {
320  return tmp<GeometricField<Type, fvPatchField, volMesh>>
321  (
322  new GeometricField<Type, fvPatchField, volMesh>
323  (
324  ddtIOobject,
325  rDeltaT
326  *(
327  alpha*rho*vf
328  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
329  )
330  )
331  );
332  }
333 }
334 
335 
336 template<class Type>
339 (
341 )
342 {
343  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
344 
345  IOobject ddtIOobject
346  (
347  "ddt("+sf.name()+')',
348  mesh().time().timeName(),
349  mesh().thisDb()
350  );
351 
353  (
355  (
356  ddtIOobject,
357  rDeltaT*(sf - sf.oldTime())
358  )
359  );
360 }
361 
362 
363 template<class Type>
366 (
368 )
369 {
370  tmp<fvMatrix<Type>> tfvm
371  (
372  new fvMatrix<Type>
373  (
374  vf,
376  )
377  );
378 
379  fvMatrix<Type>& fvm = tfvm.ref();
380 
381  scalar rDeltaT = 1.0/mesh().time().deltaTValue();
382 
383  fvm.diag() = rDeltaT*mesh().Vsc();
384 
385  if (mesh().moving())
386  {
387  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc0();
388  }
389  else
390  {
391  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc();
392  }
393 
394  return tfvm;
395 }
396 
397 
398 template<class Type>
401 (
402  const dimensionedScalar& rho,
404 )
405 {
406  tmp<fvMatrix<Type>> tfvm
407  (
408  new fvMatrix<Type>
409  (
410  vf,
411  rho.dimensions()*vf.dimensions()*dimVol/dimTime
412  )
413  );
414  fvMatrix<Type>& fvm = tfvm.ref();
415 
416  scalar rDeltaT = 1.0/mesh().time().deltaTValue();
417 
418  fvm.diag() = rDeltaT*rho.value()*mesh().Vsc();
419 
420  if (mesh().moving())
421  {
422  fvm.source() = rDeltaT
423  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc0();
424  }
425  else
426  {
427  fvm.source() = rDeltaT
428  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc();
429  }
430 
431  return tfvm;
432 }
433 
434 
435 template<class Type>
438 (
439  const volScalarField& rho,
441 )
442 {
443  tmp<fvMatrix<Type>> tfvm
444  (
445  new fvMatrix<Type>
446  (
447  vf,
448  rho.dimensions()*vf.dimensions()*dimVol/dimTime
449  )
450  );
451  fvMatrix<Type>& fvm = tfvm.ref();
452 
453  scalar rDeltaT = 1.0/mesh().time().deltaTValue();
454 
455  fvm.diag() = rDeltaT*rho.primitiveField()*mesh().Vsc();
456 
457  if (mesh().moving())
458  {
459  fvm.source() = rDeltaT
460  *rho.oldTime().primitiveField()
461  *vf.oldTime().primitiveField()*mesh().Vsc0();
462  }
463  else
464  {
465  fvm.source() = rDeltaT
466  *rho.oldTime().primitiveField()
467  *vf.oldTime().primitiveField()*mesh().Vsc();
468  }
469 
470  return tfvm;
471 }
472 
473 
474 template<class Type>
477 (
478  const volScalarField& alpha,
479  const volScalarField& rho,
481 )
482 {
483  tmp<fvMatrix<Type>> tfvm
484  (
485  new fvMatrix<Type>
486  (
487  vf,
488  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
489  )
490  );
491  fvMatrix<Type>& fvm = tfvm.ref();
492 
493  scalar rDeltaT = 1.0/mesh().time().deltaTValue();
494 
495  fvm.diag() =
496  rDeltaT*alpha.primitiveField()*rho.primitiveField()*mesh().Vsc();
497 
498  if (mesh().moving())
499  {
500  fvm.source() = rDeltaT
501  *alpha.oldTime().primitiveField()
502  *rho.oldTime().primitiveField()
503  *vf.oldTime().primitiveField()*mesh().Vsc0();
504  }
505  else
506  {
507  fvm.source() = rDeltaT
508  *alpha.oldTime().primitiveField()
509  *rho.oldTime().primitiveField()
510  *vf.oldTime().primitiveField()*mesh().Vsc();
511  }
512 
513  return tfvm;
514 }
515 
516 
517 template<class Type>
520 (
523 )
524 {
525  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
526 
527  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
528  fluxFieldType phiCorr
529  (
530  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
531  );
532 
533  return tmp<fluxFieldType>
534  (
535  new fluxFieldType
536  (
537  IOobject
538  (
539  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
540  mesh().time().timeName(),
541  mesh().thisDb()
542  ),
543  this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
544  *rDeltaT*phiCorr
545  )
546  );
547 }
548 
549 
550 template<class Type>
553 (
555  const fluxFieldType& phi
556 )
557 {
558  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
559 
560  fluxFieldType phiCorr
561  (
562  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
563  );
564 
565  return tmp<fluxFieldType>
566  (
567  new fluxFieldType
568  (
569  IOobject
570  (
571  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
572  mesh().time().timeName(),
573  mesh().thisDb()
574  ),
575  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
576  *rDeltaT*phiCorr
577  )
578  );
579 }
580 
581 
582 template<class Type>
585 (
586  const volScalarField& rho,
589 )
590 {
591  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
592 
593  if
594  (
595  U.dimensions() == dimVelocity
596  && Uf.dimensions() == rho.dimensions()*dimVelocity
597  )
598  {
600  (
601  rho.oldTime()*U.oldTime()
602  );
603 
604  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
605  fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
606 
607  return tmp<fluxFieldType>
608  (
609  new fluxFieldType
610  (
611  IOobject
612  (
613  "ddtCorr("
614  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
615  mesh().time().timeName(),
616  mesh().thisDb()
617  ),
618  this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime())
619  *rDeltaT*phiCorr
620  )
621  );
622  }
623  else if
624  (
625  U.dimensions() == rho.dimensions()*dimVelocity
626  && Uf.dimensions() == rho.dimensions()*dimVelocity
627  )
628  {
629  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
630  fluxFieldType phiCorr
631  (
632  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
633  );
634 
635  return tmp<fluxFieldType>
636  (
637  new fluxFieldType
638  (
639  IOobject
640  (
641  "ddtCorr("
642  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
643  mesh().time().timeName(),
644  mesh().thisDb()
645  ),
646  this->fvcDdtPhiCoeff
647  (
648  U.oldTime(),
649  phiUf0,
650  phiCorr,
651  rho.oldTime()
652  )*rDeltaT*phiCorr
653  )
654  );
655  }
656  else
657  {
659  << "dimensions of Uf are not correct"
660  << abort(FatalError);
661 
662  return fluxFieldType::null();
663  }
664 }
665 
666 
667 template<class Type>
670 (
671  const volScalarField& rho,
673  const fluxFieldType& phi
674 )
675 {
676  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
677 
678  if
679  (
680  U.dimensions() == dimVelocity
681  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
682  )
683  {
685  (
686  rho.oldTime()*U.oldTime()
687  );
688 
689  fluxFieldType phiCorr
690  (
691  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
692  );
693 
694  return tmp<fluxFieldType>
695  (
696  new fluxFieldType
697  (
698  IOobject
699  (
700  "ddtCorr("
701  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
702  mesh().time().timeName(),
703  mesh().thisDb()
704  ),
705  this->fvcDdtPhiCoeff
706  (
707  rhoU0,
708  phi.oldTime(),
709  phiCorr,
710  rho.oldTime()
711  )*rDeltaT*phiCorr
712  )
713  );
714  }
715  else if
716  (
717  U.dimensions() == rho.dimensions()*dimVelocity
718  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
719  )
720  {
721  fluxFieldType phiCorr
722  (
723  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
724  );
725 
726  return tmp<fluxFieldType>
727  (
728  new fluxFieldType
729  (
730  IOobject
731  (
732  "ddtCorr("
733  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
734  mesh().time().timeName(),
735  mesh().thisDb()
736  ),
737  this->fvcDdtPhiCoeff
738  (
739  U.oldTime(),
740  phi.oldTime(),
741  phiCorr,
742  rho.oldTime()
743  )*rDeltaT*phiCorr
744  )
745  );
746  }
747  else
748  {
750  << "dimensions of phi are not correct"
751  << abort(FatalError);
752 
753  return fluxFieldType::null();
754  }
755 }
756 
757 
758 template<class Type>
760 (
762 )
763 {
764  return mesh().phi();
765 }
766 
767 
768 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
769 
770 } // End namespace fv
771 
772 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
773 
774 } // End namespace Foam
775 
776 // ************************************************************************* //
const Type & value() const noexcept
Return const reference to value.
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
error FatalError
Error stream (stdout output on all processes), with additional &#39;FOAM FATAL ERROR&#39; header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:598
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:195
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
Definition: tmpI.H:235
static tmp< GeometricField< typename innerProduct< vector, Type >::type, fvsPatchField, surfaceMesh > > dotInterpolate(const surfaceVectorField &Sf, const GeometricField< Type, fvPatchField, volMesh > &tvf)
Interpolate field onto faces.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
const dimensionSet dimVol(dimVolume)
Older spelling for dimVolume.
Definition: dimensionSets.H:65
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
Generic GeometricField class.
Definition: areaFieldsFwd.H:50
Generic dimensioned Type class.
dynamicFvMesh & mesh
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
labelList fv(nPoints)
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:64
autoPtr< surfaceVectorField > Uf
const dimensionSet & dimensions() const noexcept
Return const reference to dimensions.
errorManip< error > abort(error &err)
Definition: errorManip.H:139
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Calculate the divergence of the given field.
const word & name() const noexcept
Return const reference to name.
Field< Type > & source() noexcept
Definition: fvMatrix.H:533
U
Definition: pEqn.H:72
static const word & calculatedType() noexcept
The type name for calculated patch fields.
Definition: fvPatchField.H:204
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
A special matrix type and solver, designed for finite volume solutions of scalar equations.
A class for managing temporary objects.
Definition: HashPtrTable.H:50
scalarField & diag()
Definition: lduMatrix.C:197
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Defines the attributes of an object for which implicit objectRegistry management is supported...
Definition: IOobject.H:172
const Boundary & boundaryField() const noexcept
Return const-reference to the boundary field.
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Namespace for OpenFOAM.
const dimensionSet & dimensions() const noexcept
Return dimensions.
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:127
const dimensionSet dimVelocity