EulerDdtScheme.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-2023 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 "EulerDdtScheme.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<VolField<Type>>
46 (
47  const dimensioned<Type>& dt
48 )
49 {
50  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
51 
52  const word ddtName("ddt("+dt.name()+')');
53 
54  if (mesh().moving())
55  {
56  tmp<VolField<Type>> tdtdt
57  (
59  (
60  ddtName,
61  mesh(),
63  (
64  "0",
65  dt.dimensions()/dimTime,
66  Zero
67  )
68  )
69  );
70 
71  tdtdt.ref().primitiveFieldRef() =
72  rDeltaT.value()*dt.value()*(1.0 - mesh().Vsc0()/mesh().Vsc());
73 
74  return tdtdt;
75  }
76  else
77  {
78  return VolField<Type>::New
79  (
80  ddtName,
81  mesh(),
83  (
84  "0",
85  dt.dimensions()/dimTime,
86  Zero
87  ),
89  );
90  }
91 }
92 
93 
94 template<class Type>
97 (
98  const VolField<Type>& vf
99 )
100 {
101  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
102 
103  const word ddtName("ddt("+vf.name()+')');
104 
105  if (mesh().moving())
106  {
107  return VolField<Type>::New
108  (
109  ddtName,
110  rDeltaT*
111  (
112  vf()
113  - vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
114  ),
115  rDeltaT.value()*
116  (
117  vf.boundaryField() - vf.oldTime().boundaryField()
118  )
119  );
120  }
121  else
122  {
123  return VolField<Type>::New
124  (
125  ddtName,
126  rDeltaT*(vf - vf.oldTime())
127  );
128  }
129 }
130 
131 
132 template<class Type>
135 (
136  const dimensionedScalar& rho,
137  const VolField<Type>& vf
138 )
139 {
140  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
141 
142  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
143 
144  if (mesh().moving())
145  {
146  return VolField<Type>::New
147  (
148  ddtName,
149  rDeltaT*rho*
150  (
151  vf()
152  - vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
153  ),
154  rDeltaT.value()*rho.value()*
155  (
156  vf.boundaryField() - vf.oldTime().boundaryField()
157  )
158  );
159  }
160  else
161  {
162  return VolField<Type>::New
163  (
164  ddtName,
165  rDeltaT*rho*(vf - vf.oldTime())
166  );
167  }
168 }
169 
170 
171 template<class Type>
174 (
175  const volScalarField& rho,
176  const VolField<Type>& vf
177 )
178 {
179  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
180 
181  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
182 
183  if (mesh().moving())
184  {
185  return VolField<Type>::New
186  (
187  ddtName,
188  rDeltaT*
189  (
190  rho()*vf()
191  - rho.oldTime()()
192  *vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
193  ),
194  rDeltaT.value()*
195  (
196  rho.boundaryField()*vf.boundaryField()
197  - rho.oldTime().boundaryField()
198  *vf.oldTime().boundaryField()
199  )
200  );
201  }
202  else
203  {
204  return VolField<Type>::New
205  (
206  ddtName,
207  rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
208  );
209  }
210 }
211 
212 
213 template<class Type>
216 (
217  const volScalarField& alpha,
218  const volScalarField& rho,
219  const VolField<Type>& vf
220 )
221 {
222  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
223 
224  const word ddtName("ddt("+alpha.name()+','+rho.name()+','+vf.name()+')');
225 
226  if (mesh().moving())
227  {
228  return VolField<Type>::New
229  (
230  ddtName,
231  rDeltaT*
232  (
233  alpha()
234  *rho()
235  *vf()
236 
237  - alpha.oldTime()()
238  *rho.oldTime()()
239  *vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
240  ),
241  rDeltaT.value()*
242  (
243  alpha.boundaryField()
244  *rho.boundaryField()
245  *vf.boundaryField()
246 
247  - alpha.oldTime().boundaryField()
248  *rho.oldTime().boundaryField()
249  *vf.oldTime().boundaryField()
250  )
251  );
252  }
253  else
254  {
255  return VolField<Type>::New
256  (
257  ddtName,
258  rDeltaT
259  *(
260  alpha*rho*vf
261  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
262  )
263  );
264  }
265 }
266 
267 
268 template<class Type>
271 (
272  const SurfaceField<Type>& sf
273 )
274 {
275  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
276 
277  const word ddtName("ddt("+sf.name()+')');
278 
280  (
281  ddtName,
282  rDeltaT*(sf - sf.oldTime())
283  );
284 }
285 
286 
287 template<class Type>
290 (
291  const VolField<Type>& vf
292 )
293 {
294  tmp<fvMatrix<Type>> tfvm
295  (
296  new fvMatrix<Type>
297  (
298  vf,
300  )
301  );
302 
303  fvMatrix<Type>& fvm = tfvm.ref();
304 
305  const scalar rDeltaT = 1.0/mesh().time().deltaTValue();
306 
307  fvm.diag() = rDeltaT*mesh().Vsc();
308 
309  if (mesh().moving())
310  {
311  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc0();
312  }
313  else
314  {
315  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc();
316  }
317 
318  return tfvm;
319 }
320 
321 
322 template<class Type>
325 (
326  const dimensionedScalar& rho,
327  const VolField<Type>& vf
328 )
329 {
330  tmp<fvMatrix<Type>> tfvm
331  (
332  new fvMatrix<Type>
333  (
334  vf,
335  rho.dimensions()*vf.dimensions()*dimVol/dimTime
336  )
337  );
338  fvMatrix<Type>& fvm = tfvm.ref();
339 
340  const scalar rDeltaT = 1.0/mesh().time().deltaTValue();
341 
342  fvm.diag() = rDeltaT*rho.value()*mesh().Vsc();
343 
344  if (mesh().moving())
345  {
346  fvm.source() = rDeltaT
347  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc0();
348  }
349  else
350  {
351  fvm.source() = rDeltaT
352  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc();
353  }
354 
355  return tfvm;
356 }
357 
358 
359 template<class Type>
362 (
363  const volScalarField& rho,
364  const VolField<Type>& vf
365 )
366 {
367  tmp<fvMatrix<Type>> tfvm
368  (
369  new fvMatrix<Type>
370  (
371  vf,
372  rho.dimensions()*vf.dimensions()*dimVol/dimTime
373  )
374  );
375  fvMatrix<Type>& fvm = tfvm.ref();
376 
377  const scalar rDeltaT = 1.0/mesh().time().deltaTValue();
378 
379  fvm.diag() = rDeltaT*rho.primitiveField()*mesh().Vsc();
380 
381  if (mesh().moving())
382  {
383  fvm.source() = rDeltaT
384  *rho.oldTime().primitiveField()
385  *vf.oldTime().primitiveField()*mesh().Vsc0();
386  }
387  else
388  {
389  fvm.source() = rDeltaT
390  *rho.oldTime().primitiveField()
391  *vf.oldTime().primitiveField()*mesh().Vsc();
392  }
393 
394  return tfvm;
395 }
396 
397 
398 template<class Type>
401 (
402  const volScalarField& alpha,
403  const volScalarField& rho,
404  const VolField<Type>& vf
405 )
406 {
407  tmp<fvMatrix<Type>> tfvm
408  (
409  new fvMatrix<Type>
410  (
411  vf,
412  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
413  )
414  );
415  fvMatrix<Type>& fvm = tfvm.ref();
416 
417  const scalar rDeltaT = 1.0/mesh().time().deltaTValue();
418 
419  fvm.diag() =
420  rDeltaT*alpha.primitiveField()*rho.primitiveField()*mesh().Vsc();
421 
422  if (mesh().moving())
423  {
424  fvm.source() = rDeltaT
425  *alpha.oldTime().primitiveField()
426  *rho.oldTime().primitiveField()
427  *vf.oldTime().primitiveField()*mesh().Vsc0();
428  }
429  else
430  {
431  fvm.source() = rDeltaT
432  *alpha.oldTime().primitiveField()
433  *rho.oldTime().primitiveField()
434  *vf.oldTime().primitiveField()*mesh().Vsc();
435  }
436 
437  return tfvm;
438 }
439 
440 
441 template<class Type>
444 (
445  const VolField<Type>& U,
446  const SurfaceField<Type>& Uf
447 )
448 {
449  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
450 
451  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
452  fluxFieldType phiCorr
453  (
454  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
455  );
456 
457  return fluxFieldType::New
458  (
459  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
460  this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr) *rDeltaT*phiCorr
461  );
462 }
463 
464 
465 template<class Type>
468 (
469  const VolField<Type>& U,
470  const fluxFieldType& phi
471 )
472 {
473  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
474 
475  fluxFieldType phiCorr
476  (
477  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
478  );
479 
480  return fluxFieldType::New
481  (
482  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
483  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
484  *rDeltaT*phiCorr
485  );
486 }
487 
488 
489 template<class Type>
492 (
493  const volScalarField& rho,
494  const VolField<Type>& U,
495  const SurfaceField<Type>& rhoUf
496 )
497 {
498  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
499 
500  if
501  (
502  U.dimensions() == dimVelocity
503  && rhoUf.dimensions() == rho.dimensions()*dimVelocity
504  )
505  {
506  VolField<Type> rhoU0(rho.oldTime()*U.oldTime());
507  fluxFieldType phiUf0(mesh().Sf() & rhoUf.oldTime());
508  fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
509 
510  return fluxFieldType::New
511  (
512  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + rhoUf.name() + ')',
513  this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime())
514  *rDeltaT*phiCorr
515  );
516  }
517  else if
518  (
519  U.dimensions() == rho.dimensions()*dimVelocity
520  && rhoUf.dimensions() == rho.dimensions()*dimVelocity
521  )
522  {
523  fluxFieldType phiUf0(mesh().Sf() & rhoUf.oldTime());
524  fluxFieldType phiCorr
525  (
526  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
527  );
528 
529  return fluxFieldType::New
530  (
531  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + rhoUf.name() + ')',
532  this->fvcDdtPhiCoeff
533  (
534  U.oldTime(),
535  phiUf0,
536  phiCorr,
537  rho.oldTime()
538  )*rDeltaT*phiCorr
539  );
540  }
541  else
542  {
544  << "dimensions of rhoUf are not correct"
545  << abort(FatalError);
546 
547  return fluxFieldType::null();
548  }
549 }
550 
551 
552 template<class Type>
555 (
556  const volScalarField& rho,
557  const VolField<Type>& U,
558  const fluxFieldType& phi
559 )
560 {
561  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
562 
563  if
564  (
565  U.dimensions() == dimVelocity
566  && phi.dimensions() == rho.dimensions()*dimFlux
567  )
568  {
569  VolField<Type> rhoU0
570  (
571  rho.oldTime()*U.oldTime()
572  );
573 
574  fluxFieldType phiCorr
575  (
576  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
577  );
578 
579  return fluxFieldType::New
580  (
581  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
582  this->fvcDdtPhiCoeff
583  (
584  rhoU0,
585  phi.oldTime(),
586  phiCorr,
587  rho.oldTime()
588  )*rDeltaT*phiCorr
589  );
590  }
591  else if
592  (
593  U.dimensions() == rho.dimensions()*dimVelocity
594  && phi.dimensions() == rho.dimensions()*dimFlux
595  )
596  {
597  fluxFieldType phiCorr
598  (
599  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
600  );
601 
602  return fluxFieldType::New
603  (
604  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
605  this->fvcDdtPhiCoeff
606  (
607  U.oldTime(),
608  phi.oldTime(),
609  phiCorr,
610  rho.oldTime()
611  )*rDeltaT*phiCorr
612  );
613  }
614  else
615  {
617  << "dimensions of phi are not correct"
618  << abort(FatalError);
619 
620  return fluxFieldType::null();
621  }
622 }
623 
624 
625 template<class Type>
628 (
629  const volScalarField& alpha,
630  const volScalarField& rho,
631  const VolField<Type>& U,
632  const SurfaceField<Type>& Uf
633 )
634 {
635  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
636 
637  if (U.dimensions() == dimVelocity && Uf.dimensions() == dimVelocity)
638  {
639  const volScalarField alphaRho0(alpha.oldTime()*rho.oldTime());
640 
641  return fluxFieldType::New
642  (
643  "ddtCorr("
644  + alpha.name() + rho.name() + ',' + U.name() + ',' + Uf.name()
645  + ')',
646  this->fvcDdtPhiCoeff(U.oldTime(), mesh().Sf() & Uf.oldTime())
647  *rDeltaT
648  *(
649  mesh().Sf()
650  & (
651  fvc::interpolate(alphaRho0)*Uf.oldTime()
652  - fvc::interpolate(alphaRho0*U.oldTime())
653  )
654  )
655  );
656  }
657  else
658  {
660  << "dimensions of Uf are not correct"
661  << abort(FatalError);
662 
663  return fluxFieldType::null();
664  }
665 }
666 
667 
668 template<class Type>
671 (
672  const volScalarField& alpha,
673  const volScalarField& rho,
674  const VolField<Type>& U,
675  const fluxFieldType& phi
676 )
677 {
678  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
679 
680  if (U.dimensions() == dimVelocity && phi.dimensions() == dimFlux)
681  {
682  const volScalarField alphaRho0(alpha.oldTime()*rho.oldTime());
683 
684  return fluxFieldType::New
685  (
686  "ddtCorr("
687  + alpha.name() + rho.name() + ',' + U.name() + ',' + phi.name()
688  + ')',
689  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
690  *rDeltaT
691  *(
692  fvc::interpolate(alphaRho0)*phi.oldTime()
693  -fvc::dotInterpolate(mesh().Sf(), alphaRho0*U.oldTime())
694  )
695  );
696  }
697  else
698  {
700  << "dimensions of phi are not correct"
701  << abort(FatalError);
702 
703  return fluxFieldType::null();
704  }
705 }
706 
707 
708 template<class Type>
710 (
711  const VolField<Type>&
712 )
713 {
714  return mesh().phi();
715 }
716 
717 
718 template<class Type>
720 (
721  const VolField<Type>&,
722  const label patchi
723 )
724 {
725  return mesh().phi().boundaryField()[patchi];
726 }
727 
728 
729 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
730 
731 } // End namespace fv
732 
733 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
734 
735 } // End namespace Foam
736 
737 // ************************************************************************* //
const dimensionSet & dimensions() const
Return dimensions.
Generic GeometricField class.
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
static tmp< GeometricField< Type, PatchField, GeoMesh > > New(const word &name, const Internal &, const PtrList< PatchField< Type >> &)
Return a temporary field constructed from name,.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
const word & name() const
Return name.
Definition: IOobject.H:310
This boundary condition is not designed to be evaluated; it is assumed that the value is assigned via...
Generic dimensioned Type class.
const dimensionSet & dimensions() const
Return const reference to dimensions.
const Type & value() const
Return const reference to value.
const word & name() const
Return const reference to name.
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:118
Field< Type > & source()
Definition: fvMatrix.H:307
virtual tmp< fluxFieldType > fvcDdtPhiCorr(const VolField< Type > &U, const fluxFieldType &phi)
virtual tmp< fvMatrix< Type > > fvmDdt(const VolField< Type > &)
virtual tmp< surfaceScalarField > meshPhi(const VolField< Type > &)
virtual tmp< fluxFieldType > fvcDdtUfCorr(const VolField< Type > &U, const SurfaceField< Type > &Uf)
virtual tmp< VolField< Type > > fvcDdt(const dimensioned< Type > &)
scalarField & diag()
Definition: lduMatrix.C:186
A class for managing temporary objects.
Definition: tmp.H:55
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
A class for handling words, derived from string.
Definition: word.H:62
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
A special matrix type and solver, designed for finite volume solutions of scalar equations.
Calculate the divergence of the given field.
label patchi
volScalarField sf(fieldObject, mesh)
U
Definition: pEqn.H:72
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
static tmp< SurfaceField< Type > > interpolate(const VolField< Type > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
static tmp< SurfaceField< typename innerProduct< vector, Type >::type > > dotInterpolate(const surfaceVectorField &Sf, const VolField< Type > &tvf)
Interpolate field onto faces.
Namespace for OpenFOAM.
static const zero Zero
Definition: zero.H:97
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const dimensionSet dimTime
const dimensionSet dimVol
const dimensionSet dimVelocity
error FatalError
const dimensionSet dimFlux
labelList fv(nPoints)