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 "facDiv.H"
30 #include "faMatrices.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 namespace fa
40 {
41 
42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
43 
44 template<class Type>
45 scalar backwardFaDdtScheme<Type>::deltaT_() const
46 {
47  return mesh().time().deltaT().value();
48 }
49 
50 
51 template<class Type>
52 scalar backwardFaDdtScheme<Type>::deltaT0_() const
53 {
54  return mesh().time().deltaT0().value();
55 }
56 
57 
58 template<class Type>
59 template<class GeoField>
60 scalar backwardFaDdtScheme<Type>::deltaT0_(const GeoField& vf) const
61 {
62  if (vf.oldTime().timeIndex() == vf.oldTime().oldTime().timeIndex())
63  {
64  return GREAT;
65  }
66  else
67  {
68  return deltaT0_();
69  }
70 }
71 
72 
73 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
74 
75 template<class Type>
76 tmp<GeometricField<Type, faPatchField, areaMesh>>
78 (
79  const dimensioned<Type> dt
80 )
81 {
82  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
83 
84  IOobject ddtIOobject
85  (
86  "ddt("+dt.name()+')',
87  mesh()().time().timeName(),
88  mesh()(),
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  {
103  (
105  (
106  ddtIOobject,
107  mesh(),
109  )
110  );
111 
112  tdtdt.ref().primitiveFieldRef() = rDeltaT.value()*dt.value()*
113  (
114  coefft - (coefft0*mesh().S0() - coefft00*mesh().S00())/mesh().S()
115  );
116 
117  return tdtdt;
118  }
119 
121  (
122  ddtIOobject,
123  mesh(),
126  );
127 }
128 
129 
130 template<class Type>
133 (
134  const dimensioned<Type> dt
135 )
136 {
137  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
138 
139  IOobject ddtIOobject
140  (
141  "ddt("+dt.name()+')',
142  mesh()().time().timeName(),
143  mesh()(),
146  );
147 
148  scalar deltaT = deltaT_();
149  scalar deltaT0 = deltaT0_();
150 
151  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
152  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
153  scalar coefft0 = coefft + coefft00;
154 
156  (
158  (
159  ddtIOobject,
160  mesh(),
161  -rDeltaT*(coefft0 - coefft00)*dt
162  )
163  );
164 
165  if (mesh().moving())
166  {
167  tdtdt0.ref().primitiveFieldRef() = (-rDeltaT.value()*dt.value())*
168  (
169  (coefft0*mesh().S0() - coefft00*mesh().S00())/mesh().S()
170  );
171  }
172 
173  return tdtdt0;
174 }
175 
176 
177 template<class Type>
180 (
182 )
183 {
184  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
185 
186  IOobject ddtIOobject
187  (
188  "ddt("+vf.name()+')',
189  mesh()().time().timeName(),
190  mesh()(),
193  );
194 
195  scalar deltaT = deltaT_();
196  scalar deltaT0 = deltaT0_(vf);
197 
198  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
199  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
200  scalar coefft0 = coefft + coefft00;
201 
202  if (mesh().moving())
203  {
205  (
207  (
208  ddtIOobject,
209  mesh(),
210  rDeltaT.dimensions()*vf.dimensions(),
211  rDeltaT.value()*
212  (
213  coefft*vf() -
214  (
215  coefft0*vf.oldTime()()*mesh().S0()
216  - coefft00*vf.oldTime().oldTime()()
217  *mesh().S00()
218  )/mesh().S()
219  ),
220  rDeltaT.value()*
221  (
222  coefft*vf.boundaryField() -
223  (
224  coefft0*vf.oldTime().boundaryField()
225  - coefft00*vf.oldTime().oldTime().boundaryField()
226  )
227  )
228  )
229  );
230  }
231  else
232  {
234  (
236  (
237  ddtIOobject,
238  rDeltaT*
239  (
240  coefft*vf
241  - coefft0*vf.oldTime()
242  + coefft00*vf.oldTime().oldTime()
243  )
244  )
245  );
246  }
247 }
248 
249 
250 template<class Type>
253 (
255 )
256 {
257  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
258 
259  IOobject ddtIOobject
260  (
261  "ddt0("+vf.name()+')',
262  mesh()().time().timeName(),
263  mesh()(),
266  );
267 
268  scalar deltaT = deltaT_();
269  scalar deltaT0 = deltaT0_(vf);
270 
271  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
272  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
273  scalar coefft0 = coefft + coefft00;
274 
275  if (mesh().moving())
276  {
278  (
280  (
281  ddtIOobject,
282  mesh(),
283  rDeltaT.dimensions()*vf.dimensions(),
284  rDeltaT.value()*
285  (
286  - (
287  coefft0*vf.oldTime()()*mesh().S0()
288  - coefft00*vf.oldTime().oldTime()()
289  *mesh().S00()
290  )/mesh().S()
291  ),
292  rDeltaT.value()*
293  (
294  - (
295  coefft0*vf.oldTime().boundaryField()
296  - coefft00*vf.oldTime().oldTime().boundaryField()
297  )
298  )
299  )
300  );
301  }
302  else
303  {
305  (
307  (
308  ddtIOobject,
309  rDeltaT*
310  (
311  - coefft0*vf.oldTime()
312  + coefft00*vf.oldTime().oldTime()
313  )
314  )
315  );
316  }
317 }
318 
319 
320 template<class Type>
323 (
324  const dimensionedScalar& rho,
326 )
327 {
328  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
329 
330  IOobject ddtIOobject
331  (
332  "ddt("+rho.name()+','+vf.name()+')',
333  mesh()().time().timeName(),
334  mesh()(),
337  );
338 
339  scalar deltaT = deltaT_();
340  scalar deltaT0 = deltaT0_(vf);
341 
342  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
343  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
344  scalar coefft0 = coefft + coefft00;
345 
346  if (mesh().moving())
347  {
349  (
351  (
352  ddtIOobject,
353  mesh(),
354  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
355  rDeltaT.value()*rho.value()*
356  (
357  coefft*vf.internalField() -
358  (
359  coefft0*vf.oldTime()()*mesh().S0()
360  - coefft00*vf.oldTime().oldTime()()
361  *mesh().S00()
362  )/mesh().S()
363  ),
364  rDeltaT.value()*rho.value()*
365  (
366  coefft*vf.boundaryField() -
367  (
368  coefft0*vf.oldTime().boundaryField()
369  - coefft00*vf.oldTime().oldTime().boundaryField()
370  )
371  )
372  )
373  );
374  }
375  else
376  {
378  (
380  (
381  ddtIOobject,
382  rDeltaT*rho*
383  (
384  coefft*vf
385  - coefft0*vf.oldTime()
386  + coefft00*vf.oldTime().oldTime()
387  )
388  )
389  );
390  }
391 }
392 
393 template<class Type>
396 (
397  const dimensionedScalar& rho,
399 )
400 {
401  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
402 
403  IOobject ddtIOobject
404  (
405  "ddt0("+rho.name()+','+vf.name()+')',
406  mesh()().time().timeName(),
407  mesh()(),
410  );
411 
412  scalar deltaT = deltaT_();
413  scalar deltaT0 = deltaT0_(vf);
414 
415  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
416  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
417  scalar coefft0 = coefft + coefft00;
418 
419  if (mesh().moving())
420  {
422  (
424  (
425  ddtIOobject,
426  mesh(),
427  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
428  rDeltaT.value()*rho.value()*
429  (
430  -(
431  coefft0*vf.oldTime()()*mesh().S0()
432  - coefft00*vf.oldTime().oldTime()()
433  *mesh().S00()
434  )/mesh().S()
435  ),
436  rDeltaT.value()*rho.value()*
437  (
438  -(
439  coefft0*vf.oldTime().boundaryField()
440  - coefft00*vf.oldTime().oldTime().boundaryField()
441  )
442  )
443  )
444  );
445  }
446  else
447  {
449  (
451  (
452  ddtIOobject,
453  rDeltaT*rho*
454  (
455  - coefft0*vf.oldTime()
456  + coefft00*vf.oldTime().oldTime()
457  )
458  )
459  );
460  }
461 }
462 
463 
464 template<class Type>
467 (
468  const areaScalarField& rho,
470 )
471 {
472  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
473 
474  IOobject ddtIOobject
475  (
476  "ddt("+rho.name()+','+vf.name()+')',
477  mesh()().time().timeName(),
478  mesh()(),
481  );
482 
483  scalar deltaT = deltaT_();
484  scalar deltaT0 = deltaT0_(vf);
485 
486  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
487  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
488  scalar coefft0 = coefft + coefft00;
489 
490  if (mesh().moving())
491  {
493  (
495  (
496  ddtIOobject,
497  mesh(),
498  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
499  rDeltaT.value()*
500  (
501  coefft*rho.internalField()*vf.internalField() -
502  (
503  coefft0*rho.oldTime()()
504  *vf.oldTime()()*mesh().S0()
505  - coefft00*rho.oldTime().oldTime()()
506  *vf.oldTime().oldTime()()*mesh().S00()
507  )/mesh().S()
508  ),
509  rDeltaT.value()*
510  (
511  coefft*vf.boundaryField() -
512  (
513  coefft0*rho.oldTime().boundaryField()
514  *vf.oldTime().boundaryField()
515  - coefft00*rho.oldTime().oldTime().boundaryField()
516  *vf.oldTime().oldTime().boundaryField()
517  )
518  )
519  )
520  );
521  }
522  else
523  {
525  (
527  (
528  ddtIOobject,
529  rDeltaT*
530  (
531  coefft*rho*vf
532  - coefft0*rho.oldTime()*vf.oldTime()
533  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
534  )
535  )
536  );
537  }
538 }
539 
540 
541 template<class Type>
544 (
545  const areaScalarField& rho,
547 )
548 {
549  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
550 
551  IOobject ddtIOobject
552  (
553  "ddt0("+rho.name()+','+vf.name()+')',
554  mesh()().time().timeName(),
555  mesh()(),
558  );
559 
560  scalar deltaT = deltaT_();
561  scalar deltaT0 = deltaT0_(vf);
562 
563  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
564  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
565  scalar coefft0 = coefft + coefft00;
566 
567  if (mesh().moving())
568  {
570  (
572  (
573  ddtIOobject,
574  mesh(),
575  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
576  rDeltaT.value()*
577  (
578  - (
579  coefft0*rho.oldTime()()
580  *vf.oldTime()()*mesh().S0()
581  - coefft00*rho.oldTime().oldTime()()
582  *vf.oldTime().oldTime()()*mesh().S00()
583  )/mesh().S()
584  ),
585  rDeltaT.value()*
586  (
587  - (
588  coefft0*rho.oldTime().boundaryField()
589  *vf.oldTime().boundaryField()
590  - coefft00*rho.oldTime().oldTime().boundaryField()
591  *vf.oldTime().oldTime().boundaryField()
592  )
593  )
594  )
595  );
596  }
597  else
598  {
600  (
602  (
603  ddtIOobject,
604  rDeltaT*
605  (
606  - coefft0*rho.oldTime()*vf.oldTime()
607  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
608  )
609  )
610  );
611  }
612 }
613 
614 
615 template<class Type>
618 (
620 )
621 {
622  tmp<faMatrix<Type>> tfam
623  (
624  new faMatrix<Type>
625  (
626  vf,
627  vf.dimensions()*dimArea/dimTime
628  )
629  );
630 
631  faMatrix<Type>& fam = tfam.ref();
632 
633  scalar rDeltaT = 1.0/deltaT_();
634 
635  scalar deltaT = deltaT_();
636  scalar deltaT0 = deltaT0_(vf);
637 
638  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
639  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
640  scalar coefft0 = coefft + coefft00;
641 
642  fam.diag() = (coefft*rDeltaT)*mesh().S();
643 
644  if (mesh().moving())
645  {
646  fam.source() = rDeltaT*
647  (
648  coefft0*vf.oldTime().primitiveField()*mesh().S0()
649  - coefft00*vf.oldTime().oldTime().primitiveField()
650  *mesh().S00()
651  );
652  }
653  else
654  {
655  fam.source() = rDeltaT*mesh().S()*
656  (
657  coefft0*vf.oldTime().primitiveField()
658  - coefft00*vf.oldTime().oldTime().primitiveField()
659  );
660  }
661 
662  return tfam;
663 }
664 
665 
666 template<class Type>
669 (
670  const dimensionedScalar& rho,
672 )
673 {
674  tmp<faMatrix<Type>> tfam
675  (
676  new faMatrix<Type>
677  (
678  vf,
679  rho.dimensions()*vf.dimensions()*dimArea/dimTime
680  )
681  );
682  faMatrix<Type>& fam = tfam.ref();
683 
684  scalar rDeltaT = 1.0/deltaT_();
685 
686  scalar deltaT = deltaT_();
687  scalar deltaT0 = deltaT0_(vf);
688 
689  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
690  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
691  scalar coefft0 = coefft + coefft00;
692 
693  fam.diag() = (coefft*rDeltaT*rho.value())*mesh().S();
694 
695  if (mesh().moving())
696  {
697  fam.source() = rDeltaT*rho.value()*
698  (
699  coefft0*vf.oldTime().primitiveField()*mesh().S0()
700  - coefft00*vf.oldTime().oldTime().primitiveField()
701  *mesh().S00()
702  );
703  }
704  else
705  {
706  fam.source() = rDeltaT*mesh().S()*rho.value()*
707  (
708  coefft0*vf.oldTime().primitiveField()
709  - coefft00*vf.oldTime().oldTime().primitiveField()
710  );
711  }
712 
713  return tfam;
714 }
715 
716 
717 template<class Type>
720 (
721  const areaScalarField& rho,
723 )
724 {
725  tmp<faMatrix<Type>> tfam
726  (
727  new faMatrix<Type>
728  (
729  vf,
730  rho.dimensions()*vf.dimensions()*dimArea/dimTime
731  )
732  );
733  faMatrix<Type>& fam = tfam.ref();
734 
735  scalar rDeltaT = 1.0/deltaT_();
736 
737  scalar deltaT = deltaT_();
738  scalar deltaT0 = deltaT0_(vf);
739 
740  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
741  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
742  scalar coefft0 = coefft + coefft00;
743 
744  fam.diag() = (coefft*rDeltaT)*rho.primitiveField()*mesh().S();
745 
746  if (mesh().moving())
747  {
748  fam.source() = rDeltaT*
749  (
750  coefft0*rho.oldTime().primitiveField()
751  *vf.oldTime().primitiveField()*mesh().S0()
752  - coefft00*rho.oldTime().oldTime().primitiveField()
753  *vf.oldTime().oldTime().primitiveField()*mesh().S00()
754  );
755  }
756  else
757  {
758  fam.source() = rDeltaT*mesh().S()*
759  (
760  coefft0*rho.oldTime().primitiveField()
761  *vf.oldTime().primitiveField()
762  - coefft00*rho.oldTime().oldTime().primitiveField()
763  *vf.oldTime().oldTime().primitiveField()
764  );
765  }
766 
767  return tfam;
768 }
769 
770 
771 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
772 
773 } // End namespace fa
774 
775 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
776 
777 } // End namespace Foam
778 
779 // ************************************************************************* //
Foam::IOobject::NO_WRITE
Definition: IOobject.H:195
Foam::fa::backwardFaDdtScheme::facDdt0
tmp< GeometricField< Type, faPatchField, areaMesh > > facDdt0(const dimensioned< Type >)
Definition: backwardFaDdtScheme.C:133
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
Foam::faMatrix
A special matrix type and solver, designed for finite area solutions of scalar equations....
Definition: faMatricesFwd.H:43
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::fa::backwardFaDdtScheme::famDdt
tmp< faMatrix< Type > > famDdt(const GeometricField< Type, faPatchField, areaMesh > &)
Definition: backwardFaDdtScheme.C:618
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:850
Foam::dimensioned::name
const word & name() const
Return const reference to name.
Definition: dimensionedType.C:406
faMatrices.H
facDiv.H
Calculate the divergence of the given field.
Foam::dimensioned::value
const Type & value() const
Return const reference to value.
Definition: dimensionedType.C:434
Foam::GeometricField::internalField
const Internal & internalField() const
Return a const-reference to the dimensioned internal field.
Definition: GeometricFieldI.H:43
rho
rho
Definition: readInitialConditions.H:88
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:59
Foam::calculatedFaPatchField
Author Zeljko Tukovic, FMENA Hrvoje Jasak, Wikki Ltd.
Definition: calculatedFaPatchField.H:54
Foam::fa::backwardFaDdtScheme::facDdt
tmp< GeometricField< Type, faPatchField, areaMesh > > facDdt(const dimensioned< Type >)
Definition: backwardFaDdtScheme.C:78
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedScalarFwd.H:42
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::New
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
Definition: DimensionedFieldReuseFunctions.H:105
fam
Calculate the matrix for the second temporal derivative.
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const
Return const reference to dimensions.
Definition: dimensionedType.C:420
Foam::GeometricField
Generic GeometricField class.
Definition: areaFieldsFwd.H:53
Foam::IOobject::NO_READ
Definition: IOobject.H:188
backwardFaDdtScheme.H
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:62