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  (
51  (
52  "CorDeltaT",
53  mesh(),
54  dimensionedScalar(cofrDeltaT.dimensions(), 0),
55  extrapolatedCalculatedFvPatchScalarField::typeName
56  )
57  );
58 
59  volScalarField& corDeltaT = tcorDeltaT.ref();
60 
61  const labelUList& owner = mesh().owner();
62  const labelUList& neighbour = mesh().neighbour();
63 
64  forAll(owner, facei)
65  {
66  corDeltaT[owner[facei]] =
67  max(corDeltaT[owner[facei]], cofrDeltaT[facei]);
68 
69  corDeltaT[neighbour[facei]] =
70  max(corDeltaT[neighbour[facei]], cofrDeltaT[facei]);
71  }
72 
73  const surfaceScalarField::Boundary& cofrDeltaTbf =
74  cofrDeltaT.boundaryField();
75 
76  forAll(cofrDeltaTbf, patchi)
77  {
78  const fvsPatchScalarField& pcofrDeltaT = cofrDeltaTbf[patchi];
79  const fvPatch& p = pcofrDeltaT.patch();
80  const labelUList& faceCells = p.patch().faceCells();
81 
82  forAll(pcofrDeltaT, patchFacei)
83  {
84  corDeltaT[faceCells[patchFacei]] = max
85  (
86  corDeltaT[faceCells[patchFacei]],
87  pcofrDeltaT[patchFacei]
88  );
89  }
90  }
91 
92  corDeltaT.correctBoundaryConditions();
93 
94  return tcorDeltaT;
95 }
96 
97 
98 template<class Type>
99 tmp<surfaceScalarField> CoEulerDdtScheme<Type>::CofrDeltaT() const
100 {
101  const dimensionedScalar& deltaT = mesh().time().deltaT();
102 
103  const surfaceScalarField& phi =
104  static_cast<const objectRegistry&>(mesh())
105  .lookupObject<surfaceScalarField>(phiName_);
106 
107  if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
108  {
110  (
112  *(mag(phi)/mesh().magSf())
113  *deltaT
114  );
115 
116  return max(Co/maxCo_, scalar(1))/deltaT;
117  }
118  else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
119  {
120  const volScalarField& rho =
121  static_cast<const objectRegistry&>(mesh())
122  .lookupObject<volScalarField>(rhoName_).oldTime();
123 
125  (
127  *(mag(phi)/(fvc::interpolate(rho)*mesh().magSf()))
128  *deltaT
129  );
130 
131  return max(Co/maxCo_, scalar(1))/deltaT;
132  }
133  else
134  {
136  << "Incorrect dimensions of phi: " << phi.dimensions()
137  << abort(FatalError);
138 
139  return tmp<surfaceScalarField>(nullptr);
140  }
141 }
142 
143 
144 template<class Type>
145 tmp<GeometricField<Type, fvPatchField, volMesh>>
147 (
148  const dimensioned<Type>& dt
149 )
150 {
151  const volScalarField rDeltaT(CorDeltaT());
152 
153  IOobject ddtIOobject
154  (
155  "ddt("+dt.name()+')',
156  mesh().time().timeName(),
157  mesh()
158  );
159 
160  if (mesh().moving())
161  {
163  (
165  (
166  ddtIOobject,
167  mesh(),
169  (
170  "0",
171  dt.dimensions()/dimTime,
172  Zero
173  )
174  )
175  );
176 
177  tdtdt.ref().primitiveFieldRef() =
178  rDeltaT.primitiveField()*dt.value()
179  *(1.0 - mesh().Vsc0()/mesh().Vsc());
180 
181  return tdtdt;
182  }
183  else
184  {
186  (
188  (
189  ddtIOobject,
190  mesh(),
192  (
193  "0",
194  dt.dimensions()/dimTime,
195  Zero
196  ),
198  )
199  );
200  }
201 }
202 
203 
204 template<class Type>
207 (
209 )
210 {
211  const volScalarField rDeltaT(CorDeltaT());
212 
213  IOobject ddtIOobject
214  (
215  "ddt("+vf.name()+')',
216  mesh().time().timeName(),
217  mesh()
218  );
219 
220  if (mesh().moving())
221  {
223  (
225  (
226  ddtIOobject,
227  mesh(),
228  rDeltaT.dimensions()*vf.dimensions(),
229  rDeltaT.primitiveField()*
230  (
231  vf.primitiveField()
232  - vf.oldTime().primitiveField()*mesh().Vsc0()/mesh().Vsc()
233  ),
234  rDeltaT.boundaryField()*
235  (
236  vf.boundaryField() - vf.oldTime().boundaryField()
237  )
238  )
239  );
240  }
241  else
242  {
244  (
246  (
247  ddtIOobject,
248  rDeltaT*(vf - vf.oldTime())
249  )
250  );
251  }
252 }
253 
254 
255 template<class Type>
258 (
259  const dimensionedScalar& rho,
261 )
262 {
263  const volScalarField rDeltaT(CorDeltaT());
264 
265  IOobject ddtIOobject
266  (
267  "ddt("+rho.name()+','+vf.name()+')',
268  mesh().time().timeName(),
269  mesh()
270  );
271 
272  if (mesh().moving())
273  {
275  (
277  (
278  ddtIOobject,
279  mesh(),
280  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
281  rDeltaT.primitiveField()*rho.value()*
282  (
283  vf.primitiveField()
284  - vf.oldTime().primitiveField()*mesh().Vsc0()/mesh().Vsc()
285  ),
286  rDeltaT.boundaryField()*rho.value()*
287  (
288  vf.boundaryField() - vf.oldTime().boundaryField()
289  )
290  )
291  );
292  }
293  else
294  {
296  (
298  (
299  ddtIOobject,
300  rDeltaT*rho*(vf - vf.oldTime())
301  )
302  );
303  }
304 }
305 
306 
307 template<class Type>
310 (
311  const volScalarField& rho,
313 )
314 {
315  const volScalarField rDeltaT(CorDeltaT());
316 
317  IOobject ddtIOobject
318  (
319  "ddt("+rho.name()+','+vf.name()+')',
320  mesh().time().timeName(),
321  mesh()
322  );
323 
324  if (mesh().moving())
325  {
327  (
329  (
330  ddtIOobject,
331  mesh(),
332  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
333  rDeltaT.primitiveField()*
334  (
335  rho.primitiveField()*vf.primitiveField()
336  - rho.oldTime().primitiveField()
337  *vf.oldTime().primitiveField()*mesh().Vsc0()/mesh().Vsc()
338  ),
339  rDeltaT.boundaryField()*
340  (
341  rho.boundaryField()*vf.boundaryField()
342  - rho.oldTime().boundaryField()
343  *vf.oldTime().boundaryField()
344  )
345  )
346  );
347  }
348  else
349  {
351  (
353  (
354  ddtIOobject,
355  rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
356  )
357  );
358  }
359 }
360 
361 
362 template<class Type>
365 (
366  const volScalarField& alpha,
367  const volScalarField& rho,
369 )
370 {
371  const volScalarField rDeltaT(CorDeltaT());
372 
373  IOobject ddtIOobject
374  (
375  "ddt("+alpha.name()+','+rho.name()+','+vf.name()+')',
376  mesh().time().timeName(),
377  mesh()
378  );
379 
380  if (mesh().moving())
381  {
383  (
385  (
386  ddtIOobject,
387  mesh(),
388  rDeltaT.dimensions()
389  *alpha.dimensions()*rho.dimensions()*vf.dimensions(),
390  rDeltaT.primitiveField()*
391  (
392  alpha.primitiveField()
393  *rho.primitiveField()
394  *vf.primitiveField()
395 
396  - alpha.oldTime().primitiveField()
397  *rho.oldTime().primitiveField()
398  *vf.oldTime().primitiveField()*mesh().Vsc0()/mesh().Vsc()
399  ),
400  rDeltaT.boundaryField()*
401  (
402  alpha.boundaryField()
403  *rho.boundaryField()
404  *vf.boundaryField()
405 
406  - alpha.oldTime().boundaryField()
407  *rho.oldTime().boundaryField()
408  *vf.oldTime().boundaryField()
409  )
410  )
411  );
412  }
413  else
414  {
416  (
418  (
419  ddtIOobject,
420  rDeltaT
421  *(
422  alpha*rho*vf
423  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
424  )
425  )
426  );
427  }
428 }
429 
430 
431 template<class Type>
434 (
436 )
437 {
438  tmp<fvMatrix<Type>> tfvm
439  (
440  new fvMatrix<Type>
441  (
442  vf,
444  )
445  );
446 
447  fvMatrix<Type>& fvm = tfvm.ref();
448 
449  scalarField rDeltaT(CorDeltaT()().primitiveField());
450 
451  fvm.diag() = rDeltaT*mesh().Vsc();
452 
453  if (mesh().moving())
454  {
455  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc0();
456  }
457  else
458  {
459  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc();
460  }
461 
462  return tfvm;
463 }
464 
465 
466 template<class Type>
469 (
470  const dimensionedScalar& rho,
472 )
473 {
474  tmp<fvMatrix<Type>> tfvm
475  (
476  new fvMatrix<Type>
477  (
478  vf,
480  )
481  );
482  fvMatrix<Type>& fvm = tfvm.ref();
483 
484  scalarField rDeltaT(CorDeltaT()().primitiveField());
485 
486  fvm.diag() = rDeltaT*rho.value()*mesh().Vsc();
487 
488  if (mesh().moving())
489  {
490  fvm.source() = rDeltaT
491  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc0();
492  }
493  else
494  {
495  fvm.source() = rDeltaT
496  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc();
497  }
498 
499  return tfvm;
500 }
501 
502 
503 template<class Type>
506 (
507  const volScalarField& rho,
509 )
510 {
511  tmp<fvMatrix<Type>> tfvm
512  (
513  new fvMatrix<Type>
514  (
515  vf,
517  )
518  );
519  fvMatrix<Type>& fvm = tfvm.ref();
520 
521  scalarField rDeltaT(CorDeltaT()().primitiveField());
522 
523  fvm.diag() = rDeltaT*rho.primitiveField()*mesh().Vsc();
524 
525  if (mesh().moving())
526  {
527  fvm.source() = rDeltaT
528  *rho.oldTime().primitiveField()
529  *vf.oldTime().primitiveField()*mesh().Vsc0();
530  }
531  else
532  {
533  fvm.source() = rDeltaT
534  *rho.oldTime().primitiveField()
535  *vf.oldTime().primitiveField()*mesh().Vsc();
536  }
537 
538  return tfvm;
539 }
540 
541 
542 template<class Type>
545 (
546  const volScalarField& alpha,
547  const volScalarField& rho,
549 )
550 {
551  tmp<fvMatrix<Type>> tfvm
552  (
553  new fvMatrix<Type>
554  (
555  vf,
556  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
557  )
558  );
559  fvMatrix<Type>& fvm = tfvm.ref();
560 
561  scalarField rDeltaT(CorDeltaT()().primitiveField());
562 
563  fvm.diag() =
564  rDeltaT*alpha.primitiveField()*rho.primitiveField()*mesh().Vsc();
565 
566  if (mesh().moving())
567  {
568  fvm.source() = rDeltaT
569  *alpha.oldTime().primitiveField()
570  *rho.oldTime().primitiveField()
571  *vf.oldTime().primitiveField()*mesh().Vsc0();
572  }
573  else
574  {
575  fvm.source() = rDeltaT
576  *alpha.oldTime().primitiveField()
577  *rho.oldTime().primitiveField()
578  *vf.oldTime().primitiveField()*mesh().Vsc();
579  }
580 
581  return tfvm;
582 }
583 
584 
585 template<class Type>
588 (
591 )
592 {
593  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
594 
595  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
596  fluxFieldType phiCorr
597  (
598  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
599  );
600 
601  return tmp<fluxFieldType>
602  (
603  new fluxFieldType
604  (
605  IOobject
606  (
607  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
608  mesh().time().timeName(),
609  mesh()
610  ),
611  this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
612  *rDeltaT*phiCorr
613  )
614  );
615 }
616 
617 
618 template<class Type>
621 (
623  const fluxFieldType& phi
624 )
625 {
626  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
627 
628  fluxFieldType phiCorr
629  (
630  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
631  );
632 
633  return tmp<fluxFieldType>
634  (
635  new fluxFieldType
636  (
637  IOobject
638  (
639  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
640  mesh().time().timeName(),
641  mesh()
642  ),
643  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
644  *rDeltaT*phiCorr
645  )
646  );
647 }
648 
649 
650 template<class Type>
653 (
654  const volScalarField& rho,
657 )
658 {
659  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
660 
661  if
662  (
663  U.dimensions() == dimVelocity
665  )
666  {
668  (
669  rho.oldTime()*U.oldTime()
670  );
671 
672  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
673  fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
674 
675  return tmp<fluxFieldType>
676  (
677  new fluxFieldType
678  (
679  IOobject
680  (
681  "ddtCorr("
682  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
683  mesh().time().timeName(),
684  mesh()
685  ),
686  this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime())
687  *rDeltaT*phiCorr
688  )
689  );
690  }
691  else if
692  (
695  )
696  {
697  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
698  fluxFieldType phiCorr
699  (
700  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
701  );
702 
703  return tmp<fluxFieldType>
704  (
705  new fluxFieldType
706  (
707  IOobject
708  (
709  "ddtCorr("
710  + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
711  mesh().time().timeName(),
712  mesh()
713  ),
714  this->fvcDdtPhiCoeff
715  (
716  U.oldTime(),
717  phiUf0,
718  phiCorr,
719  rho.oldTime()
720  )*rDeltaT*phiCorr
721  )
722  );
723  }
724  else
725  {
727  << "dimensions of Uf are not correct"
728  << abort(FatalError);
729 
730  return fluxFieldType::null();
731  }
732 }
733 
734 
735 template<class Type>
738 (
739  const volScalarField& rho,
741  const fluxFieldType& phi
742 )
743 {
744  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
745 
746  if
747  (
748  U.dimensions() == dimVelocity
749  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
750  )
751  {
753  (
754  rho.oldTime()*U.oldTime()
755  );
756 
757  fluxFieldType phiCorr
758  (
759  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
760  );
761 
762  return tmp<fluxFieldType>
763  (
764  new fluxFieldType
765  (
766  IOobject
767  (
768  "ddtCorr("
769  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
770  mesh().time().timeName(),
771  mesh()
772  ),
773  this->fvcDdtPhiCoeff
774  (
775  rhoU0,
776  phi.oldTime(),
777  phiCorr,
778  rho.oldTime()
779  )*rDeltaT*phiCorr
780  )
781  );
782  }
783  else if
784  (
785  U.dimensions() == rho.dimensions()*dimVelocity
786  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
787  )
788  {
789  fluxFieldType phiCorr
790  (
791  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
792  );
793 
794  return tmp<fluxFieldType>
795  (
796  new fluxFieldType
797  (
798  IOobject
799  (
800  "ddtCorr("
801  + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
802  mesh().time().timeName(),
803  mesh()
804  ),
805  this->fvcDdtPhiCoeff
806  (
807  U.oldTime(),
808  phi.oldTime(),
809  phiCorr,
810  rho.oldTime()
811  )*rDeltaT*phiCorr
812  )
813  );
814  }
815  else
816  {
818  << "dimensions of phi are not correct"
819  << abort(FatalError);
820 
821  return fluxFieldType::null();
822  }
823 }
824 
825 
826 template<class Type>
828 (
830 )
831 {
833  (
834  "meshPhi",
835  mesh(),
837  );
838 }
839 
840 
841 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
842 
843 } // End namespace fv
844 
845 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
846 
847 } // End namespace Foam
848 
849 // ************************************************************************* //
fvsPatchField< scalar > fvsPatchScalarField
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
virtual 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: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.
dynamicFvMesh & mesh
labelList fv(nPoints)
virtual tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
virtual 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
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: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
Internal & ref()
Return a reference to the dimensioned internal field.
label patchi
virtual 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 > &)
virtual tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
mesh Sf()
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