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-------------------------------------------------------------------------------
10License
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
34namespace Foam
35{
36
37// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38
39namespace fa
40{
41
42// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
43
44template<class Type>
45scalar backwardFaDdtScheme<Type>::deltaT_() const
46{
47 return mesh().time().deltaT().value();
48}
49
50
51template<class Type>
52scalar backwardFaDdtScheme<Type>::deltaT0_() const
53{
54 return mesh().time().deltaT0().value();
55}
56
57
58template<class Type>
59template<class GeoField>
60scalar 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
75template<class Type>
76tmp<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
130template<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
177template<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
250template<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
320template<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
393template<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
464template<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
541template<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
615template<class Type>
618(
620)
621{
623 (
625 (
626 vf,
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
666template<class Type>
669(
670 const dimensionedScalar& rho,
672)
673{
675 (
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
717template<class Type>
720(
721 const areaScalarField& rho,
723)
724{
726 (
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// ************************************************************************* //
const dimensionSet & dimensions() const
Return dimensions.
Generic GeometricField class.
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
const Internal & internalField() const
Return a const-reference to the dimensioned internal field.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
Author Zeljko Tukovic, FMENA Hrvoje Jasak, Wikki Ltd.
Generic dimensioned Type class.
const dimensionSet & dimensions() const
Return const reference to dimensions.
const Type & value() const
Return const reference to value.
const word & name() const
Return const reference to name.
A special matrix type and solver, designed for finite area solutions of scalar equations....
Definition: faMatrix.H:76
tmp< GeometricField< Type, faPatchField, areaMesh > > facDdt(const dimensioned< Type >)
tmp< faMatrix< Type > > famDdt(const GeometricField< Type, faPatchField, areaMesh > &)
tmp< GeometricField< Type, faPatchField, areaMesh > > facDdt0(const dimensioned< Type >)
A class for managing temporary objects.
Definition: tmp.H:65
T & ref() const
Definition: tmpI.H:227
dynamicFvMesh & mesh
Calculate the divergence of the given field.
word timeName
Definition: getTimeIndex.H:3
Namespace for OpenFOAM.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:59
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
Calculate the matrix for the second temporal derivative.