FreeFOAM The Cross-Platform CFD Toolkit
SLTSDdtScheme.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 "SLTSDdtScheme.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace fv
38 {
39 
40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41 
42 template<class Type>
43 void SLTSDdtScheme<Type>::relaxedDiag
44 (
45  scalarField& rD,
46  const surfaceScalarField& phi
47 ) const
48 {
49  const unallocLabelList& owner = mesh().owner();
50  const unallocLabelList& neighbour = mesh().neighbour();
51  scalarField diag(rD.size(), 0.0);
52 
53  forAll(owner, faceI)
54  {
55  if (phi[faceI] > 0.0)
56  {
57  diag[owner[faceI]] += phi[faceI];
58  rD[neighbour[faceI]] += phi[faceI];
59  }
60  else
61  {
62  diag[neighbour[faceI]] -= phi[faceI];
63  rD[owner[faceI]] -= phi[faceI];
64  }
65  }
66 
67  forAll(phi.boundaryField(), patchi)
68  {
69  const fvsPatchScalarField& pphi = phi.boundaryField()[patchi];
70  const unallocLabelList& faceCells = pphi.patch().patch().faceCells();
71 
72  forAll(pphi, patchFacei)
73  {
74  if (pphi[patchFacei] > 0.0)
75  {
76  diag[faceCells[patchFacei]] += pphi[patchFacei];
77  }
78  else
79  {
80  rD[faceCells[patchFacei]] -= pphi[patchFacei];
81  }
82  }
83  }
84 
85  rD += (1.0/alpha_ - 2.0)*diag;
86 }
87 
88 
89 template<class Type>
90 tmp<volScalarField> SLTSDdtScheme<Type>::SLrDeltaT() const
91 {
92  const surfaceScalarField& phi =
93  mesh().objectRegistry::lookupObject<surfaceScalarField>(phiName_);
94 
95  const dimensionedScalar& deltaT = mesh().time().deltaT();
96 
97  tmp<volScalarField> trDeltaT
98  (
99  new volScalarField
100  (
101  IOobject
102  (
103  "rDeltaT",
104  phi.instance(),
105  mesh()
106  ),
107  mesh(),
108  dimensionedScalar("rDeltaT", dimless/dimTime, 0.0),
109  zeroGradientFvPatchScalarField::typeName
110  )
111  );
112 
113  volScalarField& rDeltaT = trDeltaT();
114 
115  relaxedDiag(rDeltaT, phi);
116 
117  if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
118  {
119  rDeltaT.internalField() = max
120  (
121  rDeltaT.internalField()/mesh().V(),
122  scalar(1)/deltaT.value()
123  );
124  }
125  else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
126  {
127  const volScalarField& rho =
128  mesh().objectRegistry::lookupObject<volScalarField>(rhoName_)
129  .oldTime();
130 
131  rDeltaT.internalField() = max
132  (
133  rDeltaT.internalField()/(rho.internalField()*mesh().V()),
134  scalar(1)/deltaT.value()
135  );
136  }
137  else
138  {
139  FatalErrorIn("SLTSDdtScheme<Type>::CofrDeltaT() const")
140  << "Incorrect dimensions of phi: " << phi.dimensions()
141  << abort(FatalError);
142  }
143 
144  rDeltaT.correctBoundaryConditions();
145 
146  return trDeltaT;
147 }
148 
149 
150 template<class Type>
151 tmp<GeometricField<Type, fvPatchField, volMesh> >
153 (
154  const dimensioned<Type>& dt
155 )
156 {
157  volScalarField rDeltaT = SLrDeltaT();
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 = SLrDeltaT();
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 = SLrDeltaT();
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 = SLrDeltaT();
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 = SLrDeltaT()().internalField();
386 
387  Info<< "max/min rDeltaT " << max(rDeltaT) << " " << min(rDeltaT) << endl;
388 
389  fvm.diag() = rDeltaT*mesh().V();
390 
391  if (mesh().moving())
392  {
393  fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
394  }
395  else
396  {
397  fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V();
398  }
399 
400  return tfvm;
401 }
402 
403 
404 template<class Type>
407 (
408  const dimensionedScalar& rho,
410 )
411 {
412  tmp<fvMatrix<Type> > tfvm
413  (
414  new fvMatrix<Type>
415  (
416  vf,
418  )
419  );
420  fvMatrix<Type>& fvm = tfvm();
421 
422  scalarField rDeltaT = SLrDeltaT()().internalField();
423 
424  fvm.diag() = rDeltaT*rho.value()*mesh().V();
425 
426  if (mesh().moving())
427  {
428  fvm.source() = rDeltaT
429  *rho.value()*vf.oldTime().internalField()*mesh().V0();
430  }
431  else
432  {
433  fvm.source() = rDeltaT
434  *rho.value()*vf.oldTime().internalField()*mesh().V();
435  }
436 
437  return tfvm;
438 }
439 
440 
441 template<class Type>
444 (
445  const volScalarField& rho,
447 )
448 {
449  tmp<fvMatrix<Type> > tfvm
450  (
451  new fvMatrix<Type>
452  (
453  vf,
455  )
456  );
457  fvMatrix<Type>& fvm = tfvm();
458 
459  scalarField rDeltaT = SLrDeltaT()().internalField();
460 
461  fvm.diag() = rDeltaT*rho.internalField()*mesh().V();
462 
463  if (mesh().moving())
464  {
465  fvm.source() = rDeltaT
466  *rho.oldTime().internalField()
467  *vf.oldTime().internalField()*mesh().V0();
468  }
469  else
470  {
471  fvm.source() = rDeltaT
472  *rho.oldTime().internalField()
473  *vf.oldTime().internalField()*mesh().V();
474  }
475 
476  return tfvm;
477 }
478 
479 
480 template<class Type>
483 (
484  const volScalarField& rA,
486  const fluxFieldType& phi
487 )
488 {
489  IOobject ddtIOobject
490  (
491  "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
492  mesh().time().timeName(),
493  mesh()
494  );
495 
496  if (mesh().moving())
497  {
498  return tmp<fluxFieldType>
499  (
500  new fluxFieldType
501  (
502  ddtIOobject,
503  mesh(),
504  dimensioned<typename flux<Type>::type>
505  (
506  "0",
507  rA.dimensions()*phi.dimensions()/dimTime,
509  )
510  )
511  );
512  }
513  else
514  {
515  volScalarField rDeltaT = SLrDeltaT();
516 
517  return tmp<fluxFieldType>
518  (
519  new fluxFieldType
520  (
521  ddtIOobject,
522  fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
523  (
524  fvc::interpolate(rDeltaT*rA)*phi.oldTime()
525  - (fvc::interpolate(rDeltaT*rA*U.oldTime()) & mesh().Sf())
526  )
527  )
528  );
529  }
530 }
531 
532 
533 template<class Type>
536 (
537  const volScalarField& rA,
538  const volScalarField& rho,
540  const fluxFieldType& phi
541 )
542 {
543  IOobject ddtIOobject
544  (
545  "ddtPhiCorr("
546  + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
547  mesh().time().timeName(),
548  mesh()
549  );
550 
551  if (mesh().moving())
552  {
553  return tmp<fluxFieldType>
554  (
555  new fluxFieldType
556  (
557  ddtIOobject,
558  mesh(),
559  dimensioned<typename flux<Type>::type>
560  (
561  "0",
562  rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
564  )
565  )
566  );
567  }
568  else
569  {
570  volScalarField rDeltaT = SLrDeltaT();
571 
572  if
573  (
574  U.dimensions() == dimVelocity
575  && phi.dimensions() == dimVelocity*dimArea
576  )
577  {
578  return tmp<fluxFieldType>
579  (
580  new fluxFieldType
581  (
582  ddtIOobject,
583  fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
584  *(
585  fvc::interpolate(rDeltaT*rA*rho.oldTime())*phi.oldTime()
586  - (fvc::interpolate(rDeltaT*rA*rho.oldTime()*U.oldTime())
587  & mesh().Sf())
588  )
589  )
590  );
591  }
592  else if
593  (
594  U.dimensions() == dimVelocity
596  )
597  {
598  return tmp<fluxFieldType>
599  (
600  new fluxFieldType
601  (
602  ddtIOobject,
603  fvcDdtPhiCoeff
604  (
605  U.oldTime(),
606  phi.oldTime()/fvc::interpolate(rho.oldTime())
607  )
608  *(
609  fvc::interpolate(rDeltaT*rA*rho.oldTime())
610  *phi.oldTime()/fvc::interpolate(rho.oldTime())
611  - (
613  (
614  rDeltaT*rA*rho.oldTime()*U.oldTime()
615  ) & mesh().Sf()
616  )
617  )
618  )
619  );
620  }
621  else if
622  (
625  )
626  {
627  return tmp<fluxFieldType>
628  (
629  new fluxFieldType
630  (
631  ddtIOobject,
632  fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())
633  *(
634  fvc::interpolate(rDeltaT*rA)*phi.oldTime()
635  - (
636  fvc::interpolate(rDeltaT*rA*U.oldTime())&mesh().Sf()
637  )
638  )
639  )
640  );
641  }
642  else
643  {
645  (
646  "SLTSDdtScheme<Type>::fvcDdtPhiCorr"
647  ) << "dimensions of phi are not correct"
648  << abort(FatalError);
649 
650  return fluxFieldType::null();
651  }
652  }
653 }
654 
655 
656 template<class Type>
658 (
660 )
661 {
663  (
665  (
666  IOobject
667  (
668  "meshPhi",
669  mesh().time().timeName(),
670  mesh()
671  ),
672  mesh(),
674  )
675  );
676 }
677 
678 
679 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
680 
681 } // End namespace fv
682 
683 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
684 
685 } // End namespace Foam
686 
687 // ************************ vim: set sw=4 sts=4 et: ************************ //