SLTSDdtScheme.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-2018 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  (
101  (
102  "rDeltaT",
103  mesh(),
105  extrapolatedCalculatedFvPatchScalarField::typeName
106  )
107  );
108 
109  volScalarField& rDeltaT = trDeltaT.ref();
110 
111  relaxedDiag(rDeltaT, phi);
112 
113  if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
114  {
115  rDeltaT.primitiveFieldRef() = max
116  (
117  rDeltaT.primitiveField()/mesh().V(),
118  scalar(1)/deltaT.value()
119  );
120  }
121  else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
122  {
123  const volScalarField& rho =
124  mesh().objectRegistry::template lookupObject<volScalarField>
125  (
126  rhoName_
127  ).oldTime();
128 
129  rDeltaT.primitiveFieldRef() = max
130  (
131  rDeltaT.primitiveField()/(rho.primitiveField()*mesh().V()),
132  scalar(1)/deltaT.value()
133  );
134  }
135  else
136  {
138  << "Incorrect dimensions of phi: " << phi.dimensions()
139  << abort(FatalError);
140  }
141 
142  rDeltaT.correctBoundaryConditions();
143 
144  return trDeltaT;
145 }
146 
147 
148 template<class Type>
149 tmp<GeometricField<Type, fvPatchField, volMesh>>
151 (
152  const dimensioned<Type>& dt
153 )
154 {
155  const volScalarField rDeltaT(SLrDeltaT());
156 
157  IOobject ddtIOobject
158  (
159  "ddt("+dt.name()+')',
160  mesh().time().timeName(),
161  mesh()
162  );
163 
164  if (mesh().moving())
165  {
167  (
169  (
170  ddtIOobject,
171  mesh(),
173  (
174  "0",
175  dt.dimensions()/dimTime,
176  Zero
177  )
178  )
179  );
180 
181  tdtdt.ref().primitiveFieldRef() =
182  rDeltaT.primitiveField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
183 
184  return tdtdt;
185  }
186  else
187  {
189  (
191  (
192  ddtIOobject,
193  mesh(),
195  (
196  "0",
197  dt.dimensions()/dimTime,
198  Zero
199  ),
201  )
202  );
203  }
204 }
205 
206 
207 template<class Type>
210 (
212 )
213 {
214  const volScalarField rDeltaT(SLrDeltaT());
215 
216  IOobject ddtIOobject
217  (
218  "ddt("+vf.name()+')',
219  mesh().time().timeName(),
220  mesh()
221  );
222 
223  if (mesh().moving())
224  {
226  (
228  (
229  ddtIOobject,
230  mesh(),
231  rDeltaT.dimensions()*vf.dimensions(),
232  rDeltaT.primitiveField()*
233  (
234  vf.primitiveField()
235  - vf.oldTime().primitiveField()*mesh().V0()/mesh().V()
236  ),
237  rDeltaT.boundaryField()*
238  (
239  vf.boundaryField() - vf.oldTime().boundaryField()
240  )
241  )
242  );
243  }
244  else
245  {
247  (
249  (
250  ddtIOobject,
251  rDeltaT*(vf - vf.oldTime())
252  )
253  );
254  }
255 }
256 
257 
258 template<class Type>
261 (
262  const dimensionedScalar& rho,
264 )
265 {
266  const volScalarField rDeltaT(SLrDeltaT());
267 
268  IOobject ddtIOobject
269  (
270  "ddt("+rho.name()+','+vf.name()+')',
271  mesh().time().timeName(),
272  mesh()
273  );
274 
275  if (mesh().moving())
276  {
278  (
280  (
281  ddtIOobject,
282  mesh(),
283  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
284  rDeltaT.primitiveField()*rho.value()*
285  (
286  vf.primitiveField()
287  - vf.oldTime().primitiveField()*mesh().V0()/mesh().V()
288  ),
289  rDeltaT.boundaryField()*rho.value()*
290  (
291  vf.boundaryField() - vf.oldTime().boundaryField()
292  )
293  )
294  );
295  }
296  else
297  {
299  (
301  (
302  ddtIOobject,
303  rDeltaT*rho*(vf - vf.oldTime())
304  )
305  );
306  }
307 }
308 
309 
310 template<class Type>
313 (
314  const volScalarField& rho,
316 )
317 {
318  const volScalarField rDeltaT(SLrDeltaT());
319 
320  IOobject ddtIOobject
321  (
322  "ddt("+rho.name()+','+vf.name()+')',
323  mesh().time().timeName(),
324  mesh()
325  );
326 
327  if (mesh().moving())
328  {
330  (
332  (
333  ddtIOobject,
334  mesh(),
335  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
336  rDeltaT.primitiveField()*
337  (
338  rho.primitiveField()*vf.primitiveField()
339  - rho.oldTime().primitiveField()
340  *vf.oldTime().primitiveField()*mesh().V0()/mesh().V()
341  ),
342  rDeltaT.boundaryField()*
343  (
344  rho.boundaryField()*vf.boundaryField()
345  - rho.oldTime().boundaryField()
346  *vf.oldTime().boundaryField()
347  )
348  )
349  );
350  }
351  else
352  {
354  (
356  (
357  ddtIOobject,
358  rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
359  )
360  );
361  }
362 }
363 
364 
365 template<class Type>
368 (
369  const volScalarField& alpha,
370  const volScalarField& rho,
372 )
373 {
374  const volScalarField rDeltaT(SLrDeltaT());
375 
376  IOobject ddtIOobject
377  (
378  "ddt("+alpha.name()+','+rho.name()+','+vf.name()+')',
379  mesh().time().timeName(),
380  mesh()
381  );
382 
383  if (mesh().moving())
384  {
386  (
388  (
389  ddtIOobject,
390  mesh(),
391  rDeltaT.dimensions()
392  *alpha.dimensions()*rho.dimensions()*vf.dimensions(),
393  rDeltaT.primitiveField()*
394  (
395  alpha.primitiveField()
396  *rho.primitiveField()
397  *vf.primitiveField()
398 
399  - alpha.oldTime().primitiveField()
400  *rho.oldTime().primitiveField()
401  *vf.oldTime().primitiveField()*mesh().Vsc0()/mesh().Vsc()
402  ),
403  rDeltaT.boundaryField()*
404  (
405  alpha.boundaryField()
406  *rho.boundaryField()
407  *vf.boundaryField()
408 
409  - alpha.oldTime().boundaryField()
410  *rho.oldTime().boundaryField()
411  *vf.oldTime().boundaryField()
412  )
413  )
414  );
415  }
416  else
417  {
419  (
421  (
422  ddtIOobject,
423  rDeltaT
424  *(
425  alpha*rho*vf
426  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
427  )
428  )
429  );
430  }
431 }
432 
433 
434 template<class Type>
437 (
439 )
440 {
441  tmp<fvMatrix<Type>> tfvm
442  (
443  new fvMatrix<Type>
444  (
445  vf,
447  )
448  );
449 
450  fvMatrix<Type>& fvm = tfvm.ref();
451 
452  scalarField rDeltaT(SLrDeltaT()().primitiveField());
453 
454  fvm.diag() = rDeltaT*mesh().V();
455 
456  if (mesh().moving())
457  {
458  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().V0();
459  }
460  else
461  {
462  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().V();
463  }
464 
465  return tfvm;
466 }
467 
468 
469 template<class Type>
472 (
473  const dimensionedScalar& rho,
475 )
476 {
477  tmp<fvMatrix<Type>> tfvm
478  (
479  new fvMatrix<Type>
480  (
481  vf,
483  )
484  );
485  fvMatrix<Type>& fvm = tfvm.ref();
486 
487  scalarField rDeltaT(SLrDeltaT()().primitiveField());
488 
489  fvm.diag() = rDeltaT*rho.value()*mesh().V();
490 
491  if (mesh().moving())
492  {
493  fvm.source() = rDeltaT
494  *rho.value()*vf.oldTime().primitiveField()*mesh().V0();
495  }
496  else
497  {
498  fvm.source() = rDeltaT
499  *rho.value()*vf.oldTime().primitiveField()*mesh().V();
500  }
501 
502  return tfvm;
503 }
504 
505 
506 template<class Type>
509 (
510  const volScalarField& rho,
512 )
513 {
514  tmp<fvMatrix<Type>> tfvm
515  (
516  new fvMatrix<Type>
517  (
518  vf,
520  )
521  );
522  fvMatrix<Type>& fvm = tfvm.ref();
523 
524  scalarField rDeltaT(SLrDeltaT()().primitiveField());
525 
526  fvm.diag() = rDeltaT*rho.primitiveField()*mesh().V();
527 
528  if (mesh().moving())
529  {
530  fvm.source() = rDeltaT
531  *rho.oldTime().primitiveField()
532  *vf.oldTime().primitiveField()*mesh().V0();
533  }
534  else
535  {
536  fvm.source() = rDeltaT
537  *rho.oldTime().primitiveField()
538  *vf.oldTime().primitiveField()*mesh().V();
539  }
540 
541  return tfvm;
542 }
543 
544 
545 template<class Type>
548 (
549  const volScalarField& alpha,
550  const volScalarField& rho,
552 )
553 {
554  tmp<fvMatrix<Type>> tfvm
555  (
556  new fvMatrix<Type>
557  (
558  vf,
559  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
560  )
561  );
562  fvMatrix<Type>& fvm = tfvm.ref();
563 
564  scalarField rDeltaT(SLrDeltaT()().primitiveField());
565 
566  fvm.diag() =
567  rDeltaT*alpha.primitiveField()*rho.primitiveField()*mesh().Vsc();
568 
569  if (mesh().moving())
570  {
571  fvm.source() = rDeltaT
572  *alpha.oldTime().primitiveField()
573  *rho.oldTime().primitiveField()
574  *vf.oldTime().primitiveField()*mesh().Vsc0();
575  }
576  else
577  {
578  fvm.source() = rDeltaT
579  *alpha.oldTime().primitiveField()
580  *rho.oldTime().primitiveField()
581  *vf.oldTime().primitiveField()*mesh().Vsc();
582  }
583 
584  return tfvm;
585 }
586 
587 
588 template<class Type>
591 (
594 )
595 {
596  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
597 
598  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
599  fluxFieldType phiCorr
600  (
601  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
602  );
603 
604  return tmp<fluxFieldType>
605  (
606  new fluxFieldType
607  (
608  IOobject
609  (
610  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
611  mesh().time().timeName(),
612  mesh()
613  ),
614  this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
615  *rDeltaT*phiCorr
616  )
617  );
618 }
619 
620 
621 template<class Type>
624 (
626  const fluxFieldType& phi
627 )
628 {
629  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
630 
631  fluxFieldType phiCorr
632  (
633  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
634  );
635 
636  return tmp<fluxFieldType>
637  (
638  new fluxFieldType
639  (
640  IOobject
641  (
642  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
643  mesh().time().timeName(),
644  mesh()
645  ),
646  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
647  *rDeltaT*phiCorr
648  )
649  );
650 }
651 
652 
653 template<class Type>
656 (
657  const volScalarField& rho,
660 )
661 {
662  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
663 
664  if
665  (
666  U.dimensions() == dimVelocity
668  )
669  {
671  (
672  rho.oldTime()*U.oldTime()
673  );
674 
675  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
676  fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
677 
678  return tmp<fluxFieldType>
679  (
680  new fluxFieldType
681  (
682  IOobject
683  (
684  "ddtCorr("
685  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
686  mesh().time().timeName(),
687  mesh()
688  ),
689  this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime())
690  *rDeltaT*phiCorr
691  )
692  );
693  }
694  else if
695  (
698  )
699  {
700  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
701  fluxFieldType phiCorr
702  (
703  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
704  );
705 
706  return tmp<fluxFieldType>
707  (
708  new fluxFieldType
709  (
710  IOobject
711  (
712  "ddtCorr("
713  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
714  mesh().time().timeName(),
715  mesh()
716  ),
717  this->fvcDdtPhiCoeff
718  (
719  U.oldTime(),
720  phiUf0,
721  phiCorr,
722  rho.oldTime()
723  )*rDeltaT*phiCorr
724  )
725  );
726  }
727  else
728  {
730  << "dimensions of Uf are not correct"
731  << abort(FatalError);
732 
733  return fluxFieldType::null();
734  }
735 }
736 
737 
738 template<class Type>
741 (
742  const volScalarField& rho,
744  const fluxFieldType& phi
745 )
746 {
747  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
748 
749  if
750  (
751  U.dimensions() == dimVelocity
752  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
753  )
754  {
756  (
757  rho.oldTime()*U.oldTime()
758  );
759 
760  fluxFieldType phiCorr
761  (
762  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
763  );
764 
765  return tmp<fluxFieldType>
766  (
767  new fluxFieldType
768  (
769  IOobject
770  (
771  "ddtCorr("
772  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
773  mesh().time().timeName(),
774  mesh()
775  ),
776  this->fvcDdtPhiCoeff
777  (
778  rhoU0,
779  phi.oldTime(),
780  phiCorr,
781  rho.oldTime()
782  )*rDeltaT*phiCorr
783  )
784  );
785  }
786  else if
787  (
788  U.dimensions() == rho.dimensions()*dimVelocity
789  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
790  )
791  {
792  fluxFieldType phiCorr
793  (
794  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
795  );
796 
797  return tmp<fluxFieldType>
798  (
799  new fluxFieldType
800  (
801  IOobject
802  (
803  "ddtCorr("
804  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
805  mesh().time().timeName(),
806  mesh()
807  ),
808  this->fvcDdtPhiCoeff
809  (
810  U.oldTime(),
811  phi.oldTime(),
812  phiCorr,
813  rho.oldTime()
814  )*rDeltaT*phiCorr
815  )
816  );
817  }
818  else
819  {
821  << "dimensions of phi are not correct"
822  << abort(FatalError);
823 
824  return fluxFieldType::null();
825  }
826 }
827 
828 
829 template<class Type>
831 (
833 )
834 {
836  (
837  "meshPhi",
838  mesh(),
840  );
841 }
842 
843 
844 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
845 
846 } // End namespace fv
847 
848 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
849 
850 } // End namespace Foam
851 
852 // ************************************************************************* //
fvsPatchField< scalar > fvsPatchScalarField
rho oldTime()
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
const word & name() const
Return name.
Definition: IOobject.H:295
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:65
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Mesh &, const dimensionSet &, const word &patchFieldType=fvPatchField< scalar >::calculatedType())
Return a temporary field constructed from name, mesh, dimensionSet.
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.
virtual 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
static const zero Zero
Definition: zero.H:97
virtual 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:292
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.
virtual 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.
mesh Sf()
A class for managing temporary objects.
Definition: PtrList.H:53
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
scalarField & diag()
Definition: lduMatrix.C:186
virtual 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
virtual tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Namespace for OpenFOAM.
const dimensionSet dimVelocity