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