backwardDdtScheme.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) 2017,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 "backwardDdtScheme.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>
47 scalar backwardDdtScheme<Type>::deltaT_() const
48 {
49  return mesh().time().deltaTValue();
50 }
51 
52 
53 template<class Type>
54 scalar backwardDdtScheme<Type>::deltaT0_() const
55 {
56  return mesh().time().deltaT0Value();
57 }
58 
59 
60 template<class Type>
61 template<class GeoField>
62 scalar backwardDdtScheme<Type>::deltaT0_(const GeoField& vf) const
63 {
64  if (mesh().time().timeIndex() < 2)
65  {
66  return GREAT;
67  }
68  else
69  {
70  return deltaT0_();
71  }
72 }
73 
74 
75 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
76 
77 template<class Type>
78 tmp<GeometricField<Type, fvPatchField, volMesh>>
79 backwardDdtScheme<Type>::fvcDdt
80 (
81  const dimensioned<Type>& dt
82 )
83 {
84  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
85 
86  IOobject ddtIOobject
87  (
88  "ddt("+dt.name()+')',
89  mesh().time().timeName(),
90  mesh().thisDb()
91  );
92 
93  scalar deltaT = deltaT_();
94  scalar deltaT0 = deltaT0_();
95 
96  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
97  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
98  scalar coefft0 = coefft + coefft00;
99 
100  if (mesh().moving())
101  {
102  tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
103  (
104  new GeometricField<Type, fvPatchField, volMesh>
105  (
106  ddtIOobject,
107  mesh(),
108  dimensioned<Type>(dt.dimensions()/dimTime, Zero)
109  )
110  );
111 
112  tdtdt.ref().primitiveFieldRef() = rDeltaT.value()*dt.value()*
113  (
114  coefft - (coefft0*mesh().V0() - coefft00*mesh().V00())/mesh().V()
115  );
116 
117  return tdtdt;
118  }
119 
120  return tmp<GeometricField<Type, fvPatchField, volMesh>>::New
121  (
122  ddtIOobject,
123  mesh(),
124  dimensioned<Type>(dt.dimensions()/dimTime, Zero),
125  fvPatchFieldBase::calculatedType()
126  );
127 }
128 
129 
130 template<class Type>
131 tmp<GeometricField<Type, fvPatchField, volMesh>>
132 backwardDdtScheme<Type>::fvcDdt
133 (
134  const GeometricField<Type, fvPatchField, volMesh>& vf
135 )
136 {
137  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
138 
139  IOobject ddtIOobject
140  (
141  "ddt("+vf.name()+')',
142  mesh().time().timeName(),
143  mesh().thisDb()
144  );
145 
146  scalar deltaT = deltaT_();
147  scalar deltaT0 = deltaT0_(vf);
148 
149  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
150  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
151  scalar coefft0 = coefft + coefft00;
152 
153  if (mesh().moving())
154  {
155  tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
156  (
157  new GeometricField<Type, fvPatchField, volMesh>
158  (
159  ddtIOobject,
160  mesh(),
161  rDeltaT.dimensions()*vf.dimensions(),
162  rDeltaT.value()*
163  (
164  coefft*vf.primitiveField() -
165  (
166  coefft0*vf.oldTime().primitiveField()*mesh().V0()
167  - coefft00*vf.oldTime().oldTime().primitiveField()
168  *mesh().V00()
169  )/mesh().V()
170  ),
171  rDeltaT.value()*
172  (
173  coefft*vf.boundaryField() -
174  (
175  coefft0*vf.oldTime().boundaryField()
176  - coefft00*vf.oldTime().oldTime().boundaryField()
177  )
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*
197  (
198  coefft*vf
199  - coefft0*vf.oldTime()
200  + coefft00*vf.oldTime().oldTime()
201  )
202  )
203  );
204  }
205 }
206 
207 
208 template<class Type>
209 tmp<GeometricField<Type, fvPatchField, volMesh>>
210 backwardDdtScheme<Type>::fvcDdt
211 (
212  const dimensionedScalar& rho,
213  const GeometricField<Type, fvPatchField, volMesh>& vf
214 )
215 {
216  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
217 
218  IOobject ddtIOobject
219  (
220  "ddt("+rho.name()+','+vf.name()+')',
221  mesh().time().timeName(),
222  mesh().thisDb()
223  );
224 
225  scalar deltaT = deltaT_();
226  scalar deltaT0 = deltaT0_(vf);
227 
228  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
229  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
230  scalar coefft0 = coefft + coefft00;
231 
232  if (mesh().moving())
233  {
234  tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
235  (
236  new GeometricField<Type, fvPatchField, volMesh>
237  (
238  ddtIOobject,
239  mesh(),
240  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
241  rDeltaT.value()*rho.value()*
242  (
243  coefft*vf.primitiveField() -
244  (
245  coefft0*vf.oldTime().primitiveField()*mesh().V0()
246  - coefft00*vf.oldTime().oldTime().primitiveField()
247  *mesh().V00()
248  )/mesh().V()
249  ),
250  rDeltaT.value()*rho.value()*
251  (
252  coefft*vf.boundaryField() -
253  (
254  coefft0*vf.oldTime().boundaryField()
255  - coefft00*vf.oldTime().oldTime().boundaryField()
256  )
257  )
258  )
259  );
260 
261  // Different operation on boundary v.s. internal so re-evaluate
262  // coupled boundaries
263  tdtdt.ref().boundaryFieldRef().
264  template evaluateCoupled<coupledFvPatch>();
265 
266  return tdtdt;
267  }
268  else
269  {
270  return tmp<GeometricField<Type, fvPatchField, volMesh>>
271  (
272  new GeometricField<Type, fvPatchField, volMesh>
273  (
274  ddtIOobject,
275  rDeltaT*rho*
276  (
277  coefft*vf
278  - coefft0*vf.oldTime()
279  + coefft00*vf.oldTime().oldTime()
280  )
281  )
282  );
283  }
284 }
285 
286 
287 template<class Type>
288 tmp<GeometricField<Type, fvPatchField, volMesh>>
289 backwardDdtScheme<Type>::fvcDdt
290 (
291  const volScalarField& rho,
292  const GeometricField<Type, fvPatchField, volMesh>& vf
293 )
294 {
295  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
296 
297  IOobject ddtIOobject
298  (
299  "ddt("+rho.name()+','+vf.name()+')',
300  mesh().time().timeName(),
301  mesh().thisDb()
302  );
303 
304  scalar deltaT = deltaT_();
305  scalar deltaT0 = deltaT0_(vf);
306 
307  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
308  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
309  scalar coefft0 = coefft + coefft00;
310 
311  if (mesh().moving())
312  {
313  tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
314  (
315  new GeometricField<Type, fvPatchField, volMesh>
316  (
317  ddtIOobject,
318  mesh(),
319  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
320  rDeltaT.value()*
321  (
322  coefft*rho.primitiveField()*vf.primitiveField() -
323  (
324  coefft0*rho.oldTime().primitiveField()
325  *vf.oldTime().primitiveField()*mesh().V0()
326  - coefft00*rho.oldTime().oldTime().primitiveField()
327  *vf.oldTime().oldTime().primitiveField()*mesh().V00()
328  )/mesh().V()
329  ),
330  rDeltaT.value()*
331  (
332  coefft*rho.boundaryField()*vf.boundaryField() -
333  (
334  coefft0*rho.oldTime().boundaryField()
335  *vf.oldTime().boundaryField()
336  - coefft00*rho.oldTime().oldTime().boundaryField()
337  *vf.oldTime().oldTime().boundaryField()
338  )
339  )
340  )
341  );
342 
343  // Different operation on boundary v.s. internal so re-evaluate
344  // coupled boundaries
345  tdtdt.ref().boundaryFieldRef().
346  template evaluateCoupled<coupledFvPatch>();
347 
348  return tdtdt;
349  }
350  else
351  {
352  return tmp<GeometricField<Type, fvPatchField, volMesh>>
353  (
354  new GeometricField<Type, fvPatchField, volMesh>
355  (
356  ddtIOobject,
357  rDeltaT*
358  (
359  coefft*rho*vf
360  - coefft0*rho.oldTime()*vf.oldTime()
361  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
362  )
363  )
364  );
365  }
366 }
367 
368 
369 template<class Type>
370 tmp<GeometricField<Type, fvPatchField, volMesh>>
371 backwardDdtScheme<Type>::fvcDdt
372 (
373  const volScalarField& alpha,
374  const volScalarField& rho,
375  const GeometricField<Type, fvPatchField, volMesh>& vf
376 )
377 {
378  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
379 
380  IOobject ddtIOobject
381  (
382  "ddt("+alpha.name()+','+rho.name()+','+vf.name()+')',
383  mesh().time().timeName(),
384  mesh().thisDb()
385  );
386 
387  scalar deltaT = deltaT_();
388  scalar deltaT0 = deltaT0_(vf);
389 
390  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
391  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
392  scalar coefft0 = coefft + coefft00;
393 
394  if (mesh().moving())
395  {
396  tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
397  (
398  new GeometricField<Type, fvPatchField, volMesh>
399  (
400  ddtIOobject,
401  mesh(),
402  rDeltaT.dimensions()
403  *alpha.dimensions()*rho.dimensions()*vf.dimensions(),
404  rDeltaT.value()*
405  (
406  coefft
407  *alpha.primitiveField()
408  *rho.primitiveField()
409  *vf.primitiveField() -
410  (
411  coefft0
412  *alpha.oldTime().primitiveField()
413  *rho.oldTime().primitiveField()
414  *vf.oldTime().primitiveField()*mesh().V0()
415 
416  - coefft00
417  *alpha.oldTime().oldTime().primitiveField()
418  *rho.oldTime().oldTime().primitiveField()
419  *vf.oldTime().oldTime().primitiveField()*mesh().V00()
420  )/mesh().V()
421  ),
422  rDeltaT.value()*
423  (
424  coefft
425  *alpha.boundaryField()
426  *rho.boundaryField()
427  *vf.boundaryField() -
428  (
429  coefft0
430  *alpha.oldTime().boundaryField()
431  *rho.oldTime().boundaryField()
432  *vf.oldTime().boundaryField()
433 
434  - coefft00
435  *alpha.oldTime().oldTime().boundaryField()
436  *rho.oldTime().oldTime().boundaryField()
437  *vf.oldTime().oldTime().boundaryField()
438  )
439  )
440  )
441  );
442 
443  // Different operation on boundary v.s. internal so re-evaluate
444  // coupled boundaries
445  tdtdt.ref().boundaryFieldRef().
446  template evaluateCoupled<coupledFvPatch>();
447 
448  return tdtdt;
449  }
450  else
451  {
452  return tmp<GeometricField<Type, fvPatchField, volMesh>>
453  (
454  new GeometricField<Type, fvPatchField, volMesh>
455  (
456  ddtIOobject,
457  rDeltaT*
458  (
459  coefft*alpha*rho*vf
460  - coefft0*alpha.oldTime()*rho.oldTime()*vf.oldTime()
461  + coefft00*alpha.oldTime().oldTime()
462  *rho.oldTime().oldTime()*vf.oldTime().oldTime()
463  )
464  )
465  );
466  }
467 }
468 
469 
470 template<class Type>
471 tmp<fvMatrix<Type>>
472 backwardDdtScheme<Type>::fvmDdt
473 (
474  const GeometricField<Type, fvPatchField, volMesh>& vf
475 )
476 {
477  tmp<fvMatrix<Type>> tfvm
478  (
479  new fvMatrix<Type>
480  (
481  vf,
482  vf.dimensions()*dimVol/dimTime
483  )
484  );
485 
486  fvMatrix<Type>& fvm = tfvm.ref();
487 
488  scalar rDeltaT = 1.0/deltaT_();
489 
490  scalar deltaT = deltaT_();
491  scalar deltaT0 = deltaT0_(vf);
492 
493  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
494  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
495  scalar coefft0 = coefft + coefft00;
496 
497  fvm.diag() = (coefft*rDeltaT)*mesh().V();
498 
499  if (mesh().moving())
500  {
501  fvm.source() = rDeltaT*
502  (
503  coefft0*vf.oldTime().primitiveField()*mesh().V0()
504  - coefft00*vf.oldTime().oldTime().primitiveField()
505  *mesh().V00()
506  );
507  }
508  else
509  {
510  fvm.source() = rDeltaT*mesh().V()*
511  (
512  coefft0*vf.oldTime().primitiveField()
513  - coefft00*vf.oldTime().oldTime().primitiveField()
514  );
515  }
516 
517  return tfvm;
518 }
519 
520 
521 template<class Type>
522 tmp<fvMatrix<Type>>
523 backwardDdtScheme<Type>::fvmDdt
524 (
525  const dimensionedScalar& rho,
526  const GeometricField<Type, fvPatchField, volMesh>& vf
527 )
528 {
529  tmp<fvMatrix<Type>> tfvm
530  (
531  new fvMatrix<Type>
532  (
533  vf,
534  rho.dimensions()*vf.dimensions()*dimVol/dimTime
535  )
536  );
537  fvMatrix<Type>& fvm = tfvm.ref();
538 
539  scalar rDeltaT = 1.0/deltaT_();
540 
541  scalar deltaT = deltaT_();
542  scalar deltaT0 = deltaT0_(vf);
543 
544  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
545  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
546  scalar coefft0 = coefft + coefft00;
547 
548  fvm.diag() = (coefft*rDeltaT*rho.value())*mesh().V();
549 
550  if (mesh().moving())
551  {
552  fvm.source() = rDeltaT*rho.value()*
553  (
554  coefft0*vf.oldTime().primitiveField()*mesh().V0()
555  - coefft00*vf.oldTime().oldTime().primitiveField()
556  *mesh().V00()
557  );
558  }
559  else
560  {
561  fvm.source() = rDeltaT*mesh().V()*rho.value()*
562  (
563  coefft0*vf.oldTime().primitiveField()
564  - coefft00*vf.oldTime().oldTime().primitiveField()
565  );
566  }
567 
568  return tfvm;
569 }
570 
571 
572 template<class Type>
573 tmp<fvMatrix<Type>>
574 backwardDdtScheme<Type>::fvmDdt
575 (
576  const volScalarField& rho,
577  const GeometricField<Type, fvPatchField, volMesh>& vf
578 )
579 {
580  tmp<fvMatrix<Type>> tfvm
581  (
582  new fvMatrix<Type>
583  (
584  vf,
585  rho.dimensions()*vf.dimensions()*dimVol/dimTime
586  )
587  );
588  fvMatrix<Type>& fvm = tfvm.ref();
589 
590  scalar rDeltaT = 1.0/deltaT_();
591 
592  scalar deltaT = deltaT_();
593  scalar deltaT0 = deltaT0_(vf);
594 
595  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
596  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
597  scalar coefft0 = coefft + coefft00;
598 
599  fvm.diag() = (coefft*rDeltaT)*rho.primitiveField()*mesh().V();
600 
601  if (mesh().moving())
602  {
603  fvm.source() = rDeltaT*
604  (
605  coefft0*rho.oldTime().primitiveField()
606  *vf.oldTime().primitiveField()*mesh().V0()
607  - coefft00*rho.oldTime().oldTime().primitiveField()
608  *vf.oldTime().oldTime().primitiveField()*mesh().V00()
609  );
610  }
611  else
612  {
613  fvm.source() = rDeltaT*mesh().V()*
614  (
615  coefft0*rho.oldTime().primitiveField()
616  *vf.oldTime().primitiveField()
617  - coefft00*rho.oldTime().oldTime().primitiveField()
618  *vf.oldTime().oldTime().primitiveField()
619  );
620  }
621 
622  return tfvm;
623 }
624 
625 
626 template<class Type>
627 tmp<fvMatrix<Type>>
628 backwardDdtScheme<Type>::fvmDdt
629 (
630  const volScalarField& alpha,
631  const volScalarField& rho,
632  const GeometricField<Type, fvPatchField, volMesh>& vf
633 )
634 {
635  tmp<fvMatrix<Type>> tfvm
636  (
637  new fvMatrix<Type>
638  (
639  vf,
640  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
641  )
642  );
643  fvMatrix<Type>& fvm = tfvm.ref();
644 
645  scalar rDeltaT = 1.0/deltaT_();
646 
647  scalar deltaT = deltaT_();
648  scalar deltaT0 = deltaT0_(vf);
649 
650  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
651  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
652  scalar coefft0 = coefft + coefft00;
653 
654  fvm.diag() =
655  (coefft*rDeltaT)*alpha.primitiveField()*rho.primitiveField()*mesh().V();
656 
657  if (mesh().moving())
658  {
659  fvm.source() = rDeltaT*
660  (
661  coefft0
662  *alpha.oldTime().primitiveField()
663  *rho.oldTime().primitiveField()
664  *vf.oldTime().primitiveField()*mesh().V0()
665 
666  - coefft00
667  *alpha.oldTime().oldTime().primitiveField()
668  *rho.oldTime().oldTime().primitiveField()
669  *vf.oldTime().oldTime().primitiveField()*mesh().V00()
670  );
671  }
672  else
673  {
674  fvm.source() = rDeltaT*mesh().V()*
675  (
676  coefft0
677  *alpha.oldTime().primitiveField()
678  *rho.oldTime().primitiveField()
679  *vf.oldTime().primitiveField()
680 
681  - coefft00
682  *alpha.oldTime().oldTime().primitiveField()
683  *rho.oldTime().oldTime().primitiveField()
684  *vf.oldTime().oldTime().primitiveField()
685  );
686  }
687 
688  return tfvm;
689 }
690 
691 
692 template<class Type>
693 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
694 backwardDdtScheme<Type>::fvcDdtUfCorr
695 (
696  const GeometricField<Type, fvPatchField, volMesh>& U,
697  const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
698 )
699 {
700  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
701 
702  scalar deltaT = deltaT_();
703  scalar deltaT0 = deltaT0_(U);
704 
705  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
706  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
707  scalar coefft0 = coefft + coefft00;
708 
709  return tmp<fluxFieldType>
710  (
711  new fluxFieldType
712  (
713  IOobject
714  (
715  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
716  mesh().time().timeName(),
717  mesh().thisDb()
718  ),
719  this->fvcDdtPhiCoeff(U.oldTime(), (mesh().Sf() & Uf.oldTime()))
720  *rDeltaT
721  *(
722  mesh().Sf()
723  & (
724  (coefft0*Uf.oldTime() - coefft00*Uf.oldTime().oldTime())
725  - fvc::interpolate
726  (
727  coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
728  )
729  )
730  )
731  )
732  );
733 }
734 
735 
736 template<class Type>
737 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
738 backwardDdtScheme<Type>::fvcDdtPhiCorr
739 (
740  const GeometricField<Type, fvPatchField, volMesh>& U,
741  const fluxFieldType& phi
742 )
743 {
744  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
745 
746  scalar deltaT = deltaT_();
747  scalar deltaT0 = deltaT0_(U);
748 
749  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
750  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
751  scalar coefft0 = coefft + coefft00;
752 
753  return tmp<fluxFieldType>
754  (
755  new fluxFieldType
756  (
757  IOobject
758  (
759  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
760  mesh().time().timeName(),
761  mesh().thisDb()
762  ),
763  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
764  *rDeltaT
765  *(
766  (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
767  - fvc::dotInterpolate
768  (
769  mesh().Sf(),
770  coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
771  )
772  )
773  )
774  );
775 }
776 
777 
778 template<class Type>
779 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
780 backwardDdtScheme<Type>::fvcDdtUfCorr
781 (
782  const volScalarField& rho,
783  const GeometricField<Type, fvPatchField, volMesh>& U,
784  const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
785 )
786 {
787  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
788 
789  scalar deltaT = deltaT_();
790  scalar deltaT0 = deltaT0_(U);
791 
792  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
793  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
794  scalar coefft0 = coefft + coefft00;
795 
796  if
797  (
798  U.dimensions() == dimVelocity
799  && Uf.dimensions() == rho.dimensions()*dimVelocity
800  )
801  {
802  GeometricField<Type, fvPatchField, volMesh> rhoU0
803  (
804  rho.oldTime()*U.oldTime()
805  );
806 
807  GeometricField<Type, fvPatchField, volMesh> rhoU00
808  (
809  rho.oldTime().oldTime()*U.oldTime().oldTime()
810  );
811 
812  return tmp<fluxFieldType>
813  (
814  new fluxFieldType
815  (
816  IOobject
817  (
818  "ddtCorr("
819  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
820  mesh().time().timeName(),
821  mesh().thisDb()
822  ),
823  this->fvcDdtPhiCoeff
824  (
825  rhoU0,
826  mesh().Sf() & Uf.oldTime(),
827  rho.oldTime()
828  )
829  *rDeltaT
830  *(
831  mesh().Sf()
832  & (
833  (coefft0*Uf.oldTime() - coefft00*Uf.oldTime().oldTime())
834  - fvc::interpolate(coefft0*rhoU0 - coefft00*rhoU00)
835  )
836  )
837  )
838  );
839  }
840  else if
841  (
842  U.dimensions() == rho.dimensions()*dimVelocity
843  && Uf.dimensions() == rho.dimensions()*dimVelocity
844  )
845  {
846  return tmp<fluxFieldType>
847  (
848  new fluxFieldType
849  (
850  IOobject
851  (
852  "ddtCorr("
853  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
854  mesh().time().timeName(),
855  mesh().thisDb()
856  ),
857  this->fvcDdtPhiCoeff
858  (
859  U.oldTime(),
860  mesh().Sf() & Uf.oldTime(),
861  rho.oldTime()
862  )
863  *rDeltaT
864  *(
865  mesh().Sf()
866  & (
867  (coefft0*Uf.oldTime() - coefft00*Uf.oldTime().oldTime())
868  - fvc::interpolate
869  (
870  coefft0*U.oldTime()
871  - coefft00*U.oldTime().oldTime()
872  )
873  )
874  )
875  )
876  );
877  }
878  else
879  {
880  FatalErrorInFunction
881  << "dimensions of phi are not correct"
882  << abort(FatalError);
883 
884  return fluxFieldType::null();
885  }
886 }
887 
888 
889 template<class Type>
890 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
891 backwardDdtScheme<Type>::fvcDdtPhiCorr
892 (
893  const volScalarField& rho,
894  const GeometricField<Type, fvPatchField, volMesh>& U,
895  const fluxFieldType& phi
896 )
897 {
898  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
899 
900  scalar deltaT = deltaT_();
901  scalar deltaT0 = deltaT0_(U);
902 
903  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
904  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
905  scalar coefft0 = coefft + coefft00;
906 
907  if
908  (
909  U.dimensions() == dimVelocity
910  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
911  )
912  {
913  GeometricField<Type, fvPatchField, volMesh> rhoU0
914  (
915  rho.oldTime()*U.oldTime()
916  );
917 
918  GeometricField<Type, fvPatchField, volMesh> rhoU00
919  (
920  rho.oldTime().oldTime()*U.oldTime().oldTime()
921  );
922 
923  return tmp<fluxFieldType>
924  (
925  new fluxFieldType
926  (
927  IOobject
928  (
929  "ddtCorr("
930  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
931  mesh().time().timeName(),
932  mesh().thisDb()
933  ),
934  this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), rho.oldTime())
935  *rDeltaT
936  *(
937  (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
938  - fvc::dotInterpolate
939  (
940  mesh().Sf(),
941  coefft0*rhoU0 - coefft00*rhoU00
942  )
943  )
944  )
945  );
946  }
947  else if
948  (
949  U.dimensions() == rho.dimensions()*dimVelocity
950  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
951  )
952  {
953  return tmp<fluxFieldType>
954  (
955  new fluxFieldType
956  (
957  IOobject
958  (
959  "ddtCorr("
960  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
961  mesh().time().timeName(),
962  mesh().thisDb()
963  ),
964  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), rho.oldTime())
965  *rDeltaT
966  *(
967  (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
968  - fvc::dotInterpolate
969  (
970  mesh().Sf(),
971  coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
972  )
973  )
974  )
975  );
976  }
977  else
978  {
979  FatalErrorInFunction
980  << "dimensions of phi are not correct"
981  << abort(FatalError);
982 
983  return fluxFieldType::null();
984  }
985 }
986 
987 
988 template<class Type>
989 tmp<surfaceScalarField> backwardDdtScheme<Type>::meshPhi
990 (
991  const GeometricField<Type, fvPatchField, volMesh>& vf
992 )
993 {
994  scalar deltaT = deltaT_();
995  scalar deltaT0 = deltaT0_(vf);
996 
997  // Coefficient for t-3/2 (between times 0 and 00)
998  scalar coefft0_00 = deltaT/(deltaT + deltaT0);
999 
1000  // Coefficient for t-1/2 (between times n and 0)
1001  scalar coefftn_0 = 1 + coefft0_00;
1002 
1003  return tmp<surfaceScalarField>
1004  (
1005  new surfaceScalarField
1006  (
1007  IOobject
1008  (
1009  mesh().phi().name(),
1010  mesh().time().timeName(),
1011  mesh().thisDb(),
1012  IOobject::NO_READ,
1013  IOobject::NO_WRITE,
1014  IOobject::NO_REGISTER
1015  ),
1016  coefftn_0*mesh().phi() - coefft0_00*mesh().phi().oldTime()
1017  )
1018  );
1019 }
1020 
1021 
1022 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1023 
1024 } // End namespace fv
1025 
1026 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1027 
1028 } // End namespace Foam
1029 
1030 // ************************************************************************* //
Foam::dimensioned::value
const Type & value() const noexcept
Return const reference to value.
Definition: dimensionedType.H:366
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:1006
oldTime
Info<< "Creating field kinetic energy K\"<< endl;volScalarField K("K", 0.5 *magSqr(U));if(U.nOldTimes()){ volVectorField *Uold=&U.oldTime();volScalarField *Kold=&K.oldTime();*Kold==0.5 *magSqr(*Uold);while(Uold->nOldTimes()) { Uold=&Uold-> oldTime()
Definition: createK.H:12
Foam::fv::backwardDdtScheme::fvmDdt
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: backwardDdtScheme.C:466
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const noexcept
Return a const-reference to the internal field values.
Definition: GeometricFieldI.H:44
Foam::FatalError
error FatalError
Error stream (stdout output on all processes), with additional &#39;FOAM FATAL ERROR&#39; header text and sta...
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:598
Foam::IOobject::name
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:195
Foam::tmp::ref
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
Definition: tmpI.H:235
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::New
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.
Definition: DimensionedFieldReuseFunctions.H:124
Foam::fv::backwardDdtScheme::fvcDdt
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
Definition: backwardDdtScheme.C:73
Foam::dimVol
const dimensionSet dimVol(dimVolume)
Older spelling for dimVolume.
Definition: dimensionSets.H:65
Foam::GeometricField
Generic GeometricField class.
Definition: areaFieldsFwd.H:50
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedScalarFwd.H:33
Foam::IOobjectOption::NO_WRITE
Ignore writing from objectRegistry::writeObject()
Definition: IOobjectOption.H:86
Foam::fv::backwardDdtScheme::meshPhi
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: backwardDdtScheme.C:983
timeName
word timeName
Definition: getTimeIndex.H:3
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
Definition: exprTraits.C:127
fv
labelList fv(nPoints)
Foam::fv::backwardDdtScheme::fvcDdtPhiCorr
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
Definition: backwardDdtScheme.C:732
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:64
Uf
autoPtr< surfaceVectorField > Uf
Definition: createUfIfPresent.H:27
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const noexcept
Return const reference to dimensions.
Definition: dimensionedType.H:356
Foam::fv::backwardDdtScheme::fvcDdtUfCorr
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition: backwardDdtScheme.C:688
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:139
fvcDiv.H
Calculate the divergence of the given field.
Foam::fvc::interpolate
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Foam::dimensioned::name
const word & name() const noexcept
Return const reference to name.
Definition: dimensionedType.H:346
Foam::fvMatrix::source
Field< Type > & source() noexcept
Definition: fvMatrix.H:533
U
U
Definition: pEqn.H:72
Foam::IOobjectOption::NO_READ
Nothing to be read.
Definition: IOobjectOption.H:60
Foam::fvPatchFieldBase::calculatedType
static const word & calculatedType() noexcept
The type name for calculated patch fields.
Definition: fvPatchField.H:204
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: HashPtrTable.H:50
Foam::surfaceScalarField
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Definition: surfaceFieldsFwd.H:64
Foam::lduMatrix::diag
scalarField & diag()
Definition: lduMatrix.C:197
Foam::constant::atomic::alpha
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Definition: readThermalProperties.H:212
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported...
Definition: IOobject.H:172
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const noexcept
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:53
Foam::IOobjectOption::NO_REGISTER
Do not request registration (bool: false)
Definition: IOobjectOption.H:99
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:26
timeIndex
label timeIndex
Definition: getTimeIndex.H:24
Foam::DimensionedField::dimensions
const dimensionSet & dimensions() const noexcept
Return dimensions.
Definition: DimensionedFieldI.H:42
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:127
Foam::dimVelocity
const dimensionSet dimVelocity