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