backwardFaDdtScheme.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) 2016-2017 Wikki Ltd
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "backwardFaDdtScheme.H"
29 #include "faMatrices.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 namespace fa
39 {
40 
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 
43 template<class Type>
44 scalar backwardFaDdtScheme<Type>::deltaT_() const
45 {
46  return mesh().time().deltaTValue();
47 }
48 
49 
50 template<class Type>
51 scalar backwardFaDdtScheme<Type>::deltaT0_() const
52 {
53  return mesh().time().deltaT0Value();
54 }
55 
56 
57 template<class Type>
58 template<class GeoField>
59 scalar backwardFaDdtScheme<Type>::deltaT0_(const GeoField& vf) const
60 {
61  if (vf.oldTime().timeIndex() == vf.oldTime().oldTime().timeIndex())
62  {
63  return GREAT;
64  }
65  else
66  {
67  return mesh().time().deltaT0Value();
68  }
69 }
70 
71 
72 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
73 
74 template<class Type>
77 (
78  const dimensioned<Type> dt
79 )
80 {
81  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
82 
83  const IOobject ddtIOobject
84  (
85  mesh().thisDb().newIOobject
86  (
87  "ddt("+dt.name()+')',
89  )
90  );
91 
92  scalar deltaT = deltaT_();
93  scalar deltaT0 = deltaT0_();
94 
95  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
96  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
97  scalar coefft0 = coefft + coefft00;
98 
99  if (mesh().moving())
100  {
102  (
104  (
105  ddtIOobject,
106  mesh(),
108  )
109  );
110 
111  tdtdt.ref().primitiveFieldRef() = rDeltaT.value()*dt.value()*
112  (
113  coefft - (coefft0*mesh().S0() - coefft00*mesh().S00())/mesh().S()
114  );
115 
116  return tdtdt;
117  }
118 
120  (
121  ddtIOobject,
122  mesh(),
125  );
126 }
127 
128 
129 template<class Type>
132 (
133  const dimensioned<Type> dt
134 )
135 {
136  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
137 
138  const IOobject ddtIOobject
139  (
140  mesh().thisDb().newIOobject
141  (
142  "ddt("+dt.name()+')',
144  )
145  );
146 
147  scalar deltaT = deltaT_();
148  scalar deltaT0 = deltaT0_();
149 
150  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
151  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
152  scalar coefft0 = coefft + coefft00;
153 
155  (
157  (
158  ddtIOobject,
159  mesh(),
160  -rDeltaT*(coefft0 - coefft00)*dt
161  )
162  );
163 
164  if (mesh().moving())
165  {
166  tdtdt0.ref().primitiveFieldRef() = (-rDeltaT.value()*dt.value())*
167  (
168  (coefft0*mesh().S0() - coefft00*mesh().S00())/mesh().S()
169  );
170  }
171 
172  return tdtdt0;
173 }
174 
175 
176 template<class Type>
179 (
181 )
182 {
183  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
184 
185  const IOobject ddtIOobject
186  (
187  mesh().thisDb().newIOobject
188  (
189  "ddt("+vf.name()+')',
191  )
192  );
193 
194  scalar deltaT = deltaT_();
195  scalar deltaT0 = deltaT0_(vf);
196 
197  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
198  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
199  scalar coefft0 = coefft + coefft00;
200 
201  if (mesh().moving())
202  {
204  (
206  (
207  ddtIOobject,
208  mesh(),
209  rDeltaT.dimensions()*vf.dimensions(),
210  rDeltaT.value()*
211  (
212  coefft*vf() -
213  (
214  coefft0*vf.oldTime()()*mesh().S0()
215  - coefft00*vf.oldTime().oldTime()()
216  *mesh().S00()
217  )/mesh().S()
218  ),
219  rDeltaT.value()*
220  (
221  coefft*vf.boundaryField() -
222  (
223  coefft0*vf.oldTime().boundaryField()
224  - coefft00*vf.oldTime().oldTime().boundaryField()
225  )
226  )
227  )
228  );
229  }
230  else
231  {
232  return tmp<GeometricField<Type, faPatchField, areaMesh>>
233  (
234  new GeometricField<Type, faPatchField, areaMesh>
235  (
236  ddtIOobject,
237  rDeltaT*
238  (
239  coefft*vf
240  - coefft0*vf.oldTime()
241  + coefft00*vf.oldTime().oldTime()
242  )
243  )
244  );
245  }
246 }
247 
248 
249 template<class Type>
252 (
254 )
255 {
256  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
257 
258  const IOobject ddtIOobject
259  (
260  mesh().thisDb().newIOobject
261  (
262  "ddt0("+vf.name()+')',
264  )
265  );
266 
267  scalar deltaT = deltaT_();
268  scalar deltaT0 = deltaT0_(vf);
269 
270  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
271  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
272  scalar coefft0 = coefft + coefft00;
273 
274  if (mesh().moving())
275  {
276  return tmp<GeometricField<Type, faPatchField, areaMesh>>
277  (
278  new GeometricField<Type, faPatchField, areaMesh>
279  (
280  ddtIOobject,
281  mesh(),
282  rDeltaT.dimensions()*vf.dimensions(),
283  rDeltaT.value()*
284  (
285  - (
286  coefft0*vf.oldTime()()*mesh().S0()
287  - coefft00*vf.oldTime().oldTime()()
288  *mesh().S00()
289  )/mesh().S()
290  ),
291  rDeltaT.value()*
292  (
293  - (
294  coefft0*vf.oldTime().boundaryField()
295  - coefft00*vf.oldTime().oldTime().boundaryField()
296  )
297  )
298  )
299  );
300  }
301  else
302  {
303  return tmp<GeometricField<Type, faPatchField, areaMesh>>
304  (
305  new GeometricField<Type, faPatchField, areaMesh>
306  (
307  ddtIOobject,
308  rDeltaT*
309  (
310  - coefft0*vf.oldTime()
311  + coefft00*vf.oldTime().oldTime()
312  )
313  )
314  );
315  }
316 }
317 
318 
319 template<class Type>
322 (
323  const dimensionedScalar& rho,
325 )
326 {
327  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
328 
329  const IOobject ddtIOobject
330  (
331  mesh().thisDb().newIOobject
332  (
333  "ddt("+rho.name()+','+vf.name()+')',
335  )
336  );
337 
338  scalar deltaT = deltaT_();
339  scalar deltaT0 = deltaT0_(vf);
340 
341  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
342  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
343  scalar coefft0 = coefft + coefft00;
344 
345  if (mesh().moving())
346  {
347  return tmp<GeometricField<Type, faPatchField, areaMesh>>
348  (
349  new GeometricField<Type, faPatchField, areaMesh>
350  (
351  ddtIOobject,
352  mesh(),
353  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
354  rDeltaT.value()*rho.value()*
355  (
356  coefft*vf.internalField() -
357  (
358  coefft0*vf.oldTime()()*mesh().S0()
359  - coefft00*vf.oldTime().oldTime()()
360  *mesh().S00()
361  )/mesh().S()
362  ),
363  rDeltaT.value()*rho.value()*
364  (
365  coefft*vf.boundaryField() -
366  (
367  coefft0*vf.oldTime().boundaryField()
368  - coefft00*vf.oldTime().oldTime().boundaryField()
369  )
370  )
371  )
372  );
373  }
374  else
375  {
376  return tmp<GeometricField<Type, faPatchField, areaMesh>>
377  (
378  new GeometricField<Type, faPatchField, areaMesh>
379  (
380  ddtIOobject,
381  rDeltaT*rho*
382  (
383  coefft*vf
384  - coefft0*vf.oldTime()
385  + coefft00*vf.oldTime().oldTime()
386  )
387  )
388  );
389  }
390 }
391 
392 template<class Type>
395 (
396  const dimensionedScalar& rho,
398 )
399 {
400  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
401 
402  const IOobject ddtIOobject
403  (
404  mesh().thisDb().newIOobject
405  (
406  "ddt0("+rho.name()+','+vf.name()+')',
408  )
409  );
410 
411  scalar deltaT = deltaT_();
412  scalar deltaT0 = deltaT0_(vf);
413 
414  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
415  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
416  scalar coefft0 = coefft + coefft00;
417 
418  if (mesh().moving())
419  {
420  return tmp<GeometricField<Type, faPatchField, areaMesh>>
421  (
422  new GeometricField<Type, faPatchField, areaMesh>
423  (
424  ddtIOobject,
425  mesh(),
426  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
427  rDeltaT.value()*rho.value()*
428  (
429  -(
430  coefft0*vf.oldTime()()*mesh().S0()
431  - coefft00*vf.oldTime().oldTime()()
432  *mesh().S00()
433  )/mesh().S()
434  ),
435  rDeltaT.value()*rho.value()*
436  (
437  -(
438  coefft0*vf.oldTime().boundaryField()
439  - coefft00*vf.oldTime().oldTime().boundaryField()
440  )
441  )
442  )
443  );
444  }
445  else
446  {
447  return tmp<GeometricField<Type, faPatchField, areaMesh>>
448  (
449  new GeometricField<Type, faPatchField, areaMesh>
450  (
451  ddtIOobject,
452  rDeltaT*rho*
453  (
454  - coefft0*vf.oldTime()
455  + coefft00*vf.oldTime().oldTime()
456  )
457  )
458  );
459  }
460 }
461 
462 
463 template<class Type>
466 (
467  const areaScalarField& rho,
469 )
470 {
471  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
472 
473  const IOobject ddtIOobject
474  (
475  mesh().thisDb().newIOobject
476  (
477  "ddt("+rho.name()+','+vf.name()+')',
479  )
480  );
481 
482  scalar deltaT = deltaT_();
483  scalar deltaT0 = deltaT0_(vf);
484 
485  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
486  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
487  scalar coefft0 = coefft + coefft00;
488 
489  if (mesh().moving())
490  {
491  return tmp<GeometricField<Type, faPatchField, areaMesh>>
492  (
493  new GeometricField<Type, faPatchField, areaMesh>
494  (
495  ddtIOobject,
496  mesh(),
497  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
498  rDeltaT.value()*
499  (
500  coefft*rho.internalField()*vf.internalField() -
501  (
502  coefft0*rho.oldTime()()
503  *vf.oldTime()()*mesh().S0()
504  - coefft00*rho.oldTime().oldTime()()
505  *vf.oldTime().oldTime()()*mesh().S00()
506  )/mesh().S()
507  ),
508  rDeltaT.value()*
509  (
510  coefft*vf.boundaryField() -
511  (
512  coefft0*rho.oldTime().boundaryField()
513  *vf.oldTime().boundaryField()
514  - coefft00*rho.oldTime().oldTime().boundaryField()
515  *vf.oldTime().oldTime().boundaryField()
516  )
517  )
518  )
519  );
520  }
521  else
522  {
523  return tmp<GeometricField<Type, faPatchField, areaMesh>>
524  (
525  new GeometricField<Type, faPatchField, areaMesh>
526  (
527  ddtIOobject,
528  rDeltaT*
529  (
530  coefft*rho*vf
531  - coefft0*rho.oldTime()*vf.oldTime()
532  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
533  )
534  )
535  );
536  }
537 }
538 
539 
540 template<class Type>
543 (
544  const areaScalarField& rho,
546 )
547 {
548  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
549 
550  const IOobject ddtIOobject
551  (
552  mesh().thisDb().newIOobject
553  (
554  "ddt0("+rho.name()+','+vf.name()+')',
556  )
557  );
558 
559  scalar deltaT = deltaT_();
560  scalar deltaT0 = deltaT0_(vf);
561 
562  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
563  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
564  scalar coefft0 = coefft + coefft00;
565 
566  if (mesh().moving())
567  {
568  return tmp<GeometricField<Type, faPatchField, areaMesh>>
569  (
570  new GeometricField<Type, faPatchField, areaMesh>
571  (
572  ddtIOobject,
573  mesh(),
574  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
575  rDeltaT.value()*
576  (
577  - (
578  coefft0*rho.oldTime()()
579  *vf.oldTime()()*mesh().S0()
580  - coefft00*rho.oldTime().oldTime()()
581  *vf.oldTime().oldTime()()*mesh().S00()
582  )/mesh().S()
583  ),
584  rDeltaT.value()*
585  (
586  - (
587  coefft0*rho.oldTime().boundaryField()
588  *vf.oldTime().boundaryField()
589  - coefft00*rho.oldTime().oldTime().boundaryField()
590  *vf.oldTime().oldTime().boundaryField()
591  )
592  )
593  )
594  );
595  }
596  else
597  {
598  return tmp<GeometricField<Type, faPatchField, areaMesh>>
599  (
600  new GeometricField<Type, faPatchField, areaMesh>
601  (
602  ddtIOobject,
603  rDeltaT*
604  (
605  - coefft0*rho.oldTime()*vf.oldTime()
606  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
607  )
608  )
609  );
610  }
611 }
612 
613 
614 template<class Type>
617 (
619 )
620 {
621  tmp<faMatrix<Type>> tfam
622  (
623  new faMatrix<Type>
624  (
625  vf,
627  )
628  );
629 
630  faMatrix<Type>& fam = tfam.ref();
631 
632  scalar rDeltaT = 1.0/deltaT_();
633 
634  scalar deltaT = deltaT_();
635  scalar deltaT0 = deltaT0_(vf);
636 
637  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
638  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
639  scalar coefft0 = coefft + coefft00;
640 
641  fam.diag() = (coefft*rDeltaT)*mesh().S();
642 
643  if (mesh().moving())
644  {
645  fam.source() = rDeltaT*
646  (
647  coefft0*vf.oldTime().primitiveField()*mesh().S0()
648  - coefft00*vf.oldTime().oldTime().primitiveField()
649  *mesh().S00()
650  );
651  }
652  else
653  {
654  fam.source() = rDeltaT*mesh().S()*
655  (
656  coefft0*vf.oldTime().primitiveField()
657  - coefft00*vf.oldTime().oldTime().primitiveField()
658  );
659  }
660 
661  return tfam;
662 }
663 
664 
665 template<class Type>
668 (
669  const dimensionedScalar& rho,
671 )
672 {
673  tmp<faMatrix<Type>> tfam
674  (
675  new faMatrix<Type>
676  (
677  vf,
678  rho.dimensions()*vf.dimensions()*dimArea/dimTime
679  )
680  );
681  faMatrix<Type>& fam = tfam.ref();
682 
683  scalar rDeltaT = 1.0/deltaT_();
684 
685  scalar deltaT = deltaT_();
686  scalar deltaT0 = deltaT0_(vf);
687 
688  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
689  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
690  scalar coefft0 = coefft + coefft00;
691 
692  fam.diag() = (coefft*rDeltaT*rho.value())*mesh().S();
693 
694  if (mesh().moving())
695  {
696  fam.source() = rDeltaT*rho.value()*
697  (
698  coefft0*vf.oldTime().primitiveField()*mesh().S0()
699  - coefft00*vf.oldTime().oldTime().primitiveField()
700  *mesh().S00()
701  );
702  }
703  else
704  {
705  fam.source() = rDeltaT*mesh().S()*rho.value()*
706  (
707  coefft0*vf.oldTime().primitiveField()
708  - coefft00*vf.oldTime().oldTime().primitiveField()
709  );
710  }
711 
712  return tfam;
713 }
714 
715 
716 template<class Type>
719 (
720  const areaScalarField& rho,
722 )
723 {
724  tmp<faMatrix<Type>> tfam
725  (
726  new faMatrix<Type>
727  (
728  vf,
729  rho.dimensions()*vf.dimensions()*dimArea/dimTime
730  )
731  );
732  faMatrix<Type>& fam = tfam.ref();
733 
734  scalar rDeltaT = 1.0/deltaT_();
735 
736  scalar deltaT = deltaT_();
737  scalar deltaT0 = deltaT0_(vf);
738 
739  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
740  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
741  scalar coefft0 = coefft + coefft00;
742 
743  fam.diag() = (coefft*rDeltaT)*rho.primitiveField()*mesh().S();
744 
745  if (mesh().moving())
746  {
747  fam.source() = rDeltaT*
748  (
749  coefft0*rho.oldTime().primitiveField()
750  *vf.oldTime().primitiveField()*mesh().S0()
751  - coefft00*rho.oldTime().oldTime().primitiveField()
752  *vf.oldTime().oldTime().primitiveField()*mesh().S00()
753  );
754  }
755  else
756  {
757  fam.source() = rDeltaT*mesh().S()*
758  (
759  coefft0*rho.oldTime().primitiveField()
760  *vf.oldTime().primitiveField()
761  - coefft00*rho.oldTime().oldTime().primitiveField()
762  *vf.oldTime().oldTime().primitiveField()
763  );
764  }
765 
766  return tfam;
767 }
768 
769 
770 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
771 
772 } // End namespace fa
773 
774 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
775 
776 } // End namespace Foam
777 
778 // ************************************************************************* //
const Type & value() const noexcept
Return const reference to value.
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
tmp< GeometricField< Type, faPatchField, areaMesh > > facDdt0(const dimensioned< Type >)
tmp< faMatrix< Type > > famDdt(const GeometricField< Type, faPatchField, areaMesh > &)
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
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.
Generic GeometricField class.
Generic dimensioned Type class.
dynamicFvMesh & mesh
const dimensionSet & dimensions() const noexcept
Return const reference to dimensions.
static const word & calculatedType() noexcept
The type name for calculated patch fields.
Definition: faPatchField.H:174
tmp< GeometricField< Type, faPatchField, areaMesh > > facDdt(const dimensioned< Type >)
const word & name() const noexcept
Return const reference to name.
A special matrix type and solver, designed for finite area solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: faMatricesFwd.H:37
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
A class for managing temporary objects.
Definition: HashPtrTable.H:50
Defines the attributes of an object for which implicit objectRegistry management is supported...
Definition: IOobject.H:180
Calculate the matrix for the second temporal derivative.
Request registration (bool: true)
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