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 | Copyright (C) 2011-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "CrankNicolsonDdtScheme.H"
27 #include "surfaceInterpolate.H"
28 #include "fvcDiv.H"
29 #include "fvMatrices.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 namespace fv
39 {
40 
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 
43 template<class Type>
44 template<class GeoField>
45 CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::DDt0Field
46 (
47  const IOobject& io,
48  const fvMesh& mesh
49 )
50 :
51  GeoField(io, mesh),
52  startTimeIndex_(-2) // This field is for a restart and thus correct so set
53  // the start time-index to correspond to a previous run
54 {
55  // Set the time-index to the beginning of the run to ensure the field
56  // is updated during the first time-step
57  this->timeIndex() = mesh.time().startTimeIndex();
58 }
59 
60 
61 template<class Type>
62 template<class GeoField>
63 CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::DDt0Field
64 (
65  const IOobject& io,
66  const fvMesh& mesh,
67  const dimensioned<typename GeoField::value_type>& dimType
68 )
69 :
70  GeoField(io, mesh, dimType),
71  startTimeIndex_(mesh.time().timeIndex())
72 {}
73 
74 
75 template<class Type>
76 template<class GeoField>
77 label CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::
78 startTimeIndex() const
79 {
80  return startTimeIndex_;
81 }
82 
83 
84 template<class Type>
85 template<class GeoField>
86 GeoField& CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::
87 operator()()
88 {
89  return *this;
90 }
91 
92 
93 template<class Type>
94 template<class GeoField>
95 void CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::
96 operator=(const GeoField& gf)
97 {
98  GeoField::operator=(gf);
99 }
100 
101 
102 template<class Type>
103 template<class GeoField>
104 CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>&
105 CrankNicolsonDdtScheme<Type>::ddt0_
106 (
107  const word& name,
108  const dimensionSet& dims
109 )
110 {
111  if (!mesh().objectRegistry::template foundObject<GeoField>(name))
112  {
113  const Time& runTime = mesh().time();
114  word startTimeName = runTime.timeName(runTime.startTime().value());
115 
116  if
117  (
118  (
119  runTime.timeIndex() == runTime.startTimeIndex()
120  || runTime.timeIndex() == runTime.startTimeIndex() + 1
121  )
122  && IOobject
123  (
124  name,
125  startTimeName,
126  mesh()
127  ).headerOk()
128  )
129  {
131  (
132  new DDt0Field<GeoField>
133  (
134  IOobject
135  (
136  name,
137  startTimeName,
138  mesh(),
141  ),
142  mesh()
143  )
144  );
145  }
146  else
147  {
149  (
150  new DDt0Field<GeoField>
151  (
152  IOobject
153  (
154  name,
155  mesh().time().timeName(),
156  mesh(),
159  ),
160  mesh(),
161  dimensioned<typename GeoField::value_type>
162  (
163  "0",
164  dims/dimTime,
165  Zero
166  )
167  )
168  );
169  }
170  }
171 
172  DDt0Field<GeoField>& ddt0 = static_cast<DDt0Field<GeoField>&>
173  (
174  const_cast<GeoField&>
175  (
176  mesh().objectRegistry::template lookupObject<GeoField>(name)
177  )
178  );
179 
180  return ddt0;
181 }
182 
183 
184 template<class Type>
185 template<class GeoField>
186 bool CrankNicolsonDdtScheme<Type>::evaluate
187 (
188  const DDt0Field<GeoField>& ddt0
189 ) const
190 {
191  return ddt0.timeIndex() != mesh().time().timeIndex();
192 }
193 
194 template<class Type>
195 template<class GeoField>
196 scalar CrankNicolsonDdtScheme<Type>::coef_
197 (
198  const DDt0Field<GeoField>& ddt0
199 ) const
200 {
201  if (mesh().time().timeIndex() - ddt0.startTimeIndex() > 0)
202  {
203  return 1.0 + ocCoeff_;
204  }
205  else
206  {
207  return 1.0;
208  }
209 }
210 
211 
212 template<class Type>
213 template<class GeoField>
214 scalar CrankNicolsonDdtScheme<Type>::coef0_
215 (
216  const DDt0Field<GeoField>& ddt0
217 ) const
218 {
219  if (mesh().time().timeIndex() - ddt0.startTimeIndex() > 1)
220  {
221  return 1.0 + ocCoeff_;
222  }
223  else
224  {
225  return 1.0;
226  }
227 }
228 
229 
230 template<class Type>
231 template<class GeoField>
232 dimensionedScalar CrankNicolsonDdtScheme<Type>::rDtCoef_
233 (
234  const DDt0Field<GeoField>& ddt0
235 ) const
236 {
237  return coef_(ddt0)/mesh().time().deltaT();
238 }
239 
240 
241 template<class Type>
242 template<class GeoField>
243 dimensionedScalar CrankNicolsonDdtScheme<Type>::rDtCoef0_
244 (
245  const DDt0Field<GeoField>& ddt0
246 ) const
247 {
248  return coef0_(ddt0)/mesh().time().deltaT0();
249 }
250 
251 
252 template<class Type>
253 template<class GeoField>
254 tmp<GeoField> CrankNicolsonDdtScheme<Type>::offCentre_
255 (
256  const GeoField& ddt0
257 ) const
258 {
259  if (ocCoeff_ < 1.0)
260  {
261  return ocCoeff_*ddt0;
262  }
263  else
264  {
265  return ddt0;
266  }
267 }
268 
269 
270 template<class Type>
271 const FieldField<fvPatchField, Type>& ff
272 (
274 )
275 {
276  return bf;
277 }
278 
279 
280 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
281 
282 template<class Type>
285 (
286  const dimensioned<Type>& dt
287 )
288 {
289  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
290  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
291  (
292  "ddt0(" + dt.name() + ')',
293  dt.dimensions()
294  );
295 
296  IOobject ddtIOobject
297  (
298  "ddt(" + dt.name() + ')',
299  mesh().time().timeName(),
300  mesh()
301  );
302 
304  (
306  (
307  ddtIOobject,
308  mesh(),
310  (
311  "0",
312  dt.dimensions()/dimTime,
313  Zero
314  )
315  )
316  );
317 
318  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
319 
320  if (mesh().moving())
321  {
322  if (evaluate(ddt0))
323  {
324  dimensionedScalar rDtCoef0 = rDtCoef0_(ddt0);
325 
326  ddt0.ref() =
327  (
328  (rDtCoef0*dt)*(mesh().V0() - mesh().V00())
329  - mesh().V00()*offCentre_(ddt0.internalField())
330  )/mesh().V0();
331  }
332 
333  tdtdt.ref().ref() =
334  (
335  (rDtCoef*dt)*(mesh().V() - mesh().V0())
336  - mesh().V0()*offCentre_(ddt0.internalField())
337  )/mesh().V();
338  }
339 
340  return tdtdt;
341 }
342 
343 
344 template<class Type>
347 (
349 )
350 {
351  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
352  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
353  (
354  "ddt0(" + vf.name() + ')',
355  vf.dimensions()
356  );
357 
358  IOobject ddtIOobject
359  (
360  "ddt(" + vf.name() + ')',
361  mesh().time().timeName(),
362  mesh()
363  );
364 
365  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
366 
367  if (mesh().moving())
368  {
369  if (evaluate(ddt0))
370  {
371  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
372 
373  ddt0.primitiveFieldRef() =
374  (
375  rDtCoef0*
376  (
377  mesh().V0()*vf.oldTime().primitiveField()
378  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
379  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
380  )/mesh().V0();
381 
382  ddt0.boundaryFieldRef() =
383  (
384  rDtCoef0*
385  (
386  vf.oldTime().boundaryField()
387  - vf.oldTime().oldTime().boundaryField()
388  ) - offCentre_(ff(ddt0.boundaryField()))
389  );
390  }
391 
393  (
395  (
396  ddtIOobject,
397  (
398  rDtCoef*
399  (
400  mesh().V()*vf
401  - mesh().V0()*vf.oldTime()
402  ) - mesh().V0()*offCentre_(ddt0()())
403  )/mesh().V(),
404  rDtCoef.value()*
405  (
406  vf.boundaryField() - vf.oldTime().boundaryField()
407  ) - offCentre_(ff(ddt0.boundaryField()))
408  )
409  );
410  }
411  else
412  {
413  if (evaluate(ddt0))
414  {
415  ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime())
416  - offCentre_(ddt0());
417  }
418 
420  (
422  (
423  ddtIOobject,
424  rDtCoef*(vf - vf.oldTime()) - offCentre_(ddt0())
425  )
426  );
427  }
428 }
429 
430 
431 template<class Type>
434 (
435  const dimensionedScalar& rho,
437 )
438 {
439  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
440  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
441  (
442  "ddt0(" + rho.name() + ',' + vf.name() + ')',
443  rho.dimensions()*vf.dimensions()
444  );
445 
446  IOobject ddtIOobject
447  (
448  "ddt(" + rho.name() + ',' + vf.name() + ')',
449  mesh().time().timeName(),
450  mesh()
451  );
452 
453  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
454 
455  if (mesh().moving())
456  {
457  if (evaluate(ddt0))
458  {
459  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
460 
461  ddt0.primitiveFieldRef() =
462  (
463  rDtCoef0*rho.value()*
464  (
465  mesh().V0()*vf.oldTime().primitiveField()
466  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
467  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
468  )/mesh().V0();
469 
470  ddt0.boundaryFieldRef() =
471  (
472  rDtCoef0*rho.value()*
473  (
474  vf.oldTime().boundaryField()
475  - vf.oldTime().oldTime().boundaryField()
476  ) - offCentre_(ff(ddt0.boundaryField()))
477  );
478  }
479 
481  (
483  (
484  ddtIOobject,
485  mesh(),
486  rDtCoef.dimensions()*rho.dimensions()*vf.dimensions(),
487  (
488  rDtCoef.value()*rho.value()*
489  (
490  mesh().V()*vf.primitiveField()
491  - mesh().V0()*vf.oldTime().primitiveField()
492  ) - mesh().V0()*offCentre_(ddt0.primitiveField())
493  )/mesh().V(),
494  rDtCoef.value()*rho.value()*
495  (
496  vf.boundaryField() - vf.oldTime().boundaryField()
497  ) - offCentre_(ff(ddt0.boundaryField()))
498  )
499  );
500  }
501  else
502  {
503  if (evaluate(ddt0))
504  {
505  ddt0 = rDtCoef0_(ddt0)*rho*(vf.oldTime() - vf.oldTime().oldTime())
506  - offCentre_(ddt0());
507  }
508 
510  (
512  (
513  ddtIOobject,
514  rDtCoef*rho*(vf - vf.oldTime()) - offCentre_(ddt0())
515  )
516  );
517  }
518 }
519 
520 
521 template<class Type>
524 (
525  const volScalarField& rho,
527 )
528 {
529  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
530  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
531  (
532  "ddt0(" + rho.name() + ',' + vf.name() + ')',
533  rho.dimensions()*vf.dimensions()
534  );
535 
536  IOobject ddtIOobject
537  (
538  "ddt(" + rho.name() + ',' + vf.name() + ')',
539  mesh().time().timeName(),
540  mesh()
541  );
542 
543  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
544 
545  if (mesh().moving())
546  {
547  if (evaluate(ddt0))
548  {
549  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
550 
551  ddt0.primitiveFieldRef() =
552  (
553  rDtCoef0*
554  (
555  mesh().V0()*rho.oldTime().primitiveField()
556  *vf.oldTime().primitiveField()
557  - mesh().V00()*rho.oldTime().oldTime().primitiveField()
558  *vf.oldTime().oldTime().primitiveField()
559  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
560  )/mesh().V0();
561 
562  ddt0.boundaryFieldRef() =
563  (
564  rDtCoef0*
565  (
566  rho.oldTime().boundaryField()
567  *vf.oldTime().boundaryField()
568  - rho.oldTime().oldTime().boundaryField()
569  *vf.oldTime().oldTime().boundaryField()
570  ) - offCentre_(ff(ddt0.boundaryField()))
571  );
572  }
573 
575  (
577  (
578  ddtIOobject,
579  mesh(),
580  rDtCoef.dimensions()*rho.dimensions()*vf.dimensions(),
581  (
582  rDtCoef.value()*
583  (
584  mesh().V()*rho.primitiveField()*vf.primitiveField()
585  - mesh().V0()*rho.oldTime().primitiveField()
586  *vf.oldTime().primitiveField()
587  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
588  )/mesh().V(),
589  rDtCoef.value()*
590  (
591  rho.boundaryField()*vf.boundaryField()
592  - rho.oldTime().boundaryField()*vf.oldTime().boundaryField()
593  ) - offCentre_(ff(ddt0.boundaryField()))
594  )
595  );
596  }
597  else
598  {
599  if (evaluate(ddt0))
600  {
601  ddt0 = rDtCoef0_(ddt0)*
602  (
603  rho.oldTime()*vf.oldTime()
604  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
605  ) - offCentre_(ddt0());
606  }
607 
609  (
611  (
612  ddtIOobject,
613  rDtCoef*(rho*vf - rho.oldTime()*vf.oldTime())
614  - offCentre_(ddt0())
615  )
616  );
617  }
618 }
619 
620 
621 template<class Type>
624 (
625  const volScalarField& alpha,
626  const volScalarField& rho,
628 )
629 {
630  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
631  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
632  (
633  "ddt0(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
634  alpha.dimensions()*rho.dimensions()*vf.dimensions()
635  );
636 
637  IOobject ddtIOobject
638  (
639  "ddt(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
640  mesh().time().timeName(),
641  mesh()
642  );
643 
644  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
645 
646  if (mesh().moving())
647  {
648  if (evaluate(ddt0))
649  {
650  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
651 
652  ddt0.primitiveFieldRef() =
653  (
654  rDtCoef0*
655  (
656  mesh().V0()
657  *alpha.oldTime().primitiveField()
658  *rho.oldTime().primitiveField()
659  *vf.oldTime().primitiveField()
660 
661  - mesh().V00()
662  *alpha.oldTime().oldTime().primitiveField()
663  *rho.oldTime().oldTime().primitiveField()
664  *vf.oldTime().oldTime().primitiveField()
665  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
666  )/mesh().V0();
667 
668  ddt0.boundaryFieldRef() =
669  (
670  rDtCoef0*
671  (
672  alpha.oldTime().boundaryField()
673  *rho.oldTime().boundaryField()
674  *vf.oldTime().boundaryField()
675 
676  - alpha.oldTime().oldTime().boundaryField()
677  *rho.oldTime().oldTime().boundaryField()
678  *vf.oldTime().oldTime().boundaryField()
679  ) - offCentre_(ff(ddt0.boundaryField()))
680  );
681  }
682 
684  (
686  (
687  ddtIOobject,
688  mesh(),
689  rDtCoef.dimensions()
690  *alpha.dimensions()*rho.dimensions()*vf.dimensions(),
691  (
692  rDtCoef.value()*
693  (
694  mesh().V()
695  *alpha.primitiveField()
696  *rho.primitiveField()
697  *vf.primitiveField()
698 
699  - mesh().V0()
700  *alpha.oldTime().primitiveField()
701  *rho.oldTime().primitiveField()
702  *vf.oldTime().primitiveField()
703  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
704  )/mesh().V(),
705  rDtCoef.value()*
706  (
707  alpha.boundaryField()
708  *rho.boundaryField()
709  *vf.boundaryField()
710 
711  - alpha.oldTime().boundaryField()
712  *rho.oldTime().boundaryField()
713  *vf.oldTime().boundaryField()
714  ) - offCentre_(ff(ddt0.boundaryField()))
715  )
716  );
717  }
718  else
719  {
720  if (evaluate(ddt0))
721  {
722  ddt0 = rDtCoef0_(ddt0)*
723  (
724  alpha.oldTime()
725  *rho.oldTime()
726  *vf.oldTime()
727 
728  - alpha.oldTime().oldTime()
729  *rho.oldTime().oldTime()
730  *vf.oldTime().oldTime()
731  ) - offCentre_(ddt0());
732  }
733 
735  (
737  (
738  ddtIOobject,
739  rDtCoef
740  *(
741  alpha*rho*vf
742  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
743  )
744  - offCentre_(ddt0())
745  )
746  );
747  }
748 }
749 
750 
751 template<class Type>
754 (
756 )
757 {
758  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
759  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
760  (
761  "ddt0(" + vf.name() + ')',
762  vf.dimensions()
763  );
764 
765  tmp<fvMatrix<Type>> tfvm
766  (
767  new fvMatrix<Type>
768  (
769  vf,
771  )
772  );
773 
774  fvMatrix<Type>& fvm = tfvm.ref();
775 
776  scalar rDtCoef = rDtCoef_(ddt0).value();
777  fvm.diag() = rDtCoef*mesh().V();
778 
779  vf.oldTime().oldTime();
780 
781  if (mesh().moving())
782  {
783  if (evaluate(ddt0))
784  {
785  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
786 
787  ddt0.primitiveFieldRef() =
788  (
789  rDtCoef0*
790  (
791  mesh().V0()*vf.oldTime().primitiveField()
792  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
793  )
794  - mesh().V00()*offCentre_(ddt0.primitiveField())
795  )/mesh().V0();
796 
797  ddt0.boundaryFieldRef() =
798  (
799  rDtCoef0*
800  (
801  vf.oldTime().boundaryField()
802  - vf.oldTime().oldTime().boundaryField()
803  )
804  - offCentre_(ff(ddt0.boundaryField()))
805  );
806  }
807 
808  fvm.source() =
809  (
810  rDtCoef*vf.oldTime().primitiveField()
811  + offCentre_(ddt0.primitiveField())
812  )*mesh().V0();
813  }
814  else
815  {
816  if (evaluate(ddt0))
817  {
818  ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime())
819  - offCentre_(ddt0());
820  }
821 
822  fvm.source() =
823  (
824  rDtCoef*vf.oldTime().primitiveField()
825  + offCentre_(ddt0.primitiveField())
826  )*mesh().V();
827  }
828 
829  return tfvm;
830 }
831 
832 
833 template<class Type>
836 (
837  const dimensionedScalar& rho,
839 )
840 {
841  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
842  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
843  (
844  "ddt0(" + rho.name() + ',' + vf.name() + ')',
845  rho.dimensions()*vf.dimensions()
846  );
847 
848  tmp<fvMatrix<Type>> tfvm
849  (
850  new fvMatrix<Type>
851  (
852  vf,
854  )
855  );
856  fvMatrix<Type>& fvm = tfvm.ref();
857 
858  scalar rDtCoef = rDtCoef_(ddt0).value();
859  fvm.diag() = rDtCoef*rho.value()*mesh().V();
860 
861  vf.oldTime().oldTime();
862 
863  if (mesh().moving())
864  {
865  if (evaluate(ddt0))
866  {
867  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
868 
869  ddt0.primitiveFieldRef() =
870  (
871  rDtCoef0*rho.value()*
872  (
873  mesh().V0()*vf.oldTime().primitiveField()
874  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
875  )
876  - mesh().V00()*offCentre_(ddt0.primitiveField())
877  )/mesh().V0();
878 
879  ddt0.boundaryFieldRef() =
880  (
881  rDtCoef0*rho.value()*
882  (
883  vf.oldTime().boundaryField()
884  - vf.oldTime().oldTime().boundaryField()
885  )
886  - offCentre_(ff(ddt0.boundaryField()))
887  );
888  }
889 
890  fvm.source() =
891  (
892  rDtCoef*rho.value()*vf.oldTime().primitiveField()
893  + offCentre_(ddt0.primitiveField())
894  )*mesh().V0();
895  }
896  else
897  {
898  if (evaluate(ddt0))
899  {
900  ddt0 = rDtCoef0_(ddt0)*rho*(vf.oldTime() - vf.oldTime().oldTime())
901  - offCentre_(ddt0());
902  }
903 
904  fvm.source() =
905  (
906  rDtCoef*rho.value()*vf.oldTime().primitiveField()
907  + offCentre_(ddt0.primitiveField())
908  )*mesh().V();
909  }
910 
911  return tfvm;
912 }
913 
914 
915 template<class Type>
918 (
919  const volScalarField& rho,
921 )
922 {
923  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
924  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
925  (
926  "ddt0(" + rho.name() + ',' + vf.name() + ')',
927  rho.dimensions()*vf.dimensions()
928  );
929 
930  tmp<fvMatrix<Type>> tfvm
931  (
932  new fvMatrix<Type>
933  (
934  vf,
936  )
937  );
938  fvMatrix<Type>& fvm = tfvm.ref();
939 
940  scalar rDtCoef = rDtCoef_(ddt0).value();
941  fvm.diag() = rDtCoef*rho.primitiveField()*mesh().V();
942 
943  vf.oldTime().oldTime();
944  rho.oldTime().oldTime();
945 
946  if (mesh().moving())
947  {
948  if (evaluate(ddt0))
949  {
950  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
951 
952  ddt0.primitiveFieldRef() =
953  (
954  rDtCoef0*
955  (
956  mesh().V0()*rho.oldTime().primitiveField()
957  *vf.oldTime().primitiveField()
958  - mesh().V00()*rho.oldTime().oldTime().primitiveField()
959  *vf.oldTime().oldTime().primitiveField()
960  )
961  - mesh().V00()*offCentre_(ddt0.primitiveField())
962  )/mesh().V0();
963 
964  ddt0.boundaryFieldRef() =
965  (
966  rDtCoef0*
967  (
968  rho.oldTime().boundaryField()
969  *vf.oldTime().boundaryField()
970  - rho.oldTime().oldTime().boundaryField()
971  *vf.oldTime().oldTime().boundaryField()
972  )
973  - offCentre_(ff(ddt0.boundaryField()))
974  );
975  }
976 
977  fvm.source() =
978  (
979  rDtCoef*rho.oldTime().primitiveField()*vf.oldTime().primitiveField()
980  + offCentre_(ddt0.primitiveField())
981  )*mesh().V0();
982  }
983  else
984  {
985  if (evaluate(ddt0))
986  {
987  ddt0 = rDtCoef0_(ddt0)*
988  (
989  rho.oldTime()*vf.oldTime()
990  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
991  ) - offCentre_(ddt0());
992  }
993 
994  fvm.source() =
995  (
996  rDtCoef*rho.oldTime().primitiveField()*vf.oldTime().primitiveField()
997  + offCentre_(ddt0.primitiveField())
998  )*mesh().V();
999  }
1000 
1001  return tfvm;
1002 }
1003 
1004 
1005 template<class Type>
1009  const volScalarField& alpha,
1010  const volScalarField& rho,
1012 )
1013 {
1014  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1015  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1016  (
1017  "ddt0(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
1018  alpha.dimensions()*rho.dimensions()*vf.dimensions()
1019  );
1020 
1021  tmp<fvMatrix<Type>> tfvm
1022  (
1023  new fvMatrix<Type>
1024  (
1025  vf,
1026  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
1027  )
1028  );
1029  fvMatrix<Type>& fvm = tfvm.ref();
1030 
1031  scalar rDtCoef = rDtCoef_(ddt0).value();
1032  fvm.diag() = rDtCoef*alpha.primitiveField()*rho.primitiveField()*mesh().V();
1033 
1034  vf.oldTime().oldTime();
1035  alpha.oldTime().oldTime();
1036  rho.oldTime().oldTime();
1037 
1038  if (mesh().moving())
1039  {
1040  if (evaluate(ddt0))
1041  {
1042  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
1043 
1044  ddt0.primitiveFieldRef() =
1045  (
1046  rDtCoef0*
1047  (
1048  mesh().V0()
1049  *alpha.oldTime().primitiveField()
1050  *rho.oldTime().primitiveField()
1051  *vf.oldTime().primitiveField()
1052 
1053  - mesh().V00()
1054  *alpha.oldTime().oldTime().primitiveField()
1055  *rho.oldTime().oldTime().primitiveField()
1056  *vf.oldTime().oldTime().primitiveField()
1057  )
1058  - mesh().V00()*offCentre_(ddt0.primitiveField())
1059  )/mesh().V0();
1060 
1061  ddt0.boundaryFieldRef() =
1062  (
1063  rDtCoef0*
1064  (
1065  alpha.oldTime().boundaryField()
1066  *rho.oldTime().boundaryField()
1067  *vf.oldTime().boundaryField()
1068 
1069  - alpha.oldTime().oldTime().boundaryField()
1070  *rho.oldTime().oldTime().boundaryField()
1071  *vf.oldTime().oldTime().boundaryField()
1072  )
1073  - offCentre_(ff(ddt0.boundaryField()))
1074  );
1075  }
1076 
1077  fvm.source() =
1078  (
1079  rDtCoef
1080  *alpha.oldTime().primitiveField()
1081  *rho.oldTime().primitiveField()
1082  *vf.oldTime().primitiveField()
1083  + offCentre_(ddt0.primitiveField())
1084  )*mesh().V0();
1085  }
1086  else
1087  {
1088  if (evaluate(ddt0))
1089  {
1090  ddt0 = rDtCoef0_(ddt0)*
1091  (
1092  alpha.oldTime()
1093  *rho.oldTime()
1094  *vf.oldTime()
1095 
1096  - alpha.oldTime().oldTime()
1097  *rho.oldTime().oldTime()
1098  *vf.oldTime().oldTime()
1099  ) - offCentre_(ddt0());
1100  }
1101 
1102  fvm.source() =
1103  (
1104  rDtCoef
1105  *alpha.oldTime().primitiveField()
1106  *rho.oldTime().primitiveField()
1107  *vf.oldTime().primitiveField()
1108  + offCentre_(ddt0.primitiveField())
1109  )*mesh().V();
1110  }
1111 
1112  return tfvm;
1113 }
1114 
1115 
1116 template<class Type>
1122 )
1123 {
1124  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1125  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1126  (
1127  "ddt0(" + U.name() + ')',
1128  U.dimensions()
1129  );
1130 
1131  DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 =
1132  ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>>
1133  (
1134  "ddt0(" + Uf.name() + ')',
1135  Uf.dimensions()
1136  );
1137 
1138  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1139 
1140  if (evaluate(ddt0))
1141  {
1142  ddt0 =
1143  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1144  - offCentre_(ddt0());
1145  }
1146 
1147  if (evaluate(dUfdt0))
1148  {
1149  dUfdt0 =
1150  rDtCoef0_(dUfdt0)*(Uf.oldTime() - Uf.oldTime().oldTime())
1151  - offCentre_(dUfdt0());
1152  }
1153 
1154  return tmp<fluxFieldType>
1155  (
1156  new fluxFieldType
1157  (
1158  IOobject
1159  (
1160  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
1161  mesh().time().timeName(),
1162  mesh()
1163  ),
1164  this->fvcDdtPhiCoeff(U.oldTime(), mesh().Sf() & Uf.oldTime())
1165  *(
1166  mesh().Sf()
1167  & (
1168  (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1169  - fvc::interpolate(rDtCoef*U.oldTime() + offCentre_(ddt0()))
1170  )
1171  )
1172  )
1173  );
1174 }
1175 
1176 
1177 template<class Type>
1182  const fluxFieldType& phi
1183 )
1184 {
1185  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1186  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1187  (
1188  "ddt0(" + U.name() + ')',
1189  U.dimensions()
1190  );
1191 
1192  DDt0Field<fluxFieldType>& dphidt0 =
1193  ddt0_<fluxFieldType>
1194  (
1195  "ddt0(" + phi.name() + ')',
1196  phi.dimensions()
1197  );
1198 
1199  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1200 
1201  if (evaluate(ddt0))
1202  {
1203  ddt0 =
1204  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1205  - offCentre_(ddt0());
1206  }
1207 
1208  if (evaluate(dphidt0))
1209  {
1210  dphidt0 =
1211  rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime())
1212  - offCentre_(dphidt0());
1213  }
1214 
1215  return tmp<fluxFieldType>
1216  (
1217  new fluxFieldType
1218  (
1219  IOobject
1220  (
1221  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
1222  mesh().time().timeName(),
1223  mesh()
1224  ),
1225  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
1226  *(
1227  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1229  (
1230  mesh().Sf(),
1231  rDtCoef*U.oldTime() + offCentre_(ddt0())
1232  )
1233  )
1234  )
1235  );
1236 }
1237 
1238 
1239 template<class Type>
1243  const volScalarField& rho,
1246 )
1247 {
1248  if
1249  (
1250  U.dimensions() == dimVelocity
1251  && Uf.dimensions() == rho.dimensions()*dimVelocity
1252  )
1253  {
1254  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1255  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1256  (
1257  "ddt0(" + rho.name() + ',' + U.name() + ')',
1258  U.dimensions()
1259  );
1260 
1261  DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 =
1262  ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>>
1263  (
1264  "ddt0(" + Uf.name() + ')',
1265  Uf.dimensions()
1266  );
1267 
1268  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1269 
1271  (
1272  rho.oldTime()*U.oldTime()
1273  );
1274 
1275  if (evaluate(ddt0))
1276  {
1277  ddt0 =
1278  rDtCoef0_(ddt0)
1279  *(rhoU0 - rho.oldTime().oldTime()*U.oldTime().oldTime())
1280  - offCentre_(ddt0());
1281  }
1282 
1283  if (evaluate(dUfdt0))
1284  {
1285  dUfdt0 =
1286  rDtCoef0_(dUfdt0)
1287  *(Uf.oldTime() - Uf.oldTime().oldTime())
1288  - offCentre_(dUfdt0());
1289  }
1290 
1292  (
1293  new fluxFieldType
1294  (
1295  IOobject
1296  (
1297  "ddtCorr("
1298  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
1299  mesh().time().timeName(),
1300  mesh()
1301  ),
1302  this->fvcDdtPhiCoeff(rhoU0, mesh().Sf() & Uf.oldTime())
1303  *(
1304  mesh().Sf()
1305  & (
1306  (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1307  - fvc::interpolate(rDtCoef*rhoU0 + offCentre_(ddt0()))
1308  )
1309  )
1310  )
1311  );
1312 
1313  return ddtCorr;
1314  }
1315  else if
1316  (
1317  U.dimensions() == rho.dimensions()*dimVelocity
1318  && Uf.dimensions() == rho.dimensions()*dimVelocity
1319  )
1320  {
1321  return fvcDdtUfCorr(U, Uf);
1322  }
1323  else
1324  {
1326  << "dimensions of Uf are not correct"
1327  << abort(FatalError);
1328 
1329  return fluxFieldType::null();
1330  }
1331 }
1332 
1333 
1334 template<class Type>
1338  const volScalarField& rho,
1340  const fluxFieldType& phi
1341 )
1342 {
1343  if
1344  (
1345  U.dimensions() == dimVelocity
1346  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1347  )
1348  {
1349  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1350  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1351  (
1352  "ddt0(" + rho.name() + ',' + U.name() + ')',
1353  U.dimensions()
1354  );
1355 
1356  DDt0Field<fluxFieldType>& dphidt0 =
1357  ddt0_<fluxFieldType>
1358  (
1359  "ddt0(" + phi.name() + ')',
1360  phi.dimensions()
1361  );
1362 
1363  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1364 
1366  (
1367  rho.oldTime()*U.oldTime()
1368  );
1369 
1370  if (evaluate(ddt0))
1371  {
1372  ddt0 =
1373  rDtCoef0_(ddt0)
1374  *(rhoU0 - rho.oldTime().oldTime()*U.oldTime().oldTime())
1375  - offCentre_(ddt0());
1376  }
1377 
1378  if (evaluate(dphidt0))
1379  {
1380  dphidt0 =
1381  rDtCoef0_(dphidt0)
1382  *(phi.oldTime() - phi.oldTime().oldTime())
1383  - offCentre_(dphidt0());
1384  }
1385 
1387  (
1388  new fluxFieldType
1389  (
1390  IOobject
1391  (
1392  "ddtCorr("
1393  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
1394  mesh().time().timeName(),
1395  mesh()
1396  ),
1397  this->fvcDdtPhiCoeff(rhoU0, phi.oldTime())
1398  *(
1399  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1401  (
1402  mesh().Sf(),
1403  rDtCoef*rhoU0 + offCentre_(ddt0())
1404  )
1405  )
1406  )
1407  );
1408 
1409  return ddtCorr;
1410  }
1411  else if
1412  (
1413  U.dimensions() == rho.dimensions()*dimVelocity
1414  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1415  )
1416  {
1417  return fvcDdtPhiCorr(U, phi);
1418  }
1419  else
1420  {
1422  << "dimensions of phi are not correct"
1423  << abort(FatalError);
1424 
1425  return fluxFieldType::null();
1426  }
1427 }
1428 
1429 
1430 template<class Type>
1434 )
1435 {
1436  DDt0Field<surfaceScalarField>& meshPhi0 = ddt0_<surfaceScalarField>
1437  (
1438  "meshPhiCN_0",
1439  dimVolume
1440  );
1441 
1442  if (evaluate(meshPhi0))
1443  {
1444  meshPhi0 =
1445  coef0_(meshPhi0)*mesh().phi().oldTime() - offCentre_(meshPhi0());
1446  }
1447 
1449  (
1450  new surfaceScalarField
1451  (
1452  IOobject
1453  (
1454  mesh().phi().name(),
1455  mesh().time().timeName(),
1456  mesh(),
1459  false
1460  ),
1461  coef_(meshPhi0)*mesh().phi() - offCentre_(meshPhi0())
1462  )
1463  );
1464 }
1465 
1466 
1467 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1468 
1469 } // End namespace fv
1470 
1471 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1472 
1473 } // End namespace Foam
1474 
1475 // ************************************************************************* //
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
label timeIndex() const
Return current time index.
Definition: TimeStateI.H:35
surfaceScalarField & phi
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const surfaceVectorField & Sf() const
Return cell face area vectors.
dimensionedScalar deltaT0() const
Return old time step.
Definition: TimeStateI.H:59
const dimensionSet & dimensions() const
Return const reference to dimensions.
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const surfaceScalarField & phi() const
Return cell face motion fluxes.
static tmp< GeometricField< typename innerProduct< vector, Type >::type, fvsPatchField, surfaceMesh > > dotInterpolate(const surfaceVectorField &Sf, const GeometricField< Type, fvPatchField, volMesh > &tvf)
Interpolate field onto faces.
tmp< GeometricField< typename flux< Type >::type, fvsPatchField, surfaceMesh > > ddtCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition: fvcDdt.C:155
const DimensionedField< scalar, volMesh > & V0() const
Return old-time cell volumes.
const DimensionedField< scalar, volMesh > & V00() const
Return old-old-time cell volumes.
const dimensionSet dimVol(dimVolume)
Definition: dimensionSets.H:59
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
const Type & value() const
Return const reference to value.
Generic GeometricField class.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:715
Generic dimensioned Type class.
Generic field type.
Definition: FieldField.H:51
const word & name() const
Return const reference to name.
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
ddtScheme< Type >::fluxFieldType fluxFieldType
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
labelList fv(nPoints)
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:71
word timeName
Definition: getTimeIndex.H:3
static const zero Zero
Definition: zero.H:91
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const dimensionSet & dimensions() const
Return dimensions.
Calculate the divergence of the given field.
void store()
Transfer ownership of this object to its registry.
Definition: regIOobjectI.H:34
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
Field< Type > & source()
Definition: fvMatrix.H:291
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
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.
static const GeometricField< Type, PatchField, GeoMesh > & null()
Return a null geometric field.
tmp< surfaceScalarField > fvcDdtPhiCoeff(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi, const fluxFieldType &phiCorr)
Definition: ddtScheme.C:137
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
dimensionedScalar deltaT() const
Return time step.
Definition: TimeStateI.H:53
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label timeIndex
Definition: getTimeIndex.H:4
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
A special matrix type and solver, designed for finite volume solutions of scalar equations.
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
const fvMesh & mesh() const
Return mesh reference.
A class for managing temporary objects.
Definition: PtrList.H:54
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
scalarField & diag()
Definition: lduMatrix.C:183
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
const word & name() const
Return name.
Definition: IOobject.H:260
const FieldField< fvPatchField, Type > & ff(const FieldField< fvPatchField, Type > &bf)
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Namespace for OpenFOAM.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
const dimensionSet dimVelocity