Loading...
Searching...
No Matches
CrankNicolsonDdtScheme.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) 2011-2018 OpenFOAM Foundation
9 Copyright (C) 2020-2024 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "CrankNicolsonDdtScheme.H"
30#include "surfaceInterpolate.H"
31#include "fvcDiv.H"
32#include "fvMatrices.H"
33#include "Constant.H"
34
35// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36
37namespace Foam
38{
39
40// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41
42namespace fv
43{
44
45// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
46
47template<class Type>
48template<class GeoField>
49CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::DDt0Field
50(
51 const IOobject& io,
52 const fvMesh& mesh
53)
54:
55 GeoField(io, mesh),
56 startTimeIndex_(-2) // This field is for a restart and thus correct so set
57 // the start time-index to correspond to a previous run
58{
59 // Set the time-index to the beginning of the run to ensure the field
60 // is updated during the first time-step
61 this->timeIndex() = mesh.time().startTimeIndex();
62}
63
64
65template<class Type>
66template<class GeoField>
67CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::DDt0Field
68(
69 const IOobject& io,
70 const fvMesh& mesh,
71 const typename GeoField::value_type& value,
72 const dimensionSet& dims
73)
74:
75 GeoField(io, mesh, value, dims),
76 startTimeIndex_(mesh.time().timeIndex())
77{}
78
79
80template<class Type>
81template<class GeoField>
82GeoField& CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::operator()()
83{
84 return *this;
85}
86
87
88template<class Type>
89template<class GeoField>
90void CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::operator=
91(
92 const GeoField& gf
93)
94{
95 GeoField::operator=(gf);
96}
97
98
99template<class Type>
100template<class GeoField>
101typename CrankNicolsonDdtScheme<Type>::template DDt0Field<GeoField>&
102CrankNicolsonDdtScheme<Type>::ddt0_
103(
104 const word& name,
105 const dimensionSet& dims
106)
107{
108 if (!mesh().objectRegistry::template foundObject<GeoField>(name))
109 {
110 const Time& runTime = mesh().time();
111 word startTimeName = runTime.timeName(runTime.startTime().value());
112
113 if
114 (
115 (
116 runTime.timeIndex() == runTime.startTimeIndex()
117 || runTime.timeIndex() == runTime.startTimeIndex() + 1
118 )
119 && IOobject
120 (
121 name,
122 startTimeName,
123 mesh().thisDb()
124 ).template typeHeaderOk<DDt0Field<GeoField>>(true)
125 )
126 {
127 regIOobject::store
128 (
129 new DDt0Field<GeoField>
130 (
131 IOobject
132 (
133 name,
134 startTimeName,
135 mesh().thisDb(),
136 IOobject::MUST_READ,
137 IOobject::AUTO_WRITE,
138 IOobject::REGISTER
139 ),
140 mesh()
141 )
142 );
143 }
144 else
145 {
146 regIOobject::store
147 (
148 new DDt0Field<GeoField>
149 (
150 IOobject
151 (
152 name,
153 mesh().time().timeName(),
154 mesh().thisDb(),
155 IOobject::NO_READ,
156 IOobject::AUTO_WRITE,
157 IOobject::REGISTER
158 ),
159 mesh(),
160 Foam::zero{}, // value
161 dims/dimTime
162 )
163 );
164 }
165 }
166
167 return static_cast<DDt0Field<GeoField>&>
168 (
169 mesh().objectRegistry::template lookupObjectRef<GeoField>(name)
170 );
171}
172
173
174template<class Type>
175template<class GeoField>
176bool CrankNicolsonDdtScheme<Type>::evaluate
177(
178 DDt0Field<GeoField>& ddt0
179)
180{
181 bool evaluated = (ddt0.timeIndex() != mesh().time().timeIndex());
182 ddt0.timeIndex() = mesh().time().timeIndex();
183 return evaluated;
184}
185
186
187template<class Type>
188template<class GeoField>
189scalar CrankNicolsonDdtScheme<Type>::coef_
190(
191 const DDt0Field<GeoField>& ddt0
192) const
193{
194 if (mesh().time().timeIndex() > ddt0.startTimeIndex())
195 {
196 return 1 + ocCoeff();
197 }
198 else
199 {
200 return 1;
201 }
202}
203
204
205template<class Type>
206template<class GeoField>
207scalar CrankNicolsonDdtScheme<Type>::coef0_
208(
209 const DDt0Field<GeoField>& ddt0
210) const
211{
212 if (mesh().time().timeIndex() > ddt0.startTimeIndex() + 1)
213 {
214 return 1 + ocCoeff();
215 }
216 else
217 {
218 return 1;
219 }
220}
221
222
223template<class Type>
224template<class GeoField>
225dimensionedScalar CrankNicolsonDdtScheme<Type>::rDtCoef_
226(
227 const DDt0Field<GeoField>& ddt0
228) const
229{
230 return coef_(ddt0)/mesh().time().deltaT();
231}
232
233
234template<class Type>
235template<class GeoField>
236dimensionedScalar CrankNicolsonDdtScheme<Type>::rDtCoef0_
237(
238 const DDt0Field<GeoField>& ddt0
239) const
240{
241 return coef0_(ddt0)/mesh().time().deltaT0();
242}
243
244
245template<class Type>
246template<class GeoField>
247tmp<GeoField> CrankNicolsonDdtScheme<Type>::offCentre_
248(
249 const GeoField& ddt0
250) const
251{
252 if (ocCoeff() < 1)
253 {
254 return ocCoeff()*ddt0;
255 }
256 else
257 {
258 return ddt0;
259 }
260}
261
262
263template<class Type>
264const FieldField<fvPatchField, Type>& ff
265(
266 const FieldField<fvPatchField, Type>& bf
267)
268{
269 return bf;
270}
271
272
273// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
274
275template<class Type>
276CrankNicolsonDdtScheme<Type>::CrankNicolsonDdtScheme(const fvMesh& mesh)
277:
278 ddtScheme<Type>(mesh),
279 ocCoeff_(new Function1Types::Constant<scalar>("ocCoeff", 1))
280{
281 // Ensure the old-old-time cell volumes are available
282 // for moving meshes
283 if (mesh.moving())
284 {
285 mesh.V00();
286 }
287}
288
289
290template<class Type>
291CrankNicolsonDdtScheme<Type>::CrankNicolsonDdtScheme
292(
293 const fvMesh& mesh,
294 Istream& is
295)
296:
297 ddtScheme<Type>(mesh, is)
298{
299 token firstToken(is);
300
301 if (firstToken.isNumber())
302 {
303 const scalar ocCoeff = firstToken.number();
304 if (ocCoeff < 0 || ocCoeff > 1)
305 {
306 FatalIOErrorInFunction(is)
307 << "Off-centreing coefficient = " << ocCoeff
308 << " should be >= 0 and <= 1"
309 << exit(FatalIOError);
310 }
311
312 ocCoeff_.reset
313 (
314 new Function1Types::Constant<scalar>("ocCoeff", ocCoeff)
315 );
316 }
317 else
318 {
319 is.putBack(firstToken);
320 dictionary dict(is);
321 ocCoeff_ = Function1<scalar>::New("ocCoeff", dict, &mesh);
322 }
323
324 // Ensure the old-old-time cell volumes are available
325 // for moving meshes
326 if (mesh.moving())
327 {
328 mesh.V00();
329 }
330}
331
332
333// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
334
335template<class Type>
336tmp<GeometricField<Type, fvPatchField, volMesh>>
337CrankNicolsonDdtScheme<Type>::fvcDdt
338(
339 const dimensioned<Type>& dt
340)
341{
342 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
343 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
344 (
345 "ddt0(" + dt.name() + ')',
346 dt.dimensions()
347 );
348
349 IOobject ddtIOobject
350 (
351 "ddt(" + dt.name() + ')',
352 mesh().time().timeName(),
353 mesh().thisDb()
354 );
355
356 tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
357 (
358 new GeometricField<Type, fvPatchField, volMesh>
359 (
360 ddtIOobject,
361 mesh(),
362 Foam::zero{}, // value
363 (dt.dimensions()/dimTime)
364 )
365 );
366
367 dimensionedScalar rDtCoef = rDtCoef_(ddt0);
368
369 if (mesh().moving())
370 {
371 if (evaluate(ddt0))
372 {
373 dimensionedScalar rDtCoef0 = rDtCoef0_(ddt0);
374
375 ddt0.internalFieldRef() =
376 (
377 (rDtCoef0*dt)*(mesh().V0() - mesh().V00())
378 - mesh().V00()*offCentre_(ddt0.internalField())
379 )/mesh().V0();
380 }
381
382 tdtdt.ref().internalFieldRef() =
383 (
384 (rDtCoef*dt)*(mesh().V() - mesh().V0())
385 - mesh().V0()*offCentre_(ddt0.internalField())
386 )/mesh().V();
387
388 // Different operation on boundary v.s. internal so re-evaluate
389 // coupled boundaries
390 tdtdt.ref().boundaryFieldRef().
391 template evaluateCoupled<coupledFvPatch>();
392 }
393
394 return tdtdt;
395}
396
397
398template<class Type>
399tmp<GeometricField<Type, fvPatchField, volMesh>>
400CrankNicolsonDdtScheme<Type>::fvcDdt
401(
402 const GeometricField<Type, fvPatchField, volMesh>& vf
403)
404{
405 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
406 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
407 (
408 "ddt0(" + vf.name() + ')',
409 vf.dimensions()
410 );
411
412 // Bypass DynamicField level to avoid problems with Field level calcs
413 const auto& ddt0p = static_cast<const Field<Type>&>(ddt0.primitiveField());
414
415 IOobject ddtIOobject
416 (
417 "ddt(" + vf.name() + ')',
418 mesh().time().timeName(),
419 mesh().thisDb()
420 );
421
422 dimensionedScalar rDtCoef = rDtCoef_(ddt0);
423
424 if (mesh().moving())
425 {
426 const scalarField& V0 = mesh().V0();
427 const scalarField& V00 = mesh().V00();
428
429 if (evaluate(ddt0))
430 {
431 const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
432
433 ddt0.primitiveFieldRef() =
434 (
435 rDtCoef0*
436 (
437 V0*vf.oldTime().primitiveField()
438 - V00*vf.oldTime().oldTime().primitiveField()
439 ) - V00*offCentre_(ddt0p)
440 )/V0;
441
442 ddt0.boundaryFieldRef() =
443 (
444 rDtCoef0*
445 (
446 vf.oldTime().boundaryField()
447 - vf.oldTime().oldTime().boundaryField()
448 ) - offCentre_(ff(ddt0.boundaryField()))
449 );
450 }
451
452 tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
453 (
454 new GeometricField<Type, fvPatchField, volMesh>
455 (
456 ddtIOobject,
457 (
458 rDtCoef*
459 (
460 mesh().V()*vf
461 - mesh().V0()*vf.oldTime()
462 ) - mesh().V0()*offCentre_(ddt0()())
463 )/mesh().V(),
464 rDtCoef.value()*
465 (
466 vf.boundaryField() - vf.oldTime().boundaryField()
467 ) - offCentre_(ff(ddt0.boundaryField()))
468 )
469 );
470
471 // Different operation on boundary v.s. internal so re-evaluate
472 // coupled boundaries
473 tdtdt.ref().boundaryFieldRef().
474 template evaluateCoupled<coupledFvPatch>();
475
476 return tdtdt;
477 }
478 else
479 {
480 if (evaluate(ddt0))
481 {
482 ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime())
483 - offCentre_(ddt0());
484 }
485
486 return tmp<GeometricField<Type, fvPatchField, volMesh>>
487 (
488 new GeometricField<Type, fvPatchField, volMesh>
489 (
490 ddtIOobject,
491 rDtCoef*(vf - vf.oldTime()) - offCentre_(ddt0())
492 )
493 );
494 }
495}
496
497
498template<class Type>
499tmp<GeometricField<Type, fvPatchField, volMesh>>
500CrankNicolsonDdtScheme<Type>::fvcDdt
501(
502 const dimensionedScalar& rho,
503 const GeometricField<Type, fvPatchField, volMesh>& vf
504)
505{
506 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
507 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
508 (
509 "ddt0(" + rho.name() + ',' + vf.name() + ')',
510 rho.dimensions()*vf.dimensions()
511 );
512
513 // Bypass DynamicField level to avoid problems with Field level calcs
514 const auto& ddt0p = static_cast<const Field<Type>&>(ddt0.primitiveField());
515
516 IOobject ddtIOobject
517 (
518 "ddt(" + rho.name() + ',' + vf.name() + ')',
519 mesh().time().timeName(),
520 mesh().thisDb()
521 );
522
523 dimensionedScalar rDtCoef = rDtCoef_(ddt0);
524
525 if (mesh().moving())
526 {
527 const scalarField& V0 = mesh().V0();
528 const scalarField& V00 = mesh().V00();
529
530 if (evaluate(ddt0))
531 {
532 const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
533
534 ddt0.primitiveFieldRef() =
535 (
536 rDtCoef0*rho.value()*
537 (
538 V0*vf.oldTime().primitiveField()
539 - V00*vf.oldTime().oldTime().primitiveField()
540 ) - V00*offCentre_(ddt0p)
541 )/V0;
542
543 ddt0.boundaryFieldRef() =
544 (
545 rDtCoef0*rho.value()*
546 (
547 vf.oldTime().boundaryField()
548 - vf.oldTime().oldTime().boundaryField()
549 ) - offCentre_(ff(ddt0.boundaryField()))
550 );
551 }
552
553 tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
554 (
555 new GeometricField<Type, fvPatchField, volMesh>
556 (
557 ddtIOobject,
558 mesh(),
559 rDtCoef.dimensions()*rho.dimensions()*vf.dimensions(),
560 (
561 rDtCoef.value()*rho.value()*
562 (
563 mesh().V()*vf.primitiveField()
564 - V0*vf.oldTime().primitiveField()
565 ) - V0*offCentre_(ddt0p)
566 )/mesh().V(),
567 rDtCoef.value()*rho.value()*
568 (
569 vf.boundaryField() - vf.oldTime().boundaryField()
570 ) - offCentre_(ff(ddt0.boundaryField()))
571 )
572 );
573
574 // Different operation on boundary v.s. internal so re-evaluate
575 // coupled boundaries
576 tdtdt.ref().boundaryFieldRef().
577 template evaluateCoupled<coupledFvPatch>();
578
579 return tdtdt;
580 }
581 else
582 {
583 if (evaluate(ddt0))
584 {
585 ddt0 = rDtCoef0_(ddt0)*rho*(vf.oldTime() - vf.oldTime().oldTime())
586 - offCentre_(ddt0());
587 }
588
589 return tmp<GeometricField<Type, fvPatchField, volMesh>>
590 (
591 new GeometricField<Type, fvPatchField, volMesh>
592 (
593 ddtIOobject,
594 rDtCoef*rho*(vf - vf.oldTime()) - offCentre_(ddt0())
595 )
596 );
597 }
598}
599
600
601template<class Type>
602tmp<GeometricField<Type, fvPatchField, volMesh>>
603CrankNicolsonDdtScheme<Type>::fvcDdt
604(
605 const volScalarField& rho,
606 const GeometricField<Type, fvPatchField, volMesh>& vf
607)
608{
609 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
610 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
611 (
612 "ddt0(" + rho.name() + ',' + vf.name() + ')',
613 rho.dimensions()*vf.dimensions()
614 );
615
616 // Bypass DynamicField level to avoid problems with Field level calcs
617 const auto& ddt0p = static_cast<const Field<Type>&>(ddt0.primitiveField());
618
619 IOobject ddtIOobject
620 (
621 "ddt(" + rho.name() + ',' + vf.name() + ')',
622 mesh().time().timeName(),
623 mesh().thisDb()
624 );
625
626 dimensionedScalar rDtCoef = rDtCoef_(ddt0);
627
628 if (mesh().moving())
629 {
630 const scalarField& V0 = mesh().V0();
631 const scalarField& V00 = mesh().V00();
632
633 if (evaluate(ddt0))
634 {
635 const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
636
637 ddt0.primitiveFieldRef() =
638 (
639 rDtCoef0*
640 (
641 V0*rho.oldTime().primitiveField()
642 *vf.oldTime().primitiveField()
643 - V00*rho.oldTime().oldTime().primitiveField()
644 *vf.oldTime().oldTime().primitiveField()
645 ) - V00*offCentre_(ddt0p)
646 )/V0;
647
648 ddt0.boundaryFieldRef() =
649 (
650 rDtCoef0*
651 (
652 rho.oldTime().boundaryField()
653 *vf.oldTime().boundaryField()
654 - rho.oldTime().oldTime().boundaryField()
655 *vf.oldTime().oldTime().boundaryField()
656 ) - offCentre_(ff(ddt0.boundaryField()))
657 );
658 }
659
660 tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
661 (
662 new GeometricField<Type, fvPatchField, volMesh>
663 (
664 ddtIOobject,
665 mesh(),
666 rDtCoef.dimensions()*rho.dimensions()*vf.dimensions(),
667 (
668 rDtCoef.value()*
669 (
670 mesh().V()*rho.primitiveField()*vf.primitiveField()
671 - V0*rho.oldTime().primitiveField()
672 *vf.oldTime().primitiveField()
673 ) - V00*offCentre_(ddt0p)
674 )/mesh().V(),
675 rDtCoef.value()*
676 (
677 rho.boundaryField()*vf.boundaryField()
678 - rho.oldTime().boundaryField()*vf.oldTime().boundaryField()
679 ) - offCentre_(ff(ddt0.boundaryField()))
680 )
681 );
682
683 // Different operation on boundary v.s. internal so re-evaluate
684 // coupled boundaries
685 tdtdt.ref().boundaryFieldRef().
686 template evaluateCoupled<coupledFvPatch>();
687
688 return tdtdt;
689 }
690 else
691 {
692 if (evaluate(ddt0))
693 {
694 ddt0 = rDtCoef0_(ddt0)*
695 (
696 rho.oldTime()*vf.oldTime()
697 - rho.oldTime().oldTime()*vf.oldTime().oldTime()
698 ) - offCentre_(ddt0());
699 }
700
701 return tmp<GeometricField<Type, fvPatchField, volMesh>>
702 (
703 new GeometricField<Type, fvPatchField, volMesh>
704 (
705 ddtIOobject,
706 rDtCoef*(rho*vf - rho.oldTime()*vf.oldTime())
707 - offCentre_(ddt0())
708 )
709 );
710 }
711}
712
713
714template<class Type>
715tmp<GeometricField<Type, fvPatchField, volMesh>>
716CrankNicolsonDdtScheme<Type>::fvcDdt
717(
718 const volScalarField& alpha,
719 const volScalarField& rho,
720 const GeometricField<Type, fvPatchField, volMesh>& vf
721)
722{
723 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
724 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
725 (
726 "ddt0(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
727 alpha.dimensions()*rho.dimensions()*vf.dimensions()
728 );
729
730 // Bypass DynamicField level to avoid problems with Field level calcs
731 const auto& ddt0p = static_cast<const Field<Type>&>(ddt0.primitiveField());
732
733 IOobject ddtIOobject
734 (
735 "ddt(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
736 mesh().time().timeName(),
737 mesh().thisDb()
738 );
739
740 dimensionedScalar rDtCoef = rDtCoef_(ddt0);
741
742 if (mesh().moving())
743 {
744 const scalarField& V0 = mesh().V0();
745 const scalarField& V00 = mesh().V00();
746
747 if (evaluate(ddt0))
748 {
749 const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
750
751 ddt0.primitiveFieldRef() =
752 (
753 rDtCoef0*
754 (
755 V0
756 *alpha.oldTime().primitiveField()
757 *rho.oldTime().primitiveField()
758 *vf.oldTime().primitiveField()
759
760 - V00
761 *alpha.oldTime().oldTime().primitiveField()
762 *rho.oldTime().oldTime().primitiveField()
763 *vf.oldTime().oldTime().primitiveField()
764 ) - V00*offCentre_(ddt0p)
765 )/V0;
766
767 ddt0.boundaryFieldRef() =
768 (
769 rDtCoef0*
770 (
771 alpha.oldTime().boundaryField()
772 *rho.oldTime().boundaryField()
773 *vf.oldTime().boundaryField()
774
775 - alpha.oldTime().oldTime().boundaryField()
776 *rho.oldTime().oldTime().boundaryField()
777 *vf.oldTime().oldTime().boundaryField()
778 ) - offCentre_(ff(ddt0.boundaryField()))
779 );
780 }
781
782 tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
783 (
784 new GeometricField<Type, fvPatchField, volMesh>
785 (
786 ddtIOobject,
787 mesh(),
788 rDtCoef.dimensions()
789 *alpha.dimensions()*rho.dimensions()*vf.dimensions(),
790 (
791 rDtCoef.value()*
792 (
793 mesh().V()
794 *alpha.primitiveField()
795 *rho.primitiveField()
796 *vf.primitiveField()
797
798 - V0
799 *alpha.oldTime().primitiveField()
800 *rho.oldTime().primitiveField()
801 *vf.oldTime().primitiveField()
802 ) - V00*offCentre_(ddt0p)
803 )/mesh().V(),
804 rDtCoef.value()*
805 (
806 alpha.boundaryField()
807 *rho.boundaryField()
808 *vf.boundaryField()
809
810 - alpha.oldTime().boundaryField()
811 *rho.oldTime().boundaryField()
812 *vf.oldTime().boundaryField()
813 ) - offCentre_(ff(ddt0.boundaryField()))
814 )
815 );
816
817 // Different operation on boundary v.s. internal so re-evaluate
818 // coupled boundaries
819 tdtdt.ref().boundaryFieldRef().
820 template evaluateCoupled<coupledFvPatch>();
821
822 return tdtdt;
823 }
824 else
825 {
826 if (evaluate(ddt0))
827 {
828 ddt0 = rDtCoef0_(ddt0)*
829 (
830 alpha.oldTime()
831 *rho.oldTime()
832 *vf.oldTime()
833
834 - alpha.oldTime().oldTime()
835 *rho.oldTime().oldTime()
836 *vf.oldTime().oldTime()
837 ) - offCentre_(ddt0());
838 }
839
840 return tmp<GeometricField<Type, fvPatchField, volMesh>>
841 (
842 new GeometricField<Type, fvPatchField, volMesh>
843 (
844 ddtIOobject,
845 rDtCoef
846 *(
847 alpha*rho*vf
848 - alpha.oldTime()*rho.oldTime()*vf.oldTime()
849 )
850 - offCentre_(ddt0())
851 )
852 );
853 }
854}
855
856
857template<class Type>
858tmp<fvMatrix<Type>>
859CrankNicolsonDdtScheme<Type>::fvmDdt
860(
861 const GeometricField<Type, fvPatchField, volMesh>& vf
862)
863{
864 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
865 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
866 (
867 "ddt0(" + vf.name() + ')',
868 vf.dimensions()
869 );
870
871 // Bypass DynamicField level to avoid problems with Field level calcs
872 const auto& ddt0p = static_cast<const Field<Type>&>(ddt0.primitiveField());
873
874 tmp<fvMatrix<Type>> tfvm
875 (
876 new fvMatrix<Type>
877 (
878 vf,
879 vf.dimensions()*dimVol/dimTime
880 )
881 );
882
883 fvMatrix<Type>& fvm = tfvm.ref();
884
885 const scalar rDtCoef = rDtCoef_(ddt0).value();
886 fvm.diag() = rDtCoef*mesh().V();
887
888 vf.oldTime().oldTime();
889
890 if (mesh().moving())
891 {
892 const scalarField& V0 = mesh().V0();
893 const scalarField& V00 = mesh().V00();
894
895 if (evaluate(ddt0))
896 {
897 const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
898
899 ddt0.primitiveFieldRef() =
900 (
901 rDtCoef0*
902 (
903 V0*vf.oldTime().primitiveField()
904 - V00*vf.oldTime().oldTime().primitiveField()
905 )
906 - V00*offCentre_(ddt0p)
907 )/V0;
908
909 ddt0.boundaryFieldRef() =
910 (
911 rDtCoef0*
912 (
913 vf.oldTime().boundaryField()
914 - vf.oldTime().oldTime().boundaryField()
915 )
916 - offCentre_(ff(ddt0.boundaryField()))
917 );
918 }
919
920 fvm.source() =
921 (
922 rDtCoef*vf.oldTime().primitiveField()
923 + offCentre_(ddt0p)
924 )*V0;
925 }
926 else
927 {
928 if (evaluate(ddt0))
929 {
930 ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime())
931 - offCentre_(ddt0());
932
933 }
934
935 fvm.source() =
936 (
937 rDtCoef*vf.oldTime().primitiveField()
938 + offCentre_(ddt0p)
939 )*mesh().V();
940 }
941
942 return tfvm;
943}
944
945
946template<class Type>
947tmp<fvMatrix<Type>>
948CrankNicolsonDdtScheme<Type>::fvmDdt
949(
950 const dimensionedScalar& rho,
951 const GeometricField<Type, fvPatchField, volMesh>& vf
952)
953{
954 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
955 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
956 (
957 "ddt0(" + rho.name() + ',' + vf.name() + ')',
958 rho.dimensions()*vf.dimensions()
959 );
960
961 // Bypass DynamicField level to avoid problems with Field level calcs
962 const auto& ddt0p = static_cast<const Field<Type>&>(ddt0.primitiveField());
963
964 tmp<fvMatrix<Type>> tfvm
965 (
966 new fvMatrix<Type>
967 (
968 vf,
969 rho.dimensions()*vf.dimensions()*dimVol/dimTime
970 )
971 );
972 fvMatrix<Type>& fvm = tfvm.ref();
973
974 const scalar rDtCoef = rDtCoef_(ddt0).value();
975 fvm.diag() = rDtCoef*rho.value()*mesh().V();
976
977 vf.oldTime().oldTime();
978
979 if (mesh().moving())
980 {
981 const scalarField& V0 = mesh().V0();
982 const scalarField& V00 = mesh().V00();
983
984 if (evaluate(ddt0))
985 {
986 const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
987
988 ddt0.primitiveFieldRef() =
989 (
990 rDtCoef0*rho.value()*
991 (
992 V0*vf.oldTime().primitiveField()
993 - V00*vf.oldTime().oldTime().primitiveField()
994 )
995 - V00*offCentre_(ddt0p)
996 )/V0;
997
998 ddt0.boundaryFieldRef() =
999 (
1000 rDtCoef0*rho.value()*
1001 (
1002 vf.oldTime().boundaryField()
1003 - vf.oldTime().oldTime().boundaryField()
1004 )
1005 - offCentre_(ff(ddt0.boundaryField()))
1006 );
1007 }
1008
1009 fvm.source() =
1010 (
1011 rDtCoef*rho.value()*vf.oldTime().primitiveField()
1012 + offCentre_(ddt0p)
1013 )*V0;
1014 }
1015 else
1016 {
1017 if (evaluate(ddt0))
1018 {
1019 ddt0 = rDtCoef0_(ddt0)*rho*(vf.oldTime() - vf.oldTime().oldTime())
1020 - offCentre_(ddt0());
1021 }
1022
1023 fvm.source() =
1024 (
1025 rDtCoef*rho.value()*vf.oldTime().primitiveField()
1026 + offCentre_(ddt0p)
1027 )*mesh().V();
1028 }
1029
1030 return tfvm;
1031}
1032
1033
1034template<class Type>
1035tmp<fvMatrix<Type>>
1036CrankNicolsonDdtScheme<Type>::fvmDdt
1037(
1038 const volScalarField& rho,
1039 const GeometricField<Type, fvPatchField, volMesh>& vf
1040)
1041{
1042 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1043 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1044 (
1045 "ddt0(" + rho.name() + ',' + vf.name() + ')',
1046 rho.dimensions()*vf.dimensions()
1047 );
1048
1049 // Bypass DynamicField level to avoid problems with Field level calcs
1050 const auto& ddt0p = static_cast<const Field<Type>&>(ddt0.primitiveField());
1051
1052 tmp<fvMatrix<Type>> tfvm
1053 (
1054 new fvMatrix<Type>
1055 (
1056 vf,
1057 rho.dimensions()*vf.dimensions()*dimVol/dimTime
1058 )
1059 );
1060 fvMatrix<Type>& fvm = tfvm.ref();
1061
1062 const scalar rDtCoef = rDtCoef_(ddt0).value();
1063 fvm.diag() = rDtCoef*rho.primitiveField()*mesh().V();
1064
1065 vf.oldTime().oldTime();
1066 rho.oldTime().oldTime();
1067
1068 if (mesh().moving())
1069 {
1070 const scalarField& V0 = mesh().V0();
1071 const scalarField& V00 = mesh().V00();
1072
1073 if (evaluate(ddt0))
1074 {
1075 const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
1076
1077 ddt0.primitiveFieldRef() =
1078 (
1079 rDtCoef0*
1080 (
1081 V0*rho.oldTime().primitiveField()
1082 *vf.oldTime().primitiveField()
1083 - V00*rho.oldTime().oldTime().primitiveField()
1084 *vf.oldTime().oldTime().primitiveField()
1085 )
1086 - V00*offCentre_(ddt0p)
1087 )/V0;
1088
1089 ddt0.boundaryFieldRef() =
1090 (
1091 rDtCoef0*
1092 (
1093 rho.oldTime().boundaryField()
1094 *vf.oldTime().boundaryField()
1095 - rho.oldTime().oldTime().boundaryField()
1096 *vf.oldTime().oldTime().boundaryField()
1097 )
1098 - offCentre_(ff(ddt0.boundaryField()))
1099 );
1100 }
1101
1102 fvm.source() =
1103 (
1104 rDtCoef*rho.oldTime().primitiveField()*vf.oldTime().primitiveField()
1105 + offCentre_(ddt0p)
1106 )*V0;
1107 }
1108 else
1109 {
1110 if (evaluate(ddt0))
1111 {
1112 ddt0 = rDtCoef0_(ddt0)*
1113 (
1114 rho.oldTime()*vf.oldTime()
1115 - rho.oldTime().oldTime()*vf.oldTime().oldTime()
1116 ) - offCentre_(ddt0());
1117 }
1118
1119 fvm.source() =
1120 (
1121 rDtCoef*rho.oldTime().primitiveField()*vf.oldTime().primitiveField()
1122 + offCentre_(ddt0p)
1123 )*mesh().V();
1124 }
1125
1126 return tfvm;
1127}
1128
1129
1130template<class Type>
1131tmp<fvMatrix<Type>>
1132CrankNicolsonDdtScheme<Type>::fvmDdt
1133(
1134 const volScalarField& alpha,
1135 const volScalarField& rho,
1136 const GeometricField<Type, fvPatchField, volMesh>& vf
1137)
1138{
1139 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1140 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1141 (
1142 "ddt0(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
1143 alpha.dimensions()*rho.dimensions()*vf.dimensions()
1144 );
1145
1146 // Bypass DynamicField level to avoid problems with Field level calcs
1147 const auto& ddt0p = static_cast<const Field<Type>&>(ddt0.primitiveField());
1148
1149 tmp<fvMatrix<Type>> tfvm
1150 (
1151 new fvMatrix<Type>
1152 (
1153 vf,
1154 alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
1155 )
1156 );
1157 fvMatrix<Type>& fvm = tfvm.ref();
1158
1159 const scalar rDtCoef = rDtCoef_(ddt0).value();
1160 fvm.diag() = rDtCoef*alpha.primitiveField()*rho.primitiveField()*mesh().V();
1161
1162 vf.oldTime().oldTime();
1163 alpha.oldTime().oldTime();
1164 rho.oldTime().oldTime();
1165
1166 if (mesh().moving())
1167 {
1168 const scalarField& V0 = mesh().V0();
1169 const scalarField& V00 = mesh().V00();
1170
1171 if (evaluate(ddt0))
1172 {
1173 const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
1174
1175 ddt0.primitiveFieldRef() =
1176 (
1177 rDtCoef0*
1178 (
1179 V0
1180 *alpha.oldTime().primitiveField()
1181 *rho.oldTime().primitiveField()
1182 *vf.oldTime().primitiveField()
1183
1184 - V00
1185 *alpha.oldTime().oldTime().primitiveField()
1186 *rho.oldTime().oldTime().primitiveField()
1187 *vf.oldTime().oldTime().primitiveField()
1188 )
1189 - V00*offCentre_(ddt0p)
1190 )/V0;
1191
1192 ddt0.boundaryFieldRef() =
1193 (
1194 rDtCoef0*
1195 (
1196 alpha.oldTime().boundaryField()
1197 *rho.oldTime().boundaryField()
1198 *vf.oldTime().boundaryField()
1199
1200 - alpha.oldTime().oldTime().boundaryField()
1201 *rho.oldTime().oldTime().boundaryField()
1202 *vf.oldTime().oldTime().boundaryField()
1203 )
1204 - offCentre_(ff(ddt0.boundaryField()))
1205 );
1206 }
1207
1208 fvm.source() =
1209 (
1210 rDtCoef
1211 *alpha.oldTime().primitiveField()
1212 *rho.oldTime().primitiveField()
1213 *vf.oldTime().primitiveField()
1214 + offCentre_(ddt0p)
1215 )*V0;
1216 }
1217 else
1218 {
1219 if (evaluate(ddt0))
1220 {
1221 ddt0 = rDtCoef0_(ddt0)*
1222 (
1223 alpha.oldTime()
1224 *rho.oldTime()
1225 *vf.oldTime()
1226
1227 - alpha.oldTime().oldTime()
1228 *rho.oldTime().oldTime()
1229 *vf.oldTime().oldTime()
1230 ) - offCentre_(ddt0());
1231 }
1232
1233 fvm.source() =
1234 (
1235 rDtCoef
1236 *alpha.oldTime().primitiveField()
1237 *rho.oldTime().primitiveField()
1238 *vf.oldTime().primitiveField()
1239 + offCentre_(ddt0p)
1240 )*mesh().V();
1241 }
1242
1243 return tfvm;
1244}
1245
1246
1247template<class Type>
1248tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1249CrankNicolsonDdtScheme<Type>::fvcDdtUfCorr
1250(
1251 const GeometricField<Type, fvPatchField, volMesh>& U,
1252 const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
1253)
1254{
1255 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1256 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1257 (
1258 "ddtCorrDdt0(" + U.name() + ')',
1259 U.dimensions()
1260 );
1261
1262 DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 =
1263 ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>>
1264 (
1265 "ddtCorrDdt0(" + Uf.name() + ')',
1266 Uf.dimensions()
1267 );
1268
1269 dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1270
1271 if (evaluate(ddt0))
1272 {
1273 ddt0 =
1274 rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1275 - offCentre_(ddt0());
1276 }
1277
1278 if (evaluate(dUfdt0))
1279 {
1280 dUfdt0 =
1281 rDtCoef0_(dUfdt0)*(Uf.oldTime() - Uf.oldTime().oldTime())
1282 - offCentre_(dUfdt0());
1283 }
1284
1285 return tmp<fluxFieldType>
1286 (
1287 new fluxFieldType
1288 (
1289 IOobject
1290 (
1291 "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
1292 mesh().time().timeName(),
1293 mesh().thisDb()
1294 ),
1295 this->fvcDdtPhiCoeff(U.oldTime(), mesh().Sf() & Uf.oldTime())
1296 *(
1297 mesh().Sf()
1298 & (
1299 (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1300 - fvc::interpolate(rDtCoef*U.oldTime() + offCentre_(ddt0()))
1301 )
1302 )
1303 )
1304 );
1305}
1306
1307
1308template<class Type>
1309tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1310CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
1311(
1312 const GeometricField<Type, fvPatchField, volMesh>& U,
1313 const fluxFieldType& phi
1314)
1315{
1316 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1317 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1318 (
1319 "ddtCorrDdt0(" + U.name() + ')',
1320 U.dimensions()
1321 );
1322
1323 DDt0Field<fluxFieldType>& dphidt0 =
1324 ddt0_<fluxFieldType>
1325 (
1326 "ddtCorrDdt0(" + phi.name() + ')',
1327 phi.dimensions()
1328 );
1329 dphidt0.setOriented();
1330
1331 dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1332
1333 if (evaluate(ddt0))
1334 {
1335 ddt0 =
1336 rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1337 - offCentre_(ddt0());
1338 }
1339
1340 if (evaluate(dphidt0))
1341 {
1342 dphidt0 =
1343 rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime())
1344 - offCentre_(dphidt0());
1345 }
1346
1347 return tmp<fluxFieldType>
1348 (
1349 new fluxFieldType
1350 (
1351 IOobject
1352 (
1353 "ddtCorr(" + U.name() + ',' + phi.name() + ')',
1354 mesh().time().timeName(),
1355 mesh().thisDb()
1356 ),
1357 this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
1358 *(
1359 (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1360 - fvc::dotInterpolate
1361 (
1362 mesh().Sf(),
1363 rDtCoef*U.oldTime() + offCentre_(ddt0())
1364 )
1365 )
1366 )
1367 );
1368}
1369
1370
1371template<class Type>
1372tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1373CrankNicolsonDdtScheme<Type>::fvcDdtUfCorr
1374(
1375 const volScalarField& rho,
1376 const GeometricField<Type, fvPatchField, volMesh>& U,
1377 const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
1378)
1379{
1380 if
1381 (
1382 U.dimensions() == dimVelocity
1383 && Uf.dimensions() == rho.dimensions()*dimVelocity
1384 )
1385 {
1386 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1387 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1388 (
1389 "ddtCorrDdt0(" + rho.name() + ',' + U.name() + ')',
1390 rho.dimensions()*U.dimensions()
1391 );
1392
1393 DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 =
1394 ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>>
1395 (
1396 "ddtCorrDdt0(" + Uf.name() + ')',
1397 Uf.dimensions()
1398 );
1399
1400 dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1401
1402 GeometricField<Type, fvPatchField, volMesh> rhoU0
1403 (
1404 rho.oldTime()*U.oldTime()
1405 );
1406
1407 if (evaluate(ddt0))
1408 {
1409 ddt0 =
1410 rDtCoef0_(ddt0)
1411 *(rhoU0 - rho.oldTime().oldTime()*U.oldTime().oldTime())
1412 - offCentre_(ddt0());
1413 }
1414
1415 if (evaluate(dUfdt0))
1416 {
1417 dUfdt0 =
1418 rDtCoef0_(dUfdt0)
1419 *(Uf.oldTime() - Uf.oldTime().oldTime())
1420 - offCentre_(dUfdt0());
1421 }
1422
1423 tmp<fluxFieldType> ddtCorr
1424 (
1425 new fluxFieldType
1426 (
1427 IOobject
1428 (
1429 "ddtCorr("
1430 + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
1431 mesh().time().timeName(),
1432 mesh().thisDb()
1433 ),
1434 this->fvcDdtPhiCoeff
1435 (
1436 rhoU0,
1437 mesh().Sf() & Uf.oldTime(),
1438 rho.oldTime()
1439 )
1440 *(
1441 mesh().Sf()
1442 & (
1443 (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1444 - fvc::interpolate(rDtCoef*rhoU0 + offCentre_(ddt0()))
1445 )
1446 )
1447 )
1448 );
1449
1450 return ddtCorr;
1451 }
1452 else if
1453 (
1454 U.dimensions() == rho.dimensions()*dimVelocity
1455 && Uf.dimensions() == rho.dimensions()*dimVelocity
1456 )
1457 {
1458 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1459 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1460 (
1461 "ddtCorrDdt0(" + U.name() + ')',
1462 U.dimensions()
1463 );
1464
1465 DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 =
1466 ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>>
1467 (
1468 "ddtCorrDdt0(" + Uf.name() + ')',
1469 Uf.dimensions()
1470 );
1471
1472 dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1473
1474 if (evaluate(ddt0))
1475 {
1476 ddt0 =
1477 rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1478 - offCentre_(ddt0());
1479 }
1480
1481 if (evaluate(dUfdt0))
1482 {
1483 dUfdt0 =
1484 rDtCoef0_(dUfdt0)*(Uf.oldTime() - Uf.oldTime().oldTime())
1485 - offCentre_(dUfdt0());
1486 }
1487
1488 return tmp<fluxFieldType>
1489 (
1490 new fluxFieldType
1491 (
1492 IOobject
1493 (
1494 "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
1495 mesh().time().timeName(),
1496 mesh().thisDb()
1497 ),
1498 this->fvcDdtPhiCoeff
1499 (
1500 U.oldTime(),
1501 mesh().Sf() & Uf.oldTime(),
1502 rho.oldTime()
1503 )
1504 *(
1505 mesh().Sf()
1506 & (
1507 (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1508 - fvc::interpolate
1509 (
1510 rDtCoef*U.oldTime() + offCentre_(ddt0())
1511 )
1512 )
1513 )
1514 )
1515 );
1516 }
1517 else
1518 {
1519 FatalErrorInFunction
1520 << "dimensions of Uf are not correct"
1521 << abort(FatalError);
1522
1523 return fluxFieldType::null();
1524 }
1525}
1526
1527
1528template<class Type>
1529tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1530CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
1531(
1532 const volScalarField& rho,
1533 const GeometricField<Type, fvPatchField, volMesh>& U,
1534 const fluxFieldType& phi
1535)
1536{
1537 if
1538 (
1539 U.dimensions() == dimVelocity
1540 && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1541 )
1542 {
1543 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1544 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1545 (
1546 "ddtCorrDdt0(" + rho.name() + ',' + U.name() + ')',
1547 rho.dimensions()*U.dimensions()
1548 );
1549
1550 DDt0Field<fluxFieldType>& dphidt0 =
1551 ddt0_<fluxFieldType>
1552 (
1553 "ddtCorrDdt0(" + phi.name() + ')',
1554 phi.dimensions()
1555 );
1556
1557 dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1558
1559 GeometricField<Type, fvPatchField, volMesh> rhoU0
1560 (
1561 rho.oldTime()*U.oldTime()
1562 );
1563
1564 if (evaluate(ddt0))
1565 {
1566 ddt0 =
1567 rDtCoef0_(ddt0)
1568 *(rhoU0 - rho.oldTime().oldTime()*U.oldTime().oldTime())
1569 - offCentre_(ddt0());
1570 }
1571
1572 if (evaluate(dphidt0))
1573 {
1574 dphidt0 =
1575 rDtCoef0_(dphidt0)
1576 *(phi.oldTime() - phi.oldTime().oldTime())
1577 - offCentre_(dphidt0());
1578 }
1579
1580 tmp<fluxFieldType> ddtCorr
1581 (
1582 new fluxFieldType
1583 (
1584 IOobject
1585 (
1586 "ddtCorr("
1587 + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
1588 mesh().time().timeName(),
1589 mesh().thisDb()
1590 ),
1591 this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), rho.oldTime())
1592 *(
1593 (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1594 - fvc::dotInterpolate
1595 (
1596 mesh().Sf(),
1597 rDtCoef*rhoU0 + offCentre_(ddt0())
1598 )
1599 )
1600 )
1601 );
1602
1603 return ddtCorr;
1604 }
1605 else if
1606 (
1607 U.dimensions() == rho.dimensions()*dimVelocity
1608 && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1609 )
1610 {
1611 DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1612 ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1613 (
1614 "ddtCorrDdt0(" + U.name() + ')',
1615 U.dimensions()
1616 );
1617
1618 DDt0Field<fluxFieldType>& dphidt0 =
1619 ddt0_<fluxFieldType>
1620 (
1621 "ddtCorrDdt0(" + phi.name() + ')',
1622 phi.dimensions()
1623 );
1624
1625 dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1626
1627 if (evaluate(ddt0))
1628 {
1629 ddt0 =
1630 rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1631 - offCentre_(ddt0());
1632 }
1633
1634 if (evaluate(dphidt0))
1635 {
1636 dphidt0 =
1637 rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime())
1638 - offCentre_(dphidt0());
1639 }
1640
1641 return tmp<fluxFieldType>
1642 (
1643 new fluxFieldType
1644 (
1645 IOobject
1646 (
1647 "ddtCorr(" + U.name() + ',' + phi.name() + ')',
1648 mesh().time().timeName(),
1649 mesh().thisDb()
1650 ),
1651 this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), rho.oldTime())
1652 *(
1653 (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1654 - fvc::dotInterpolate
1655 (
1656 mesh().Sf(),
1657 rDtCoef*U.oldTime() + offCentre_(ddt0())
1658 )
1659 )
1660 )
1661 );
1662 }
1663 else
1664 {
1665 FatalErrorInFunction
1666 << "dimensions of phi are not correct"
1667 << abort(FatalError);
1668
1669 return fluxFieldType::null();
1670 }
1671}
1672
1673
1674template<class Type>
1675tmp<surfaceScalarField> CrankNicolsonDdtScheme<Type>::meshPhi
1676(
1677 const GeometricField<Type, fvPatchField, volMesh>& vf
1678)
1679{
1680 DDt0Field<surfaceScalarField>& meshPhi0 = ddt0_<surfaceScalarField>
1681 (
1682 "meshPhiCN_0",
1683 dimVolume
1684 );
1685
1686 meshPhi0.setOriented();
1687
1688 if (evaluate(meshPhi0))
1689 {
1690 meshPhi0 =
1691 coef0_(meshPhi0)*mesh().phi().oldTime() - offCentre_(meshPhi0());
1692 }
1693
1694 return tmp<surfaceScalarField>
1695 (
1696 new surfaceScalarField
1697 (
1698 IOobject
1699 (
1700 mesh().phi().name(),
1701 mesh().time().timeName(),
1702 mesh().thisDb(),
1703 IOobject::NO_READ,
1704 IOobject::NO_WRITE,
1705 IOobject::NO_REGISTER
1706 ),
1707 coef_(meshPhi0)*mesh().phi() - offCentre_(meshPhi0())
1708 )
1709 );
1710}
1711
1712
1713// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1714
1715} // End namespace fv
1716
1717// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1718
1719} // End namespace Foam
1720
1721// ************************************************************************* //
ocCoeff
scalar ocCoeff
Definition alphaEqn.H:6
Foam::DimensionedField::dimensions
const dimensionSet & dimensions() const noexcept
Return dimensions.
Definition DimensionedField.H:693
Foam::FieldField
A field of fields is a PtrList of fields with reference counting.
Definition FieldField.H:77
Foam::Field
Generic templated field type that is much like a Foam::List except that it is expected to hold numeri...
Definition Field.H:172
Foam::Function1Types::Constant
Templated function that returns a constant value.
Definition Constant.H:71
Foam::GeometricField
Generic GeometricField class.
Definition GeometricField.H:75
Foam::GeometricField::null
static const this_type & null() noexcept
Return a null GeometricField (reference to a nullObject).
Definition GeometricField.H:193
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition GeometricField.C:961
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const noexcept
Return a const-reference to the internal field values.
Definition GeometricFieldI.H:36
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const noexcept
Return const-reference to the boundary field.
Definition GeometricFieldI.H:45
Foam::IOobjectOption::NO_REGISTER
@ NO_REGISTER
Do not request registration (bool: false).
Definition IOobjectOption.H:99
Foam::IOobjectOption::REGISTER
@ REGISTER
Request registration (bool: true).
Definition IOobjectOption.H:102
Foam::IOobjectOption::NO_READ
@ NO_READ
Nothing to be read.
Definition IOobjectOption.H:60
Foam::IOobjectOption::MUST_READ
@ MUST_READ
Reading required.
Definition IOobjectOption.H:63
Foam::IOobjectOption::NO_WRITE
@ NO_WRITE
Ignore writing from objectRegistry::writeObject().
Definition IOobjectOption.H:86
Foam::IOobjectOption::AUTO_WRITE
@ AUTO_WRITE
Automatically write from objectRegistry::writeObject().
Definition IOobjectOption.H:89
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition IOobject.H:191
Foam::IOobject::name
const word & name() const noexcept
Return the object name.
Definition IOobjectI.H:205
Foam::Istream
An Istream is an abstract base class for all input systems (streams, files, token lists etc)....
Definition Istream.H:60
Foam::Istream::putBack
void putBack(const token &tok)
Put back a token (copy). Only a single put back is permitted.
Definition Istream.C:71
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition Time.H:75
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition dictionary.H:133
Foam::dimensionSet
Dimension set for the base types, which can be used to implement rigorous dimension checking for alge...
Definition dimensionSet.H:106
Foam::dimensioned
Generic dimensioned Type class.
Definition dimensionedType.H:70
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const noexcept
Return const reference to dimensions.
Definition dimensionedType.H:362
Foam::dimensioned::name
const word & name() const noexcept
Return const reference to name.
Definition dimensionedType.H:352
Foam::dimensioned::value
const Type & value() const noexcept
Return const reference to value.
Definition dimensionedType.H:372
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() noexcept
Definition fvMatrix.H:535
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition fvMesh.H:85
Foam::fv::CrankNicolsonDdtScheme
Second-oder Crank-Nicolson implicit ddt using the current and previous time-step fields as well as th...
Definition CrankNicolsonDdtScheme.H:115
Foam::fv::CrankNicolsonDdtScheme::fvcDdt
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
Definition CrankNicolsonDdtScheme.C:331
Foam::fv::CrankNicolsonDdtScheme::fvcDdtUfCorr
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition CrankNicolsonDdtScheme.C:1243
Foam::fv::CrankNicolsonDdtScheme::ocCoeff
scalar ocCoeff() const
Return the current off-centreing coefficient.
Definition CrankNicolsonDdtScheme.H:276
Foam::fv::CrankNicolsonDdtScheme::fvmDdt
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
Definition CrankNicolsonDdtScheme.C:853
Foam::fv::CrankNicolsonDdtScheme::fluxFieldType
ddtScheme< Type >::fluxFieldType fluxFieldType
Definition CrankNicolsonDdtScheme.H:334
Foam::fv::CrankNicolsonDdtScheme::fvcDdtPhiCorr
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
Definition CrankNicolsonDdtScheme.C:1304
Foam::fv::CrankNicolsonDdtScheme::mesh
const fvMesh & mesh() const
Return mesh reference.
Definition CrankNicolsonDdtScheme.H:268
Foam::fv::CrankNicolsonDdtScheme::meshPhi
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
Definition CrankNicolsonDdtScheme.C:1669
Foam::fv::ddtScheme
Abstract base class for ddt schemes.
Definition ddtScheme.H:85
Foam::fv::ddtScheme::ddtScheme
ddtScheme(const ddtScheme &)=delete
No copy construct.
Foam::fv::ddtScheme::fvcDdtPhiCoeff
tmp< surfaceScalarField > fvcDdtPhiCoeff(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi, const fluxFieldType &phiCorr)
Definition ddtScheme.C:122
Foam::regIOobject::store
bool store()
Register object with its registry and transfer ownership to the registry.
Definition regIOobjectI.H:36
Foam::tmp
A class for managing temporary objects.
Definition tmp.H:75
Foam::tmp::ref
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
Definition tmpI.H:235
Foam::token
A token holds an item read from Istream.
Definition token.H:70
Foam::token::isNumber
bool isNumber() const noexcept
Token is (signed/unsigned) integer type, FLOAT or DOUBLE.
Definition tokenI.H:992
Foam::token::number
scalar number() const
Return label, float or double value.
Definition tokenI.H:998
Foam::word
A class for handling words, derived from Foam::string.
Definition word.H:66
Foam::zero
A class representing the concept of 0 (zero) that can be used to avoid manipulating objects known to ...
Definition zero.H:58
U
U
Definition pEqn.H:72
mesh
dynamicFvMesh & mesh
Definition createDynamicFvMesh.H:6
runTime
engineTime & runTime
Definition createEngineTime.H:13
oldTime
Info<< "Creating field kinetic energy K\n"<< endl;volScalarField K("K", 0.5 *magSqr(U));if(U.nOldTimes()){ volVectorField *Uold=&U.oldTime();volScalarField *Kold=&K.oldTime(); *Kold==0.5 *magSqr(*Uold);while(Uold->nOldTimes()) { Uold=&Uold-> oldTime()
Uf
autoPtr< surfaceVectorField > Uf
Definition createUfIfPresent.H:27
FatalIOErrorInFunction
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition error.H:629
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition error.H:600
io
const auto & io
Definition findMeshDefinitionDict.H:139
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.
name
auto & name
Definition getAllFaRegionOptions.H:175
timeName
word timeName
Definition getTimeIndex.H:3
Foam::fv
Namespace for finite-volume.
Definition atmAmbientTurbSource.C:30
Foam::fv::ff
const FieldField< fvPatchField, Type > & ff(const FieldField< fvPatchField, Type > &bf)
Definition CrankNicolsonDdtScheme.C:258
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.
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::fvc::ddtCorr
tmp< GeometricField< typename flux< Type >::type, fvsPatchField, surfaceMesh > > ddtCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition fvcDdt.C:165
Foam::fvm
Namespace of functions to calculate implicit derivatives returning a matrix.
Foam::stringOps::evaluate
string evaluate(label fieldWidth, const std::string &s, size_t pos=0, size_t len=std::string::npos)
String evaluation with specified (positive, non-zero) field width.
Foam
Namespace for OpenFOAM.
Definition atmBoundaryLayer.C:27
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition dimensionSets.H:51
Foam::volScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition volFieldsFwd.H:72
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition dimensionSets.H:57
Foam::dimVelocity
const dimensionSet dimVelocity
Foam::surfaceScalarField
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Definition surfaceFieldsFwd.H:62
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition primitiveFieldsFwd.H:46
Foam::abort
errorManip< error > abort(error &err)
Definition errorManip.H:139
Foam::FatalIOError
IOerror FatalIOError
Error stream (stdout output on all processes), with additional 'FOAM FATAL IO ERROR' header text and ...
Foam::FatalError
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
Foam::dimVolume
const dimensionSet dimVolume(pow3(dimLength))
Definition dimensionSets.H:58
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
Definition exprTraits.C:127
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition dimensionedScalarFwd.H:35
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition errorManip.H:125
Foam::dimVol
const dimensionSet dimVol(dimVolume)
Older spelling for dimVolume.
Definition dimensionSets.H:65
timeIndex
label timeIndex
Definition getTimeIndex.H:24
alpha
volScalarField & alpha
Definition readThermalProperties.H:212
dict
dictionary dict
Definition searchingEngine.H:11