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