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 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(),
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()
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.ref() =
385  (
386  (rDtCoef0*dt)*(mesh().V0() - mesh().V00())
387  - mesh().V00()*offCentre_(ddt0.internalField())
388  )/mesh().V0();
389  }
390 
391  tdtdt.ref().ref() =
392  (
393  (rDtCoef*dt)*(mesh().V() - mesh().V0())
394  - mesh().V0()*offCentre_(ddt0.internalField())
395  )/mesh().V();
396  }
397 
398  return tdtdt;
399 }
400 
401 
402 template<class Type>
403 tmp<GeometricField<Type, fvPatchField, volMesh>>
404 CrankNicolsonDdtScheme<Type>::fvcDdt
405 (
406  const GeometricField<Type, fvPatchField, volMesh>& vf
407 )
408 {
409  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
410  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
411  (
412  "ddt0(" + vf.name() + ')',
413  vf.dimensions()
414  );
415 
416  IOobject ddtIOobject
417  (
418  "ddt(" + vf.name() + ')',
419  mesh().time().timeName(),
420  mesh()
421  );
422 
423  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
424 
425  if (mesh().moving())
426  {
427  if (evaluate(ddt0))
428  {
429  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
430 
431  ddt0.primitiveFieldRef() =
432  (
433  rDtCoef0*
434  (
435  mesh().V0()*vf.oldTime().primitiveField()
436  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
437  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
438  )/mesh().V0();
439 
440  ddt0.boundaryFieldRef() =
441  (
442  rDtCoef0*
443  (
444  vf.oldTime().boundaryField()
445  - vf.oldTime().oldTime().boundaryField()
446  ) - offCentre_(ff(ddt0.boundaryField()))
447  );
448  }
449 
450  return tmp<GeometricField<Type, fvPatchField, volMesh>>
451  (
452  new GeometricField<Type, fvPatchField, volMesh>
453  (
454  ddtIOobject,
455  (
456  rDtCoef*
457  (
458  mesh().V()*vf
459  - mesh().V0()*vf.oldTime()
460  ) - mesh().V0()*offCentre_(ddt0()())
461  )/mesh().V(),
462  rDtCoef.value()*
463  (
464  vf.boundaryField() - vf.oldTime().boundaryField()
465  ) - offCentre_(ff(ddt0.boundaryField()))
466  )
467  );
468  }
469  else
470  {
471  if (evaluate(ddt0))
472  {
473  ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime())
474  - offCentre_(ddt0());
475  }
476 
477  return tmp<GeometricField<Type, fvPatchField, volMesh>>
478  (
479  new GeometricField<Type, fvPatchField, volMesh>
480  (
481  ddtIOobject,
482  rDtCoef*(vf - vf.oldTime()) - offCentre_(ddt0())
483  )
484  );
485  }
486 }
487 
488 
489 template<class Type>
490 tmp<GeometricField<Type, fvPatchField, volMesh>>
491 CrankNicolsonDdtScheme<Type>::fvcDdt
492 (
493  const dimensionedScalar& rho,
494  const GeometricField<Type, fvPatchField, volMesh>& vf
495 )
496 {
497  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
498  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
499  (
500  "ddt0(" + rho.name() + ',' + vf.name() + ')',
501  rho.dimensions()*vf.dimensions()
502  );
503 
504  IOobject ddtIOobject
505  (
506  "ddt(" + rho.name() + ',' + vf.name() + ')',
507  mesh().time().timeName(),
508  mesh()
509  );
510 
511  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
512 
513  if (mesh().moving())
514  {
515  if (evaluate(ddt0))
516  {
517  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
518 
519  ddt0.primitiveFieldRef() =
520  (
521  rDtCoef0*rho.value()*
522  (
523  mesh().V0()*vf.oldTime().primitiveField()
524  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
525  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
526  )/mesh().V0();
527 
528  ddt0.boundaryFieldRef() =
529  (
530  rDtCoef0*rho.value()*
531  (
532  vf.oldTime().boundaryField()
533  - vf.oldTime().oldTime().boundaryField()
534  ) - offCentre_(ff(ddt0.boundaryField()))
535  );
536  }
537 
538  return tmp<GeometricField<Type, fvPatchField, volMesh>>
539  (
540  new GeometricField<Type, fvPatchField, volMesh>
541  (
542  ddtIOobject,
543  mesh(),
544  rDtCoef.dimensions()*rho.dimensions()*vf.dimensions(),
545  (
546  rDtCoef.value()*rho.value()*
547  (
548  mesh().V()*vf.primitiveField()
549  - mesh().V0()*vf.oldTime().primitiveField()
550  ) - mesh().V0()*offCentre_(ddt0.primitiveField())
551  )/mesh().V(),
552  rDtCoef.value()*rho.value()*
553  (
554  vf.boundaryField() - vf.oldTime().boundaryField()
555  ) - offCentre_(ff(ddt0.boundaryField()))
556  )
557  );
558  }
559  else
560  {
561  if (evaluate(ddt0))
562  {
563  ddt0 = rDtCoef0_(ddt0)*rho*(vf.oldTime() - vf.oldTime().oldTime())
564  - offCentre_(ddt0());
565  }
566 
567  return tmp<GeometricField<Type, fvPatchField, volMesh>>
568  (
569  new GeometricField<Type, fvPatchField, volMesh>
570  (
571  ddtIOobject,
572  rDtCoef*rho*(vf - vf.oldTime()) - offCentre_(ddt0())
573  )
574  );
575  }
576 }
577 
578 
579 template<class Type>
580 tmp<GeometricField<Type, fvPatchField, volMesh>>
581 CrankNicolsonDdtScheme<Type>::fvcDdt
582 (
583  const volScalarField& rho,
584  const GeometricField<Type, fvPatchField, volMesh>& vf
585 )
586 {
587  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
588  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
589  (
590  "ddt0(" + rho.name() + ',' + vf.name() + ')',
591  rho.dimensions()*vf.dimensions()
592  );
593 
594  IOobject ddtIOobject
595  (
596  "ddt(" + rho.name() + ',' + vf.name() + ')',
597  mesh().time().timeName(),
598  mesh()
599  );
600 
601  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
602 
603  if (mesh().moving())
604  {
605  if (evaluate(ddt0))
606  {
607  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
608 
609  ddt0.primitiveFieldRef() =
610  (
611  rDtCoef0*
612  (
613  mesh().V0()*rho.oldTime().primitiveField()
614  *vf.oldTime().primitiveField()
615  - mesh().V00()*rho.oldTime().oldTime().primitiveField()
616  *vf.oldTime().oldTime().primitiveField()
617  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
618  )/mesh().V0();
619 
620  ddt0.boundaryFieldRef() =
621  (
622  rDtCoef0*
623  (
624  rho.oldTime().boundaryField()
625  *vf.oldTime().boundaryField()
626  - rho.oldTime().oldTime().boundaryField()
627  *vf.oldTime().oldTime().boundaryField()
628  ) - offCentre_(ff(ddt0.boundaryField()))
629  );
630  }
631 
632  return tmp<GeometricField<Type, fvPatchField, volMesh>>
633  (
634  new GeometricField<Type, fvPatchField, volMesh>
635  (
636  ddtIOobject,
637  mesh(),
638  rDtCoef.dimensions()*rho.dimensions()*vf.dimensions(),
639  (
640  rDtCoef.value()*
641  (
642  mesh().V()*rho.primitiveField()*vf.primitiveField()
643  - mesh().V0()*rho.oldTime().primitiveField()
644  *vf.oldTime().primitiveField()
645  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
646  )/mesh().V(),
647  rDtCoef.value()*
648  (
649  rho.boundaryField()*vf.boundaryField()
650  - rho.oldTime().boundaryField()*vf.oldTime().boundaryField()
651  ) - offCentre_(ff(ddt0.boundaryField()))
652  )
653  );
654  }
655  else
656  {
657  if (evaluate(ddt0))
658  {
659  ddt0 = rDtCoef0_(ddt0)*
660  (
661  rho.oldTime()*vf.oldTime()
662  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
663  ) - offCentre_(ddt0());
664  }
665 
666  return tmp<GeometricField<Type, fvPatchField, volMesh>>
667  (
668  new GeometricField<Type, fvPatchField, volMesh>
669  (
670  ddtIOobject,
671  rDtCoef*(rho*vf - rho.oldTime()*vf.oldTime())
672  - offCentre_(ddt0())
673  )
674  );
675  }
676 }
677 
678 
679 template<class Type>
680 tmp<GeometricField<Type, fvPatchField, volMesh>>
681 CrankNicolsonDdtScheme<Type>::fvcDdt
682 (
683  const volScalarField& alpha,
684  const volScalarField& rho,
685  const GeometricField<Type, fvPatchField, volMesh>& vf
686 )
687 {
688  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
689  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
690  (
691  "ddt0(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
692  alpha.dimensions()*rho.dimensions()*vf.dimensions()
693  );
694 
695  IOobject ddtIOobject
696  (
697  "ddt(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
698  mesh().time().timeName(),
699  mesh()
700  );
701 
702  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
703 
704  if (mesh().moving())
705  {
706  if (evaluate(ddt0))
707  {
708  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
709 
710  ddt0.primitiveFieldRef() =
711  (
712  rDtCoef0*
713  (
714  mesh().V0()
715  *alpha.oldTime().primitiveField()
716  *rho.oldTime().primitiveField()
717  *vf.oldTime().primitiveField()
718 
719  - mesh().V00()
720  *alpha.oldTime().oldTime().primitiveField()
721  *rho.oldTime().oldTime().primitiveField()
722  *vf.oldTime().oldTime().primitiveField()
723  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
724  )/mesh().V0();
725 
726  ddt0.boundaryFieldRef() =
727  (
728  rDtCoef0*
729  (
730  alpha.oldTime().boundaryField()
731  *rho.oldTime().boundaryField()
732  *vf.oldTime().boundaryField()
733 
734  - alpha.oldTime().oldTime().boundaryField()
735  *rho.oldTime().oldTime().boundaryField()
736  *vf.oldTime().oldTime().boundaryField()
737  ) - offCentre_(ff(ddt0.boundaryField()))
738  );
739  }
740 
741  return tmp<GeometricField<Type, fvPatchField, volMesh>>
742  (
743  new GeometricField<Type, fvPatchField, volMesh>
744  (
745  ddtIOobject,
746  mesh(),
747  rDtCoef.dimensions()
748  *alpha.dimensions()*rho.dimensions()*vf.dimensions(),
749  (
750  rDtCoef.value()*
751  (
752  mesh().V()
753  *alpha.primitiveField()
754  *rho.primitiveField()
755  *vf.primitiveField()
756 
757  - mesh().V0()
758  *alpha.oldTime().primitiveField()
759  *rho.oldTime().primitiveField()
760  *vf.oldTime().primitiveField()
761  ) - mesh().V00()*offCentre_(ddt0.primitiveField())
762  )/mesh().V(),
763  rDtCoef.value()*
764  (
765  alpha.boundaryField()
766  *rho.boundaryField()
767  *vf.boundaryField()
768 
769  - alpha.oldTime().boundaryField()
770  *rho.oldTime().boundaryField()
771  *vf.oldTime().boundaryField()
772  ) - offCentre_(ff(ddt0.boundaryField()))
773  )
774  );
775  }
776  else
777  {
778  if (evaluate(ddt0))
779  {
780  ddt0 = rDtCoef0_(ddt0)*
781  (
782  alpha.oldTime()
783  *rho.oldTime()
784  *vf.oldTime()
785 
786  - alpha.oldTime().oldTime()
787  *rho.oldTime().oldTime()
788  *vf.oldTime().oldTime()
789  ) - offCentre_(ddt0());
790  }
791 
792  return tmp<GeometricField<Type, fvPatchField, volMesh>>
793  (
794  new GeometricField<Type, fvPatchField, volMesh>
795  (
796  ddtIOobject,
797  rDtCoef
798  *(
799  alpha*rho*vf
800  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
801  )
802  - offCentre_(ddt0())
803  )
804  );
805  }
806 }
807 
808 
809 template<class Type>
810 tmp<fvMatrix<Type>>
811 CrankNicolsonDdtScheme<Type>::fvmDdt
812 (
813  const GeometricField<Type, fvPatchField, volMesh>& vf
814 )
815 {
816  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
817  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
818  (
819  "ddt0(" + vf.name() + ')',
820  vf.dimensions()
821  );
822 
823  tmp<fvMatrix<Type>> tfvm
824  (
825  new fvMatrix<Type>
826  (
827  vf,
828  vf.dimensions()*dimVol/dimTime
829  )
830  );
831 
832  fvMatrix<Type>& fvm = tfvm.ref();
833 
834  const scalar rDtCoef = rDtCoef_(ddt0).value();
835  fvm.diag() = rDtCoef*mesh().V();
836 
837  vf.oldTime().oldTime();
838 
839  if (mesh().moving())
840  {
841  if (evaluate(ddt0))
842  {
843  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
844 
845  ddt0.primitiveFieldRef() =
846  (
847  rDtCoef0*
848  (
849  mesh().V0()*vf.oldTime().primitiveField()
850  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
851  )
852  - mesh().V00()*offCentre_(ddt0.primitiveField())
853  )/mesh().V0();
854 
855  ddt0.boundaryFieldRef() =
856  (
857  rDtCoef0*
858  (
859  vf.oldTime().boundaryField()
860  - vf.oldTime().oldTime().boundaryField()
861  )
862  - offCentre_(ff(ddt0.boundaryField()))
863  );
864  }
865 
866  fvm.source() =
867  (
868  rDtCoef*vf.oldTime().primitiveField()
869  + offCentre_(ddt0.primitiveField())
870  )*mesh().V0();
871  }
872  else
873  {
874  if (evaluate(ddt0))
875  {
876  ddt0 = rDtCoef0_(ddt0)*(vf.oldTime() - vf.oldTime().oldTime())
877  - offCentre_(ddt0());
878 
879  }
880 
881  fvm.source() =
882  (
883  rDtCoef*vf.oldTime().primitiveField()
884  + offCentre_(ddt0.primitiveField())
885  )*mesh().V();
886  }
887 
888  return tfvm;
889 }
890 
891 
892 template<class Type>
893 tmp<fvMatrix<Type>>
894 CrankNicolsonDdtScheme<Type>::fvmDdt
895 (
896  const dimensionedScalar& rho,
897  const GeometricField<Type, fvPatchField, volMesh>& vf
898 )
899 {
900  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
901  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
902  (
903  "ddt0(" + rho.name() + ',' + vf.name() + ')',
904  rho.dimensions()*vf.dimensions()
905  );
906 
907  tmp<fvMatrix<Type>> tfvm
908  (
909  new fvMatrix<Type>
910  (
911  vf,
912  rho.dimensions()*vf.dimensions()*dimVol/dimTime
913  )
914  );
915  fvMatrix<Type>& fvm = tfvm.ref();
916 
917  const scalar rDtCoef = rDtCoef_(ddt0).value();
918  fvm.diag() = rDtCoef*rho.value()*mesh().V();
919 
920  vf.oldTime().oldTime();
921 
922  if (mesh().moving())
923  {
924  if (evaluate(ddt0))
925  {
926  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
927 
928  ddt0.primitiveFieldRef() =
929  (
930  rDtCoef0*rho.value()*
931  (
932  mesh().V0()*vf.oldTime().primitiveField()
933  - mesh().V00()*vf.oldTime().oldTime().primitiveField()
934  )
935  - mesh().V00()*offCentre_(ddt0.primitiveField())
936  )/mesh().V0();
937 
938  ddt0.boundaryFieldRef() =
939  (
940  rDtCoef0*rho.value()*
941  (
942  vf.oldTime().boundaryField()
943  - vf.oldTime().oldTime().boundaryField()
944  )
945  - offCentre_(ff(ddt0.boundaryField()))
946  );
947  }
948 
949  fvm.source() =
950  (
951  rDtCoef*rho.value()*vf.oldTime().primitiveField()
952  + offCentre_(ddt0.primitiveField())
953  )*mesh().V0();
954  }
955  else
956  {
957  if (evaluate(ddt0))
958  {
959  ddt0 = rDtCoef0_(ddt0)*rho*(vf.oldTime() - vf.oldTime().oldTime())
960  - offCentre_(ddt0());
961  }
962 
963  fvm.source() =
964  (
965  rDtCoef*rho.value()*vf.oldTime().primitiveField()
966  + offCentre_(ddt0.primitiveField())
967  )*mesh().V();
968  }
969 
970  return tfvm;
971 }
972 
973 
974 template<class Type>
975 tmp<fvMatrix<Type>>
976 CrankNicolsonDdtScheme<Type>::fvmDdt
977 (
978  const volScalarField& rho,
979  const GeometricField<Type, fvPatchField, volMesh>& vf
980 )
981 {
982  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
983  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
984  (
985  "ddt0(" + rho.name() + ',' + vf.name() + ')',
986  rho.dimensions()*vf.dimensions()
987  );
988 
989  tmp<fvMatrix<Type>> tfvm
990  (
991  new fvMatrix<Type>
992  (
993  vf,
994  rho.dimensions()*vf.dimensions()*dimVol/dimTime
995  )
996  );
997  fvMatrix<Type>& fvm = tfvm.ref();
998 
999  const scalar rDtCoef = rDtCoef_(ddt0).value();
1000  fvm.diag() = rDtCoef*rho.primitiveField()*mesh().V();
1001 
1002  vf.oldTime().oldTime();
1003  rho.oldTime().oldTime();
1004 
1005  if (mesh().moving())
1006  {
1007  if (evaluate(ddt0))
1008  {
1009  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
1010 
1011  ddt0.primitiveFieldRef() =
1012  (
1013  rDtCoef0*
1014  (
1015  mesh().V0()*rho.oldTime().primitiveField()
1016  *vf.oldTime().primitiveField()
1017  - mesh().V00()*rho.oldTime().oldTime().primitiveField()
1018  *vf.oldTime().oldTime().primitiveField()
1019  )
1020  - mesh().V00()*offCentre_(ddt0.primitiveField())
1021  )/mesh().V0();
1022 
1023  ddt0.boundaryFieldRef() =
1024  (
1025  rDtCoef0*
1026  (
1027  rho.oldTime().boundaryField()
1028  *vf.oldTime().boundaryField()
1029  - rho.oldTime().oldTime().boundaryField()
1030  *vf.oldTime().oldTime().boundaryField()
1031  )
1032  - offCentre_(ff(ddt0.boundaryField()))
1033  );
1034  }
1035 
1036  fvm.source() =
1037  (
1038  rDtCoef*rho.oldTime().primitiveField()*vf.oldTime().primitiveField()
1039  + offCentre_(ddt0.primitiveField())
1040  )*mesh().V0();
1041  }
1042  else
1043  {
1044  if (evaluate(ddt0))
1045  {
1046  ddt0 = rDtCoef0_(ddt0)*
1047  (
1048  rho.oldTime()*vf.oldTime()
1049  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
1050  ) - offCentre_(ddt0());
1051  }
1052 
1053  fvm.source() =
1054  (
1055  rDtCoef*rho.oldTime().primitiveField()*vf.oldTime().primitiveField()
1056  + offCentre_(ddt0.primitiveField())
1057  )*mesh().V();
1058  }
1059 
1060  return tfvm;
1061 }
1062 
1063 
1064 template<class Type>
1065 tmp<fvMatrix<Type>>
1066 CrankNicolsonDdtScheme<Type>::fvmDdt
1067 (
1068  const volScalarField& alpha,
1069  const volScalarField& rho,
1070  const GeometricField<Type, fvPatchField, volMesh>& vf
1071 )
1072 {
1073  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1074  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1075  (
1076  "ddt0(" + alpha.name() + ',' + rho.name() + ',' + vf.name() + ')',
1077  alpha.dimensions()*rho.dimensions()*vf.dimensions()
1078  );
1079 
1080  tmp<fvMatrix<Type>> tfvm
1081  (
1082  new fvMatrix<Type>
1083  (
1084  vf,
1085  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
1086  )
1087  );
1088  fvMatrix<Type>& fvm = tfvm.ref();
1089 
1090  const scalar rDtCoef = rDtCoef_(ddt0).value();
1091  fvm.diag() = rDtCoef*alpha.primitiveField()*rho.primitiveField()*mesh().V();
1092 
1093  vf.oldTime().oldTime();
1094  alpha.oldTime().oldTime();
1095  rho.oldTime().oldTime();
1096 
1097  if (mesh().moving())
1098  {
1099  if (evaluate(ddt0))
1100  {
1101  const scalar rDtCoef0 = rDtCoef0_(ddt0).value();
1102 
1103  ddt0.primitiveFieldRef() =
1104  (
1105  rDtCoef0*
1106  (
1107  mesh().V0()
1108  *alpha.oldTime().primitiveField()
1109  *rho.oldTime().primitiveField()
1110  *vf.oldTime().primitiveField()
1111 
1112  - mesh().V00()
1113  *alpha.oldTime().oldTime().primitiveField()
1114  *rho.oldTime().oldTime().primitiveField()
1115  *vf.oldTime().oldTime().primitiveField()
1116  )
1117  - mesh().V00()*offCentre_(ddt0.primitiveField())
1118  )/mesh().V0();
1119 
1120  ddt0.boundaryFieldRef() =
1121  (
1122  rDtCoef0*
1123  (
1124  alpha.oldTime().boundaryField()
1125  *rho.oldTime().boundaryField()
1126  *vf.oldTime().boundaryField()
1127 
1128  - alpha.oldTime().oldTime().boundaryField()
1129  *rho.oldTime().oldTime().boundaryField()
1130  *vf.oldTime().oldTime().boundaryField()
1131  )
1132  - offCentre_(ff(ddt0.boundaryField()))
1133  );
1134  }
1135 
1136  fvm.source() =
1137  (
1138  rDtCoef
1139  *alpha.oldTime().primitiveField()
1140  *rho.oldTime().primitiveField()
1141  *vf.oldTime().primitiveField()
1142  + offCentre_(ddt0.primitiveField())
1143  )*mesh().V0();
1144  }
1145  else
1146  {
1147  if (evaluate(ddt0))
1148  {
1149  ddt0 = rDtCoef0_(ddt0)*
1150  (
1151  alpha.oldTime()
1152  *rho.oldTime()
1153  *vf.oldTime()
1154 
1155  - alpha.oldTime().oldTime()
1156  *rho.oldTime().oldTime()
1157  *vf.oldTime().oldTime()
1158  ) - offCentre_(ddt0());
1159  }
1160 
1161  fvm.source() =
1162  (
1163  rDtCoef
1164  *alpha.oldTime().primitiveField()
1165  *rho.oldTime().primitiveField()
1166  *vf.oldTime().primitiveField()
1167  + offCentre_(ddt0.primitiveField())
1168  )*mesh().V();
1169  }
1170 
1171  return tfvm;
1172 }
1173 
1174 
1175 template<class Type>
1176 tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1177 CrankNicolsonDdtScheme<Type>::fvcDdtUfCorr
1178 (
1179  const GeometricField<Type, fvPatchField, volMesh>& U,
1180  const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
1181 )
1182 {
1183  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1184  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1185  (
1186  "ddtCorrDdt0(" + U.name() + ')',
1187  U.dimensions()
1188  );
1189 
1190  DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 =
1191  ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>>
1192  (
1193  "ddtCorrDdt0(" + Uf.name() + ')',
1194  Uf.dimensions()
1195  );
1196 
1197  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1198 
1199  if (evaluate(ddt0))
1200  {
1201  ddt0 =
1202  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1203  - offCentre_(ddt0());
1204  }
1205 
1206  if (evaluate(dUfdt0))
1207  {
1208  dUfdt0 =
1209  rDtCoef0_(dUfdt0)*(Uf.oldTime() - Uf.oldTime().oldTime())
1210  - offCentre_(dUfdt0());
1211  }
1212 
1213  return tmp<fluxFieldType>
1214  (
1215  new fluxFieldType
1216  (
1217  IOobject
1218  (
1219  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
1220  mesh().time().timeName(),
1221  mesh()
1222  ),
1223  this->fvcDdtPhiCoeff(U.oldTime(), mesh().Sf() & Uf.oldTime())
1224  *(
1225  mesh().Sf()
1226  & (
1227  (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1228  - fvc::interpolate(rDtCoef*U.oldTime() + offCentre_(ddt0()))
1229  )
1230  )
1231  )
1232  );
1233 }
1234 
1235 
1236 template<class Type>
1237 tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1238 CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
1239 (
1240  const GeometricField<Type, fvPatchField, volMesh>& U,
1241  const fluxFieldType& phi
1242 )
1243 {
1244  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1245  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1246  (
1247  "ddtCorrDdt0(" + U.name() + ')',
1248  U.dimensions()
1249  );
1250 
1251  DDt0Field<fluxFieldType>& dphidt0 =
1252  ddt0_<fluxFieldType>
1253  (
1254  "ddtCorrDdt0(" + phi.name() + ')',
1255  phi.dimensions()
1256  );
1257  dphidt0.setOriented();
1258 
1259  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1260 
1261  if (evaluate(ddt0))
1262  {
1263  ddt0 =
1264  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1265  - offCentre_(ddt0());
1266  }
1267 
1268  if (evaluate(dphidt0))
1269  {
1270  dphidt0 =
1271  rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime())
1272  - offCentre_(dphidt0());
1273  }
1274 
1275  return tmp<fluxFieldType>
1276  (
1277  new fluxFieldType
1278  (
1279  IOobject
1280  (
1281  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
1282  mesh().time().timeName(),
1283  mesh()
1284  ),
1285  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
1286  *(
1287  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1288  - fvc::dotInterpolate
1289  (
1290  mesh().Sf(),
1291  rDtCoef*U.oldTime() + offCentre_(ddt0())
1292  )
1293  )
1294  )
1295  );
1296 }
1297 
1298 
1299 template<class Type>
1300 tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1301 CrankNicolsonDdtScheme<Type>::fvcDdtUfCorr
1302 (
1303  const volScalarField& rho,
1304  const GeometricField<Type, fvPatchField, volMesh>& U,
1305  const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
1306 )
1307 {
1308  if
1309  (
1310  U.dimensions() == dimVelocity
1311  && Uf.dimensions() == rho.dimensions()*dimVelocity
1312  )
1313  {
1314  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1315  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1316  (
1317  "ddtCorrDdt0(" + rho.name() + ',' + U.name() + ')',
1318  rho.dimensions()*U.dimensions()
1319  );
1320 
1321  DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 =
1322  ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>>
1323  (
1324  "ddtCorrDdt0(" + Uf.name() + ')',
1325  Uf.dimensions()
1326  );
1327 
1328  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1329 
1330  GeometricField<Type, fvPatchField, volMesh> rhoU0
1331  (
1332  rho.oldTime()*U.oldTime()
1333  );
1334 
1335  if (evaluate(ddt0))
1336  {
1337  ddt0 =
1338  rDtCoef0_(ddt0)
1339  *(rhoU0 - rho.oldTime().oldTime()*U.oldTime().oldTime())
1340  - offCentre_(ddt0());
1341  }
1342 
1343  if (evaluate(dUfdt0))
1344  {
1345  dUfdt0 =
1346  rDtCoef0_(dUfdt0)
1347  *(Uf.oldTime() - Uf.oldTime().oldTime())
1348  - offCentre_(dUfdt0());
1349  }
1350 
1351  tmp<fluxFieldType> ddtCorr
1352  (
1353  new fluxFieldType
1354  (
1355  IOobject
1356  (
1357  "ddtCorr("
1358  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
1359  mesh().time().timeName(),
1360  mesh()
1361  ),
1362  this->fvcDdtPhiCoeff
1363  (
1364  rhoU0,
1365  mesh().Sf() & Uf.oldTime(),
1366  rho.oldTime()
1367  )
1368  *(
1369  mesh().Sf()
1370  & (
1371  (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1372  - fvc::interpolate(rDtCoef*rhoU0 + offCentre_(ddt0()))
1373  )
1374  )
1375  )
1376  );
1377 
1378  return ddtCorr;
1379  }
1380  else if
1381  (
1382  U.dimensions() == rho.dimensions()*dimVelocity
1383  && Uf.dimensions() == rho.dimensions()*dimVelocity
1384  )
1385  {
1386  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1387  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1388  (
1389  "ddtCorrDdt0(" + U.name() + ')',
1390  U.dimensions()
1391  );
1392 
1393  DDt0Field<GeometricField<Type, fvsPatchField, surfaceMesh>>& dUfdt0 =
1394  ddt0_<GeometricField<Type, fvsPatchField, surfaceMesh>>
1395  (
1396  "ddtCorrDdt0(" + Uf.name() + ')',
1397  Uf.dimensions()
1398  );
1399 
1400  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1401 
1402  if (evaluate(ddt0))
1403  {
1404  ddt0 =
1405  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1406  - offCentre_(ddt0());
1407  }
1408 
1409  if (evaluate(dUfdt0))
1410  {
1411  dUfdt0 =
1412  rDtCoef0_(dUfdt0)*(Uf.oldTime() - Uf.oldTime().oldTime())
1413  - offCentre_(dUfdt0());
1414  }
1415 
1416  return tmp<fluxFieldType>
1417  (
1418  new fluxFieldType
1419  (
1420  IOobject
1421  (
1422  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
1423  mesh().time().timeName(),
1424  mesh()
1425  ),
1426  this->fvcDdtPhiCoeff
1427  (
1428  U.oldTime(),
1429  mesh().Sf() & Uf.oldTime(),
1430  rho.oldTime()
1431  )
1432  *(
1433  mesh().Sf()
1434  & (
1435  (rDtCoef*Uf.oldTime() + offCentre_(dUfdt0()))
1436  - fvc::interpolate
1437  (
1438  rDtCoef*U.oldTime() + offCentre_(ddt0())
1439  )
1440  )
1441  )
1442  )
1443  );
1444  }
1445  else
1446  {
1447  FatalErrorInFunction
1448  << "dimensions of Uf are not correct"
1449  << abort(FatalError);
1450 
1451  return fluxFieldType::null();
1452  }
1453 }
1454 
1455 
1456 template<class Type>
1457 tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
1458 CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
1459 (
1460  const volScalarField& rho,
1461  const GeometricField<Type, fvPatchField, volMesh>& U,
1462  const fluxFieldType& phi
1463 )
1464 {
1465  if
1466  (
1467  U.dimensions() == dimVelocity
1468  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1469  )
1470  {
1471  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1472  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1473  (
1474  "ddtCorrDdt0(" + rho.name() + ',' + U.name() + ')',
1475  rho.dimensions()*U.dimensions()
1476  );
1477 
1478  DDt0Field<fluxFieldType>& dphidt0 =
1479  ddt0_<fluxFieldType>
1480  (
1481  "ddtCorrDdt0(" + phi.name() + ')',
1482  phi.dimensions()
1483  );
1484 
1485  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1486 
1487  GeometricField<Type, fvPatchField, volMesh> rhoU0
1488  (
1489  rho.oldTime()*U.oldTime()
1490  );
1491 
1492  if (evaluate(ddt0))
1493  {
1494  ddt0 =
1495  rDtCoef0_(ddt0)
1496  *(rhoU0 - rho.oldTime().oldTime()*U.oldTime().oldTime())
1497  - offCentre_(ddt0());
1498  }
1499 
1500  if (evaluate(dphidt0))
1501  {
1502  dphidt0 =
1503  rDtCoef0_(dphidt0)
1504  *(phi.oldTime() - phi.oldTime().oldTime())
1505  - offCentre_(dphidt0());
1506  }
1507 
1508  tmp<fluxFieldType> ddtCorr
1509  (
1510  new fluxFieldType
1511  (
1512  IOobject
1513  (
1514  "ddtCorr("
1515  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
1516  mesh().time().timeName(),
1517  mesh()
1518  ),
1519  this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), rho.oldTime())
1520  *(
1521  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1522  - fvc::dotInterpolate
1523  (
1524  mesh().Sf(),
1525  rDtCoef*rhoU0 + offCentre_(ddt0())
1526  )
1527  )
1528  )
1529  );
1530 
1531  return ddtCorr;
1532  }
1533  else if
1534  (
1535  U.dimensions() == rho.dimensions()*dimVelocity
1536  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
1537  )
1538  {
1539  DDt0Field<GeometricField<Type, fvPatchField, volMesh>>& ddt0 =
1540  ddt0_<GeometricField<Type, fvPatchField, volMesh>>
1541  (
1542  "ddtCorrDdt0(" + U.name() + ')',
1543  U.dimensions()
1544  );
1545 
1546  DDt0Field<fluxFieldType>& dphidt0 =
1547  ddt0_<fluxFieldType>
1548  (
1549  "ddtCorrDdt0(" + phi.name() + ')',
1550  phi.dimensions()
1551  );
1552 
1553  dimensionedScalar rDtCoef = rDtCoef_(ddt0);
1554 
1555  if (evaluate(ddt0))
1556  {
1557  ddt0 =
1558  rDtCoef0_(ddt0)*(U.oldTime() - U.oldTime().oldTime())
1559  - offCentre_(ddt0());
1560  }
1561 
1562  if (evaluate(dphidt0))
1563  {
1564  dphidt0 =
1565  rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime())
1566  - offCentre_(dphidt0());
1567  }
1568 
1569  return tmp<fluxFieldType>
1570  (
1571  new fluxFieldType
1572  (
1573  IOobject
1574  (
1575  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
1576  mesh().time().timeName(),
1577  mesh()
1578  ),
1579  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), rho.oldTime())
1580  *(
1581  (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
1582  - fvc::dotInterpolate
1583  (
1584  mesh().Sf(),
1585  rDtCoef*U.oldTime() + offCentre_(ddt0())
1586  )
1587  )
1588  )
1589  );
1590  }
1591  else
1592  {
1593  FatalErrorInFunction
1594  << "dimensions of phi are not correct"
1595  << abort(FatalError);
1596 
1597  return fluxFieldType::null();
1598  }
1599 }
1600 
1601 
1602 template<class Type>
1603 tmp<surfaceScalarField> CrankNicolsonDdtScheme<Type>::meshPhi
1604 (
1605  const GeometricField<Type, fvPatchField, volMesh>& vf
1606 )
1607 {
1608  DDt0Field<surfaceScalarField>& meshPhi0 = ddt0_<surfaceScalarField>
1609  (
1610  "meshPhiCN_0",
1611  dimVolume
1612  );
1613 
1614  meshPhi0.setOriented();
1615 
1616  if (evaluate(meshPhi0))
1617  {
1618  meshPhi0 =
1619  coef0_(meshPhi0)*mesh().phi().oldTime() - offCentre_(meshPhi0());
1620  }
1621 
1622  return tmp<surfaceScalarField>
1623  (
1624  new surfaceScalarField
1625  (
1626  IOobject
1627  (
1628  mesh().phi().name(),
1629  mesh().time().timeName(),
1630  mesh(),
1631  IOobject::NO_READ,
1632  IOobject::NO_WRITE,
1633  false
1634  ),
1635  coef_(meshPhi0)*mesh().phi() - offCentre_(meshPhi0())
1636  )
1637  );
1638 }
1639 
1640 
1641 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1642 
1643 } // End namespace fv
1644 
1645 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1646 
1647 } // End namespace Foam
1648 
1649 // ************************************************************************* //
Foam::IOobject::NO_WRITE
@ NO_WRITE
Definition: IOobject.H:191
Foam::fv::CrankNicolsonDdtScheme::mesh
const fvMesh & mesh() const
Definition: CrankNicolsonDdtScheme.H:224
runTime
engineTime & runTime
Definition: createEngineTime.H:13
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:165
Foam::IOobject::AUTO_WRITE
@ AUTO_WRITE
Definition: IOobject.H:190
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:69
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:63
Foam::FieldField
A field of fields is a PtrList of fields with reference counting.
Definition: FieldField.H:49
Foam::fvc::dotInterpolate
static tmp< GeometricField< typename innerProduct< vector, Type >::type, fvsPatchField, surfaceMesh > > dotInterpolate(const surfaceVectorField &Sf, const GeometricField< Type, fvPatchField, volMesh > &tvf)
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:57
Foam::Zero
static constexpr const zero Zero
Definition: zero.H:131
Foam::dimVelocity
const dimensionSet dimVelocity
Foam::constant::atomic::alpha
const dimensionedScalar alpha
Definition: readThermalProperties.H:212
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const
Definition: GeometricFieldI.H:46
Foam::fv::CrankNicolsonDdtScheme::fvcDdtUfCorr
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition: CrankNicolsonDdtScheme.C:1171
Uf
autoPtr< surfaceVectorField > Uf
Definition: createUfIfPresent.H:27
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Definition: GeometricField.C:843
fvcDiv.H
Calculate the divergence of the given field.
Foam::dimensioned::name
const word & name() const
Definition: dimensionedType.C:399
Foam::FatalIOError
IOerror FatalIOError
Foam::token
A token holds an item read from Istream.
Definition: token.H:64
Foam::dimensioned::value
const Type & value() const
Definition: dimensionedType.C:427
Foam::dimensionSet
Dimension set for the base types, which can be used to implement rigorous dimension checking for alge...
Definition: dimensionSet.H:104
Foam::fvMesh::V00
const DimensionedField< scalar, volMesh > & V00() const
Definition: fvMeshGeometry.C:226
Foam::Function1
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
Definition: propellerInfo.H:287
Foam::regIOobject::store
bool store()
Definition: regIOobjectI.H:30
fvMatrices.H
A special matrix type and solver, designed for finite volume solutions of scalar equations.
Foam::token::isNumber
bool isNumber() const noexcept
Definition: tokenI.H:580
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
Foam::fv::CrankNicolsonDdtScheme::ocCoeff
scalar ocCoeff() const
Definition: CrankNicolsonDdtScheme.H:230
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:220
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Foam::token::number
scalar number() const
Definition: tokenI.H:586
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::Function1Types::Constant
Templated function that returns a constant value.
Definition: Constant.H:68
Foam::Istream
An Istream is an abstract base class for all input systems (streams, files, token lists etc)....
Definition: Istream.H:57
rho
rho
Definition: readInitialConditions.H:88
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:36
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::fv::CrankNicolsonDdtScheme::fvmDdt
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: CrankNicolsonDdtScheme.C:805
ddtCorr
ddtCorr
Definition: readControls.H:9
Foam::fv::CrankNicolsonDdtScheme::fvcDdt
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
Definition: CrankNicolsonDdtScheme.C:341
timeName
word timeName
Definition: getTimeIndex.H:3
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:119
CrankNicolsonDdtScheme.H
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::fv::CrankNicolsonDdtScheme::fvcDdtPhiCorr
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
Definition: CrankNicolsonDdtScheme.C:1232
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedScalarFwd.H:36
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:81
Foam
Definition: atmBoundaryLayer.C:26
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:139
Foam::fvc::interpolate
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Foam::fv::CrankNicolsonDdtScheme::meshPhi
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: CrankNicolsonDdtScheme.C:1597
ocCoeff
scalar ocCoeff
Definition: alphaEqn.H:6
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:125
fv
labelList fv(nPoints)
U
U
Definition: pEqn.H:72
Foam::fv::ff
const FieldField< fvPatchField, Type > & ff(const FieldField< fvPatchField, Type > &bf)
Definition: CrankNicolsonDdtScheme.C:268
FatalErrorInFunction
#define FatalErrorInFunction
Definition: error.H:465
Foam::surfaceScalarField
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Definition: surfaceFieldsFwd.H:50
Foam::Istream::putBack
void putBack(const token &tok)
Definition: Istream.C:63
Foam::fvMatrix::source
Field< Type > & source()
Definition: fvMatrix.H:441
Foam::polyMesh::moving
bool moving() const noexcept
Definition: polyMesh.H:518
Foam::fvMatrix< Type >
Constant.H
Foam::fv::ddtScheme
Abstract base class for ddt schemes.
Definition: ddtScheme.H:79
Foam::name
word name(const expressions::valueTypeCode typeCode)
Definition: exprTraits.C:52
Foam::dimVol
const dimensionSet dimVol(dimVolume)
Definition: dimensionSets.H:65
FatalIOErrorInFunction
#define FatalIOErrorInFunction(ios)
Definition: error.H:494
timeIndex
label timeIndex
Definition: getTimeIndex.H:24
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const
Definition: dimensionedType.C:413
Foam::dimVolume
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
Foam::GeometricField
Generic GeometricField class.
Definition: areaFieldsFwd.H:49
Foam::IOobject::NO_READ
@ NO_READ
Definition: IOobject.H:184
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Definition: GeometricFieldI.H:55
Foam::IOobject::MUST_READ
@ MUST_READ
Definition: IOobject.H:181
Foam::stringOps::evaluate
string evaluate(label fieldWidth, const std::string &s, size_t pos=0, size_t len=std::string::npos)
Definition: stringOpsEvaluate.C:30

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.17