FreeFOAM The Cross-Platform CFD Toolkit
CrankNicholsonDdtScheme.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) 1991-2010 OpenCFD Ltd.
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 
28 #include <finiteVolume/fvcDiv.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 CrankNicholsonDdtScheme<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 corespond 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 CrankNicholsonDdtScheme<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 CrankNicholsonDdtScheme<Type>::DDt0Field<GeoField>::
78 startTimeIndex() const
79 {
80  return startTimeIndex_;
81 }
82 
83 
84 template<class Type>
85 template<class GeoField>
86 GeoField& CrankNicholsonDdtScheme<Type>::DDt0Field<GeoField>::
87 operator()()
88 {
89  return *this;
90 }
91 
92 
93 template<class Type>
94 template<class GeoField>
95 void CrankNicholsonDdtScheme<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 CrankNicholsonDdtScheme<Type>::DDt0Field<GeoField>&
105 CrankNicholsonDdtScheme<Type>::ddt0_
106 (
107  const word& name,
108  const dimensionSet& dims
109 )
110 {
111  if (!mesh().objectRegistry::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::lookupObject<GeoField>(name)
177  )
178  );
179 
180  return ddt0;
181 }
182 
183 
184 template<class Type>
185 template<class GeoField>
186 bool CrankNicholsonDdtScheme<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 CrankNicholsonDdtScheme<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 CrankNicholsonDdtScheme<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 CrankNicholsonDdtScheme<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 CrankNicholsonDdtScheme<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> CrankNicholsonDdtScheme<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 (
628 )
629 {
630  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
631  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
632  (
633  "ddt0(" + vf.name() + ')',
634  vf.dimensions()
635  );
636 
637  tmp<fvMatrix<Type> > tfvm
638  (
639  new fvMatrix<Type>
640  (
641  vf,
643  )
644  );
645 
646  fvMatrix<Type>& fvm = tfvm();
647 
648  scalar rDtCoef = rDtCoef_(ddt0).value();
649 
650  fvm.diag() = rDtCoef*mesh().V();
651 
652  vf.oldTime().oldTime();
653 
654  if (mesh().moving())
655  {
656  if (evaluate(ddt0))
657  {
658  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
659 
660  ddt0.internalField() =
661  (
662  rDtCoef0*
663  (
664  mesh().V0()*vf.oldTime().internalField()
665  - mesh().V00()*vf.oldTime().oldTime().internalField()
666  )
667  - mesh().V00()*offCentre_(ddt0.internalField())
668  )/mesh().V0();
669 
670  ddt0.boundaryField() =
671  (
672  rDtCoef0*
673  (
674  vf.oldTime().boundaryField()
675  - vf.oldTime().oldTime().boundaryField()
676  )
677  - offCentre_(ff(ddt0.boundaryField()))
678  );
679  }
680 
681  fvm.source() =
682  (
683  rDtCoef*vf.oldTime().internalField()
684  + offCentre_(ddt0.internalField())
685  )*mesh().V0();
686  }
687  else
688  {
689  if (evaluate(ddt0))
690  {
691  ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime())
692  - offCentre_(ddt0());
693  }
694 
695  fvm.source() =
696  (
697  rDtCoef*vf.oldTime().internalField()
698  + offCentre_(ddt0.internalField())
699  )*mesh().V();
700  }
701 
702  return tfvm;
703 }
704 
705 
706 template<class Type>
709 (
710  const dimensionedScalar& rho,
712 )
713 {
714  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
715  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
716  (
717  "ddt0(" + rho.name() + ',' + vf.name() + ')',
718  rho.dimensions()*vf.dimensions()
719  );
720 
721  tmp<fvMatrix<Type> > tfvm
722  (
723  new fvMatrix<Type>
724  (
725  vf,
727  )
728  );
729  fvMatrix<Type>& fvm = tfvm();
730 
731  scalar rDtCoef = rDtCoef_(ddt0).value();
732  fvm.diag() = rDtCoef*rho.value()*mesh().V();
733 
734  vf.oldTime().oldTime();
735 
736  if (mesh().moving())
737  {
738  if (evaluate(ddt0))
739  {
740  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
741 
742  ddt0.internalField() =
743  (
744  rDtCoef0*rho.value()*
745  (
746  mesh().V0()*vf.oldTime().internalField()
747  - mesh().V00()*vf.oldTime().oldTime().internalField()
748  )
749  - mesh().V00()*offCentre_(ddt0.internalField())
750  )/mesh().V0();
751 
752  ddt0.boundaryField() =
753  (
754  rDtCoef0*rho.value()*
755  (
756  vf.oldTime().boundaryField()
757  - vf.oldTime().oldTime().boundaryField()
758  )
759  - offCentre_(ff(ddt0.boundaryField()))
760  );
761  }
762 
763  fvm.source() =
764  (
765  rDtCoef*rho.value()*vf.oldTime().internalField()
766  + offCentre_(ddt0.internalField())
767  )*mesh().V0();
768  }
769  else
770  {
771  if (evaluate(ddt0))
772  {
773  ddt0 = rDtCoef0_(ddt0)*rho*(vf.oldTime() - vf.oldTime().oldTime())
774  - offCentre_(ddt0());
775  }
776 
777  fvm.source() =
778  (
779  rDtCoef*rho.value()*vf.oldTime().internalField()
780  + offCentre_(ddt0.internalField())
781  )*mesh().V();
782  }
783 
784  return tfvm;
785 }
786 
787 
788 template<class Type>
791 (
792  const volScalarField& rho,
794 )
795 {
796  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
797  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
798  (
799  "ddt0(" + rho.name() + ',' + vf.name() + ')',
800  rho.dimensions()*vf.dimensions()
801  );
802 
803  tmp<fvMatrix<Type> > tfvm
804  (
805  new fvMatrix<Type>
806  (
807  vf,
809  )
810  );
811  fvMatrix<Type>& fvm = tfvm();
812 
813  scalar rDtCoef = rDtCoef_(ddt0).value();
814  fvm.diag() = rDtCoef*rho.internalField()*mesh().V();
815 
816  vf.oldTime().oldTime();
817  rho.oldTime().oldTime();
818 
819  if (mesh().moving())
820  {
821  if (evaluate(ddt0))
822  {
823  scalar rDtCoef0 = rDtCoef0_(ddt0).value();
824 
825  ddt0.internalField() =
826  (
827  rDtCoef0*
828  (
829  mesh().V0()*rho.oldTime().internalField()
830  *vf.oldTime().internalField()
831  - mesh().V00()*rho.oldTime().oldTime().internalField()
832  *vf.oldTime().oldTime().internalField()
833  )
834  - mesh().V00()*offCentre_(ddt0.internalField())
835  )/mesh().V0();
836 
837  ddt0.boundaryField() =
838  (
839  rDtCoef0*
840  (
841  rho.oldTime().boundaryField()
842  *vf.oldTime().boundaryField()
843  - rho.oldTime().oldTime().boundaryField()
844  *vf.oldTime().oldTime().boundaryField()
845  )
846  - offCentre_(ff(ddt0.boundaryField()))
847  );
848  }
849 
850  fvm.source() =
851  (
852  rDtCoef*rho.internalField()*vf.oldTime().internalField()
853  + offCentre_(ddt0.internalField())
854  )*mesh().V0();
855  }
856  else
857  {
858  if (evaluate(ddt0))
859  {
860  ddt0 = rDtCoef0_(ddt0)*
861  (
862  rho.oldTime()*vf.oldTime()
863  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
864  ) - offCentre_(ddt0());
865  }
866 
867  fvm.source() =
868  (
869  rDtCoef*rho.oldTime().internalField()*vf.oldTime().internalField()
870  + offCentre_(ddt0.internalField())
871  )*mesh().V();
872  }
873 
874  return tfvm;
875 }
876 
877 
878 
879 template<class Type>
882 (
883  const volScalarField& rA,
885  const fluxFieldType& phi
886 )
887 {
888  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& dUdt0 =
889  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
890  (
891  "ddt0(" + U.name() + ')',
892  U.dimensions()
893  );
894 
895  DDt0Field<fluxFieldType>& dphidt0 =
896  ddt0_<fluxFieldType>
897  (
898  "ddt0(" + phi.name() + ')',
899  phi.dimensions()
900  );
901 
902  IOobject ddtIOobject
903  (
904  "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
905  mesh().time().timeName(),
906  mesh()
907  );
908 
909  dimensionedScalar rDtCoef = rDtCoef_(dUdt0);
910 
911  if (mesh().moving())
912  {
913  return tmp<fluxFieldType>
914  (
915  new fluxFieldType
916  (
917  ddtIOobject,
918  mesh(),
919  dimensioned<typename flux<Type>::type>
920  (
921  "0",
922  rDtCoef.dimensions()*rA.dimensions()*phi.dimensions(),
924  )
925  )
926  );
927  }
928  else
929  {
930  if (evaluate(dUdt0))
931  {
932  dUdt0 =
933  rDtCoef0_(dUdt0)*(U.oldTime() - U.oldTime().oldTime())
934  - offCentre_(dUdt0());
935  }
936 
937  if (evaluate(dphidt0))
938  {
939  dphidt0 =
940  rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime())
941  - offCentre_(dphidt0());
942  }
943 
944  return tmp<fluxFieldType>
945  (
946  new fluxFieldType
947  (
948  ddtIOobject,
949  fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
950  *fvc::interpolate(rA)
951  *(
952  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
953  - (
955  (
956  (rDtCoef*U.oldTime() + offCentre_(dUdt0()))
957  ) & mesh().Sf()
958  )
959  )
960  )
961  );
962  }
963 }
964 
965 
966 template<class Type>
969 (
970  const volScalarField& rA,
971  const volScalarField& rho,
973  const fluxFieldType& phi
974 )
975 {
976  DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& dUdt0 =
977  ddt0_<GeometricField<Type, fvPatchField, volMesh> >
978  (
979  "ddt0(" + U.name() + ')',
980  U.dimensions()
981  );
982 
983  DDt0Field<fluxFieldType>& dphidt0 =
984  ddt0_<fluxFieldType>
985  (
986  "ddt0(" + phi.name() + ')',
987  U.dimensions()*dimArea
988  );
989 
990  IOobject ddtIOobject
991  (
992  "ddtPhiCorr("
993  + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
994  mesh().time().timeName(),
995  mesh()
996  );
997 
998  dimensionedScalar rDtCoef = rDtCoef_(dUdt0);
999 
1000  if (mesh().moving())
1001  {
1002  return tmp<fluxFieldType>
1003  (
1004  new fluxFieldType
1005  (
1006  ddtIOobject,
1007  mesh(),
1008  dimensioned<typename flux<Type>::type>
1009  (
1010  "0",
1011  rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
1013  )
1014  )
1015  );
1016  }
1017  else
1018  {
1019  if
1020  (
1021  U.dimensions() == dimVelocity
1022  && phi.dimensions() == dimVelocity*dimArea
1023  )
1024  {
1025  if (evaluate(dUdt0))
1026  {
1027  dUdt0 = rDtCoef0_(dUdt0)*
1028  (
1029  U.oldTime() - U.oldTime().oldTime()
1030  ) - offCentre_(dUdt0());
1031  }
1032 
1033  if (evaluate(dphidt0))
1034  {
1035  dphidt0 = rDtCoef0_(dphidt0)*
1036  (
1037  phi.oldTime()
1038  - fvc::interpolate(rho.oldTime().oldTime()/rho.oldTime())
1039  *phi.oldTime().oldTime()
1040  ) - offCentre_(dphidt0());
1041  }
1042 
1043  return tmp<fluxFieldType>
1044  (
1045  new fluxFieldType
1046  (
1047  ddtIOobject,
1048  fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
1049  (
1050  fvc::interpolate(rA*rho.oldTime())
1051  *(rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1052  - (
1054  (
1055  rA*rho.oldTime()
1056  *(rDtCoef*U.oldTime() + offCentre_(dUdt0()))
1057  ) & mesh().Sf()
1058  )
1059  )
1060  )
1061  );
1062  }
1063  else if
1064  (
1065  U.dimensions() == dimVelocity
1066  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1067  )
1068  {
1069  if (evaluate(dUdt0))
1070  {
1071  dUdt0 = rDtCoef0_(dUdt0)*
1072  (
1073  U.oldTime() - U.oldTime().oldTime()
1074  ) - offCentre_(dUdt0());
1075  }
1076 
1077  if (evaluate(dphidt0))
1078  {
1079  dphidt0 = rDtCoef0_(dphidt0)*
1080  (
1081  phi.oldTime()
1082  /fvc::interpolate(rho.oldTime())
1083  - phi.oldTime().oldTime()
1084  /fvc::interpolate(rho.oldTime().oldTime())
1085  ) - offCentre_(dphidt0());
1086  }
1087 
1088  return tmp<fluxFieldType>
1089  (
1090  new fluxFieldType
1091  (
1092  ddtIOobject,
1093  fvcDdtPhiCoeff
1094  (
1095  U.oldTime(),
1096  phi.oldTime()/fvc::interpolate(rho.oldTime())
1097  )
1098  *(
1099  fvc::interpolate(rA*rho.oldTime())
1100  *(
1101  rDtCoef*phi.oldTime()/fvc::interpolate(rho.oldTime())
1102  + offCentre_(dphidt0())
1103  )
1104  - (
1106  (
1107  rA*rho.oldTime()
1108  *(rDtCoef*U.oldTime() + offCentre_(dUdt0()))
1109  ) & mesh().Sf()
1110  )
1111  )
1112  )
1113  );
1114  }
1115  else if
1116  (
1117  U.dimensions() == rho.dimensions()*dimVelocity
1118  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1119  )
1120  {
1121  if (evaluate(dUdt0))
1122  {
1123  dUdt0 = rDtCoef0_(dUdt0)*
1124  (
1125  U.oldTime() - U.oldTime().oldTime()
1126  ) - offCentre_(dUdt0());
1127  }
1128 
1129  if (evaluate(dphidt0))
1130  {
1131  dphidt0 = rDtCoef0_(dphidt0)*
1132  (
1133  phi.oldTime() - phi.oldTime().oldTime()
1134  ) - offCentre_(dphidt0());
1135  }
1136 
1137  return tmp<fluxFieldType>
1138  (
1139  new fluxFieldType
1140  (
1141  ddtIOobject,
1142  fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())*
1143  (
1144  fvc::interpolate(rA)
1145  *(rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1146  - (
1148  (
1149  rA*(rDtCoef*U.oldTime() + offCentre_(dUdt0()))
1150  ) & mesh().Sf()
1151  )
1152  )
1153  )
1154  );
1155  }
1156  else
1157  {
1158  FatalErrorIn
1159  (
1160  "CrankNicholsonDdtScheme<Type>::fvcDdtPhiCorr"
1161  ) << "dimensions of phi are not correct"
1162  << abort(FatalError);
1163 
1164  return fluxFieldType::null();
1165  }
1166  }
1167 }
1168 
1169 
1170 template<class Type>
1174 )
1175 {
1176  DDt0Field<surfaceScalarField>& meshPhi0 = ddt0_<surfaceScalarField>
1177  (
1178  "meshPhiCN_0",
1179  dimVolume
1180  );
1181 
1182  if (evaluate(meshPhi0))
1183  {
1184  meshPhi0 =
1185  coef0_(meshPhi0)*mesh().phi().oldTime() - offCentre_(meshPhi0());
1186  }
1187 
1188  return coef_(meshPhi0)*mesh().phi() - offCentre_(meshPhi0());
1189 }
1190 
1191 
1192 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1193 
1194 } // End namespace fv
1195 
1196 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1197 
1198 } // End namespace Foam
1199 
1200 // ************************ vim: set sw=4 sts=4 et: ************************ //