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-2022 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<GeometricField<Type, fvPatchField, volMesh>>
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<GeometricField<Type, fvPatchField, volMesh>> tdtdt
95  (
96  GeometricField<Type, fvPatchField, volMesh>::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<GeometricField<Type, fvPatchField, volMesh>>
119  (
120  GeometricField<Type, fvPatchField, volMesh>::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<GeometricField<Type, fvPatchField, volMesh>>
139 backwardDdtScheme<Type>::fvcDdt
140 (
141  const GeometricField<Type, fvPatchField, volMesh>& 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<GeometricField<Type, fvPatchField, volMesh>>
158  (
159  GeometricField<Type, fvPatchField, volMesh>::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<GeometricField<Type, fvPatchField, volMesh>>
185  (
186  GeometricField<Type, fvPatchField, volMesh>::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<GeometricField<Type, fvPatchField, volMesh>>
203 backwardDdtScheme<Type>::fvcDdt
204 (
205  const dimensionedScalar& rho,
206  const GeometricField<Type, fvPatchField, volMesh>& 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<GeometricField<Type, fvPatchField, volMesh>>
223  (
224  GeometricField<Type, fvPatchField, volMesh>::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<GeometricField<Type, fvPatchField, volMesh>>
250  (
251  GeometricField<Type, fvPatchField, volMesh>::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<GeometricField<Type, fvPatchField, volMesh>>
268 backwardDdtScheme<Type>::fvcDdt
269 (
270  const volScalarField& rho,
271  const GeometricField<Type, fvPatchField, volMesh>& 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<GeometricField<Type, fvPatchField, volMesh>>
288  (
289  GeometricField<Type, fvPatchField, volMesh>::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<GeometricField<Type, fvPatchField, volMesh>>
318  (
319  GeometricField<Type, fvPatchField, volMesh>::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<GeometricField<Type, fvPatchField, volMesh>>
336 backwardDdtScheme<Type>::fvcDdt
337 (
338  const volScalarField& alpha,
339  const volScalarField& rho,
340  const GeometricField<Type, fvPatchField, volMesh>& 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<GeometricField<Type, fvPatchField, volMesh>>
357  (
358  GeometricField<Type, fvPatchField, volMesh>::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<GeometricField<Type, fvPatchField, volMesh>>
397  (
398  GeometricField<Type, fvPatchField, volMesh>::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 GeometricField<Type, fvPatchField, volMesh>& 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 GeometricField<Type, fvPatchField, volMesh>& 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 GeometricField<Type, fvPatchField, volMesh>& 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 GeometricField<Type, fvPatchField, volMesh>& 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 GeometricField<Type, fvPatchField, volMesh>& U,
641  const GeometricField<Type, fvsPatchField, surfaceMesh>& 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 GeometricField<Type, fvPatchField, volMesh>& 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 GeometricField<Type, fvPatchField, volMesh>& U,
712  const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
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  && Uf.dimensions() == rho.dimensions()*dimVelocity
728  )
729  {
730  const GeometricField<Type, fvPatchField, volMesh> rhoU0
731  (
732  rho.oldTime()*U.oldTime()
733  );
734 
735  const GeometricField<Type, fvPatchField, volMesh> rhoU00
736  (
737  rho.oldTime().oldTime()*U.oldTime().oldTime()
738  );
739 
740  return fluxFieldType::New
741  (
742  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
743  this->fvcDdtPhiCoeff
744  (
745  rhoU0,
746  mesh().Sf() & Uf.oldTime(),
747  rho.oldTime()
748  )
749  *rDeltaT
750  *(
751  mesh().Sf()
752  & (
753  (coefft0*Uf.oldTime() - coefft00*Uf.oldTime().oldTime())
754  - fvc::interpolate(coefft0*rhoU0 - coefft00*rhoU00)
755  )
756  )
757  );
758  }
759  else if
760  (
761  U.dimensions() == rho.dimensions()*dimVelocity
762  && Uf.dimensions() == rho.dimensions()*dimVelocity
763  )
764  {
765  return fluxFieldType::New
766  (
767  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
768  this->fvcDdtPhiCoeff
769  (
770  U.oldTime(),
771  mesh().Sf() & Uf.oldTime(),
772  rho.oldTime()
773  )
774  *rDeltaT
775  *(
776  mesh().Sf()
777  & (
778  (coefft0*Uf.oldTime() - coefft00*Uf.oldTime().oldTime())
779  - fvc::interpolate
780  (
781  coefft0*U.oldTime()
782  - coefft00*U.oldTime().oldTime()
783  )
784  )
785  )
786  );
787  }
788  else
789  {
790  FatalErrorInFunction
791  << "dimensions of phi are not correct"
792  << abort(FatalError);
793 
794  return fluxFieldType::null();
795  }
796 }
797 
798 
799 template<class Type>
800 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
801 backwardDdtScheme<Type>::fvcDdtPhiCorr
802 (
803  const volScalarField& rho,
804  const GeometricField<Type, fvPatchField, volMesh>& U,
805  const fluxFieldType& phi
806 )
807 {
808  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
809 
810  const scalar deltaT = deltaT_();
811  const scalar deltaT0 = deltaT0_(U);
812 
813  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
814  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
815  const scalar coefft0 = coefft + coefft00;
816 
817  if
818  (
819  U.dimensions() == dimVelocity
820  && phi.dimensions() == rho.dimensions()*dimFlux
821  )
822  {
823  const GeometricField<Type, fvPatchField, volMesh> rhoU0
824  (
825  rho.oldTime()*U.oldTime()
826  );
827 
828  const GeometricField<Type, fvPatchField, volMesh> rhoU00
829  (
830  rho.oldTime().oldTime()*U.oldTime().oldTime()
831  );
832 
833  return fluxFieldType::New
834  (
835  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
836  this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), rho.oldTime())
837  *rDeltaT
838  *(
839  (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
840  - fvc::dotInterpolate
841  (
842  mesh().Sf(),
843  coefft0*rhoU0 - coefft00*rhoU00
844  )
845  )
846  );
847  }
848  else if
849  (
850  U.dimensions() == rho.dimensions()*dimVelocity
851  && phi.dimensions() == rho.dimensions()*dimFlux
852  )
853  {
854  return fluxFieldType::New
855  (
856  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
857  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), rho.oldTime())
858  *rDeltaT
859  *(
860  (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
861  - fvc::dotInterpolate
862  (
863  mesh().Sf(),
864  coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
865  )
866  )
867  );
868  }
869  else
870  {
871  FatalErrorInFunction
872  << "dimensions of phi are not correct"
873  << abort(FatalError);
874 
875  return fluxFieldType::null();
876  }
877 }
878 
879 
880 template<class Type>
881 tmp<surfaceScalarField> backwardDdtScheme<Type>::meshPhi
882 (
883  const GeometricField<Type, fvPatchField, volMesh>& vf
884 )
885 {
886  const scalar deltaT = deltaT_();
887  const scalar deltaT0 = deltaT0_(vf);
888 
889  // Coefficient for t-3/2 (between times 0 and 00)
890  const scalar coefft0_00 = deltaT/(deltaT + deltaT0);
891 
892  // Coefficient for t-1/2 (between times n and 0)
893  const scalar coefftn_0 = 1 + coefft0_00;
894 
895  return surfaceScalarField::New
896  (
897  mesh().phi().name(),
898  coefftn_0*mesh().phi() - coefft0_00*mesh().phi().oldTime()
899  );
900 }
901 
902 
903 template<class Type>
904 tmp<scalarField> backwardDdtScheme<Type>::meshPhi
905 (
906  const GeometricField<Type, fvPatchField, volMesh>& vf,
907  const label patchi
908 )
909 {
910  const scalar deltaT = deltaT_();
911  const scalar deltaT0 = deltaT0_(vf);
912 
913  // Coefficient for t-3/2 (between times 0 and 00)
914  const scalar coefft0_00 = deltaT/(deltaT + deltaT0);
915 
916  // Coefficient for t-1/2 (between times n and 0)
917  const scalar coefftn_0 = 1 + coefft0_00;
918 
919  return
920  (
921  coefftn_0*mesh().phi().boundaryField()[patchi]
922  - coefft0_00*mesh().phi().oldTime().boundaryField()[patchi]
923  );
924 }
925 
926 
927 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
928 
929 } // End namespace fv
930 
931 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
932 
933 } // End namespace Foam
934 
935 // ************************************************************************* //
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:1108
Foam::fv::backwardDdtScheme::fvmDdt
virtual tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: backwardDdtScheme.C:417
Foam::IOobject::name
const word & name() const
Return name.
Definition: IOobject.H:315
Foam::compressible::New
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
Definition: compressibleMomentumTransportModelTemplates.C:34
Foam::FatalError
error FatalError
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:60
Foam::tmp::ref
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
alpha
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
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::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::fv::backwardDdtScheme::fvcDdt
virtual tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
Definition: backwardDdtScheme.C:77
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Definition: GeometricFieldI.H:51
Foam::GeometricField
Generic GeometricField class.
Definition: surfaceFieldsFwd.H:52
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedScalarFwd.H:41
mesh
fvMesh & mesh
Definition: createMeshesPostProcess.H:9
Foam::fv::backwardDdtScheme::meshPhi
virtual tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: backwardDdtScheme.C:882
Foam::dimTime
const dimensionSet dimTime
Foam::DimensionedField::dimensions
const dimensionSet & dimensions() const
Return dimensions.
Definition: DimensionedFieldI.H:46
Foam::word
A class for handling words, derived from string.
Definition: word.H:59
fv
labelList fv(nPoints)
Foam::dimVol
const dimensionSet dimVol
Foam::dimFlux
const dimensionSet dimFlux
Foam::calculatedFvPatchField
This boundary condition is not designed to be evaluated; it is assumed that the value is assigned via...
Definition: calculatedFvPatchField.H:61
Foam::dimensioned::value
const Type & value() const
Return const reference to value.
Definition: dimensionedType.C:283
Foam::fv::backwardDdtScheme::fvcDdtPhiCorr
virtual tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
Definition: backwardDdtScheme.C:675
phi
phi
Definition: correctPhi.H:3
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:72
Foam::Zero
static const zero Zero
Definition: zero.H:97
Foam::fv::backwardDdtScheme::fvcDdtUfCorr
virtual tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition: backwardDdtScheme.C:639
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:131
fvcDiv.H
Calculate the divergence of the given field.
Foam::dimVelocity
const dimensionSet dimVelocity
Foam::fvMatrix::source
Field< Type > & source()
Definition: fvMatrix.H:300
Foam::dimensioned::name
const word & name() const
Return const reference to name.
Definition: dimensionedType.C:257
Foam::name
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
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.
oldTime
rho oldTime()
patchi
label patchi
Definition: getPatchFieldScalar.H:1
fvMatrices.H
A special matrix type and solver, designed for finite volume solutions of scalar equations.
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const
Return const reference to dimensions.
Definition: dimensionedType.C:270
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:53
Foam::lduMatrix::diag
scalarField & diag()
Definition: lduMatrix.C:186
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.H:213