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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2018 OpenFOAM Foundation
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "CrankNicolsonDdtScheme.H"
29 #include "surfaceInterpolate.H"
30 #include "fvcDiv.H"
31 #include "fvMatrices.H"
32 #include "Constant.H"
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 
39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
40 
41 namespace fv
42 {
43 
44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
45 
46 template<class Type>
47 template<class GeoField>
48 CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::DDt0Field
49 (
50  const IOobject& io,
51  const fvMesh& mesh
52 )
53 :
54  GeoField(io, mesh),
55  startTimeIndex_(-2) // This field is for a restart and thus correct so set
56  // the start time-index to correspond to a previous run
57 {
58  // Set the time-index to the beginning of the run to ensure the field
59  // is updated during the first time-step
60  this->timeIndex() = mesh.time().startTimeIndex();
61 }
62 
63 
64 template<class Type>
65 template<class GeoField>
66 CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::DDt0Field
67 (
68  const IOobject& io,
69  const fvMesh& mesh,
70  const dimensioned<typename GeoField::value_type>& dimType
71 )
72 :
73  GeoField(io, mesh, dimType),
74  startTimeIndex_(mesh.time().timeIndex())
75 {}
76 
77 
78 template<class Type>
79 template<class GeoField>
80 label CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::startTimeIndex() const
81 {
82  return startTimeIndex_;
83 }
84 
85 
86 template<class Type>
87 template<class GeoField>
88 GeoField& CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::operator()()
89 {
90  return *this;
91 }
92 
93 
94 template<class Type>
95 template<class GeoField>
96 void CrankNicolsonDdtScheme<Type>::DDt0Field<GeoField>::operator=
97 (
98  const GeoField& gf
99 )
100 {
101  GeoField::operator=(gf);
102 }
103 
104 
105 template<class Type>
106 template<class GeoField>
107 typename CrankNicolsonDdtScheme<Type>::template DDt0Field<GeoField>&
108 CrankNicolsonDdtScheme<Type>::ddt0_
109 (
110  const word& name,
111  const dimensionSet& dims
112 )
113 {
114  if (!mesh().objectRegistry::template foundObject<GeoField>(name))
115  {
116  const Time& runTime = mesh().time();
117  word startTimeName = runTime.timeName(runTime.startTime().value());
118 
119  if
120  (
121  (
122  runTime.timeIndex() == runTime.startTimeIndex()
123  || runTime.timeIndex() == runTime.startTimeIndex() + 1
124  )
125  && IOobject
126  (
127  name,
128  startTimeName,
129  mesh()
130  ).template typeHeaderOk<DDt0Field<GeoField>>(true)
131  )
132  {
133  regIOobject::store
134  (
135  new DDt0Field<GeoField>
136  (
137  IOobject
138  (
139  name,
140  startTimeName,
141  mesh(),
142  IOobject::MUST_READ,
143  IOobject::AUTO_WRITE
144  ),
145  mesh()
146  )
147  );
148  }
149  else
150  {
151  regIOobject::store
152  (
153  new DDt0Field<GeoField>
154  (
155  IOobject
156  (
157  name,
158  mesh().time().timeName(),
159  mesh(),
160  IOobject::NO_READ,
161  IOobject::AUTO_WRITE
162  ),
163  mesh(),
164  dimensioned<typename GeoField::value_type>
165  (
166  dims/dimTime,
167  Zero
168  )
169  )
170  );
171  }
172  }
173 
174  DDt0Field<GeoField>& ddt0 = static_cast<DDt0Field<GeoField>&>
175  (
176  mesh().objectRegistry::template lookupObjectRef<GeoField>(name)
177  );
178 
179  return ddt0;
180 }
181 
182 
183 template<class Type>
184 template<class GeoField>
185 bool CrankNicolsonDdtScheme<Type>::evaluate
186 (
187  const DDt0Field<GeoField>& ddt0
188 ) const
189 {
190  return ddt0.timeIndex() != mesh().time().timeIndex();
191 }
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())
202  {
203  return 1 + ocCoeff();
204  }
205  else
206  {
207  return 1;
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 + ocCoeff();
222  }
223  else
224  {
225  return 1;
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)
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 (
273  const FieldField<fvPatchField, Type>& bf
274 )
275 {
276  return bf;
277 }
278 
279 
280 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
281 
282 template<class Type>
283 CrankNicolsonDdtScheme<Type>::CrankNicolsonDdtScheme(const fvMesh& mesh)
284 :
285  ddtScheme<Type>(mesh),
286  ocCoeff_(new Function1Types::Constant<scalar>("ocCoeff", 1))
287 {
288  // Ensure the old-old-time cell volumes are available
289  // for moving meshes
290  if (mesh.moving())
291  {
292  mesh.V00();
293  }
294 }
295 
296 
297 template<class Type>
298 CrankNicolsonDdtScheme<Type>::CrankNicolsonDdtScheme
299 (
300  const fvMesh& mesh,
301  Istream& is
302 )
303 :
304  ddtScheme<Type>(mesh, is)
305 {
306  token firstToken(is);
307 
308  if (firstToken.isNumber())
309  {
310  const scalar ocCoeff = firstToken.number();
311  if (ocCoeff < 0 || ocCoeff > 1)
312  {
313  FatalIOErrorInFunction(is)
314  << "Off-centreing coefficient = " << ocCoeff
315  << " should be >= 0 and <= 1"
316  << exit(FatalIOError);
317  }
318 
319  ocCoeff_.reset
320  (
321  new Function1Types::Constant<scalar>("ocCoeff", ocCoeff)
322  );
323  }
324  else
325  {
326  is.putBack(firstToken);
327  dictionary dict(is);
328  ocCoeff_ = Function1<scalar>::New("ocCoeff", dict);
329  }
330 
331  // Ensure the old-old-time cell volumes are available
332  // for moving meshes
333  if (mesh.moving())
334  {
335  mesh.V00();
336  }
337 }
338 
339 
340 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
341 
342 template<class Type>
343 tmp<GeometricField<Type, fvPatchField, volMesh>>
344 CrankNicolsonDdtScheme<Type>::fvcDdt
345 (
346  const dimensioned<Type>& dt
347 )
348 {
349  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
350  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
351  (
352  "ddt0(" + dt.name() + ')',
353  dt.dimensions()
354  );
355 
356  IOobject ddtIOobject
357  (
358  "ddt(" + dt.name() + ')',
359  mesh().time().timeName(),
360  mesh()
361  );
362 
363  tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
364  (
365  new GeometricField<Type, fvPatchField, volMesh>
366  (
367  ddtIOobject,
368  mesh(),
369  dimensioned<Type>(dt.dimensions()/dimTime, Zero)
370  )
371  );
372 
373  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
374 
375  if (mesh().moving())
376  {
377  if (evaluate(ddt0))
378  {
379  dimensionedScalar rDtCoef0 = rDtCoef0_(ddt0);
380 
381  ddt0.ref() =
382  (
383  (rDtCoef0*dt)*(mesh().V0() - mesh().V00())
384  - mesh().V00()*offCentre_(ddt0.internalField())
385  )/mesh().V0();
386  }
387 
388  tdtdt.ref().ref() =
389  (
390  (rDtCoef*dt)*(mesh().V() - mesh().V0())
391  - mesh().V0()*offCentre_(ddt0.internalField())
392  )/mesh().V();
393  }
394 
395  return tdtdt;
396 }
397 
398 
399 template<class Type>
400 tmp<GeometricField<Type, fvPatchField, volMesh>>
401 CrankNicolsonDdtScheme<Type>::fvcDdt
402 (
403  const GeometricField<Type, fvPatchField, volMesh>& vf
404 )
405 {
406  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
407  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
408  (
409  "ddt0(" + vf.name() + ')',
410  vf.dimensions()
411  );
412 
413  IOobject ddtIOobject
414  (
415  "ddt(" + vf.name() + ')',
416  mesh().time().timeName(),
417  mesh()
418  );
419 
420  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
421 
422  if (mesh().moving())
423  {
424  if (evaluate(ddt0))
425  {
426  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
427 
428  ddt0.primitiveFieldRef() =
429  (
430  rDtCoef0*
431  (
432  mesh().V0()*vf.oldTime().primitiveField()
433  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
434  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
435  )/mesh().V0();
436 
437  ddt0.boundaryFieldRef() =
438  (
439  rDtCoef0*
440  (
441  vf.oldTime().boundaryField()
442  - vf.oldTime().oldTime().boundaryField()
443  ) - offCentre_(ff(ddt0.boundaryField()))
444  );
445  }
446 
447  return tmp<GeometricField<Type, fvPatchField, volMesh>>
448  (
449  new GeometricField<Type, fvPatchField, volMesh>
450  (
451  ddtIOobject,
452  (
453  rDtCoef*
454  (
455  mesh().V()*vf
456  - mesh().V0()*vf.oldTime()
457  ) - mesh().V0()*offCentre_(ddt0()())
458  )/mesh().V(),
459  rDtCoef.value()*
460  (
461  vf.boundaryField() - vf.oldTime().boundaryField()
462  ) - offCentre_(ff(ddt0.boundaryField()))
463  )
464  );
465  }
466  else
467  {
468  if (evaluate(ddt0))
469  {
470  ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime())
471  - offCentre_(ddt0());
472  }
473 
474  return tmp<GeometricField<Type, fvPatchField, volMesh>>
475  (
476  new GeometricField<Type, fvPatchField, volMesh>
477  (
478  ddtIOobject,
479  rDtCoef*(vf - vf.oldTime()) - offCentre_(ddt0())
480  )
481  );
482  }
483 }
484 
485 
486 template<class Type>
487 tmp<GeometricField<Type, fvPatchField, volMesh>>
488 CrankNicolsonDdtScheme<Type>::fvcDdt
489 (
490  const dimensionedScalar& rho,
491  const GeometricField<Type, fvPatchField, volMesh>& vf
492 )
493 {
494  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
495  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
496  (
497  "ddt0(" + rho.name() + ',' + vf.name() + ')',
498  rho.dimensions()*vf.dimensions()
499  );
500 
501  IOobject ddtIOobject
502  (
503  "ddt(" + rho.name() + ',' + vf.name() + ')',
504  mesh().time().timeName(),
505  mesh()
506  );
507 
508  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
509 
510  if (mesh().moving())
511  {
512  if (evaluate(ddt0))
513  {
514  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
515 
516  ddt0.primitiveFieldRef() =
517  (
518  rDtCoef0*rho.value()*
519  (
520  mesh().V0()*vf.oldTime().primitiveField()
521  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
522  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
523  )/mesh().V0();
524 
525  ddt0.boundaryFieldRef() =
526  (
527  rDtCoef0*rho.value()*
528  (
529  vf.oldTime().boundaryField()
530  - vf.oldTime().oldTime().boundaryField()
531  ) - offCentre_(ff(ddt0.boundaryField()))
532  );
533  }
534 
535  return tmp<GeometricField<Type, fvPatchField, volMesh>>
536  (
537  new GeometricField<Type, fvPatchField, volMesh>
538  (
539  ddtIOobject,
540  mesh(),
541  rDtCoef.dimensions()*rho.dimensions()*vf.dimensions(),
542  (
543  rDtCoef.value()*rho.value()*
544  (
545  mesh().V()*vf.primitiveField()
546  - mesh().V0()*vf.oldTime().primitiveField()
547  ) - mesh().V0()*offCentre_(ddt0.primitiveField())
548  )/mesh().V(),
549  rDtCoef.value()*rho.value()*
550  (
551  vf.boundaryField() - vf.oldTime().boundaryField()
552  ) - offCentre_(ff(ddt0.boundaryField()))
553  )
554  );
555  }
556  else
557  {
558  if (evaluate(ddt0))
559  {
560  ddt0 = rDtCoef0_(ddt0)*rho*(vf.oldTime() - vf.oldTime().oldTime())
561  - offCentre_(ddt0());
562  }
563 
564  return tmp<GeometricField<Type, fvPatchField, volMesh>>
565  (
566  new GeometricField<Type, fvPatchField, volMesh>
567  (
568  ddtIOobject,
569  rDtCoef*rho*(vf - vf.oldTime()) - offCentre_(ddt0())
570  )
571  );
572  }
573 }
574 
575 
576 template<class Type>
577 tmp<GeometricField<Type, fvPatchField, volMesh>>
578 CrankNicolsonDdtScheme<Type>::fvcDdt
579 (
580  const volScalarField& rho,
581  const GeometricField<Type, fvPatchField, volMesh>& vf
582 )
583 {
584  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
585  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
586  (
587  "ddt0(" + rho.name() + ',' + vf.name() + ')',
588  rho.dimensions()*vf.dimensions()
589  );
590 
591  IOobject ddtIOobject
592  (
593  "ddt(" + rho.name() + ',' + vf.name() + ')',
594  mesh().time().timeName(),
595  mesh()
596  );
597 
598  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
599 
600  if (mesh().moving())
601  {
602  if (evaluate(ddt0))
603  {
604  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
605 
606  ddt0.primitiveFieldRef() =
607  (
608  rDtCoef0*
609  (
610  mesh().V0()*rho.oldTime().primitiveField()
611  *vf.oldTime().primitiveField()
612  - mesh().V00()*rho.oldTime().oldTime().primitiveField()
613  *vf.oldTime().oldTime().primitiveField()
614  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
615  )/mesh().V0();
616 
617  ddt0.boundaryFieldRef() =
618  (
619  rDtCoef0*
620  (
621  rho.oldTime().boundaryField()
622  *vf.oldTime().boundaryField()
623  - rho.oldTime().oldTime().boundaryField()
624  *vf.oldTime().oldTime().boundaryField()
625  ) - offCentre_(ff(ddt0.boundaryField()))
626  );
627  }
628 
629  return tmp<GeometricField<Type, fvPatchField, volMesh>>
630  (
631  new GeometricField<Type, fvPatchField, volMesh>
632  (
633  ddtIOobject,
634  mesh(),
635  rDtCoef.dimensions()*rho.dimensions()*vf.dimensions(),
636  (
637  rDtCoef.value()*
638  (
639  mesh().V()*rho.primitiveField()*vf.primitiveField()
640  - mesh().V0()*rho.oldTime().primitiveField()
641  *vf.oldTime().primitiveField()
642  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
643  )/mesh().V(),
644  rDtCoef.value()*
645  (
646  rho.boundaryField()*vf.boundaryField()
647  - rho.oldTime().boundaryField()*vf.oldTime().boundaryField()
648  ) - offCentre_(ff(ddt0.boundaryField()))
649  )
650  );
651  }
652  else
653  {
654  if (evaluate(ddt0))
655  {
656  ddt0 = rDtCoef0_(ddt0)*
657  (
658  rho.oldTime()*vf.oldTime()
659  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
660  ) - offCentre_(ddt0());
661  }
662 
663  return tmp<GeometricField<Type, fvPatchField, volMesh>>
664  (
665  new GeometricField<Type, fvPatchField, volMesh>
666  (
667  ddtIOobject,
668  rDtCoef*(rho*vf - rho.oldTime()*vf.oldTime())
669  - offCentre_(ddt0())
670  )
671  );
672  }
673 }
674 
675 
676 template<class Type>
677 tmp<GeometricField<Type, fvPatchField, volMesh>>
678 CrankNicolsonDdtScheme<Type>::fvcDdt
679 (
680  const volScalarField& alpha,
681  const volScalarField& rho,
682  const GeometricField<Type, fvPatchField, volMesh>& vf
683 )
684 {
685  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
686  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
687  (
688  "ddt0(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
689  alpha.dimensions()*rho.dimensions()*vf.dimensions()
690  );
691 
692  IOobject ddtIOobject
693  (
694  "ddt(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
695  mesh().time().timeName(),
696  mesh()
697  );
698 
699  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
700 
701  if (mesh().moving())
702  {
703  if (evaluate(ddt0))
704  {
705  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
706 
707  ddt0.primitiveFieldRef() =
708  (
709  rDtCoef0*
710  (
711  mesh().V0()
712  *alpha.oldTime().primitiveField()
713  *rho.oldTime().primitiveField()
714  *vf.oldTime().primitiveField()
715 
716  - mesh().V00()
717  *alpha.oldTime().oldTime().primitiveField()
718  *rho.oldTime().oldTime().primitiveField()
719  *vf.oldTime().oldTime().primitiveField()
720  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
721  )/mesh().V0();
722 
723  ddt0.boundaryFieldRef() =
724  (
725  rDtCoef0*
726  (
727  alpha.oldTime().boundaryField()
728  *rho.oldTime().boundaryField()
729  *vf.oldTime().boundaryField()
730 
731  - alpha.oldTime().oldTime().boundaryField()
732  *rho.oldTime().oldTime().boundaryField()
733  *vf.oldTime().oldTime().boundaryField()
734  ) - offCentre_(ff(ddt0.boundaryField()))
735  );
736  }
737 
738  return tmp<GeometricField<Type, fvPatchField, volMesh>>
739  (
740  new GeometricField<Type, fvPatchField, volMesh>
741  (
742  ddtIOobject,
743  mesh(),
744  rDtCoef.dimensions()
745  *alpha.dimensions()*rho.dimensions()*vf.dimensions(),
746  (
747  rDtCoef.value()*
748  (
749  mesh().V()
750  *alpha.primitiveField()
751  *rho.primitiveField()
752  *vf.primitiveField()
753 
754  - mesh().V0()
755  *alpha.oldTime().primitiveField()
756  *rho.oldTime().primitiveField()
757  *vf.oldTime().primitiveField()
758  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
759  )/mesh().V(),
760  rDtCoef.value()*
761  (
762  alpha.boundaryField()
763  *rho.boundaryField()
764  *vf.boundaryField()
765 
766  - alpha.oldTime().boundaryField()
767  *rho.oldTime().boundaryField()
768  *vf.oldTime().boundaryField()
769  ) - offCentre_(ff(ddt0.boundaryField()))
770  )
771  );
772  }
773  else
774  {
775  if (evaluate(ddt0))
776  {
777  ddt0 = rDtCoef0_(ddt0)*
778  (
779  alpha.oldTime()
780  *rho.oldTime()
781  *vf.oldTime()
782 
783  - alpha.oldTime().oldTime()
784  *rho.oldTime().oldTime()
785  *vf.oldTime().oldTime()
786  ) - offCentre_(ddt0());
787  }
788 
789  return tmp<GeometricField<Type, fvPatchField, volMesh>>
790  (
791  new GeometricField<Type, fvPatchField, volMesh>
792  (
793  ddtIOobject,
794  rDtCoef
795  *(
796  alpha*rho*vf
797  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
798  )
799  - offCentre_(ddt0())
800  )
801  );
802  }
803 }
804 
805 
806 template<class Type>
807 tmp<fvMatrix<Type>>
808 CrankNicolsonDdtScheme<Type>::fvmDdt
809 (
810  const GeometricField<Type, fvPatchField, volMesh>& vf
811 )
812 {
813  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
814  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
815  (
816  "ddt0(" + vf.name() + ')',
817  vf.dimensions()
818  );
819 
820  tmp<fvMatrix<Type>> tfvm
821  (
822  new fvMatrix<Type>
823  (
824  vf,
825  vf.dimensions()*dimVol/dimTime
826  )
827  );
828 
829  fvMatrix<Type>& fvm = tfvm.ref();
830 
831  const scalar rDtCoef = rDtCoef_(ddt0).value();
832  fvm.diag() = rDtCoef*mesh().V();
833 
834  vf.oldTime().oldTime();
835 
836  if (mesh().moving())
837  {
838  if (evaluate(ddt0))
839  {
840  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
841 
842  ddt0.primitiveFieldRef() =
843  (
844  rDtCoef0*
845  (
846  mesh().V0()*vf.oldTime().primitiveField()
847  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
848  )
849  - mesh().V00()*offCentre_(ddt0.primitiveField())
850  )/mesh().V0();
851 
852  ddt0.boundaryFieldRef() =
853  (
854  rDtCoef0*
855  (
856  vf.oldTime().boundaryField()
857  - vf.oldTime().oldTime().boundaryField()
858  )
859  - offCentre_(ff(ddt0.boundaryField()))
860  );
861  }
862 
863  fvm.source() =
864  (
865  rDtCoef*vf.oldTime().primitiveField()
866  + offCentre_(ddt0.primitiveField())
867  )*mesh().V0();
868  }
869  else
870  {
871  if (evaluate(ddt0))
872  {
873  ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime())
874  - offCentre_(ddt0());
875  }
876 
877  fvm.source() =
878  (
879  rDtCoef*vf.oldTime().primitiveField()
880  + offCentre_(ddt0.primitiveField())
881  )*mesh().V();
882  }
883 
884  return tfvm;
885 }
886 
887 
888 template<class Type>
889 tmp<fvMatrix<Type>>
890 CrankNicolsonDdtScheme<Type>::fvmDdt
891 (
892  const dimensionedScalar& rho,
893  const GeometricField<Type, fvPatchField, volMesh>& vf
894 )
895 {
896  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
897  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
898  (
899  "ddt0(" + rho.name() + ',' + vf.name() + ')',
900  rho.dimensions()*vf.dimensions()
901  );
902 
903  tmp<fvMatrix<Type>> tfvm
904  (
905  new fvMatrix<Type>
906  (
907  vf,
908  rho.dimensions()*vf.dimensions()*dimVol/dimTime
909  )
910  );
911  fvMatrix<Type>& fvm = tfvm.ref();
912 
913  const scalar rDtCoef = rDtCoef_(ddt0).value();
914  fvm.diag() = rDtCoef*rho.value()*mesh().V();
915 
916  vf.oldTime().oldTime();
917 
918  if (mesh().moving())
919  {
920  if (evaluate(ddt0))
921  {
922  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
923 
924  ddt0.primitiveFieldRef() =
925  (
926  rDtCoef0*rho.value()*
927  (
928  mesh().V0()*vf.oldTime().primitiveField()
929  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
930  )
931  - mesh().V00()*offCentre_(ddt0.primitiveField())
932  )/mesh().V0();
933 
934  ddt0.boundaryFieldRef() =
935  (
936  rDtCoef0*rho.value()*
937  (
938  vf.oldTime().boundaryField()
939  - vf.oldTime().oldTime().boundaryField()
940  )
941  - offCentre_(ff(ddt0.boundaryField()))
942  );
943  }
944 
945  fvm.source() =
946  (
947  rDtCoef*rho.value()*vf.oldTime().primitiveField()
948  + offCentre_(ddt0.primitiveField())
949  )*mesh().V0();
950  }
951  else
952  {
953  if (evaluate(ddt0))
954  {
955  ddt0 = rDtCoef0_(ddt0)*rho*(vf.oldTime() - vf.oldTime().oldTime())
956  - offCentre_(ddt0());
957  }
958 
959  fvm.source() =
960  (
961  rDtCoef*rho.value()*vf.oldTime().primitiveField()
962  + offCentre_(ddt0.primitiveField())
963  )*mesh().V();
964  }
965 
966  return tfvm;
967 }
968 
969 
970 template<class Type>
971 tmp<fvMatrix<Type>>
972 CrankNicolsonDdtScheme<Type>::fvmDdt
973 (
974  const volScalarField& rho,
975  const GeometricField<Type, fvPatchField, volMesh>& vf
976 )
977 {
978  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
979  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
980  (
981  "ddt0(" + rho.name() + ',' + vf.name() + ')',
982  rho.dimensions()*vf.dimensions()
983  );
984 
985  tmp<fvMatrix<Type>> tfvm
986  (
987  new fvMatrix<Type>
988  (
989  vf,
990  rho.dimensions()*vf.dimensions()*dimVol/dimTime
991  )
992  );
993  fvMatrix<Type>& fvm = tfvm.ref();
994 
995  const scalar rDtCoef = rDtCoef_(ddt0).value();
996  fvm.diag() = rDtCoef*rho.primitiveField()*mesh().V();
997 
998  vf.oldTime().oldTime();
999  rho.oldTime().oldTime();
1000 
1001  if (mesh().moving())
1002  {
1003  if (evaluate(ddt0))
1004  {
1005  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
1006 
1007  ddt0.primitiveFieldRef() =
1008  (
1009  rDtCoef0*
1010  (
1011  mesh().V0()*rho.oldTime().primitiveField()
1012  *vf.oldTime().primitiveField()
1013  - mesh().V00()*rho.oldTime().oldTime().primitiveField()
1014  *vf.oldTime().oldTime().primitiveField()
1015  )
1016  - mesh().V00()*offCentre_(ddt0.primitiveField())
1017  )/mesh().V0();
1018 
1019  ddt0.boundaryFieldRef() =
1020  (
1021  rDtCoef0*
1022  (
1023  rho.oldTime().boundaryField()
1024  *vf.oldTime().boundaryField()
1025  - rho.oldTime().oldTime().boundaryField()
1026  *vf.oldTime().oldTime().boundaryField()
1027  )
1028  - offCentre_(ff(ddt0.boundaryField()))
1029  );
1030  }
1031 
1032  fvm.source() =
1033  (
1034  rDtCoef*rho.oldTime().primitiveField()*vf.oldTime().primitiveField()
1035  + offCentre_(ddt0.primitiveField())
1036  )*mesh().V0();
1037  }
1038  else
1039  {
1040  if (evaluate(ddt0))
1041  {
1042  ddt0 = rDtCoef0_(ddt0)*
1043  (
1044  rho.oldTime()*vf.oldTime()
1045  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
1046  ) - offCentre_(ddt0());
1047  }
1048 
1049  fvm.source() =
1050  (
1051  rDtCoef*rho.oldTime().primitiveField()*vf.oldTime().primitiveField()
1052  + offCentre_(ddt0.primitiveField())
1053  )*mesh().V();
1054  }
1055 
1056  return tfvm;
1057 }
1058 
1059 
1060 template<class Type>
1061 tmp<fvMatrix<Type>>
1062 CrankNicolsonDdtScheme<Type>::fvmDdt
1063 (
1064  const volScalarField& alpha,
1065  const volScalarField& rho,
1066  const GeometricField<Type, fvPatchField, volMesh>& vf
1067 )
1068 {
1069  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1070  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1071  (
1072  "ddt0(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
1073  alpha.dimensions()*rho.dimensions()*vf.dimensions()
1074  );
1075 
1076  tmp<fvMatrix<Type>> tfvm
1077  (
1078  new fvMatrix<Type>
1079  (
1080  vf,
1081  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
1082  )
1083  );
1084  fvMatrix<Type>& fvm = tfvm.ref();
1085 
1086  const scalar rDtCoef = rDtCoef_(ddt0).value();
1087  fvm.diag() = rDtCoef*alpha.primitiveField()*rho.primitiveField()*mesh().V();
1088 
1089  vf.oldTime().oldTime();
1090  alpha.oldTime().oldTime();
1091  rho.oldTime().oldTime();
1092 
1093  if (mesh().moving())
1094  {
1095  if (evaluate(ddt0))
1096  {
1097  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
1098 
1099  ddt0.primitiveFieldRef() =
1100  (
1101  rDtCoef0*
1102  (
1103  mesh().V0()
1104  *alpha.oldTime().primitiveField()
1105  *rho.oldTime().primitiveField()
1106  *vf.oldTime().primitiveField()
1107 
1108  - mesh().V00()
1109  *alpha.oldTime().oldTime().primitiveField()
1110  *rho.oldTime().oldTime().primitiveField()
1111  *vf.oldTime().oldTime().primitiveField()
1112  )
1113  - mesh().V00()*offCentre_(ddt0.primitiveField())
1114  )/mesh().V0();
1115 
1116  ddt0.boundaryFieldRef() =
1117  (
1118  rDtCoef0*
1119  (
1120  alpha.oldTime().boundaryField()
1121  *rho.oldTime().boundaryField()
1122  *vf.oldTime().boundaryField()
1123 
1124  - alpha.oldTime().oldTime().boundaryField()
1125  *rho.oldTime().oldTime().boundaryField()
1126  *vf.oldTime().oldTime().boundaryField()
1127  )
1128  - offCentre_(ff(ddt0.boundaryField()))
1129  );
1130  }
1131 
1132  fvm.source() =
1133  (
1134  rDtCoef
1135  *alpha.oldTime().primitiveField()
1136  *rho.oldTime().primitiveField()
1137  *vf.oldTime().primitiveField()
1138  + offCentre_(ddt0.primitiveField())
1139  )*mesh().V0();
1140  }
1141  else
1142  {
1143  if (evaluate(ddt0))
1144  {
1145  ddt0 = rDtCoef0_(ddt0)*
1146  (
1147  alpha.oldTime()
1148  *rho.oldTime()
1149  *vf.oldTime()
1150 
1151  - alpha.oldTime().oldTime()
1152  *rho.oldTime().oldTime()
1153  *vf.oldTime().oldTime()
1154  ) - offCentre_(ddt0());
1155  }
1156 
1157  fvm.source() =
1158  (
1159  rDtCoef
1160  *alpha.oldTime().primitiveField()
1161  *rho.oldTime().primitiveField()
1162  *vf.oldTime().primitiveField()
1163  + offCentre_(ddt0.primitiveField())
1164  )*mesh().V();
1165  }
1166 
1167  return tfvm;
1168 }
1169 
1170 
1171 template<class Type>
1172 tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1173 CrankNicolsonDdtScheme<Type>::fvcDdtUfCorr
1174 (
1175  const GeometricField<Type, fvPatchField, volMesh>& U,
1176  const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
1177 )
1178 {
1179  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1180  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1181  (
1182  "ddtCorrDdt0(" + U.name() + ')',
1183  U.dimensions()
1184  );
1185 
1186  DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 =
1187  ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>>
1188  (
1189  "ddtCorrDdt0(" + Uf.name() + ')',
1190  Uf.dimensions()
1191  );
1192 
1193  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1194 
1195  if (evaluate(ddt0))
1196  {
1197  ddt0 =
1198  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1199  - offCentre_(ddt0());
1200  }
1201 
1202  if (evaluate(dUfdt0))
1203  {
1204  dUfdt0 =
1205  rDtCoef0_(dUfdt0)*(Uf.oldTime() - Uf.oldTime().oldTime())
1206  - offCentre_(dUfdt0());
1207  }
1208 
1209  return tmp<fluxFieldType>
1210  (
1211  new fluxFieldType
1212  (
1213  IOobject
1214  (
1215  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
1216  mesh().time().timeName(),
1217  mesh()
1218  ),
1219  this->fvcDdtPhiCoeff(U.oldTime(), mesh().Sf() & Uf.oldTime())
1220  *(
1221  mesh().Sf()
1222  & (
1223  (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1224  - fvc::interpolate(rDtCoef*U.oldTime() + offCentre_(ddt0()))
1225  )
1226  )
1227  )
1228  );
1229 }
1230 
1231 
1232 template<class Type>
1233 tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1234 CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
1235 (
1236  const GeometricField<Type, fvPatchField, volMesh>& U,
1237  const fluxFieldType& phi
1238 )
1239 {
1240  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1241  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1242  (
1243  "ddtCorrDdt0(" + U.name() + ')',
1244  U.dimensions()
1245  );
1246 
1247  DDt0Field<fluxFieldType>& dphidt0 =
1248  ddt0_<fluxFieldType>
1249  (
1250  "ddtCorrDdt0(" + phi.name() + ')',
1251  phi.dimensions()
1252  );
1253  dphidt0.setOriented();
1254 
1255  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1256 
1257  if (evaluate(ddt0))
1258  {
1259  ddt0 =
1260  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1261  - offCentre_(ddt0());
1262  }
1263 
1264  if (evaluate(dphidt0))
1265  {
1266  dphidt0 =
1267  rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime())
1268  - offCentre_(dphidt0());
1269  }
1270 
1271  return tmp<fluxFieldType>
1272  (
1273  new fluxFieldType
1274  (
1275  IOobject
1276  (
1277  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
1278  mesh().time().timeName(),
1279  mesh()
1280  ),
1281  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
1282  *(
1283  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1284  - fvc::dotInterpolate
1285  (
1286  mesh().Sf(),
1287  rDtCoef*U.oldTime() + offCentre_(ddt0())
1288  )
1289  )
1290  )
1291  );
1292 }
1293 
1294 
1295 template<class Type>
1296 tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1297 CrankNicolsonDdtScheme<Type>::fvcDdtUfCorr
1298 (
1299  const volScalarField& rho,
1300  const GeometricField<Type, fvPatchField, volMesh>& U,
1301  const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
1302 )
1303 {
1304  if
1305  (
1306  U.dimensions() == dimVelocity
1307  && Uf.dimensions() == rho.dimensions()*dimVelocity
1308  )
1309  {
1310  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1311  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1312  (
1313  "ddtCorrDdt0(" + rho.name() + ',' + U.name() + ')',
1314  rho.dimensions()*U.dimensions()
1315  );
1316 
1317  DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 =
1318  ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>>
1319  (
1320  "ddtCorrDdt0(" + Uf.name() + ')',
1321  Uf.dimensions()
1322  );
1323 
1324  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1325 
1326  GeometricField<Type, fvPatchField, volMesh> rhoU0
1327  (
1328  rho.oldTime()*U.oldTime()
1329  );
1330 
1331  if (evaluate(ddt0))
1332  {
1333  ddt0 =
1334  rDtCoef0_(ddt0)
1335  *(rhoU0 - rho.oldTime().oldTime()*U.oldTime().oldTime())
1336  - offCentre_(ddt0());
1337  }
1338 
1339  if (evaluate(dUfdt0))
1340  {
1341  dUfdt0 =
1342  rDtCoef0_(dUfdt0)
1343  *(Uf.oldTime() - Uf.oldTime().oldTime())
1344  - offCentre_(dUfdt0());
1345  }
1346 
1347  tmp<fluxFieldType> ddtCorr
1348  (
1349  new fluxFieldType
1350  (
1351  IOobject
1352  (
1353  "ddtCorr("
1354  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
1355  mesh().time().timeName(),
1356  mesh()
1357  ),
1358  this->fvcDdtPhiCoeff
1359  (
1360  rhoU0,
1361  mesh().Sf() & Uf.oldTime(),
1362  rho.oldTime()
1363  )
1364  *(
1365  mesh().Sf()
1366  & (
1367  (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1368  - fvc::interpolate(rDtCoef*rhoU0 + offCentre_(ddt0()))
1369  )
1370  )
1371  )
1372  );
1373 
1374  return ddtCorr;
1375  }
1376  else if
1377  (
1378  U.dimensions() == rho.dimensions()*dimVelocity
1379  && Uf.dimensions() == rho.dimensions()*dimVelocity
1380  )
1381  {
1382  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1383  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1384  (
1385  "ddtCorrDdt0(" + U.name() + ')',
1386  U.dimensions()
1387  );
1388 
1389  DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 =
1390  ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>>
1391  (
1392  "ddtCorrDdt0(" + Uf.name() + ')',
1393  Uf.dimensions()
1394  );
1395 
1396  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1397 
1398  if (evaluate(ddt0))
1399  {
1400  ddt0 =
1401  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1402  - offCentre_(ddt0());
1403  }
1404 
1405  if (evaluate(dUfdt0))
1406  {
1407  dUfdt0 =
1408  rDtCoef0_(dUfdt0)*(Uf.oldTime() - Uf.oldTime().oldTime())
1409  - offCentre_(dUfdt0());
1410  }
1411 
1412  return tmp<fluxFieldType>
1413  (
1414  new fluxFieldType
1415  (
1416  IOobject
1417  (
1418  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
1419  mesh().time().timeName(),
1420  mesh()
1421  ),
1422  this->fvcDdtPhiCoeff
1423  (
1424  U.oldTime(),
1425  mesh().Sf() & Uf.oldTime(),
1426  rho.oldTime()
1427  )
1428  *(
1429  mesh().Sf()
1430  & (
1431  (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1432  - fvc::interpolate
1433  (
1434  rDtCoef*U.oldTime() + offCentre_(ddt0())
1435  )
1436  )
1437  )
1438  )
1439  );
1440  }
1441  else
1442  {
1443  FatalErrorInFunction
1444  << "dimensions of Uf are not correct"
1445  << abort(FatalError);
1446 
1447  return fluxFieldType::null();
1448  }
1449 }
1450 
1451 
1452 template<class Type>
1453 tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1454 CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
1455 (
1456  const volScalarField& rho,
1457  const GeometricField<Type, fvPatchField, volMesh>& U,
1458  const fluxFieldType& phi
1459 )
1460 {
1461  if
1462  (
1463  U.dimensions() == dimVelocity
1464  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1465  )
1466  {
1467  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1468  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1469  (
1470  "ddtCorrDdt0(" + rho.name() + ',' + U.name() + ')',
1471  rho.dimensions()*U.dimensions()
1472  );
1473 
1474  DDt0Field<fluxFieldType>& dphidt0 =
1475  ddt0_<fluxFieldType>
1476  (
1477  "ddtCorrDdt0(" + phi.name() + ')',
1478  phi.dimensions()
1479  );
1480 
1481  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1482 
1483  GeometricField<Type, fvPatchField, volMesh> rhoU0
1484  (
1485  rho.oldTime()*U.oldTime()
1486  );
1487 
1488  if (evaluate(ddt0))
1489  {
1490  ddt0 =
1491  rDtCoef0_(ddt0)
1492  *(rhoU0 - rho.oldTime().oldTime()*U.oldTime().oldTime())
1493  - offCentre_(ddt0());
1494  }
1495 
1496  if (evaluate(dphidt0))
1497  {
1498  dphidt0 =
1499  rDtCoef0_(dphidt0)
1500  *(phi.oldTime() - phi.oldTime().oldTime())
1501  - offCentre_(dphidt0());
1502  }
1503 
1504  tmp<fluxFieldType> ddtCorr
1505  (
1506  new fluxFieldType
1507  (
1508  IOobject
1509  (
1510  "ddtCorr("
1511  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
1512  mesh().time().timeName(),
1513  mesh()
1514  ),
1515  this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), rho.oldTime())
1516  *(
1517  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1518  - fvc::dotInterpolate
1519  (
1520  mesh().Sf(),
1521  rDtCoef*rhoU0 + offCentre_(ddt0())
1522  )
1523  )
1524  )
1525  );
1526 
1527  return ddtCorr;
1528  }
1529  else if
1530  (
1531  U.dimensions() == rho.dimensions()*dimVelocity
1532  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1533  )
1534  {
1535  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1536  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1537  (
1538  "ddtCorrDdt0(" + U.name() + ')',
1539  U.dimensions()
1540  );
1541 
1542  DDt0Field<fluxFieldType>& dphidt0 =
1543  ddt0_<fluxFieldType>
1544  (
1545  "ddtCorrDdt0(" + phi.name() + ')',
1546  phi.dimensions()
1547  );
1548 
1549  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1550 
1551  if (evaluate(ddt0))
1552  {
1553  ddt0 =
1554  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1555  - offCentre_(ddt0());
1556  }
1557 
1558  if (evaluate(dphidt0))
1559  {
1560  dphidt0 =
1561  rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime())
1562  - offCentre_(dphidt0());
1563  }
1564 
1565  return tmp<fluxFieldType>
1566  (
1567  new fluxFieldType
1568  (
1569  IOobject
1570  (
1571  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
1572  mesh().time().timeName(),
1573  mesh()
1574  ),
1575  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), rho.oldTime())
1576  *(
1577  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1578  - fvc::dotInterpolate
1579  (
1580  mesh().Sf(),
1581  rDtCoef*U.oldTime() + offCentre_(ddt0())
1582  )
1583  )
1584  )
1585  );
1586  }
1587  else
1588  {
1589  FatalErrorInFunction
1590  << "dimensions of phi are not correct"
1591  << abort(FatalError);
1592 
1593  return fluxFieldType::null();
1594  }
1595 }
1596 
1597 
1598 template<class Type>
1599 tmp<surfaceScalarField> CrankNicolsonDdtScheme<Type>::meshPhi
1600 (
1601  const GeometricField<Type, fvPatchField, volMesh>& vf
1602 )
1603 {
1604  DDt0Field<surfaceScalarField>& meshPhi0 = ddt0_<surfaceScalarField>
1605  (
1606  "meshPhiCN_0",
1607  dimVolume
1608  );
1609 
1610  meshPhi0.setOriented();
1611 
1612  if (evaluate(meshPhi0))
1613  {
1614  meshPhi0 =
1615  coef0_(meshPhi0)*mesh().phi().oldTime() - offCentre_(meshPhi0());
1616  }
1617 
1618  return tmp<surfaceScalarField>
1619  (
1620  new surfaceScalarField
1621  (
1622  IOobject
1623  (
1624  mesh().phi().name(),
1625  mesh().time().timeName(),
1626  mesh(),
1627  IOobject::NO_READ,
1628  IOobject::NO_WRITE,
1629  false
1630  ),
1631  coef_(meshPhi0)*mesh().phi() - offCentre_(meshPhi0())
1632  )
1633  );
1634 }
1635 
1636 
1637 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1638 
1639 } // End namespace fv
1640 
1641 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1642 
1643 } // End namespace Foam
1644 
1645 // ************************************************************************* //
Foam::IOobject::NO_WRITE
Definition: IOobject.H:130
Foam::fv::CrankNicolsonDdtScheme::mesh
const fvMesh & mesh() const
Return mesh reference.
Definition: CrankNicolsonDdtScheme.H:227
runTime
engineTime & runTime
Definition: createEngineTime.H:13
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:104
Foam::IOobject::AUTO_WRITE
Definition: IOobject.H:129
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:73
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:62
Foam::FieldField
A field of fields is a PtrList of fields with reference counting.
Definition: FieldField.H:53
Foam::fvc::dotInterpolate
static tmp< GeometricField< typename innerProduct< vector, Type >::type, fvsPatchField, surfaceMesh > > dotInterpolate(const surfaceVectorField &Sf, const GeometricField< Type, fvPatchField, volMesh > &tvf)
Interpolate field onto faces.
Foam::token::isNumber
bool isNumber() const
Token is LABEL, FLOAT or DOUBLE.
Definition: tokenI.H:551
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:59
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::dimVelocity
const dimensionSet dimVelocity
Foam::polyMesh::moving
bool moving() const
Is mesh moving.
Definition: polyMesh.H:513
Foam::constant::atomic::alpha
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Definition: readThermalProperties.H:212
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Definition: GeometricFieldI.H:53
Foam::fv::CrankNicolsonDdtScheme::fvcDdtUfCorr
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition: CrankNicolsonDdtScheme.C:1174
Uf
autoPtr< surfaceVectorField > Uf
Definition: createUfIfPresent.H:33
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:850
fvcDiv.H
Calculate the divergence of the given field.
Foam::dimensioned::name
const word & name() const
Return const reference to name.
Definition: dimensionedType.C:406
Foam::FatalIOError
IOerror FatalIOError
Foam::token
A token holds an item read from Istream.
Definition: token.H:69
Foam::dimensionSet
Dimension set for the base types.
Definition: dimensionSet.H:65
Foam::fvMesh::V00
const DimensionedField< scalar, volMesh > & V00() const
Return old-old-time cell volumes.
Definition: fvMeshGeometry.C:233
rho
rho
Definition: readInitialConditions.H:88
Foam::regIOobject::store
bool store()
Definition: regIOobjectI.H:37
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:54
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:60
Foam::token::number
scalar number() const
Return label, float or double value.
Definition: tokenI.H:557
Foam::fv::CrankNicolsonDdtScheme
Second-oder Crank-Nicolson implicit ddt using the current and previous time-step fields as well as th...
Definition: CrankNicolsonDdtScheme.H:115
Foam::Function1Types::Constant
Templated function that returns a constant value.
Definition: Constant.H:59
Foam::Istream
An Istream is an abstract base class for all input systems (streams, files, token lists etc)....
Definition: Istream.H:61
Foam::name
word name(const complex &c)
Return string representation of complex.
Definition: complex.C:76
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:43
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::fv::CrankNicolsonDdtScheme::fvmDdt
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: CrankNicolsonDdtScheme.C:809
ddtCorr
ddtCorr
Definition: readControls.H:9
Foam::fv::CrankNicolsonDdtScheme::fvcDdt
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
Definition: CrankNicolsonDdtScheme.C:345
timeName
word timeName
Definition: getTimeIndex.H:3
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:121
CrankNicolsonDdtScheme.H
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::fv::CrankNicolsonDdtScheme::fvcDdtPhiCorr
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
Definition: CrankNicolsonDdtScheme.C:1235
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedScalarFwd.H:43
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:84
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:137
Foam::fvc::interpolate
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Foam::fv::CrankNicolsonDdtScheme::meshPhi
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: CrankNicolsonDdtScheme.C:1600
ocCoeff
scalar ocCoeff
Definition: alphaEqn.H:6
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
fv
labelList fv(nPoints)
U
U
Definition: pEqn.H:72
Foam::fv::ff
const FieldField< fvPatchField, Type > & ff(const FieldField< fvPatchField, Type > &bf)
Definition: CrankNicolsonDdtScheme.C:272
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:372
Foam::Istream::putBack
void putBack(const token &tok)
Put back token.
Definition: Istream.C:53
Foam::fvMatrix::source
Field< Type > & source()
Definition: fvMatrix.H:295
Foam::fvMatrix
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvPatchField.H:76
Constant.H
Foam::fv::ddtScheme
Abstract base class for ddt schemes.
Definition: ddtScheme.H:83
Foam::dimVol
const dimensionSet dimVol(dimVolume)
Older spelling for dimVolume.
Definition: dimensionSets.H:65
FatalIOErrorInFunction
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:392
timeIndex
label timeIndex
Definition: getTimeIndex.H:4
Foam::Function1::New
static autoPtr< Function1< Type > > New(const word &entryName, const dictionary &dict, const word &redirectType=word::null)
Selector.
Definition: Function1New.C:35
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const
Return const reference to dimensions.
Definition: dimensionedType.C:420
Foam::dimVolume
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:61
Foam::GeometricField< Type, fvPatchField, volMesh >
Foam::IOobject::NO_READ
Definition: IOobject.H:123
Foam::stringOps::evaluate
string evaluate(const std::string &s, size_t pos=0, size_t len=std::string::npos)
Definition: stringOpsEvaluate.C:37
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:62
Foam::IOobject::MUST_READ
Definition: IOobject.H:120