boundedBackwardFaDdtScheme.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 "boundedBackwardFaDdtScheme.H"
29 #include "faMatrices.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 namespace fa
39 {
40 
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 
43 scalar boundedBackwardFaDdtScheme::deltaT_() const
44 {
45  return mesh().time().deltaTValue();
46 }
47 
48 
49 scalar boundedBackwardFaDdtScheme::deltaT0_() const
50 {
51  return mesh().time().deltaT0Value();
52 }
53 
54 
55 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
56 
57 tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
58 (
59  const dimensionedScalar dt
60 )
61 {
62  // No change compared to backward
63 
64  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
65 
66  const IOobject ddtIOobject
67  (
68  mesh().thisDb().newIOobject
69  (
70  "ddt("+dt.name()+')',
71  { IOobject::REGISTER }
72  )
73  );
74 
75  scalar deltaT = deltaT_();
76  scalar deltaT0 = deltaT0_();
77 
78  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
79  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
80  scalar coefft0 = coefft + coefft00;
81 
82  if (mesh().moving())
83  {
84  tmp<areaScalarField> tdtdt
85  (
86  new areaScalarField
87  (
88  ddtIOobject,
89  mesh(),
90  dimensionedScalar(dt.dimensions()/dimTime, Zero)
91  )
92  );
93 
94  tdtdt.ref().primitiveFieldRef() = rDeltaT.value()*dt.value()*
95  (
96  coefft - (coefft0*mesh().S0() - coefft00*mesh().S00())/mesh().S()
97  );
98 
99  return tdtdt;
100  }
101 
102  return tmp<areaScalarField>::New
103  (
104  ddtIOobject,
105  mesh(),
106  dimensionedScalar(dt.dimensions()/dimTime, Zero),
107  faPatchFieldBase::calculatedType()
108  );
109 }
110 
111 
112 tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
113 (
114  const dimensionedScalar dt
115 )
116 {
117  // No change compared to backward
118 
119  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
120 
121  const IOobject ddtIOobject
122  (
123  mesh().thisDb().newIOobject
124  (
125  "ddt("+dt.name()+')',
126  { IOobject::REGISTER }
127  )
128  );
129 
130  scalar deltaT = deltaT_();
131  scalar deltaT0 = deltaT0_();
132 
133  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
134  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
135  scalar coefft0 = coefft + coefft00;
136 
137  tmp<areaScalarField> tdtdt0
138  (
139  new areaScalarField
140  (
141  ddtIOobject,
142  mesh(),
143  -rDeltaT*(coefft0 - coefft00)*dt
144  )
145  );
146 
147  if (mesh().moving())
148  {
149  tdtdt0.ref().primitiveFieldRef() = (-rDeltaT.value()*dt.value())*
150  (
151  (coefft0*mesh().S0() - coefft00*mesh().S00())/mesh().S()
152  );
153  }
154 
155  return tdtdt0;
156 }
157 
158 
159 tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
160 (
161  const areaScalarField& vf
162 )
163 {
164  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
165 
166  const IOobject ddtIOobject
167  (
168  mesh().thisDb().newIOobject
169  (
170  "ddt("+vf.name()+')',
171  { IOobject::REGISTER }
172  )
173  );
174 
175  scalar deltaT = deltaT_();
176  scalar deltaT0 = deltaT0_(vf);
177 
178  // Calculate unboundedness indicator
179  // Note: all times moved by one because access to internal field
180  // copies current field into the old-time level.
181  areaScalarField phict
182  (
183  mag
184  (
185  vf.oldTime().oldTime()
186  - vf.oldTime().oldTime().oldTime()
187  )/
188  (
189  mag
190  (
191  vf.oldTime()
192  - vf.oldTime().oldTime()
193  )
194  + dimensionedScalar("small", vf.dimensions(), SMALL)
195  )
196  );
197 
198  areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
199 
200  areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
201  areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
202  areaScalarField coefft0(coefft + coefft00);
203 
204  if (mesh().moving())
205  {
206  return tmp<areaScalarField>
207  (
208  new areaScalarField
209  (
210  ddtIOobject,
211  mesh(),
212  rDeltaT.dimensions()*vf.dimensions(),
213  rDeltaT.value()*
214  (
215  coefft*vf.primitiveField() -
216  (
217  coefft0.primitiveField()
218  *vf.oldTime().primitiveField()*mesh().S0()
219  - coefft00.primitiveField()
220  *vf.oldTime().oldTime().primitiveField()
221  *mesh().S00()
222  )/mesh().S()
223  ),
224  rDeltaT.value()*
225  (
226  coefft.boundaryField()*vf.boundaryField() -
227  (
228  coefft0.boundaryField()*
229  vf.oldTime().boundaryField()
230  - coefft00.boundaryField()*
231  vf.oldTime().oldTime().boundaryField()
232  )
233  )
234  )
235  );
236  }
237  else
238  {
239  return tmp<areaScalarField>
240  (
241  new areaScalarField
242  (
243  ddtIOobject,
244  rDeltaT*
245  (
246  coefft*vf
247  - coefft0*vf.oldTime()
248  + coefft00*vf.oldTime().oldTime()
249  )
250  )
251  );
252  }
253 }
254 
255 
256 tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
257 (
258  const areaScalarField& vf
259 )
260 {
261  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
262 
263  const IOobject ddtIOobject
264  (
265  mesh().thisDb().newIOobject
266  (
267  "ddt0("+vf.name()+')',
268  { IOobject::REGISTER }
269  )
270  );
271 
272  scalar deltaT = deltaT_();
273  scalar deltaT0 = deltaT0_(vf);
274 
275  // Calculate unboundedness indicator
276  // Note: all times moved by one because access to internal field
277  // copies current field into the old-time level.
278  areaScalarField phict
279  (
280  mag
281  (
282  vf.oldTime().oldTime()
283  - vf.oldTime().oldTime().oldTime()
284  )/
285  (
286  mag
287  (
288  vf.oldTime()
289  - vf.oldTime().oldTime()
290  )
291  + dimensionedScalar("small", vf.dimensions(), SMALL)
292  )
293  );
294 
295  areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
296 
297  areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
298  areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
299  areaScalarField coefft0(coefft + coefft00);
300 
301  if (mesh().moving())
302  {
303  return tmp<areaScalarField>
304  (
305  new areaScalarField
306  (
307  ddtIOobject,
308  mesh(),
309  rDeltaT.dimensions()*vf.dimensions(),
310  rDeltaT.value()*
311  (
312  - (
313  coefft0.primitiveField()*
314  vf.oldTime().primitiveField()*mesh().S0()
315  - coefft00.primitiveField()*
316  vf.oldTime().oldTime().primitiveField()
317  *mesh().S00()
318  )/mesh().S()
319  ),
320  rDeltaT.value()*
321  (
322  - (
323  coefft0.boundaryField()*
324  vf.oldTime().boundaryField()
325  - coefft00.boundaryField()*
326  vf.oldTime().oldTime().boundaryField()
327  )
328  )
329  )
330  );
331  }
332  else
333  {
334  return tmp<areaScalarField>
335  (
336  new areaScalarField
337  (
338  ddtIOobject,
339  rDeltaT*
340  (
341  - coefft0*vf.oldTime()
342  + coefft00*vf.oldTime().oldTime()
343  )
344  )
345  );
346  }
347 }
348 
349 
350 tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
351 (
352  const dimensionedScalar& rho,
353  const areaScalarField& vf
354 )
355 {
356  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
357 
358  const IOobject ddtIOobject
359  (
360  mesh().thisDb().newIOobject
361  (
362  "ddt("+rho.name()+','+vf.name()+')',
363  { IOobject::REGISTER }
364  )
365  );
366 
367  scalar deltaT = deltaT_();
368  scalar deltaT0 = deltaT0_(vf);
369 
370  // Calculate unboundedness indicator
371  // Note: all times moved by one because access to internal field
372  // copies current field into the old-time level.
373  areaScalarField phict
374  (
375  mag
376  (
377  vf.oldTime().oldTime()
378  - vf.oldTime().oldTime().oldTime()
379  )/
380  (
381  mag
382  (
383  vf.oldTime()
384  - vf.oldTime().oldTime()
385  )
386  + dimensionedScalar("small", vf.dimensions(), SMALL)
387  )
388  );
389 
390  areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
391 
392  areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
393  areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
394  areaScalarField coefft0(coefft + coefft00);
395 
396  if (mesh().moving())
397  {
398  return tmp<areaScalarField>
399  (
400  new areaScalarField
401  (
402  ddtIOobject,
403  mesh(),
404  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
405  rDeltaT.value()*rho.value()*
406  (
407  coefft*vf.primitiveField() -
408  (
409  coefft0.primitiveField()*
410  vf.oldTime().primitiveField()*mesh().S0()
411  - coefft00.primitiveField()*
412  vf.oldTime().oldTime().primitiveField()
413  *mesh().S00()
414  )/mesh().S()
415  ),
416  rDeltaT.value()*rho.value()*
417  (
418  coefft.boundaryField()*vf.boundaryField() -
419  (
420  coefft0.boundaryField()*
421  vf.oldTime().boundaryField()
422  - coefft00.boundaryField()*
423  vf.oldTime().oldTime().boundaryField()
424  )
425  )
426  )
427  );
428  }
429  else
430  {
431  return tmp<areaScalarField>
432  (
433  new areaScalarField
434  (
435  ddtIOobject,
436  rDeltaT*rho*
437  (
438  coefft*vf
439  - coefft0*vf.oldTime()
440  + coefft00*vf.oldTime().oldTime()
441  )
442  )
443  );
444  }
445 }
446 
447 tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
448 (
449  const dimensionedScalar& rho,
450  const areaScalarField& vf
451 )
452 {
453  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
454 
455  const IOobject ddtIOobject
456  (
457  mesh().thisDb().newIOobject
458  (
459  "ddt0("+rho.name()+','+vf.name()+')',
460  { IOobject::REGISTER }
461  )
462  );
463 
464  scalar deltaT = deltaT_();
465  scalar deltaT0 = deltaT0_(vf);
466 
467  // Calculate unboundedness indicator
468  // Note: all times moved by one because access to internal field
469  // copies current field into the old-time level.
470  areaScalarField phict
471  (
472  mag
473  (
474  vf.oldTime().oldTime()
475  - vf.oldTime().oldTime().oldTime()
476  )/
477  (
478  mag
479  (
480  vf.oldTime()
481  - vf.oldTime().oldTime()
482  )
483  + dimensionedScalar("small", vf.dimensions(), SMALL)
484  )
485  );
486 
487  areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
488 
489  areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
490  areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
491  areaScalarField coefft0(coefft + coefft00);
492 
493  if (mesh().moving())
494  {
495  return tmp<areaScalarField>
496  (
497  new areaScalarField
498  (
499  ddtIOobject,
500  mesh(),
501  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
502  rDeltaT.value()*rho.value()*
503  (
504  -(
505  coefft0.primitiveField()*
506  vf.oldTime().primitiveField()*mesh().S0()
507  - coefft00.primitiveField()*
508  vf.oldTime().oldTime().primitiveField()
509  *mesh().S00()
510  )/mesh().S()
511  ),
512  rDeltaT.value()*rho.value()*
513  (
514  -(
515  coefft0.boundaryField()*
516  vf.oldTime().boundaryField()
517  - coefft00.boundaryField()*
518  vf.oldTime().oldTime().boundaryField()
519  )
520  )
521  )
522  );
523  }
524  else
525  {
526  return tmp<areaScalarField>
527  (
528  new areaScalarField
529  (
530  ddtIOobject,
531  rDeltaT*rho*
532  (
533  - coefft0*vf.oldTime()
534  + coefft00*vf.oldTime().oldTime()
535  )
536  )
537  );
538  }
539 }
540 
541 
542 tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
543 (
544  const areaScalarField& rho,
545  const areaScalarField& vf
546 )
547 {
548  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
549 
550  const IOobject ddtIOobject
551  (
552  mesh().thisDb().newIOobject
553  (
554  "ddt("+rho.name()+','+vf.name()+')',
555  { IOobject::REGISTER }
556  )
557  );
558 
559  scalar deltaT = deltaT_();
560  scalar deltaT0 = deltaT0_(vf);
561 
562  // Calculate unboundedness indicator
563  // Note: all times moved by one because access to internal field
564  // copies current field into the old-time level.
565  areaScalarField phict
566  (
567  mag
568  (
569  rho.oldTime().oldTime()*vf.oldTime().oldTime()
570  - rho.oldTime().oldTime().oldTime()*vf.oldTime().oldTime().oldTime()
571  )/
572  (
573  mag
574  (
575  rho.oldTime()*vf.oldTime()
576  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
577  )
578  + dimensionedScalar("small", rho.dimensions()*vf.dimensions(), SMALL)
579  )
580  );
581 
582  areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
583 
584  areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
585  areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
586  areaScalarField coefft0(coefft + coefft00);
587 
588  if (mesh().moving())
589  {
590  return tmp<areaScalarField>
591  (
592  new areaScalarField
593  (
594  ddtIOobject,
595  mesh(),
596  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
597  rDeltaT.value()*
598  (
599  coefft*rho.primitiveField()*vf.primitiveField() -
600  (
601  coefft0.primitiveField()*
602  rho.oldTime().primitiveField()*
603  vf.oldTime().primitiveField()*mesh().S0()
604  - coefft00.primitiveField()*
605  rho.oldTime().oldTime().primitiveField()
606  *vf.oldTime().oldTime().primitiveField()*mesh().S00()
607  )/mesh().S()
608  ),
609  rDeltaT.value()*
610  (
611  coefft.boundaryField()*vf.boundaryField() -
612  (
613  coefft0.boundaryField()*
614  rho.oldTime().boundaryField()*
615  vf.oldTime().boundaryField()
616  - coefft00.boundaryField()*
617  rho.oldTime().oldTime().boundaryField()*
618  vf.oldTime().oldTime().boundaryField()
619  )
620  )
621  )
622  );
623  }
624  else
625  {
626  return tmp<areaScalarField>
627  (
628  new areaScalarField
629  (
630  ddtIOobject,
631  rDeltaT*
632  (
633  coefft*rho*vf
634  - coefft0*rho.oldTime()*vf.oldTime()
635  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
636  )
637  )
638  );
639  }
640 }
641 
642 
643 tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
644 (
645  const areaScalarField& rho,
646  const areaScalarField& vf
647 )
648 {
649  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
650 
651  const IOobject ddtIOobject
652  (
653  mesh().thisDb().newIOobject
654  (
655  "ddt0("+rho.name()+','+vf.name()+')',
656  { IOobject::REGISTER }
657  )
658  );
659 
660  scalar deltaT = deltaT_();
661  scalar deltaT0 = deltaT0_(vf);
662 
663  // Calculate unboundedness indicator
664  // Note: all times moved by one because access to internal field
665  // copies current field into the old-time level.
666  areaScalarField phict
667  (
668  mag
669  (
670  rho.oldTime().oldTime()*vf.oldTime().oldTime()
671  - rho.oldTime().oldTime().oldTime()*vf.oldTime().oldTime().oldTime()
672  )/
673  (
674  mag
675  (
676  rho.oldTime()*vf.oldTime()
677  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
678  )
679  + dimensionedScalar("small", rho.dimensions()*vf.dimensions(), SMALL)
680  )
681  );
682 
683  areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
684 
685  areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
686  areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
687  areaScalarField coefft0(coefft + coefft00);
688 
689  if (mesh().moving())
690  {
691  return tmp<areaScalarField>
692  (
693  new areaScalarField
694  (
695  ddtIOobject,
696  mesh(),
697  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
698  rDeltaT.value()*
699  (
700  - (
701  coefft0.primitiveField()*
702  rho.oldTime().primitiveField()*
703  vf.oldTime().primitiveField()*mesh().S0()
704  - coefft00.primitiveField()*
705  rho.oldTime().oldTime().primitiveField()*
706  vf.oldTime().oldTime().primitiveField()*mesh().S00()
707  )/mesh().S()
708  ),
709  rDeltaT.value()*
710  (
711  - (
712  coefft0.boundaryField()*
713  rho.oldTime().boundaryField()*
714  vf.oldTime().boundaryField()
715  - coefft00.boundaryField()*
716  rho.oldTime().oldTime().boundaryField()*
717  vf.oldTime().oldTime().boundaryField()
718  )
719  )
720  )
721  );
722  }
723  else
724  {
725  return tmp<areaScalarField>
726  (
727  new areaScalarField
728  (
729  ddtIOobject,
730  rDeltaT*
731  (
732  - coefft0*rho.oldTime()*vf.oldTime()
733  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
734  )
735  )
736  );
737  }
738 }
739 
740 
741 tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
742 (
743  const areaScalarField& vf
744 )
745 {
746  tmp<faScalarMatrix> tfam
747  (
748  new faScalarMatrix
749  (
750  vf,
751  vf.dimensions()*dimArea/dimTime
752  )
753  );
754 
755  faScalarMatrix& fam = tfam.ref();
756 
757  scalar rDeltaT = 1.0/deltaT_();
758 
759  scalar deltaT = deltaT_();
760  scalar deltaT0 = deltaT0_(vf);
761 
762  // Calculate unboundedness indicator
763  // Note: all times moved by one because access to internal field
764  // copies current field into the old-time level.
765  scalarField phict
766  (
767  mag
768  (
769  vf.oldTime().oldTime().internalField()
770  - vf.oldTime().oldTime().oldTime().internalField()
771  )/
772  (
773  mag
774  (
775  vf.oldTime().internalField()
776  - vf.oldTime().oldTime().internalField()
777  )
778  + SMALL
779  )
780  );
781 
782  scalarField limiter(pos(phict) - pos(phict - 1.0));
783 
784  scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
785  scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)));
786  scalarField coefft0(coefft + coefft00);
787 
788  fam.diag() = (coefft*rDeltaT)*mesh().S();
789 
790  if (mesh().moving())
791  {
792  fam.source() = rDeltaT*
793  (
794  coefft0*vf.oldTime().primitiveField()*mesh().S0()
795  - coefft00*vf.oldTime().oldTime().primitiveField()
796  *mesh().S00()
797  );
798  }
799  else
800  {
801  fam.source() = rDeltaT*mesh().S()*
802  (
803  coefft0*vf.oldTime().primitiveField()
804  - coefft00*vf.oldTime().oldTime().primitiveField()
805  );
806  }
807 
808  return tfam;
809 }
810 
811 
812 tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
813 (
814  const dimensionedScalar& rho,
815  const areaScalarField& vf
816 )
817 {
818  tmp<faScalarMatrix> tfam
819  (
820  new faScalarMatrix
821  (
822  vf,
823  rho.dimensions()*vf.dimensions()*dimArea/dimTime
824  )
825  );
826  faScalarMatrix& fam = tfam.ref();
827 
828  scalar rDeltaT = 1.0/deltaT_();
829 
830  scalar deltaT = deltaT_();
831  scalar deltaT0 = deltaT0_(vf);
832 
833  // Calculate unboundedness indicator
834  // Note: all times moved by one because access to internal field
835  // copies current field into the old-time level.
836  scalarField phict
837  (
838  mag
839  (
840  vf.oldTime().oldTime().internalField()
841  - vf.oldTime().oldTime().oldTime().internalField()
842  )/
843  (
844  mag
845  (
846  vf.oldTime().internalField()
847  - vf.oldTime().oldTime().internalField()
848  )
849  + SMALL
850  )
851  );
852 
853  scalarField limiter(pos(phict) - pos(phict - 1.0));
854 
855  scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
856  scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)));
857  scalarField coefft0(coefft + coefft00);
858 
859  fam.diag() = (coefft*rDeltaT*rho.value())*mesh().S();
860 
861  if (mesh().moving())
862  {
863  fam.source() = rDeltaT*rho.value()*
864  (
865  coefft0*vf.oldTime().primitiveField()*mesh().S0()
866  - coefft00*vf.oldTime().oldTime().primitiveField()
867  *mesh().S00()
868  );
869  }
870  else
871  {
872  fam.source() = rDeltaT*mesh().S()*rho.value()*
873  (
874  coefft0*vf.oldTime().primitiveField()
875  - coefft00*vf.oldTime().oldTime().primitiveField()
876  );
877  }
878 
879  return tfam;
880 }
881 
882 
883 tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
884 (
885  const areaScalarField& rho,
886  const areaScalarField& vf
887 )
888 {
889  tmp<faScalarMatrix> tfam
890  (
891  new faScalarMatrix
892  (
893  vf,
894  rho.dimensions()*vf.dimensions()*dimArea/dimTime
895  )
896  );
897  faScalarMatrix& fam = tfam.ref();
898 
899  scalar rDeltaT = 1.0/deltaT_();
900 
901  scalar deltaT = deltaT_();
902  scalar deltaT0 = deltaT0_(vf);
903 
904  // Calculate unboundedness indicator
905  // Note: all times moved by one because access to internal field
906  // copies current field into the old-time level.
907  scalarField phict
908  (
909  mag
910  (
911  rho.oldTime().oldTime().internalField()*
912  vf.oldTime().oldTime().internalField()
913  - rho.oldTime().oldTime().oldTime().internalField()*
914  vf.oldTime().oldTime().oldTime().internalField()
915  )/
916  (
917  mag
918  (
919  rho.oldTime().internalField()*
920  vf.oldTime().internalField()
921  - rho.oldTime().oldTime().internalField()*
922  vf.oldTime().oldTime().internalField()
923  )
924  + SMALL
925  )
926  );
927 
928  scalarField limiter(pos(phict) - pos(phict - 1.0));
929 
930  scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
931  scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)));
932  scalarField coefft0(coefft + coefft00);
933 
934  fam.diag() = (coefft*rDeltaT)*rho.primitiveField()*mesh().S();
935 
936  if (mesh().moving())
937  {
938  fam.source() = rDeltaT*
939  (
940  coefft0*rho.oldTime().primitiveField()
941  *vf.oldTime().primitiveField()*mesh().S0()
942  - coefft00*rho.oldTime().oldTime().primitiveField()
943  *vf.oldTime().oldTime().primitiveField()*mesh().S00()
944  );
945  }
946  else
947  {
948  fam.source() = rDeltaT*mesh().S()*
949  (
950  coefft0*rho.oldTime().primitiveField()
951  *vf.oldTime().primitiveField()
952  - coefft00*rho.oldTime().oldTime().primitiveField()
953  *vf.oldTime().oldTime().primitiveField()
954  );
955  }
956 
957  return tfam;
958 }
959 
960 
961 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
962 
963 defineTypeNameAndDebug(boundedBackwardFaDdtScheme, 0);
964 
965 faDdtScheme<scalar>::addIstreamConstructorToTable<boundedBackwardFaDdtScheme>
966  addboundedBackwardFaDdtSchemeIstreamConstructorToTable_;
967 
968 
969 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
970 
971 } // End namespace fa
972 
973 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
974 
975 } // End namespace Foam
976 
977 // ************************************************************************* //
Foam::dimensioned::value
const Type & value() const noexcept
Return const reference to value.
Definition: dimensionedType.H:366
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:938
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const noexcept
Return a const-reference to the internal field values.
Definition: GeometricFieldI.H:36
Foam::fa::boundedBackwardFaDdtScheme::facDdt0
tmp< areaScalarField > facDdt0(const dimensionedScalar)
Definition: boundedBackwardFaDdtScheme.C:106
Foam::TimeState::deltaTValue
scalar deltaTValue() const noexcept
Return time step value.
Definition: TimeStateI.H:49
Foam::faMesh::S00
const DimensionedField< scalar, areaMesh > & S00() const
Return old-old-time face areas.
Definition: faMesh.C:1159
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::fa::boundedBackwardFaDdtScheme::facDdt
tmp< areaScalarField > facDdt(const dimensionedScalar)
Definition: boundedBackwardFaDdtScheme.C:51
Foam::faMesh::time
const Time & time() const
Return reference to time.
Definition: faMesh.C:1022
Foam::IOobject::name
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:195
Foam::tmp::ref
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
Definition: tmpI.H:235
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:44
Foam::GeometricField::internalField
const Internal & internalField() const noexcept
Return a const-reference to the dimensioned internal field.
Definition: GeometricFieldI.H:27
Foam::fa::addboundedBackwardFaDdtSchemeIstreamConstructorToTable_
faDdtScheme< scalar >::addIstreamConstructorToTable< boundedBackwardFaDdtScheme > addboundedBackwardFaDdtSchemeIstreamConstructorToTable_
Definition: boundedBackwardFaDdtScheme.C:959
Foam::faMesh::S
const DimensionedField< scalar, areaMesh > & S() const
Return face areas.
Definition: faMesh.C:1133
Foam::limiter
tmp< areaScalarField > limiter(const areaScalarField &phi)
Definition: faNVDscheme.C:31
Foam::pos
dimensionedScalar pos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:170
Foam::fa::boundedBackwardFaDdtScheme::famDdt
tmp< faScalarMatrix > famDdt(const areaScalarField &)
Definition: boundedBackwardFaDdtScheme.C:735
Foam::fa::defineTypeNameAndDebug
defineTypeNameAndDebug(limitHeight, 0)
Foam::tmp::New
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
Definition: tmp.H:206
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const noexcept
Return const reference to dimensions.
Definition: dimensionedType.H:356
Foam::fa::boundedBackwardFaDdtScheme::mesh
const faMesh & mesh() const
Return mesh reference.
Definition: boundedBackwardFaDdtScheme.H:131
Foam::TimeState::deltaT0Value
scalar deltaT0Value() const noexcept
Return old time step value.
Definition: TimeStateI.H:55
Foam::faPatchFieldBase::calculatedType
static const word & calculatedType() noexcept
The type name for calculated patch fields.
Definition: faPatchField.H:174
Foam::dimensioned::name
const word & name() const noexcept
Return const reference to name.
Definition: dimensionedType.H:346
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:33
Foam::faMesh::moving
bool moving() const
Is mesh moving.
Definition: faMesh.H:1481
Foam::faScalarMatrix
Template specialisation for scalar faMatrix.
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
Foam::tmp
A class for managing temporary objects.
Definition: HashPtrTable.H:50
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported...
Definition: IOobject.H:180
fam
Calculate the matrix for the second temporal derivative.
Foam::IOobjectOption::REGISTER
Request registration (bool: true)
Definition: IOobjectOption.H:102
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const noexcept
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:45
Foam::TimeState::deltaT
dimensionedScalar deltaT() const
Return time step.
Definition: TimeStateI.H:61
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:26
Foam::faMesh::S0
const DimensionedField< scalar, areaMesh > & S0() const
Return old-time face areas.
Definition: faMesh.C:1145
Foam::DimensionedField::dimensions
const dimensionSet & dimensions() const noexcept
Return dimensions.
Definition: DimensionedFieldI.H:34
Foam::areaScalarField
GeometricField< scalar, faPatchField, areaMesh > areaScalarField
Definition: areaFieldsFwd.H:72
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:127