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) 2011-2016 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 "SLTSDdtScheme.H"
27 #include "surfaceInterpolate.H"
28 #include "fvMatrices.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 labelUList& owner = mesh().owner();
50  const labelUList& 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 labelUList& 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::template
94  lookupObject<surfaceScalarField>(phiName_);
95 
96  const dimensionedScalar& deltaT = mesh().time().deltaT();
97 
98  tmp<volScalarField> trDeltaT
99  (
100  new volScalarField
101  (
102  IOobject
103  (
104  "rDeltaT",
105  phi.instance(),
106  mesh()
107  ),
108  mesh(),
109  dimensionedScalar("rDeltaT", dimless/dimTime, 0.0),
110  extrapolatedCalculatedFvPatchScalarField::typeName
111  )
112  );
113 
114  volScalarField& rDeltaT = trDeltaT.ref();
115 
116  relaxedDiag(rDeltaT, phi);
117 
118  if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
119  {
120  rDeltaT.primitiveFieldRef() = max
121  (
122  rDeltaT.primitiveField()/mesh().V(),
123  scalar(1)/deltaT.value()
124  );
125  }
126  else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
127  {
128  const volScalarField& rho =
129  mesh().objectRegistry::template lookupObject<volScalarField>
130  (
131  rhoName_
132  ).oldTime();
133 
134  rDeltaT.primitiveFieldRef() = max
135  (
136  rDeltaT.primitiveField()/(rho.primitiveField()*mesh().V()),
137  scalar(1)/deltaT.value()
138  );
139  }
140  else
141  {
143  << "Incorrect dimensions of phi: " << phi.dimensions()
144  << abort(FatalError);
145  }
146 
147  rDeltaT.correctBoundaryConditions();
148 
149  return trDeltaT;
150 }
151 
152 
153 template<class Type>
154 tmp<GeometricField<Type, fvPatchField, volMesh>>
156 (
157  const dimensioned<Type>& dt
158 )
159 {
160  const volScalarField rDeltaT(SLrDeltaT());
161 
162  IOobject ddtIOobject
163  (
164  "ddt("+dt.name()+')',
165  mesh().time().timeName(),
166  mesh()
167  );
168 
169  if (mesh().moving())
170  {
172  (
174  (
175  ddtIOobject,
176  mesh(),
178  (
179  "0",
180  dt.dimensions()/dimTime,
181  Zero
182  )
183  )
184  );
185 
186  tdtdt.ref().primitiveFieldRef() =
187  rDeltaT.primitiveField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
188 
189  return tdtdt;
190  }
191  else
192  {
194  (
196  (
197  ddtIOobject,
198  mesh(),
200  (
201  "0",
202  dt.dimensions()/dimTime,
203  Zero
204  ),
206  )
207  );
208  }
209 }
210 
211 
212 template<class Type>
215 (
217 )
218 {
219  const volScalarField rDeltaT(SLrDeltaT());
220 
221  IOobject ddtIOobject
222  (
223  "ddt("+vf.name()+')',
224  mesh().time().timeName(),
225  mesh()
226  );
227 
228  if (mesh().moving())
229  {
231  (
233  (
234  ddtIOobject,
235  mesh(),
236  rDeltaT.dimensions()*vf.dimensions(),
237  rDeltaT.primitiveField()*
238  (
239  vf.primitiveField()
240  - vf.oldTime().primitiveField()*mesh().V0()/mesh().V()
241  ),
242  rDeltaT.boundaryField()*
243  (
244  vf.boundaryField() - vf.oldTime().boundaryField()
245  )
246  )
247  );
248  }
249  else
250  {
252  (
254  (
255  ddtIOobject,
256  rDeltaT*(vf - vf.oldTime())
257  )
258  );
259  }
260 }
261 
262 
263 template<class Type>
266 (
267  const dimensionedScalar& rho,
269 )
270 {
271  const volScalarField rDeltaT(SLrDeltaT());
272 
273  IOobject ddtIOobject
274  (
275  "ddt("+rho.name()+','+vf.name()+')',
276  mesh().time().timeName(),
277  mesh()
278  );
279 
280  if (mesh().moving())
281  {
283  (
285  (
286  ddtIOobject,
287  mesh(),
288  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
289  rDeltaT.primitiveField()*rho.value()*
290  (
291  vf.primitiveField()
292  - vf.oldTime().primitiveField()*mesh().V0()/mesh().V()
293  ),
294  rDeltaT.boundaryField()*rho.value()*
295  (
296  vf.boundaryField() - vf.oldTime().boundaryField()
297  )
298  )
299  );
300  }
301  else
302  {
304  (
306  (
307  ddtIOobject,
308  rDeltaT*rho*(vf - vf.oldTime())
309  )
310  );
311  }
312 }
313 
314 
315 template<class Type>
318 (
319  const volScalarField& rho,
321 )
322 {
323  const volScalarField rDeltaT(SLrDeltaT());
324 
325  IOobject ddtIOobject
326  (
327  "ddt("+rho.name()+','+vf.name()+')',
328  mesh().time().timeName(),
329  mesh()
330  );
331 
332  if (mesh().moving())
333  {
335  (
337  (
338  ddtIOobject,
339  mesh(),
340  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
341  rDeltaT.primitiveField()*
342  (
343  rho.primitiveField()*vf.primitiveField()
344  - rho.oldTime().primitiveField()
345  *vf.oldTime().primitiveField()*mesh().V0()/mesh().V()
346  ),
347  rDeltaT.boundaryField()*
348  (
349  rho.boundaryField()*vf.boundaryField()
350  - rho.oldTime().boundaryField()
351  *vf.oldTime().boundaryField()
352  )
353  )
354  );
355  }
356  else
357  {
359  (
361  (
362  ddtIOobject,
363  rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
364  )
365  );
366  }
367 }
368 
369 
370 template<class Type>
373 (
374  const volScalarField& alpha,
375  const volScalarField& rho,
377 )
378 {
379  const volScalarField rDeltaT(SLrDeltaT());
380 
381  IOobject ddtIOobject
382  (
383  "ddt("+alpha.name()+','+rho.name()+','+vf.name()+')',
384  mesh().time().timeName(),
385  mesh()
386  );
387 
388  if (mesh().moving())
389  {
391  (
393  (
394  ddtIOobject,
395  mesh(),
396  rDeltaT.dimensions()
397  *alpha.dimensions()*rho.dimensions()*vf.dimensions(),
398  rDeltaT.primitiveField()*
399  (
400  alpha.primitiveField()
401  *rho.primitiveField()
402  *vf.primitiveField()
403 
404  - alpha.oldTime().primitiveField()
405  *rho.oldTime().primitiveField()
406  *vf.oldTime().primitiveField()*mesh().Vsc0()/mesh().Vsc()
407  ),
408  rDeltaT.boundaryField()*
409  (
410  alpha.boundaryField()
411  *rho.boundaryField()
412  *vf.boundaryField()
413 
414  - alpha.oldTime().boundaryField()
415  *rho.oldTime().boundaryField()
416  *vf.oldTime().boundaryField()
417  )
418  )
419  );
420  }
421  else
422  {
424  (
426  (
427  ddtIOobject,
428  rDeltaT
429  *(
430  alpha*rho*vf
431  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
432  )
433  )
434  );
435  }
436 }
437 
438 
439 template<class Type>
442 (
444 )
445 {
446  tmp<fvMatrix<Type>> tfvm
447  (
448  new fvMatrix<Type>
449  (
450  vf,
452  )
453  );
454 
455  fvMatrix<Type>& fvm = tfvm.ref();
456 
457  scalarField rDeltaT(SLrDeltaT()().primitiveField());
458 
459  fvm.diag() = rDeltaT*mesh().V();
460 
461  if (mesh().moving())
462  {
463  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().V0();
464  }
465  else
466  {
467  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().V();
468  }
469 
470  return tfvm;
471 }
472 
473 
474 template<class Type>
477 (
478  const dimensionedScalar& rho,
480 )
481 {
482  tmp<fvMatrix<Type>> tfvm
483  (
484  new fvMatrix<Type>
485  (
486  vf,
488  )
489  );
490  fvMatrix<Type>& fvm = tfvm.ref();
491 
492  scalarField rDeltaT(SLrDeltaT()().primitiveField());
493 
494  fvm.diag() = rDeltaT*rho.value()*mesh().V();
495 
496  if (mesh().moving())
497  {
498  fvm.source() = rDeltaT
499  *rho.value()*vf.oldTime().primitiveField()*mesh().V0();
500  }
501  else
502  {
503  fvm.source() = rDeltaT
504  *rho.value()*vf.oldTime().primitiveField()*mesh().V();
505  }
506 
507  return tfvm;
508 }
509 
510 
511 template<class Type>
514 (
515  const volScalarField& rho,
517 )
518 {
519  tmp<fvMatrix<Type>> tfvm
520  (
521  new fvMatrix<Type>
522  (
523  vf,
525  )
526  );
527  fvMatrix<Type>& fvm = tfvm.ref();
528 
529  scalarField rDeltaT(SLrDeltaT()().primitiveField());
530 
531  fvm.diag() = rDeltaT*rho.primitiveField()*mesh().V();
532 
533  if (mesh().moving())
534  {
535  fvm.source() = rDeltaT
536  *rho.oldTime().primitiveField()
537  *vf.oldTime().primitiveField()*mesh().V0();
538  }
539  else
540  {
541  fvm.source() = rDeltaT
542  *rho.oldTime().primitiveField()
543  *vf.oldTime().primitiveField()*mesh().V();
544  }
545 
546  return tfvm;
547 }
548 
549 
550 template<class Type>
553 (
554  const volScalarField& alpha,
555  const volScalarField& rho,
557 )
558 {
559  tmp<fvMatrix<Type>> tfvm
560  (
561  new fvMatrix<Type>
562  (
563  vf,
564  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
565  )
566  );
567  fvMatrix<Type>& fvm = tfvm.ref();
568 
569  scalarField rDeltaT(SLrDeltaT()().primitiveField());
570 
571  fvm.diag() =
572  rDeltaT*alpha.primitiveField()*rho.primitiveField()*mesh().Vsc();
573 
574  if (mesh().moving())
575  {
576  fvm.source() = rDeltaT
577  *alpha.oldTime().primitiveField()
578  *rho.oldTime().primitiveField()
579  *vf.oldTime().primitiveField()*mesh().Vsc0();
580  }
581  else
582  {
583  fvm.source() = rDeltaT
584  *alpha.oldTime().primitiveField()
585  *rho.oldTime().primitiveField()
586  *vf.oldTime().primitiveField()*mesh().Vsc();
587  }
588 
589  return tfvm;
590 }
591 
592 
593 template<class Type>
596 (
599 )
600 {
601  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
602 
603  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
604  fluxFieldType phiCorr
605  (
606  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
607  );
608 
609  return tmp<fluxFieldType>
610  (
611  new fluxFieldType
612  (
613  IOobject
614  (
615  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
616  mesh().time().timeName(),
617  mesh()
618  ),
619  this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
620  *rDeltaT*phiCorr
621  )
622  );
623 }
624 
625 
626 template<class Type>
629 (
631  const fluxFieldType& phi
632 )
633 {
634  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
635 
636  fluxFieldType phiCorr
637  (
638  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
639  );
640 
641  return tmp<fluxFieldType>
642  (
643  new fluxFieldType
644  (
645  IOobject
646  (
647  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
648  mesh().time().timeName(),
649  mesh()
650  ),
651  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
652  *rDeltaT*phiCorr
653  )
654  );
655 }
656 
657 
658 template<class Type>
661 (
662  const volScalarField& rho,
665 )
666 {
667  if
668  (
669  U.dimensions() == dimVelocity
671  )
672  {
673  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
674 
676  (
677  rho.oldTime()*U.oldTime()
678  );
679 
680  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
681  fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
682 
683  return tmp<fluxFieldType>
684  (
685  new fluxFieldType
686  (
687  IOobject
688  (
689  "ddtCorr("
690  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
691  mesh().time().timeName(),
692  mesh()
693  ),
694  this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr)*rDeltaT*phiCorr
695  )
696  );
697  }
698  else if
699  (
702  )
703  {
704  return fvcDdtUfCorr(U, Uf);
705  }
706  else
707  {
709  << "dimensions of Uf are not correct"
710  << abort(FatalError);
711 
712  return fluxFieldType::null();
713  }
714 }
715 
716 
717 template<class Type>
720 (
721  const volScalarField& rho,
723  const fluxFieldType& phi
724 )
725 {
726  if
727  (
728  U.dimensions() == dimVelocity
729  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
730  )
731  {
732  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
733 
735  (
736  rho.oldTime()*U.oldTime()
737  );
738 
739  fluxFieldType phiCorr
740  (
741  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
742  );
743 
744  return tmp<fluxFieldType>
745  (
746  new fluxFieldType
747  (
748  IOobject
749  (
750  "ddtCorr("
751  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
752  mesh().time().timeName(),
753  mesh()
754  ),
755  this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr)
756  *rDeltaT*phiCorr
757  )
758  );
759  }
760  else if
761  (
762  U.dimensions() == rho.dimensions()*dimVelocity
763  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
764  )
765  {
766  return fvcDdtPhiCorr(U, phi);
767  }
768  else
769  {
771  << "dimensions of phi are not correct"
772  << abort(FatalError);
773 
774  return fluxFieldType::null();
775  }
776 }
777 
778 
779 template<class Type>
781 (
783 )
784 {
786  (
788  (
789  IOobject
790  (
791  "meshPhi",
792  mesh().time().timeName(),
793  mesh(),
796  false
797  ),
798  mesh(),
800  )
801  );
802 }
803 
804 
805 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
806 
807 } // End namespace fv
808 
809 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
810 
811 } // End namespace Foam
812 
813 // ************************************************************************* //
fvsPatchField< scalar > fvsPatchScalarField
rho oldTime()
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
const word & name() const
Return name.
Definition: IOobject.H:291
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const Boundary & boundaryField() const
Return const-reference to the boundary field.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
static tmp< GeometricField< typename innerProduct< vector, Type >::type, fvsPatchField, surfaceMesh > > dotInterpolate(const surfaceVectorField &Sf, const GeometricField< Type, fvPatchField, volMesh > &tvf)
Interpolate field onto faces.
const dimensionSet dimVol(dimVolume)
Definition: dimensionSets.H:59
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Generic dimensioned Type class.
UList< label > labelUList
Definition: UList.H:64
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
const dimensionSet & dimensions() const
Return dimensions.
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
dynamicFvMesh & mesh
labelList fv(nPoints)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
This boundary condition is not designed to be evaluated; it is assmued that the value is assigned via...
const Type & value() const
Return const reference to value.
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
word timeName
Definition: getTimeIndex.H:3
static const zero Zero
Definition: zero.H:91
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Field< Type > & source()
Definition: fvMatrix.H:294
const word & name() const
Return const reference to name.
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.
const dimensionSet dimDensity
label patchi
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
A special matrix type and solver, designed for finite volume solutions of scalar equations.
const dimensionSet & dimensions() const
Return const reference to dimensions.
A class for managing temporary objects.
Definition: PtrList.H:53
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
scalarField & diag()
Definition: lduMatrix.C:186
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
tmp< volScalarField > trDeltaT
Definition: createRDeltaT.H:3
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Namespace for OpenFOAM.
const dimensionSet dimVelocity