FreeFOAM The Cross-Platform CFD Toolkit
CoEulerDdtScheme.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 | Copyright (C) 1991-2010 OpenCFD Ltd.
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 "CoEulerDdtScheme.H"
28 #include <finiteVolume/fvcDiv.H>
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 namespace fv
39 {
40 
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 
43 template<class Type>
44 tmp<volScalarField> CoEulerDdtScheme<Type>::CorDeltaT() const
45 {
46  surfaceScalarField cofrDeltaT = CofrDeltaT();
47 
48  tmp<volScalarField> tcorDeltaT
49  (
50  new volScalarField
51  (
52  IOobject
53  (
54  "CorDeltaT",
55  cofrDeltaT.instance(),
56  mesh()
57  ),
58  mesh(),
59  dimensionedScalar("CorDeltaT", cofrDeltaT.dimensions(), 0.0),
60  zeroGradientFvPatchScalarField::typeName
61  )
62  );
63 
64  volScalarField& corDeltaT = tcorDeltaT();
65 
66  const unallocLabelList& owner = mesh().owner();
67  const unallocLabelList& neighbour = mesh().neighbour();
68 
69  forAll(owner, faceI)
70  {
71  corDeltaT[owner[faceI]] =
72  max(corDeltaT[owner[faceI]], cofrDeltaT[faceI]);
73 
74  corDeltaT[neighbour[faceI]] =
75  max(corDeltaT[neighbour[faceI]], cofrDeltaT[faceI]);
76  }
77 
78  forAll(corDeltaT.boundaryField(), patchi)
79  {
80  const fvsPatchScalarField& pcofrDeltaT =
81  cofrDeltaT.boundaryField()[patchi];
82 
83  const fvPatch& p = pcofrDeltaT.patch();
84  const unallocLabelList& faceCells = p.patch().faceCells();
85 
86  forAll(pcofrDeltaT, patchFacei)
87  {
88  corDeltaT[faceCells[patchFacei]] = max
89  (
90  corDeltaT[faceCells[patchFacei]],
91  pcofrDeltaT[patchFacei]
92  );
93  }
94  }
95 
96  corDeltaT.correctBoundaryConditions();
97 
98  //corDeltaT = max(corDeltaT, max(corDeltaT)/100.0);
99 
100  return tcorDeltaT;
101 }
102 
103 
104 template<class Type>
105 tmp<surfaceScalarField> CoEulerDdtScheme<Type>::CofrDeltaT() const
106 {
107  const dimensionedScalar& deltaT = mesh().time().deltaT();
108 
109  const surfaceScalarField& phi =
110  static_cast<const objectRegistry&>(mesh())
111  .lookupObject<surfaceScalarField>(phiName_);
112 
113  if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
114  {
116  (
118  *(mag(phi)/mesh().magSf())
119  *deltaT
120  );
121 
122  return max(Co/maxCo_, scalar(1))/deltaT;
123  }
124  else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
125  {
126  const volScalarField& rho =
127  static_cast<const objectRegistry&>(mesh())
128  .lookupObject<volScalarField>(rhoName_).oldTime();
129 
131  (
133  *(mag(phi)/(fvc::interpolate(rho)*mesh().magSf()))
134  *deltaT
135  );
136 
137  return max(Co/maxCo_, scalar(1))/deltaT;
138  }
139  else
140  {
141  FatalErrorIn("CoEulerDdtScheme<Type>::CofrDeltaT() const")
142  << "Incorrect dimensions of phi: " << phi.dimensions()
143  << abort(FatalError);
144 
145  return tmp<surfaceScalarField>(NULL);
146  }
147 }
148 
149 
150 template<class Type>
151 tmp<GeometricField<Type, fvPatchField, volMesh> >
153 (
154  const dimensioned<Type>& dt
155 )
156 {
157  volScalarField rDeltaT = CorDeltaT();
158 
159  IOobject ddtIOobject
160  (
161  "ddt("+dt.name()+')',
162  mesh().time().timeName(),
163  mesh()
164  );
165 
166  if (mesh().moving())
167  {
169  (
171  (
172  ddtIOobject,
173  mesh(),
175  (
176  "0",
177  dt.dimensions()/dimTime,
179  )
180  )
181  );
182 
183  tdtdt().internalField() =
184  rDeltaT.internalField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
185 
186  return tdtdt;
187  }
188  else
189  {
191  (
193  (
194  ddtIOobject,
195  mesh(),
197  (
198  "0",
199  dt.dimensions()/dimTime,
201  ),
203  )
204  );
205  }
206 }
207 
208 
209 template<class Type>
212 (
214 )
215 {
216  volScalarField rDeltaT = CorDeltaT();
217 
218  IOobject ddtIOobject
219  (
220  "ddt("+vf.name()+')',
221  mesh().time().timeName(),
222  mesh()
223  );
224 
225  if (mesh().moving())
226  {
228  (
230  (
231  ddtIOobject,
232  mesh(),
233  rDeltaT.dimensions()*vf.dimensions(),
234  rDeltaT.internalField()*
235  (
236  vf.internalField()
237  - vf.oldTime().internalField()*mesh().V0()/mesh().V()
238  ),
239  rDeltaT.boundaryField()*
240  (
241  vf.boundaryField() - vf.oldTime().boundaryField()
242  )
243  )
244  );
245  }
246  else
247  {
249  (
251  (
252  ddtIOobject,
253  rDeltaT*(vf - vf.oldTime())
254  )
255  );
256  }
257 }
258 
259 
260 template<class Type>
263 (
264  const dimensionedScalar& rho,
266 )
267 {
268  volScalarField rDeltaT = CorDeltaT();
269 
270  IOobject ddtIOobject
271  (
272  "ddt("+rho.name()+','+vf.name()+')',
273  mesh().time().timeName(),
274  mesh()
275  );
276 
277  if (mesh().moving())
278  {
280  (
282  (
283  ddtIOobject,
284  mesh(),
285  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
286  rDeltaT.internalField()*rho.value()*
287  (
288  vf.internalField()
289  - vf.oldTime().internalField()*mesh().V0()/mesh().V()
290  ),
291  rDeltaT.boundaryField()*rho.value()*
292  (
293  vf.boundaryField() - vf.oldTime().boundaryField()
294  )
295  )
296  );
297  }
298  else
299  {
301  (
303  (
304  ddtIOobject,
305  rDeltaT*rho*(vf - vf.oldTime())
306  )
307  );
308  }
309 }
310 
311 
312 template<class Type>
315 (
316  const volScalarField& rho,
318 )
319 {
320  volScalarField rDeltaT = CorDeltaT();
321 
322  IOobject ddtIOobject
323  (
324  "ddt("+rho.name()+','+vf.name()+')',
325  mesh().time().timeName(),
326  mesh()
327  );
328 
329  if (mesh().moving())
330  {
332  (
334  (
335  ddtIOobject,
336  mesh(),
337  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
338  rDeltaT.internalField()*
339  (
340  rho.internalField()*vf.internalField()
341  - rho.oldTime().internalField()
342  *vf.oldTime().internalField()*mesh().V0()/mesh().V()
343  ),
344  rDeltaT.boundaryField()*
345  (
346  rho.boundaryField()*vf.boundaryField()
347  - rho.oldTime().boundaryField()
348  *vf.oldTime().boundaryField()
349  )
350  )
351  );
352  }
353  else
354  {
356  (
358  (
359  ddtIOobject,
360  rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
361  )
362  );
363  }
364 }
365 
366 
367 template<class Type>
370 (
372 )
373 {
374  tmp<fvMatrix<Type> > tfvm
375  (
376  new fvMatrix<Type>
377  (
378  vf,
380  )
381  );
382 
383  fvMatrix<Type>& fvm = tfvm();
384 
385  scalarField rDeltaT = CorDeltaT()().internalField();
386 
387  fvm.diag() = rDeltaT*mesh().V();
388 
389  if (mesh().moving())
390  {
391  fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
392  }
393  else
394  {
395  fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V();
396  }
397 
398  return tfvm;
399 }
400 
401 
402 template<class Type>
405 (
406  const dimensionedScalar& rho,
408 )
409 {
410  tmp<fvMatrix<Type> > tfvm
411  (
412  new fvMatrix<Type>
413  (
414  vf,
416  )
417  );
418  fvMatrix<Type>& fvm = tfvm();
419 
420  scalarField rDeltaT = CorDeltaT()().internalField();
421 
422  fvm.diag() = rDeltaT*rho.value()*mesh().V();
423 
424  if (mesh().moving())
425  {
426  fvm.source() = rDeltaT
427  *rho.value()*vf.oldTime().internalField()*mesh().V0();
428  }
429  else
430  {
431  fvm.source() = rDeltaT
432  *rho.value()*vf.oldTime().internalField()*mesh().V();
433  }
434 
435  return tfvm;
436 }
437 
438 
439 template<class Type>
442 (
443  const volScalarField& rho,
445 )
446 {
447  tmp<fvMatrix<Type> > tfvm
448  (
449  new fvMatrix<Type>
450  (
451  vf,
453  )
454  );
455  fvMatrix<Type>& fvm = tfvm();
456 
457  scalarField rDeltaT = CorDeltaT()().internalField();
458 
459  fvm.diag() = rDeltaT*rho.internalField()*mesh().V();
460 
461  if (mesh().moving())
462  {
463  fvm.source() = rDeltaT
464  *rho.oldTime().internalField()
465  *vf.oldTime().internalField()*mesh().V0();
466  }
467  else
468  {
469  fvm.source() = rDeltaT
470  *rho.oldTime().internalField()
471  *vf.oldTime().internalField()*mesh().V();
472  }
473 
474  return tfvm;
475 }
476 
477 
478 template<class Type>
481 (
482  const volScalarField& rA,
484  const fluxFieldType& phi
485 )
486 {
487  IOobject ddtIOobject
488  (
489  "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
490  mesh().time().timeName(),
491  mesh()
492  );
493 
494  if (mesh().moving())
495  {
496  return tmp<fluxFieldType>
497  (
498  new fluxFieldType
499  (
500  ddtIOobject,
501  mesh(),
502  dimensioned<typename flux<Type>::type>
503  (
504  "0",
505  rA.dimensions()*phi.dimensions()/dimTime,
507  )
508  )
509  );
510  }
511  else
512  {
513  volScalarField rDeltaT = CorDeltaT();
514 
515  return tmp<fluxFieldType>
516  (
517  new fluxFieldType
518  (
519  ddtIOobject,
520  fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
521  (
522  fvc::interpolate(rDeltaT*rA)*phi.oldTime()
523  - (fvc::interpolate(rDeltaT*rA*U.oldTime()) & mesh().Sf())
524  )
525  )
526  );
527  }
528 }
529 
530 
531 template<class Type>
534 (
535  const volScalarField& rA,
536  const volScalarField& rho,
538  const fluxFieldType& phi
539 )
540 {
541  IOobject ddtIOobject
542  (
543  "ddtPhiCorr("
544  + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
545  mesh().time().timeName(),
546  mesh()
547  );
548 
549  if (mesh().moving())
550  {
551  return tmp<fluxFieldType>
552  (
553  new fluxFieldType
554  (
555  ddtIOobject,
556  mesh(),
557  dimensioned<typename flux<Type>::type>
558  (
559  "0",
560  rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
562  )
563  )
564  );
565  }
566  else
567  {
568  volScalarField rDeltaT = CorDeltaT();
569 
570  if
571  (
572  U.dimensions() == dimVelocity
573  && phi.dimensions() == dimVelocity*dimArea
574  )
575  {
576  return tmp<fluxFieldType>
577  (
578  new fluxFieldType
579  (
580  ddtIOobject,
581  fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
582  *(
583  fvc::interpolate(rDeltaT*rA*rho.oldTime())*phi.oldTime()
584  - (fvc::interpolate(rDeltaT*rA*rho.oldTime()*U.oldTime())
585  & mesh().Sf())
586  )
587  )
588  );
589  }
590  else if
591  (
592  U.dimensions() == dimVelocity
594  )
595  {
596  return tmp<fluxFieldType>
597  (
598  new fluxFieldType
599  (
600  ddtIOobject,
601  fvcDdtPhiCoeff
602  (
603  U.oldTime(),
604  phi.oldTime()/fvc::interpolate(rho.oldTime())
605  )
606  *(
607  fvc::interpolate(rDeltaT*rA*rho.oldTime())
608  *phi.oldTime()/fvc::interpolate(rho.oldTime())
609  - (
611  (
612  rDeltaT*rA*rho.oldTime()*U.oldTime()
613  ) & mesh().Sf()
614  )
615  )
616  )
617  );
618  }
619  else if
620  (
623  )
624  {
625  return tmp<fluxFieldType>
626  (
627  new fluxFieldType
628  (
629  ddtIOobject,
630  fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())
631  *(
632  fvc::interpolate(rDeltaT*rA)*phi.oldTime()
633  - (
634  fvc::interpolate(rDeltaT*rA*U.oldTime())&mesh().Sf()
635  )
636  )
637  )
638  );
639  }
640  else
641  {
643  (
644  "CoEulerDdtScheme<Type>::fvcDdtPhiCorr"
645  ) << "dimensions of phi are not correct"
646  << abort(FatalError);
647 
648  return fluxFieldType::null();
649  }
650  }
651 }
652 
653 
654 template<class Type>
656 (
658 )
659 {
661  (
663  (
664  IOobject
665  (
666  "meshPhi",
667  mesh().time().timeName(),
668  mesh()
669  ),
670  mesh(),
672  )
673  );
674 }
675 
676 
677 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
678 
679 } // End namespace fv
680 
681 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
682 
683 } // End namespace Foam
684 
685 // ************************ vim: set sw=4 sts=4 et: ************************ //