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