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-2019 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  const word ddtName("ddt("+dt.name()+')');
154 
155  if (mesh().moving())
156  {
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  {
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 (
204 )
205 {
206  const volScalarField rDeltaT(CorDeltaT());
207 
208  const word ddtName("ddt("+vf.name()+')');
209 
210  if (mesh().moving())
211  {
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  {
232  (
234  (
235  ddtName,
236  rDeltaT*(vf - vf.oldTime())
237  )
238  );
239  }
240 }
241 
242 
243 template<class Type>
246 (
247  const dimensionedScalar& rho,
249 )
250 {
251  const volScalarField rDeltaT(CorDeltaT());
252 
253  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
254 
255  if (mesh().moving())
256  {
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  {
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,
294 )
295 {
296  const volScalarField rDeltaT(CorDeltaT());
297 
298  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
299 
300  if (mesh().moving())
301  {
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  {
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,
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  {
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  {
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 (
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,
438 )
439 {
440  tmp<fvMatrix<Type>> tfvm
441  (
442  new fvMatrix<Type>
443  (
444  vf,
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,
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  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,
515 )
516 {
517  tmp<fvMatrix<Type>> tfvm
518  (
519  new fvMatrix<Type>
520  (
521  vf,
522  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
523  )
524  );
525  fvMatrix<Type>& fvm = tfvm.ref();
526 
527  const scalarField rDeltaT(CorDeltaT()().primitiveField());
528 
529  fvm.diag() =
530  rDeltaT*alpha.primitiveField()*rho.primitiveField()*mesh().Vsc();
531 
532  if (mesh().moving())
533  {
534  fvm.source() = rDeltaT
535  *alpha.oldTime().primitiveField()
536  *rho.oldTime().primitiveField()
537  *vf.oldTime().primitiveField()*mesh().Vsc0();
538  }
539  else
540  {
541  fvm.source() = rDeltaT
542  *alpha.oldTime().primitiveField()
543  *rho.oldTime().primitiveField()
544  *vf.oldTime().primitiveField()*mesh().Vsc();
545  }
546 
547  return tfvm;
548 }
549 
550 
551 template<class Type>
554 (
557 )
558 {
559  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
560 
561  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
562  fluxFieldType phiCorr
563  (
564  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
565  );
566 
567  return fluxFieldType::New
568  (
569  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
570  this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)*rDeltaT*phiCorr
571  );
572 }
573 
574 
575 template<class Type>
578 (
580  const fluxFieldType& phi
581 )
582 {
583  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
584 
585  fluxFieldType phiCorr
586  (
587  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
588  );
589 
590  return fluxFieldType::New
591  (
592  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
593  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
594  *rDeltaT*phiCorr
595  );
596 }
597 
598 
599 template<class Type>
602 (
603  const volScalarField& rho,
606 )
607 {
608  const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
609 
610  if
611  (
612  U.dimensions() == dimVelocity
614  )
615  {
617  (
618  rho.oldTime()*U.oldTime()
619  );
620 
621  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
622  fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
623 
624  return fluxFieldType::New
625  (
626  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
627  this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime())
628  *rDeltaT*phiCorr
629  );
630  }
631  else if
632  (
635  )
636  {
637  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
638  fluxFieldType phiCorr
639  (
640  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
641  );
642 
643  return fluxFieldType::New
644  (
645  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
646  this->fvcDdtPhiCoeff
647  (
648  U.oldTime(),
649  phiUf0,
650  phiCorr,
651  rho.oldTime()
652  )*rDeltaT*phiCorr
653  );
654  }
655  else
656  {
658  << "dimensions of Uf are not correct"
659  << abort(FatalError);
660 
661  return fluxFieldType::null();
662  }
663 }
664 
665 
666 template<class Type>
669 (
670  const volScalarField& rho,
672  const fluxFieldType& phi
673 )
674 {
675  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
676 
677  if
678  (
679  U.dimensions() == dimVelocity
680  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
681  )
682  {
684  (
685  rho.oldTime()*U.oldTime()
686  );
687 
688  fluxFieldType phiCorr
689  (
690  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
691  );
692 
693  return fluxFieldType::New
694  (
695  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
696  this->fvcDdtPhiCoeff
697  (
698  rhoU0,
699  phi.oldTime(),
700  phiCorr,
701  rho.oldTime()
702  )*rDeltaT*phiCorr
703  );
704  }
705  else if
706  (
707  U.dimensions() == rho.dimensions()*dimVelocity
708  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
709  )
710  {
711  fluxFieldType phiCorr
712  (
713  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
714  );
715 
716  return fluxFieldType::New
717  (
718  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
719  this->fvcDdtPhiCoeff
720  (
721  U.oldTime(),
722  phi.oldTime(),
723  phiCorr,
724  rho.oldTime()
725  )*rDeltaT*phiCorr
726  );
727  }
728  else
729  {
731  << "dimensions of phi are not correct"
732  << abort(FatalError);
733 
734  return fluxFieldType::null();
735  }
736 }
737 
738 
739 template<class Type>
741 (
743 )
744 {
746  (
747  "meshPhi",
748  mesh(),
750  );
751 }
752 
753 
754 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
755 
756 } // End namespace fv
757 
758 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
759 
760 } // End namespace Foam
761 
762 // ************************************************************************* //
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:303
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.
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< scalar >> &)
Return a temporary field constructed from name,.
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
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
const dimensionSet & dimensions() const
Return dimensions.
dynamicFvMesh & mesh
autoPtr< BasicCompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleMomentumTransportModel::transportModel &transport)
A class for handling words, derived from string.
Definition: word.H:59
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 assumed 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)
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
A class for managing temporary objects.
Definition: PtrList.H:53
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
scalarField & diag()
Definition: lduMatrix.C:186
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