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-2024 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<VolField<Type>>
147 (
148  const dimensioned<Type>& dt
149 )
150 {
151  const volScalarField rDeltaT(CorDeltaT());
152 
153  const word ddtName("ddt("+dt.name()+')');
154 
155  if (mesh().moving())
156  {
157  tmp<VolField<Type>> tdtdt
158  (
160  (
161  ddtName,
162  mesh(),
164  (
165  "0",
166  dt.dimensions()/dimTime,
167  Zero
168  )
169  )
170  );
171 
172  tdtdt.ref().primitiveFieldRef() =
173  rDeltaT.primitiveField()*dt.value()
174  *(1.0 - mesh().Vsc0()/mesh().Vsc());
175 
176  return tdtdt;
177  }
178  else
179  {
180  return tmp<VolField<Type>>
181  (
183  (
184  ddtName,
185  mesh(),
187  (
188  "0",
189  dt.dimensions()/dimTime,
190  Zero
191  ),
193  )
194  );
195  }
196 }
197 
198 
199 template<class Type>
202 (
203  const VolField<Type>& vf
204 )
205 {
206  const volScalarField rDeltaT(CorDeltaT());
207 
208  const word ddtName("ddt("+vf.name()+')');
209 
210  if (mesh().moving())
211  {
212  return tmp<VolField<Type>>
213  (
215  (
216  ddtName,
217  rDeltaT()*
218  (
219  vf()
220  - vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
221  ),
222  rDeltaT.boundaryField()*
223  (
224  vf.boundaryField() - vf.oldTime().boundaryField()
225  )
226  )
227  );
228  }
229  else
230  {
231  return tmp<VolField<Type>>
232  (
234  (
235  ddtName,
236  rDeltaT*(vf - vf.oldTime())
237  )
238  );
239  }
240 }
241 
242 
243 template<class Type>
246 (
247  const dimensionedScalar& rho,
248  const VolField<Type>& vf
249 )
250 {
251  const volScalarField rDeltaT(CorDeltaT());
252 
253  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
254 
255  if (mesh().moving())
256  {
257  return tmp<VolField<Type>>
258  (
260  (
261  ddtName,
262  rDeltaT()*rho.value()*
263  (
264  vf()
265  - vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
266  ),
267  rDeltaT.boundaryField()*rho.value()*
268  (
269  vf.boundaryField() - vf.oldTime().boundaryField()
270  )
271  )
272  );
273  }
274  else
275  {
276  return tmp<VolField<Type>>
277  (
279  (
280  ddtName,
281  rDeltaT*rho*(vf - vf.oldTime())
282  )
283  );
284  }
285 }
286 
287 
288 template<class Type>
291 (
292  const volScalarField& rho,
293  const VolField<Type>& vf
294 )
295 {
296  const volScalarField rDeltaT(CorDeltaT());
297 
298  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
299 
300  if (mesh().moving())
301  {
302  return tmp<VolField<Type>>
303  (
305  (
306  ddtName,
307  rDeltaT()*
308  (
309  rho()*vf()
310  - rho.oldTime()()
311  *vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
312  ),
313  rDeltaT.boundaryField()*
314  (
315  rho.boundaryField()*vf.boundaryField()
316  - rho.oldTime().boundaryField()
317  *vf.oldTime().boundaryField()
318  )
319  )
320  );
321  }
322  else
323  {
324  return tmp<VolField<Type>>
325  (
327  (
328  ddtName,
329  rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
330  )
331  );
332  }
333 }
334 
335 
336 template<class Type>
339 (
340  const volScalarField& alpha,
341  const volScalarField& rho,
342  const VolField<Type>& vf
343 )
344 {
345  const volScalarField rDeltaT(CorDeltaT());
346 
347  const word ddtName("ddt("+alpha.name()+','+rho.name()+','+vf.name()+')');
348 
349  if (mesh().moving())
350  {
351  return tmp<VolField<Type>>
352  (
354  (
355  ddtName,
356  rDeltaT()*
357  (
358  alpha()
359  *rho()
360  *vf()
361 
362  - alpha.oldTime()()
363  *rho.oldTime()()
364  *vf.oldTime()()*mesh().Vsc0()/mesh().Vsc()
365  ),
366  rDeltaT.boundaryField()*
367  (
368  alpha.boundaryField()
369  *rho.boundaryField()
370  *vf.boundaryField()
371 
372  - alpha.oldTime().boundaryField()
373  *rho.oldTime().boundaryField()
374  *vf.oldTime().boundaryField()
375  )
376  )
377  );
378  }
379  else
380  {
381  return tmp<VolField<Type>>
382  (
384  (
385  ddtName,
386  rDeltaT
387  *(
388  alpha*rho*vf
389  - alpha.oldTime()*rho.oldTime()*vf.oldTime()
390  )
391  )
392  );
393  }
394 }
395 
396 
397 template<class Type>
400 (
401  const VolField<Type>& vf
402 )
403 {
404  tmp<fvMatrix<Type>> tfvm
405  (
406  new fvMatrix<Type>
407  (
408  vf,
410  )
411  );
412 
413  fvMatrix<Type>& fvm = tfvm.ref();
414 
415  const scalarField rDeltaT(CorDeltaT()().primitiveField());
416 
417  fvm.diag() = rDeltaT*mesh().Vsc();
418 
419  if (mesh().moving())
420  {
421  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc0();
422  }
423  else
424  {
425  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc();
426  }
427 
428  return tfvm;
429 }
430 
431 
432 template<class Type>
435 (
436  const dimensionedScalar& rho,
437  const VolField<Type>& vf
438 )
439 {
440  tmp<fvMatrix<Type>> tfvm
441  (
442  new fvMatrix<Type>
443  (
444  vf,
445  rho.dimensions()*vf.dimensions()*dimVolume/dimTime
446  )
447  );
448  fvMatrix<Type>& fvm = tfvm.ref();
449 
450  const scalarField rDeltaT(CorDeltaT()().primitiveField());
451 
452  fvm.diag() = rDeltaT*rho.value()*mesh().Vsc();
453 
454  if (mesh().moving())
455  {
456  fvm.source() = rDeltaT
457  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc0();
458  }
459  else
460  {
461  fvm.source() = rDeltaT
462  *rho.value()*vf.oldTime().primitiveField()*mesh().Vsc();
463  }
464 
465  return tfvm;
466 }
467 
468 
469 template<class Type>
472 (
473  const volScalarField& rho,
474  const VolField<Type>& vf
475 )
476 {
477  tmp<fvMatrix<Type>> tfvm
478  (
479  new fvMatrix<Type>
480  (
481  vf,
482  rho.dimensions()*vf.dimensions()*dimVolume/dimTime
483  )
484  );
485  fvMatrix<Type>& fvm = tfvm.ref();
486 
487  const scalarField rDeltaT(CorDeltaT()().primitiveField());
488 
489  fvm.diag() = rDeltaT*rho.primitiveField()*mesh().Vsc();
490 
491  if (mesh().moving())
492  {
493  fvm.source() = rDeltaT
494  *rho.oldTime().primitiveField()
495  *vf.oldTime().primitiveField()*mesh().Vsc0();
496  }
497  else
498  {
499  fvm.source() = rDeltaT
500  *rho.oldTime().primitiveField()
501  *vf.oldTime().primitiveField()*mesh().Vsc();
502  }
503 
504  return tfvm;
505 }
506 
507 
508 template<class Type>
511 (
512  const volScalarField& alpha,
513  const volScalarField& rho,
514  const VolField<Type>& vf
515 )
516 {
517  tmp<fvMatrix<Type>> tfvm
518  (
519  new fvMatrix<Type>
520  (
521  vf,
522  alpha.dimensions()
523  *rho.dimensions()
524  *vf.dimensions()
525  *dimVolume
526  /dimTime
527  )
528  );
529  fvMatrix<Type>& fvm = tfvm.ref();
530 
531  const scalarField rDeltaT(CorDeltaT()().primitiveField());
532 
533  fvm.diag() =
534  rDeltaT*alpha.primitiveField()*rho.primitiveField()*mesh().Vsc();
535 
536  if (mesh().moving())
537  {
538  fvm.source() = rDeltaT
539  *alpha.oldTime().primitiveField()
540  *rho.oldTime().primitiveField()
541  *vf.oldTime().primitiveField()*mesh().Vsc0();
542  }
543  else
544  {
545  fvm.source() = rDeltaT
546  *alpha.oldTime().primitiveField()
547  *rho.oldTime().primitiveField()
548  *vf.oldTime().primitiveField()*mesh().Vsc();
549  }
550 
551  return tfvm;
552 }
553 
554 
555 template<class Type>
558 (
559  const VolField<Type>& U,
560  const SurfaceField<Type>& Uf
561 )
562 {
563  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
564 
565  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
566  fluxFieldType phiCorr
567  (
568  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
569  );
570 
571  return fluxFieldType::New
572  (
573  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
574  this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)*rDeltaT*phiCorr
575  );
576 }
577 
578 
579 template<class Type>
582 (
583  const VolField<Type>& U,
584  const fluxFieldType& phi
585 )
586 {
587  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
588 
589  fluxFieldType phiCorr
590  (
591  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
592  );
593 
594  return fluxFieldType::New
595  (
596  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
597  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
598  *rDeltaT*phiCorr
599  );
600 }
601 
602 
603 template<class Type>
606 (
607  const volScalarField& rho,
608  const VolField<Type>& U,
609  const SurfaceField<Type>& rhoUf
610 )
611 {
612  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
613 
614  if
615  (
616  U.dimensions() == dimVelocity
617  && rhoUf.dimensions() == dimDensity*dimVelocity
618  )
619  {
620  VolField<Type> rhoU0
621  (
622  rho.oldTime()*U.oldTime()
623  );
624 
625  fluxFieldType phiUf0(mesh().Sf() & rhoUf.oldTime());
626  fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
627 
628  return fluxFieldType::New
629  (
630  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + rhoUf.name() + ')',
631  this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime())
632  *rDeltaT*phiCorr
633  );
634  }
635  else if
636  (
637  U.dimensions() == dimDensity*dimVelocity
638  && rhoUf.dimensions() == dimDensity*dimVelocity
639  )
640  {
641  fluxFieldType phiUf0(mesh().Sf() & rhoUf.oldTime());
642  fluxFieldType phiCorr
643  (
644  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
645  );
646 
647  return fluxFieldType::New
648  (
649  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + rhoUf.name() + ')',
650  this->fvcDdtPhiCoeff
651  (
652  U.oldTime(),
653  phiUf0,
654  phiCorr,
655  rho.oldTime()
656  )*rDeltaT*phiCorr
657  );
658  }
659  else
660  {
662  << "dimensions of rhoUf are not correct"
663  << abort(FatalError);
664 
665  return fluxFieldType::null();
666  }
667 }
668 
669 
670 template<class Type>
673 (
674  const volScalarField& rho,
675  const VolField<Type>& U,
676  const fluxFieldType& phi
677 )
678 {
679  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
680 
681  if
682  (
683  U.dimensions() == dimVelocity
684  && phi.dimensions() == rho.dimensions()*dimVolumetricFlux
685  )
686  {
687  VolField<Type> rhoU0
688  (
689  rho.oldTime()*U.oldTime()
690  );
691 
692  fluxFieldType phiCorr
693  (
694  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
695  );
696 
697  return fluxFieldType::New
698  (
699  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
700  this->fvcDdtPhiCoeff
701  (
702  rhoU0,
703  phi.oldTime(),
704  phiCorr,
705  rho.oldTime()
706  )*rDeltaT*phiCorr
707  );
708  }
709  else if
710  (
711  U.dimensions() == rho.dimensions()*dimVelocity
712  && phi.dimensions() == rho.dimensions()*dimVolumetricFlux
713  )
714  {
715  fluxFieldType phiCorr
716  (
717  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
718  );
719 
720  return fluxFieldType::New
721  (
722  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
723  this->fvcDdtPhiCoeff
724  (
725  U.oldTime(),
726  phi.oldTime(),
727  phiCorr,
728  rho.oldTime()
729  )*rDeltaT*phiCorr
730  );
731  }
732  else
733  {
735  << "dimensions of phi are not correct"
736  << abort(FatalError);
737 
738  return fluxFieldType::null();
739  }
740 }
741 
742 
743 template<class Type>
745 (
746  const VolField<Type>&
747 )
748 {
750  (
751  "meshPhi",
752  mesh(),
754  );
755 }
756 
757 
758 template<class Type>
760 (
761  const VolField<Type>&,
762  const label patchi
763 )
764 {
765  return tmp<scalarField>
766  (
767  new scalarField(mesh().boundary()[patchi].size(), 0)
768  );
769 }
770 
771 
772 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
773 
774 } // End namespace fv
775 
776 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
777 
778 } // End namespace Foam
779 
780 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const dimensionSet & dimensions() const
Return dimensions.
Generic GeometricField class.
const Internal::FieldType & primitiveField() const
Return a const-reference to the primitive field.
GeometricBoundaryField< Type, PatchField, GeoMesh > Boundary
Type of the boundary field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
static tmp< GeometricField< Type, PatchField, GeoMesh > > New(const word &name, const Internal &, const PtrList< PatchField< Type >> &, const HashPtrTable< Source > &=HashPtrTable< Source >())
Return a temporary field constructed from name,.
const word & name() const
Return name.
Definition: IOobject.H:310
const FieldType & oldTime() const
Return the old time field.
Definition: OldTimeField.C:315
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 > fvcDdtUfCorr(const VolField< Type > &U, const SurfaceField< Type > &Uf)
virtual tmp< fvMatrix< Type > > fvmDdt(const VolField< Type > &)
virtual tmp< surfaceScalarField > meshPhi(const VolField< Type > &)
virtual tmp< VolField< Type > > fvcDdt(const dimensioned< Type > &)
virtual tmp< fluxFieldType > fvcDdtPhiCorr(const VolField< Type > &U, const fluxFieldType &phi)
scalarField & diag()
Definition: lduMatrix.C:186
const surfaceScalarField & deltaCoeffs() const
Return reference to cell-centre difference coefficients.
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:334
A special matrix type and solver, designed for finite volume solutions of scalar equations.
Calculate the divergence of the given field.
label patchi
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
const dimensionSet dimVolumetricFlux
errorManip< error > abort(error &err)
Definition: errorManip.H:131
SurfaceField< scalar > surfaceScalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const dimensionSet dimTime
const dimensionSet dimDensity
const dimensionSet dimVolume
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:64
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
const dimensionSet dimVelocity
error FatalError
UList< label > labelUList
Definition: UList.H:65
fvsPatchField< scalar > fvsPatchScalarField
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
faceListList boundary(nPatches)
labelList fv(nPoints)
volScalarField & p