CoEulerDdtScheme.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 "CoEulerDdtScheme.H"
27 #include "surfaceInterpolate.H"
28 #include "fvcDiv.H"
29 #include "fvMatrices.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  const 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  extrapolatedCalculatedFvPatchScalarField::typeName
61  )
62  );
63 
64  volScalarField& corDeltaT = tcorDeltaT.ref();
65 
66  const labelUList& owner = mesh().owner();
67  const labelUList& 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  const surfaceScalarField::Boundary& cofrDeltaTbf =
79  cofrDeltaT.boundaryField();
80 
81  forAll(cofrDeltaTbf, patchi)
82  {
83  const fvsPatchScalarField& pcofrDeltaT = cofrDeltaTbf[patchi];
84  const fvPatch& p = pcofrDeltaT.patch();
85  const labelUList& faceCells = p.patch().faceCells();
86 
87  forAll(pcofrDeltaT, patchFacei)
88  {
89  corDeltaT[faceCells[patchFacei]] = max
90  (
91  corDeltaT[faceCells[patchFacei]],
92  pcofrDeltaT[patchFacei]
93  );
94  }
95  }
96 
97  corDeltaT.correctBoundaryConditions();
98 
99  return tcorDeltaT;
100 }
101 
102 
103 template<class Type>
104 tmp<surfaceScalarField> CoEulerDdtScheme<Type>::CofrDeltaT() const
105 {
106  const dimensionedScalar& deltaT = mesh().time().deltaT();
107 
108  const surfaceScalarField& phi =
109  static_cast<const objectRegistry&>(mesh())
110  .lookupObject<surfaceScalarField>(phiName_);
111 
112  if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
113  {
115  (
117  *(mag(phi)/mesh().magSf())
118  *deltaT
119  );
120 
121  return max(Co/maxCo_, scalar(1))/deltaT;
122  }
123  else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
124  {
125  const volScalarField& rho =
126  static_cast<const objectRegistry&>(mesh())
127  .lookupObject<volScalarField>(rhoName_).oldTime();
128 
130  (
132  *(mag(phi)/(fvc::interpolate(rho)*mesh().magSf()))
133  *deltaT
134  );
135 
136  return max(Co/maxCo_, scalar(1))/deltaT;
137  }
138  else
139  {
141  << "Incorrect dimensions of phi: " << phi.dimensions()
142  << abort(FatalError);
143 
144  return tmp<surfaceScalarField>(nullptr);
145  }
146 }
147 
148 
149 template<class Type>
150 tmp<GeometricField<Type, fvPatchField, volMesh>>
152 (
153  const dimensioned<Type>& dt
154 )
155 {
156  const volScalarField rDeltaT(CorDeltaT());
157 
158  IOobject ddtIOobject
159  (
160  "ddt("+dt.name()+')',
161  mesh().time().timeName(),
162  mesh()
163  );
164 
165  if (mesh().moving())
166  {
168  (
170  (
171  ddtIOobject,
172  mesh(),
174  (
175  "0",
176  dt.dimensions()/dimTime,
177  Zero
178  )
179  )
180  );
181 
182  tdtdt.ref().primitiveFieldRef() =
183  rDeltaT.primitiveField()*dt.value()
184  *(1.0 - mesh().Vsc0()/mesh().Vsc());
185 
186  return tdtdt;
187  }
188  else
189  {
191  (
193  (
194  ddtIOobject,
195  mesh(),
197  (
198  "0",
199  dt.dimensions()/dimTime,
200  Zero
201  ),
203  )
204  );
205  }
206 }
207 
208 
209 template<class Type>
212 (
214 )
215 {
216  const 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.primitiveField()*
235  (
236  vf.primitiveField()
237  - vf.oldTime().primitiveField()*mesh().Vsc0()/mesh().Vsc()
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  const 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.primitiveField()*rho.value()*
287  (
288  vf.primitiveField()
289  - vf.oldTime().primitiveField()*mesh().Vsc0()/mesh().Vsc()
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  const 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.primitiveField()*
339  (
340  rho.primitiveField()*vf.primitiveField()
341  - rho.oldTime().primitiveField()
342  *vf.oldTime().primitiveField()*mesh().Vsc0()/mesh().Vsc()
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 (
371  const volScalarField& alpha,
372  const volScalarField& rho,
374 )
375 {
376  const volScalarField rDeltaT(CorDeltaT());
377 
378  IOobject ddtIOobject
379  (
380  "ddt("+alpha.name()+','+rho.name()+','+vf.name()+')',
381  mesh().time().timeName(),
382  mesh()
383  );
384 
385  if (mesh().moving())
386  {
388  (
390  (
391  ddtIOobject,
392  mesh(),
393  rDeltaT.dimensions()
394  *alpha.dimensions()*rho.dimensions()*vf.dimensions(),
395  rDeltaT.primitiveField()*
396  (
397  alpha.primitiveField()
398  *rho.primitiveField()
399  *vf.primitiveField()
400 
401  - alpha.oldTime().primitiveField()
402  *rho.oldTime().primitiveField()
403  *vf.oldTime().primitiveField()*mesh().Vsc0()/mesh().Vsc()
404  ),
405  rDeltaT.boundaryField()*
406  (
407  alpha.boundaryField()
408  *rho.boundaryField()
409  *vf.boundaryField()
410 
411  - alpha.oldTime().boundaryField()
412  *rho.oldTime().boundaryField()
413  *vf.oldTime().boundaryField()
414  )
415  )
416  );
417  }
418  else
419  {
421  (
423  (
424  ddtIOobject,
425  rDeltaT
426  *(
427  alpha*rho*vf
428  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
429  )
430  )
431  );
432  }
433 }
434 
435 
436 template<class Type>
439 (
441 )
442 {
443  tmp<fvMatrix<Type>> tfvm
444  (
445  new fvMatrix<Type>
446  (
447  vf,
449  )
450  );
451 
452  fvMatrix<Type>& fvm = tfvm.ref();
453 
454  scalarField rDeltaT(CorDeltaT()().primitiveField());
455 
456  fvm.diag() = rDeltaT*mesh().Vsc();
457 
458  if (mesh().moving())
459  {
460  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc0();
461  }
462  else
463  {
464  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc();
465  }
466 
467  return tfvm;
468 }
469 
470 
471 template<class Type>
474 (
475  const dimensionedScalar& rho,
477 )
478 {
479  tmp<fvMatrix<Type>> tfvm
480  (
481  new fvMatrix<Type>
482  (
483  vf,
485  )
486  );
487  fvMatrix<Type>& fvm = tfvm.ref();
488 
489  scalarField rDeltaT(CorDeltaT()().primitiveField());
490 
491  fvm.diag() = rDeltaT*rho.value()*mesh().Vsc();
492 
493  if (mesh().moving())
494  {
495  fvm.source() = rDeltaT
496  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc0();
497  }
498  else
499  {
500  fvm.source() = rDeltaT
501  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc();
502  }
503 
504  return tfvm;
505 }
506 
507 
508 template<class Type>
511 (
512  const volScalarField& rho,
514 )
515 {
516  tmp<fvMatrix<Type>> tfvm
517  (
518  new fvMatrix<Type>
519  (
520  vf,
522  )
523  );
524  fvMatrix<Type>& fvm = tfvm.ref();
525 
526  scalarField rDeltaT(CorDeltaT()().primitiveField());
527 
528  fvm.diag() = rDeltaT*rho.primitiveField()*mesh().Vsc();
529 
530  if (mesh().moving())
531  {
532  fvm.source() = rDeltaT
533  *rho.oldTime().primitiveField()
534  *vf.oldTime().primitiveField()*mesh().Vsc0();
535  }
536  else
537  {
538  fvm.source() = rDeltaT
539  *rho.oldTime().primitiveField()
540  *vf.oldTime().primitiveField()*mesh().Vsc();
541  }
542 
543  return tfvm;
544 }
545 
546 
547 template<class Type>
550 (
551  const volScalarField& alpha,
552  const volScalarField& rho,
554 )
555 {
556  tmp<fvMatrix<Type>> tfvm
557  (
558  new fvMatrix<Type>
559  (
560  vf,
561  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
562  )
563  );
564  fvMatrix<Type>& fvm = tfvm.ref();
565 
566  scalarField rDeltaT(CorDeltaT()().primitiveField());
567 
568  fvm.diag() =
569  rDeltaT*alpha.primitiveField()*rho.primitiveField()*mesh().Vsc();
570 
571  if (mesh().moving())
572  {
573  fvm.source() = rDeltaT
574  *alpha.oldTime().primitiveField()
575  *rho.oldTime().primitiveField()
576  *vf.oldTime().primitiveField()*mesh().Vsc0();
577  }
578  else
579  {
580  fvm.source() = rDeltaT
581  *alpha.oldTime().primitiveField()
582  *rho.oldTime().primitiveField()
583  *vf.oldTime().primitiveField()*mesh().Vsc();
584  }
585 
586  return tfvm;
587 }
588 
589 
590 template<class Type>
593 (
596 )
597 {
598  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
599 
600  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
601  fluxFieldType phiCorr
602  (
603  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
604  );
605 
606  return tmp<fluxFieldType>
607  (
608  new fluxFieldType
609  (
610  IOobject
611  (
612  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
613  mesh().time().timeName(),
614  mesh()
615  ),
616  this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
617  *rDeltaT*phiCorr
618  )
619  );
620 }
621 
622 
623 template<class Type>
626 (
628  const fluxFieldType& phi
629 )
630 {
631  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
632 
633  fluxFieldType phiCorr
634  (
635  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
636  );
637 
638  return tmp<fluxFieldType>
639  (
640  new fluxFieldType
641  (
642  IOobject
643  (
644  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
645  mesh().time().timeName(),
646  mesh()
647  ),
648  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
649  *rDeltaT*phiCorr
650  )
651  );
652 }
653 
654 
655 template<class Type>
658 (
659  const volScalarField& rho,
662 )
663 {
664  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
665 
666  if
667  (
668  U.dimensions() == dimVelocity
670  )
671  {
673  (
674  rho.oldTime()*U.oldTime()
675  );
676 
677  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
678  fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
679 
680  return tmp<fluxFieldType>
681  (
682  new fluxFieldType
683  (
684  IOobject
685  (
686  "ddtCorr("
687  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
688  mesh().time().timeName(),
689  mesh()
690  ),
691  this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime())
692  *rDeltaT*phiCorr
693  )
694  );
695  }
696  else if
697  (
700  )
701  {
702  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
703  fluxFieldType phiCorr
704  (
705  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
706  );
707 
708  return tmp<fluxFieldType>
709  (
710  new fluxFieldType
711  (
712  IOobject
713  (
714  "ddtCorr("
715  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
716  mesh().time().timeName(),
717  mesh()
718  ),
719  this->fvcDdtPhiCoeff
720  (
721  U.oldTime(),
722  phiUf0,
723  phiCorr,
724  rho.oldTime()
725  )*rDeltaT*phiCorr
726  )
727  );
728  }
729  else
730  {
732  << "dimensions of Uf are not correct"
733  << abort(FatalError);
734 
735  return fluxFieldType::null();
736  }
737 }
738 
739 
740 template<class Type>
743 (
744  const volScalarField& rho,
746  const fluxFieldType& phi
747 )
748 {
749  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
750 
751  if
752  (
753  U.dimensions() == dimVelocity
754  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
755  )
756  {
758  (
759  rho.oldTime()*U.oldTime()
760  );
761 
762  fluxFieldType phiCorr
763  (
764  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
765  );
766 
767  return tmp<fluxFieldType>
768  (
769  new fluxFieldType
770  (
771  IOobject
772  (
773  "ddtCorr("
774  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
775  mesh().time().timeName(),
776  mesh()
777  ),
778  this->fvcDdtPhiCoeff
779  (
780  rhoU0,
781  phi.oldTime(),
782  phiCorr,
783  rho.oldTime()
784  )*rDeltaT*phiCorr
785  )
786  );
787  }
788  else if
789  (
790  U.dimensions() == rho.dimensions()*dimVelocity
791  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
792  )
793  {
794  fluxFieldType phiCorr
795  (
796  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
797  );
798 
799  return tmp<fluxFieldType>
800  (
801  new fluxFieldType
802  (
803  IOobject
804  (
805  "ddtCorr("
806  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
807  mesh().time().timeName(),
808  mesh()
809  ),
810  this->fvcDdtPhiCoeff
811  (
812  U.oldTime(),
813  phi.oldTime(),
814  phiCorr,
815  rho.oldTime()
816  )*rDeltaT*phiCorr
817  )
818  );
819  }
820  else
821  {
823  << "dimensions of phi are not correct"
824  << abort(FatalError);
825 
826  return fluxFieldType::null();
827  }
828 }
829 
830 
831 template<class Type>
833 (
835 )
836 {
838  (
840  (
841  IOobject
842  (
843  "meshPhi",
844  mesh().time().timeName(),
845  mesh(),
848  ),
849  mesh(),
851  )
852  );
853 }
854 
855 
856 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
857 
858 } // End namespace fv
859 
860 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
861 
862 } // End namespace Foam
863 
864 // ************************************************************************* //
fvsPatchField< scalar > fvsPatchScalarField
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
const word & name() const
Return name.
Definition: IOobject.H:297
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.
dynamicFvMesh & mesh
labelList fv(nPoints)
tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
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:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Calculate the divergence of the given field.
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
Internal & ref()
Return a reference to the dimensioned internal field.
label patchi
tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
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.
dimensioned< scalar > mag(const dimensioned< Type > &)
tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
A class for managing temporary objects.
Definition: PtrList.H:53
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
scalarField & diag()
Definition: lduMatrix.C:186
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
const surfaceScalarField & deltaCoeffs() const
Return reference to cell-centre difference coefficients.
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Namespace for OpenFOAM.
const dimensionSet dimVelocity