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-2015 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>
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>
64 (
65  const IOobject& io,
66  const fvMesh& mesh,
68 )
69 :
70  GeoField(io, mesh, dimType),
71  startTimeIndex_(mesh.time().timeIndex())
72 {}
73 
74 
75 template<class Type>
76 template<class GeoField>
79 {
80  return startTimeIndex_;
81 }
82 
83 
84 template<class Type>
85 template<class GeoField>
88 {
89  return *this;
90 }
91 
92 
93 template<class Type>
94 template<class GeoField>
96 operator=(const GeoField& gf)
97 {
98  GeoField::operator=(gf);
99 }
100 
101 
102 template<class Type>
103 template<class GeoField>
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  (
133  (
134  IOobject
135  (
136  name,
137  startTimeName,
138  mesh(),
141  ),
142  mesh()
143  )
144  );
145  }
146  else
147  {
149  (
151  (
152  IOobject
153  (
154  name,
155  mesh().time().timeName(),
156  mesh(),
159  ),
160  mesh(),
162  (
163  "0",
164  dims/dimTime,
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>
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>
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>
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>
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>
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>
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>
272 (
274 )
275 {
276  return bf;
277 }
278 
279 
280 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
281 
282 template<class Type>
285 (
286  const dimensioned<Type>& dt
287 )
288 {
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,
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.dimensionedInternalField() =
327  (
328  (rDtCoef0*dt)*(mesh().V0() - mesh().V00())
329  - mesh().V00()*offCentre_(ddt0.dimensionedInternalField())
330  )/mesh().V0();
331  }
332 
333  tdtdt().dimensionedInternalField() =
334  (
335  (rDtCoef*dt)*(mesh().V() - mesh().V0())
336  - mesh().V0()*offCentre_(ddt0.dimensionedInternalField())
337  )/mesh().V();
338  }
339 
340  return tdtdt;
341 }
342 
343 
344 template<class Type>
347 (
349 )
350 {
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.internalField() =
374  (
375  rDtCoef0*
376  (
377  mesh().V0()*vf.oldTime().internalField()
378  - mesh().V00()*vf.oldTime().oldTime().internalField()
379  ) - mesh().V00()*offCentre_(ddt0.internalField())
380  )/mesh().V0();
381 
382  ddt0.boundaryField() =
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  mesh(),
398  rDtCoef.dimensions()*vf.dimensions(),
399  (
400  rDtCoef.value()*
401  (
402  mesh().V()*vf.internalField()
403  - mesh().V0()*vf.oldTime().internalField()
404  ) - mesh().V0()*offCentre_(ddt0.internalField())
405  )/mesh().V(),
406  rDtCoef.value()*
407  (
408  vf.boundaryField() - vf.oldTime().boundaryField()
409  ) - offCentre_(ff(ddt0.boundaryField()))
410  )
411  );
412  }
413  else
414  {
415  if (evaluate(ddt0))
416  {
417  ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime())
418  - offCentre_(ddt0());
419  }
420 
422  (
424  (
425  ddtIOobject,
426  rDtCoef*(vf - vf.oldTime()) - offCentre_(ddt0())
427  )
428  );
429  }
430 }
431 
432 
433 template<class Type>
436 (
437  const dimensionedScalar& rho,
439 )
440 {
442  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
443  (
444  "ddt0(" + rho.name() + ',' + vf.name() + ')',
445  rho.dimensions()*vf.dimensions()
446  );
447 
448  IOobject ddtIOobject
449  (
450  "ddt(" + rho.name() + ',' + vf.name() + ')',
451  mesh().time().timeName(),
452  mesh()
453  );
454 
455  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
456 
457  if (mesh().moving())
458  {
459  if (evaluate(ddt0))
460  {
461  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
462 
463  ddt0.internalField() =
464  (
465  rDtCoef0*rho.value()*
466  (
467  mesh().V0()*vf.oldTime().internalField()
468  - mesh().V00()*vf.oldTime().oldTime().internalField()
469  ) - mesh().V00()*offCentre_(ddt0.internalField())
470  )/mesh().V0();
471 
472  ddt0.boundaryField() =
473  (
474  rDtCoef0*rho.value()*
475  (
476  vf.oldTime().boundaryField()
477  - vf.oldTime().oldTime().boundaryField()
478  ) - offCentre_(ff(ddt0.boundaryField()))
479  );
480  }
481 
483  (
485  (
486  ddtIOobject,
487  mesh(),
488  rDtCoef.dimensions()*rho.dimensions()*vf.dimensions(),
489  (
490  rDtCoef.value()*rho.value()*
491  (
492  mesh().V()*vf.internalField()
493  - mesh().V0()*vf.oldTime().internalField()
494  ) - mesh().V0()*offCentre_(ddt0.internalField())
495  )/mesh().V(),
496  rDtCoef.value()*rho.value()*
497  (
498  vf.boundaryField() - vf.oldTime().boundaryField()
499  ) - offCentre_(ff(ddt0.boundaryField()))
500  )
501  );
502  }
503  else
504  {
505  if (evaluate(ddt0))
506  {
507  ddt0 = rDtCoef0_(ddt0)*rho*(vf.oldTime() - vf.oldTime().oldTime())
508  - offCentre_(ddt0());
509  }
510 
512  (
514  (
515  ddtIOobject,
516  rDtCoef*rho*(vf - vf.oldTime()) - offCentre_(ddt0())
517  )
518  );
519  }
520 }
521 
522 
523 template<class Type>
526 (
527  const volScalarField& rho,
529 )
530 {
532  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
533  (
534  "ddt0(" + rho.name() + ',' + vf.name() + ')',
535  rho.dimensions()*vf.dimensions()
536  );
537 
538  IOobject ddtIOobject
539  (
540  "ddt(" + rho.name() + ',' + vf.name() + ')',
541  mesh().time().timeName(),
542  mesh()
543  );
544 
545  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
546 
547  if (mesh().moving())
548  {
549  if (evaluate(ddt0))
550  {
551  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
552 
553  ddt0.internalField() =
554  (
555  rDtCoef0*
556  (
557  mesh().V0()*rho.oldTime().internalField()
558  *vf.oldTime().internalField()
559  - mesh().V00()*rho.oldTime().oldTime().internalField()
560  *vf.oldTime().oldTime().internalField()
561  ) - mesh().V00()*offCentre_(ddt0.internalField())
562  )/mesh().V0();
563 
564  ddt0.boundaryField() =
565  (
566  rDtCoef0*
567  (
568  rho.oldTime().boundaryField()
569  *vf.oldTime().boundaryField()
570  - rho.oldTime().oldTime().boundaryField()
571  *vf.oldTime().oldTime().boundaryField()
572  ) - offCentre_(ff(ddt0.boundaryField()))
573  );
574  }
575 
577  (
579  (
580  ddtIOobject,
581  mesh(),
582  rDtCoef.dimensions()*rho.dimensions()*vf.dimensions(),
583  (
584  rDtCoef.value()*
585  (
586  mesh().V()*rho.internalField()*vf.internalField()
587  - mesh().V0()*rho.oldTime().internalField()
588  *vf.oldTime().internalField()
589  ) - mesh().V00()*offCentre_(ddt0.internalField())
590  )/mesh().V(),
591  rDtCoef.value()*
592  (
593  rho.boundaryField()*vf.boundaryField()
594  - rho.oldTime().boundaryField()*vf.oldTime().boundaryField()
595  ) - offCentre_(ff(ddt0.boundaryField()))
596  )
597  );
598  }
599  else
600  {
601  if (evaluate(ddt0))
602  {
603  ddt0 = rDtCoef0_(ddt0)*
604  (
605  rho.oldTime()*vf.oldTime()
606  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
607  ) - offCentre_(ddt0());
608  }
609 
611  (
613  (
614  ddtIOobject,
615  rDtCoef*(rho*vf - rho.oldTime()*vf.oldTime())
616  - offCentre_(ddt0())
617  )
618  );
619  }
620 }
621 
622 
623 template<class Type>
626 (
627  const volScalarField& alpha,
628  const volScalarField& rho,
630 )
631 {
633  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
634  (
635  "ddt0(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
636  alpha.dimensions()*rho.dimensions()*vf.dimensions()
637  );
638 
639  IOobject ddtIOobject
640  (
641  "ddt(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
642  mesh().time().timeName(),
643  mesh()
644  );
645 
646  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
647 
648  if (mesh().moving())
649  {
650  if (evaluate(ddt0))
651  {
652  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
653 
654  ddt0.internalField() =
655  (
656  rDtCoef0*
657  (
658  mesh().V0()
659  *alpha.oldTime().internalField()
660  *rho.oldTime().internalField()
661  *vf.oldTime().internalField()
662 
663  - mesh().V00()
664  *alpha.oldTime().oldTime().internalField()
665  *rho.oldTime().oldTime().internalField()
666  *vf.oldTime().oldTime().internalField()
667  ) - mesh().V00()*offCentre_(ddt0.internalField())
668  )/mesh().V0();
669 
670  ddt0.boundaryField() =
671  (
672  rDtCoef0*
673  (
674  alpha.oldTime().boundaryField()
675  *rho.oldTime().boundaryField()
676  *vf.oldTime().boundaryField()
677 
678  - alpha.oldTime().oldTime().boundaryField()
679  *rho.oldTime().oldTime().boundaryField()
680  *vf.oldTime().oldTime().boundaryField()
681  ) - offCentre_(ff(ddt0.boundaryField()))
682  );
683  }
684 
686  (
688  (
689  ddtIOobject,
690  mesh(),
691  rDtCoef.dimensions()
692  *alpha.dimensions()*rho.dimensions()*vf.dimensions(),
693  (
694  rDtCoef.value()*
695  (
696  mesh().V()
697  *alpha.internalField()
698  *rho.internalField()
699  *vf.internalField()
700 
701  - mesh().V0()
702  *alpha.oldTime().internalField()
703  *rho.oldTime().internalField()
704  *vf.oldTime().internalField()
705  ) - mesh().V00()*offCentre_(ddt0.internalField())
706  )/mesh().V(),
707  rDtCoef.value()*
708  (
709  alpha.boundaryField()
710  *rho.boundaryField()
711  *vf.boundaryField()
712 
713  - alpha.oldTime().boundaryField()
714  *rho.oldTime().boundaryField()
715  *vf.oldTime().boundaryField()
716  ) - offCentre_(ff(ddt0.boundaryField()))
717  )
718  );
719  }
720  else
721  {
722  if (evaluate(ddt0))
723  {
724  ddt0 = rDtCoef0_(ddt0)*
725  (
726  alpha.oldTime()
727  *rho.oldTime()
728  *vf.oldTime()
729 
730  - alpha.oldTime().oldTime()
731  *rho.oldTime().oldTime()
732  *vf.oldTime().oldTime()
733  ) - offCentre_(ddt0());
734  }
735 
737  (
739  (
740  ddtIOobject,
741  rDtCoef
742  *(
743  alpha*rho*vf
744  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
745  )
746  - offCentre_(ddt0())
747  )
748  );
749  }
750 }
751 
752 
753 template<class Type>
756 (
758 )
759 {
761  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
762  (
763  "ddt0(" + vf.name() + ')',
764  vf.dimensions()
765  );
766 
767  tmp<fvMatrix<Type> > tfvm
768  (
769  new fvMatrix<Type>
770  (
771  vf,
772  vf.dimensions()*dimVol/dimTime
773  )
774  );
775 
776  fvMatrix<Type>& fvm = tfvm();
777 
778  scalar rDtCoef = rDtCoef_(ddt0).value();
779  fvm.diag() = rDtCoef*mesh().V();
780 
781  vf.oldTime().oldTime();
782 
783  if (mesh().moving())
784  {
785  if (evaluate(ddt0))
786  {
787  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
788 
789  ddt0.internalField() =
790  (
791  rDtCoef0*
792  (
793  mesh().V0()*vf.oldTime().internalField()
794  - mesh().V00()*vf.oldTime().oldTime().internalField()
795  )
796  - mesh().V00()*offCentre_(ddt0.internalField())
797  )/mesh().V0();
798 
799  ddt0.boundaryField() =
800  (
801  rDtCoef0*
802  (
803  vf.oldTime().boundaryField()
804  - vf.oldTime().oldTime().boundaryField()
805  )
806  - offCentre_(ff(ddt0.boundaryField()))
807  );
808  }
809 
810  fvm.source() =
811  (
812  rDtCoef*vf.oldTime().internalField()
813  + offCentre_(ddt0.internalField())
814  )*mesh().V0();
815  }
816  else
817  {
818  if (evaluate(ddt0))
819  {
820  ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime())
821  - offCentre_(ddt0());
822  }
823 
824  fvm.source() =
825  (
826  rDtCoef*vf.oldTime().internalField()
827  + offCentre_(ddt0.internalField())
828  )*mesh().V();
829  }
830 
831  return tfvm;
832 }
833 
834 
835 template<class Type>
838 (
839  const dimensionedScalar& rho,
841 )
842 {
844  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
845  (
846  "ddt0(" + rho.name() + ',' + vf.name() + ')',
847  rho.dimensions()*vf.dimensions()
848  );
849 
850  tmp<fvMatrix<Type> > tfvm
851  (
852  new fvMatrix<Type>
853  (
854  vf,
855  rho.dimensions()*vf.dimensions()*dimVol/dimTime
856  )
857  );
858  fvMatrix<Type>& fvm = tfvm();
859 
860  scalar rDtCoef = rDtCoef_(ddt0).value();
861  fvm.diag() = rDtCoef*rho.value()*mesh().V();
862 
863  vf.oldTime().oldTime();
864 
865  if (mesh().moving())
866  {
867  if (evaluate(ddt0))
868  {
869  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
870 
871  ddt0.internalField() =
872  (
873  rDtCoef0*rho.value()*
874  (
875  mesh().V0()*vf.oldTime().internalField()
876  - mesh().V00()*vf.oldTime().oldTime().internalField()
877  )
878  - mesh().V00()*offCentre_(ddt0.internalField())
879  )/mesh().V0();
880 
881  ddt0.boundaryField() =
882  (
883  rDtCoef0*rho.value()*
884  (
885  vf.oldTime().boundaryField()
886  - vf.oldTime().oldTime().boundaryField()
887  )
888  - offCentre_(ff(ddt0.boundaryField()))
889  );
890  }
891 
892  fvm.source() =
893  (
894  rDtCoef*rho.value()*vf.oldTime().internalField()
895  + offCentre_(ddt0.internalField())
896  )*mesh().V0();
897  }
898  else
899  {
900  if (evaluate(ddt0))
901  {
902  ddt0 = rDtCoef0_(ddt0)*rho*(vf.oldTime() - vf.oldTime().oldTime())
903  - offCentre_(ddt0());
904  }
905 
906  fvm.source() =
907  (
908  rDtCoef*rho.value()*vf.oldTime().internalField()
909  + offCentre_(ddt0.internalField())
910  )*mesh().V();
911  }
912 
913  return tfvm;
914 }
915 
916 
917 template<class Type>
920 (
921  const volScalarField& rho,
923 )
924 {
926  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
927  (
928  "ddt0(" + rho.name() + ',' + vf.name() + ')',
929  rho.dimensions()*vf.dimensions()
930  );
931 
932  tmp<fvMatrix<Type> > tfvm
933  (
934  new fvMatrix<Type>
935  (
936  vf,
937  rho.dimensions()*vf.dimensions()*dimVol/dimTime
938  )
939  );
940  fvMatrix<Type>& fvm = tfvm();
941 
942  scalar rDtCoef = rDtCoef_(ddt0).value();
943  fvm.diag() = rDtCoef*rho.internalField()*mesh().V();
944 
945  vf.oldTime().oldTime();
946  rho.oldTime().oldTime();
947 
948  if (mesh().moving())
949  {
950  if (evaluate(ddt0))
951  {
952  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
953 
954  ddt0.internalField() =
955  (
956  rDtCoef0*
957  (
958  mesh().V0()*rho.oldTime().internalField()
959  *vf.oldTime().internalField()
960  - mesh().V00()*rho.oldTime().oldTime().internalField()
961  *vf.oldTime().oldTime().internalField()
962  )
963  - mesh().V00()*offCentre_(ddt0.internalField())
964  )/mesh().V0();
965 
966  ddt0.boundaryField() =
967  (
968  rDtCoef0*
969  (
970  rho.oldTime().boundaryField()
971  *vf.oldTime().boundaryField()
972  - rho.oldTime().oldTime().boundaryField()
973  *vf.oldTime().oldTime().boundaryField()
974  )
975  - offCentre_(ff(ddt0.boundaryField()))
976  );
977  }
978 
979  fvm.source() =
980  (
981  rDtCoef*rho.oldTime().internalField()*vf.oldTime().internalField()
982  + offCentre_(ddt0.internalField())
983  )*mesh().V0();
984  }
985  else
986  {
987  if (evaluate(ddt0))
988  {
989  ddt0 = rDtCoef0_(ddt0)*
990  (
991  rho.oldTime()*vf.oldTime()
992  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
993  ) - offCentre_(ddt0());
994  }
995 
996  fvm.source() =
997  (
998  rDtCoef*rho.oldTime().internalField()*vf.oldTime().internalField()
999  + offCentre_(ddt0.internalField())
1000  )*mesh().V();
1001  }
1002 
1003  return tfvm;
1004 }
1005 
1006 
1007 template<class Type>
1011  const volScalarField& alpha,
1012  const volScalarField& rho,
1014 )
1015 {
1017  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
1018  (
1019  "ddt0(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
1020  alpha.dimensions()*rho.dimensions()*vf.dimensions()
1021  );
1022 
1023  tmp<fvMatrix<Type> > tfvm
1024  (
1025  new fvMatrix<Type>
1026  (
1027  vf,
1028  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
1029  )
1030  );
1031  fvMatrix<Type>& fvm = tfvm();
1032 
1033  scalar rDtCoef = rDtCoef_(ddt0).value();
1034  fvm.diag() = rDtCoef*alpha.internalField()*rho.internalField()*mesh().V();
1035 
1036  vf.oldTime().oldTime();
1037  alpha.oldTime().oldTime();
1038  rho.oldTime().oldTime();
1039 
1040  if (mesh().moving())
1041  {
1042  if (evaluate(ddt0))
1043  {
1044  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
1045 
1046  ddt0.internalField() =
1047  (
1048  rDtCoef0*
1049  (
1050  mesh().V0()
1051  *alpha.oldTime().internalField()
1052  *rho.oldTime().internalField()
1053  *vf.oldTime().internalField()
1054 
1055  - mesh().V00()
1056  *alpha.oldTime().oldTime().internalField()
1057  *rho.oldTime().oldTime().internalField()
1058  *vf.oldTime().oldTime().internalField()
1059  )
1060  - mesh().V00()*offCentre_(ddt0.internalField())
1061  )/mesh().V0();
1062 
1063  ddt0.boundaryField() =
1064  (
1065  rDtCoef0*
1066  (
1067  alpha.oldTime().boundaryField()
1068  *rho.oldTime().boundaryField()
1069  *vf.oldTime().boundaryField()
1070 
1071  - alpha.oldTime().oldTime().boundaryField()
1072  *rho.oldTime().oldTime().boundaryField()
1073  *vf.oldTime().oldTime().boundaryField()
1074  )
1075  - offCentre_(ff(ddt0.boundaryField()))
1076  );
1077  }
1078 
1079  fvm.source() =
1080  (
1081  rDtCoef
1082  *alpha.oldTime().internalField()
1083  *rho.oldTime().internalField()
1084  *vf.oldTime().internalField()
1085  + offCentre_(ddt0.internalField())
1086  )*mesh().V0();
1087  }
1088  else
1089  {
1090  if (evaluate(ddt0))
1091  {
1092  ddt0 = rDtCoef0_(ddt0)*
1093  (
1094  alpha.oldTime()
1095  *rho.oldTime()
1096  *vf.oldTime()
1097 
1098  - alpha.oldTime().oldTime()
1099  *rho.oldTime().oldTime()
1100  *vf.oldTime().oldTime()
1101  ) - offCentre_(ddt0());
1102  }
1103 
1104  fvm.source() =
1105  (
1106  rDtCoef
1107  *alpha.oldTime().internalField()
1108  *rho.oldTime().internalField()
1109  *vf.oldTime().internalField()
1110  + offCentre_(ddt0.internalField())
1111  )*mesh().V();
1112  }
1113 
1114  return tfvm;
1115 }
1116 
1117 
1118 template<class Type>
1124 )
1125 {
1127  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
1128  (
1129  "ddt0(" + U.name() + ')',
1130  U.dimensions()
1131  );
1132 
1134  ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh> >
1135  (
1136  "ddt0(" + Uf.name() + ')',
1137  Uf.dimensions()
1138  );
1139 
1140  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1141 
1142  if (evaluate(ddt0))
1143  {
1144  ddt0 =
1145  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1146  - offCentre_(ddt0());
1147  }
1148 
1149  if (evaluate(dUfdt0))
1150  {
1151  dUfdt0 =
1152  rDtCoef0_(dUfdt0)*(Uf.oldTime() - Uf.oldTime().oldTime())
1153  - offCentre_(dUfdt0());
1154  }
1155 
1156  return tmp<fluxFieldType>
1157  (
1158  new fluxFieldType
1159  (
1160  IOobject
1161  (
1162  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
1163  mesh().time().timeName(),
1164  mesh()
1165  ),
1166  this->fvcDdtPhiCoeff(U.oldTime(), mesh().Sf() & Uf.oldTime())
1167  *(
1168  mesh().Sf()
1169  & (
1170  (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1171  - fvc::interpolate(rDtCoef*U.oldTime() + offCentre_(ddt0()))
1172  )
1173  )
1174  )
1175  );
1176 }
1177 
1178 
1179 template<class Type>
1184  const fluxFieldType& phi
1185 )
1186 {
1188  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
1189  (
1190  "ddt0(" + U.name() + ')',
1191  U.dimensions()
1192  );
1193 
1194  DDt0Field<fluxFieldType>& dphidt0 =
1195  ddt0_<fluxFieldType>
1196  (
1197  "ddt0(" + phi.name() + ')',
1198  phi.dimensions()
1199  );
1200 
1201  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1202 
1203  if (evaluate(ddt0))
1204  {
1205  ddt0 =
1206  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1207  - offCentre_(ddt0());
1208  }
1209 
1210  if (evaluate(dphidt0))
1211  {
1212  dphidt0 =
1213  rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime())
1214  - offCentre_(dphidt0());
1215  }
1216 
1217  return tmp<fluxFieldType>
1218  (
1219  new fluxFieldType
1220  (
1221  IOobject
1222  (
1223  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
1224  mesh().time().timeName(),
1225  mesh()
1226  ),
1227  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
1228  *(
1229  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1230  - (
1231  mesh().Sf()
1232  & fvc::interpolate(rDtCoef*U.oldTime() + offCentre_(ddt0()))
1233  )
1234  )
1235  )
1236  );
1237 }
1238 
1239 
1240 template<class Type>
1244  const volScalarField& rho,
1247 )
1248 {
1249  if
1250  (
1251  U.dimensions() == dimVelocity
1252  && Uf.dimensions() == rho.dimensions()*dimVelocity
1253  )
1254  {
1256  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
1257  (
1258  "ddt0(" + rho.name() + ',' + U.name() + ')',
1259  U.dimensions()
1260  );
1261 
1263  ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh> >
1264  (
1265  "ddt0(" + Uf.name() + ')',
1266  Uf.dimensions()
1267  );
1268 
1269  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1270 
1272  (
1273  rho.oldTime()*U.oldTime()
1274  );
1275 
1276  if (evaluate(ddt0))
1277  {
1278  ddt0 =
1279  rDtCoef0_(ddt0)
1280  *(rhoU0 - rho.oldTime().oldTime()*U.oldTime().oldTime())
1281  - offCentre_(ddt0());
1282  }
1283 
1284  if (evaluate(dUfdt0))
1285  {
1286  dUfdt0 =
1287  rDtCoef0_(dUfdt0)
1288  *(Uf.oldTime() - Uf.oldTime().oldTime())
1289  - offCentre_(dUfdt0());
1290  }
1291 
1293  (
1294  new fluxFieldType
1295  (
1296  IOobject
1297  (
1298  "ddtCorr("
1299  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
1300  mesh().time().timeName(),
1301  mesh()
1302  ),
1303  this->fvcDdtPhiCoeff(rhoU0, mesh().Sf() & Uf.oldTime())
1304  *(
1305  mesh().Sf()
1306  & (
1307  (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1308  - fvc::interpolate(rDtCoef*rhoU0 + offCentre_(ddt0()))
1309  )
1310  )
1311  )
1312  );
1313 
1314  return ddtCorr;
1315  }
1316  else if
1317  (
1318  U.dimensions() == rho.dimensions()*dimVelocity
1319  && Uf.dimensions() == rho.dimensions()*dimVelocity
1320  )
1321  {
1322  return fvcDdtUfCorr(U, Uf);
1323  }
1324  else
1325  {
1327  << "dimensions of Uf are not correct"
1328  << abort(FatalError);
1329 
1330  return fluxFieldType::null();
1331  }
1332 }
1333 
1334 
1335 template<class Type>
1339  const volScalarField& rho,
1341  const fluxFieldType& phi
1342 )
1343 {
1344  if
1345  (
1346  U.dimensions() == dimVelocity
1347  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1348  )
1349  {
1351  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
1352  (
1353  "ddt0(" + rho.name() + ',' + U.name() + ')',
1354  U.dimensions()
1355  );
1356 
1357  DDt0Field<fluxFieldType>& dphidt0 =
1358  ddt0_<fluxFieldType>
1359  (
1360  "ddt0(" + phi.name() + ')',
1361  phi.dimensions()
1362  );
1363 
1364  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1365 
1367  (
1368  rho.oldTime()*U.oldTime()
1369  );
1370 
1371  if (evaluate(ddt0))
1372  {
1373  ddt0 =
1374  rDtCoef0_(ddt0)
1375  *(rhoU0 - rho.oldTime().oldTime()*U.oldTime().oldTime())
1376  - offCentre_(ddt0());
1377  }
1378 
1379  if (evaluate(dphidt0))
1380  {
1381  dphidt0 =
1382  rDtCoef0_(dphidt0)
1383  *(phi.oldTime() - phi.oldTime().oldTime())
1384  - offCentre_(dphidt0());
1385  }
1386 
1388  (
1389  new fluxFieldType
1390  (
1391  IOobject
1392  (
1393  "ddtCorr("
1394  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
1395  mesh().time().timeName(),
1396  mesh()
1397  ),
1398  this->fvcDdtPhiCoeff(rhoU0, phi.oldTime())
1399  *(
1400  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1401  - (
1402  mesh().Sf()
1403  & fvc::interpolate(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 // ************************************************************************* //
Foam::IOobject
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
Foam::word
A class for handling words, derived from string.
Definition: word.H:59
Foam::FieldField
Generic field type.
Definition: FieldField.H:51
Foam::IOobject::AUTO_WRITE
@ AUTO_WRITE
Definition: IOobject.H:117
Foam::fv::CrankNicolsonDdtScheme::operator=
void operator=(const CrankNicolsonDdtScheme &)
Disallow default bitwise assignment.
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:118
Foam::dimVelocity
const dimensionSet dimVelocity
Foam::fv::CrankNicolsonDdtScheme::evaluate
bool evaluate(const DDt0Field< GeoField > &ddt0) const
Check if the ddt0 needs to be evaluated for this time-step.
Definition: CrankNicolsonDdtScheme.C:187
Foam::fv::CrankNicolsonDdtScheme::rDtCoef0_
dimensionedScalar rDtCoef0_(const DDt0Field< GeoField > &) const
Return the reciprocal old time-step coefficient for Euler for the.
Definition: CrankNicolsonDdtScheme.C:244
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::constant::atomic::alpha
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Definition: readThermalProperties.H:216
Foam::fv::CrankNicolsonDdtScheme::ddt0_
DDt0Field< GeoField > & ddt0_(const word &name, const dimensionSet &dims)
Definition: CrankNicolsonDdtScheme.C:106
Foam::fv::CrankNicolsonDdtScheme::DDt0Field::startTimeIndex
label startTimeIndex() const
Return the start-time index.
Definition: CrankNicolsonDdtScheme.C:78
Foam::fvc::interpolate
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &vf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Definition: surfaceInterpolate.C:110
Foam::fv::CrankNicolsonDdtScheme::fvcDdtUfCorr
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition: CrankNicolsonDdtScheme.C:1121
Foam::fv::CrankNicolsonDdtScheme::coef_
scalar coef_(const DDt0Field< GeoField > &) const
Return the coefficient for Euler scheme for the first time-step.
Definition: CrankNicolsonDdtScheme.C:197
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:804
fvcDiv.H
Calculate the divergence of the given field.
Foam::dimensioned::name
const word & name() const
Return const reference to name.
Definition: dimensionedType.C:235
Foam::GeometricField::boundaryField
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
Definition: GeometricField.C:735
Foam::IOobject::MUST_READ
@ MUST_READ
Definition: IOobject.H:108
Foam::fv::CrankNicolsonDdtScheme::DDt0Field::DDt0Field
DDt0Field(const IOobject &io, const fvMesh &mesh)
Constructor from file for restart.
Definition: CrankNicolsonDdtScheme.C:46
Foam::dimensioned::value
const Type & value() const
Return const reference to value.
Definition: dimensionedType.C:261
Foam::dimensionSet
Dimension set for the base types.
Definition: dimensionSet.H:116
Foam::IOobject::NO_WRITE
@ NO_WRITE
Definition: IOobject.H:118
U
U
Definition: pEqn.H:46
fvMatrices.H
A special matrix type and solver, designed for finite volume solutions of scalar equations.
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Foam::label
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
Foam::fv::CrankNicolsonDdtScheme
Second-oder Crank-Nicolson implicit ddt using the current and previous time-step fields as well as th...
Definition: CrankNicolsonDdtScheme.H:94
Foam::IOobject::NO_READ
@ NO_READ
Definition: IOobject.H:111
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:155
Foam::GeometricField::internalField
InternalField & internalField()
Return internal field.
Definition: GeometricField.C:724
Foam::fv::CrankNicolsonDdtScheme::rDtCoef_
dimensionedScalar rDtCoef_(const DDt0Field< GeoField > &) const
Return the reciprocal time-step coefficient for Euler for the.
Definition: CrankNicolsonDdtScheme.C:233
Foam::fv::CrankNicolsonDdtScheme::fvmDdt
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: CrankNicolsonDdtScheme.C:756
Foam::fv::CrankNicolsonDdtScheme::DDt0Field
Class to store the ddt0 fields on the objectRegistry for use in the.
Definition: CrankNicolsonDdtScheme.H:104
Foam::fv::CrankNicolsonDdtScheme::fvcDdt
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
Definition: CrankNicolsonDdtScheme.C:285
timeIndex
label timeIndex
Definition: getTimeIndex.H:4
Foam::FatalError
error FatalError
CrankNicolsonDdtScheme.H
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:18
Foam::fv::CrankNicolsonDdtScheme::fvcDdtPhiCorr
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
Definition: CrankNicolsonDdtScheme.C:1182
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedScalarFwd.H:41
surfaceInterpolate.H
Surface Interpolation.
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
Foam
Namespace for OpenFOAM.
Definition: combustionModel.C:30
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::fv::CrankNicolsonDdtScheme::meshPhi
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: CrankNicolsonDdtScheme.C:1432
rho
rho
Definition: pEqn.H:3
fv
labelList fv(nPoints)
Foam::fv::CrankNicolsonDdtScheme::coef0_
scalar coef0_(const DDt0Field< GeoField > &) const
Return the old time-step coefficient for Euler scheme for the.
Definition: CrankNicolsonDdtScheme.C:215
Foam::fv::ff
const FieldField< fvPatchField, Type > & ff(const FieldField< fvPatchField, Type > &bf)
Definition: CrankNicolsonDdtScheme.C:272
Foam::regIOobject::store
void store()
Transfer ownership of this object to its registry.
Definition: regIOobjectI.H:34
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:318
Foam::Time::timeName
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:741
Foam::pTraits
Traits class for primitives.
Definition: pTraits.H:50
Foam::fvMatrix::source
Field< Type > & source()
Definition: fvMatrix.H:291
timeName
word timeName
Definition: getTimeIndex.H:3
Foam::fvMatrix< Type >
Foam::fv::CrankNicolsonDdtScheme::offCentre_
tmp< GeoField > offCentre_(const GeoField &ddt0) const
Return ddt0 multiplied by the off-centreing coefficient.
Definition: CrankNicolsonDdtScheme.C:255
Foam::fvMesh::time
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
Foam::dimVol
const dimensionSet dimVol(dimVolume)
Definition: dimensionSets.H:59
Foam::TimeState::timeIndex
label timeIndex() const
Return current time index.
Definition: TimeState.C:73
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const
Return const reference to dimensions.
Definition: dimensionedType.C:248
Foam::dimVolume
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
Foam::GeometricField
Generic GeometricField class.
Definition: surfaceFieldsFwd.H:52
Uf
Uf
Definition: pEqn.H:78
Foam::Time::startTime
virtual dimensionedScalar startTime() const
Return start time.
Definition: Time.C:872
Foam::Time::startTimeIndex
virtual label startTimeIndex() const
Return start time index.
Definition: Time.C:866
Foam::name
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47