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-2014 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  pTraits<typename GeoField::value_type>::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,
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 {
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.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 {
441  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
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 {
531  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
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 {
632  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
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 {
760  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
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,
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 {
843  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
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,
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 {
925  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
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,
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 {
1016  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
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 {
1126  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
1127  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
1128  (
1129  "ddt0(" + U.name() + ')',
1130  U.dimensions()
1131  );
1132 
1133  DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh> >& dUfdt0 =
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 {
1187  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
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  {
1255  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
1256  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
1257  (
1258  "ddt0(" + rho.name() + ',' + U.name() + ')',
1259  U.dimensions()
1260  );
1261 
1262  DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh> >& dUfdt0 =
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  {
1326  FatalErrorIn
1327  (
1328  "CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr"
1329  ) << "dimensions of Uf are not correct"
1330  << abort(FatalError);
1331 
1332  return fluxFieldType::null();
1333  }
1334 }
1335 
1336 
1337 template<class Type>
1341  const volScalarField& rho,
1343  const fluxFieldType& phi
1344 )
1345 {
1346  if
1347  (
1348  U.dimensions() == dimVelocity
1349  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1350  )
1351  {
1352  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
1353  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
1354  (
1355  "ddt0(" + rho.name() + ',' + U.name() + ')',
1356  U.dimensions()
1357  );
1358 
1359  DDt0Field<fluxFieldType>& dphidt0 =
1360  ddt0_<fluxFieldType>
1361  (
1362  "ddt0(" + phi.name() + ')',
1363  phi.dimensions()
1364  );
1365 
1366  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1367 
1369  (
1370  rho.oldTime()*U.oldTime()
1371  );
1372 
1373  if (evaluate(ddt0))
1374  {
1375  ddt0 =
1376  rDtCoef0_(ddt0)
1377  *(rhoU0 - rho.oldTime().oldTime()*U.oldTime().oldTime())
1378  - offCentre_(ddt0());
1379  }
1380 
1381  if (evaluate(dphidt0))
1382  {
1383  dphidt0 =
1384  rDtCoef0_(dphidt0)
1385  *(phi.oldTime() - phi.oldTime().oldTime())
1386  - offCentre_(dphidt0());
1387  }
1388 
1390  (
1391  new fluxFieldType
1392  (
1393  IOobject
1394  (
1395  "ddtCorr("
1396  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
1397  mesh().time().timeName(),
1398  mesh()
1399  ),
1400  this->fvcDdtPhiCoeff(rhoU0, phi.oldTime())
1401  *(
1402  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1403  - (
1404  mesh().Sf()
1405  & fvc::interpolate(rDtCoef*rhoU0 + offCentre_(ddt0()))
1406  )
1407  )
1408  )
1409  );
1410 
1411  return ddtCorr;
1412  }
1413  else if
1414  (
1415  U.dimensions() == rho.dimensions()*dimVelocity
1416  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1417  )
1418  {
1419  return fvcDdtPhiCorr(U, phi);
1420  }
1421  else
1422  {
1423  FatalErrorIn
1424  (
1425  "CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr"
1426  ) << "dimensions of phi are not correct"
1427  << abort(FatalError);
1428 
1429  return fluxFieldType::null();
1430  }
1431 }
1432 
1433 
1434 template<class Type>
1438 )
1439 {
1440  DDt0Field<surfaceScalarField>& meshPhi0 = ddt0_<surfaceScalarField>
1441  (
1442  "meshPhiCN_0",
1443  dimVolume
1444  );
1445 
1446  if (evaluate(meshPhi0))
1447  {
1448  meshPhi0 =
1449  coef0_(meshPhi0)*mesh().phi().oldTime() - offCentre_(meshPhi0());
1450  }
1451 
1453  (
1454  new surfaceScalarField
1455  (
1456  IOobject
1457  (
1458  mesh().phi().name(),
1459  mesh().time().timeName(),
1460  mesh(),
1463  false
1464  ),
1465  coef_(meshPhi0)*mesh().phi() - offCentre_(meshPhi0())
1466  )
1467  );
1468 }
1469 
1470 
1471 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1472 
1473 } // End namespace fv
1474 
1475 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1476 
1477 } // End namespace Foam
1478 
1479 // ************************************************************************* //
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &vf, const surfaceScalarField &faceFlux, Istream &schemeData)
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
Generic field type.
Definition: FieldField.H:51
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
const surfaceVectorField & Sf() const
Return cell face area vectors.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
const surfaceScalarField & phi() const
Return cell face motion fluxes.
const dimensionSet dimVol(dimVolume)
Definition: dimensionSets.H:59
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
label timeIndex
Definition: getTimeIndex.H:4
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
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
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:741
const DimensionedField< scalar, volMesh > & V00() const
Return old-old-time cell volumes.
InternalField & internalField()
Return internal field.
Surface Interpolation.
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:68
dimensionedScalar deltaT0() const
Return old time step.
Definition: TimeState.C:85
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
Namespace for OpenFOAM.
const dimensionSet & dimensions() const
Return const reference to dimensions.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
Generic dimensioned Type class.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
label timeIndex() const
Return current time index.
Definition: TimeState.C:73
const word & name() const
Return const reference to name.
void store()
Transfer ownership of this object to its registry.
Definition: regIOobjectI.H:34
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
const dimensionSet & dimensions() const
Return dimensions.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const word & name() const
Return name.
Definition: IOobject.H:260
scalarField & diag()
Definition: lduMatrix.C:183
const FieldField< fvPatchField, Type > & ff(const FieldField< fvPatchField, Type > &bf)
surfaceScalarField & phi
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
static const GeometricField< Type, PatchField, GeoMesh > & null()
Return a null geometric field.
Generic GeometricField class.
Calculate the divergence of the given field.
Traits class for primitives.
Definition: pTraits.H:50
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
ddtScheme< Type >::fluxFieldType fluxFieldType
error FatalError
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
labelList fv(nPoints)
const fvMesh & mesh() const
Return mesh reference.
dimensionedScalar deltaT() const
Return time step.
Definition: TimeState.C:79
Field< Type > & source()
Definition: fvMatrix.H:291
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
const dimensionSet dimVelocity
tmp< surfaceScalarField > fvcDdtPhiCoeff(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi, const fluxFieldType &phiCorr)
Definition: ddtScheme.C:141
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
const DimensionedField< scalar, volMesh > & V0() const
Return old-time cell volumes.
const Type & value() const
Return const reference to value.
word timeName
Definition: getTimeIndex.H:3
A special matrix type and solver, designed for finite volume solutions of scalar equations.