SLTSDdtScheme.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 "SLTSDdtScheme.H"
27 #include "surfaceInterpolate.H"
28 #include "fvMatrices.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace fv
38 {
39 
40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41 
42 template<class Type>
43 void SLTSDdtScheme<Type>::relaxedDiag
44 (
45  scalarField& rD,
46  const surfaceScalarField& phi
47 ) const
48 {
49  const labelUList& owner = mesh().owner();
50  const labelUList& neighbour = mesh().neighbour();
51  scalarField diag(rD.size(), 0.0);
52 
53  forAll(owner, facei)
54  {
55  if (phi[facei] > 0.0)
56  {
57  diag[owner[facei]] += phi[facei];
58  rD[neighbour[facei]] += phi[facei];
59  }
60  else
61  {
62  diag[neighbour[facei]] -= phi[facei];
63  rD[owner[facei]] -= phi[facei];
64  }
65  }
66 
67  forAll(phi.boundaryField(), patchi)
68  {
69  const fvsPatchScalarField& pphi = phi.boundaryField()[patchi];
70  const labelUList& faceCells = pphi.patch().patch().faceCells();
71 
72  forAll(pphi, patchFacei)
73  {
74  if (pphi[patchFacei] > 0.0)
75  {
76  diag[faceCells[patchFacei]] += pphi[patchFacei];
77  }
78  else
79  {
80  rD[faceCells[patchFacei]] -= pphi[patchFacei];
81  }
82  }
83  }
84 
85  rD += (1.0/alpha_ - 2.0)*diag;
86 }
87 
88 
89 template<class Type>
90 tmp<volScalarField> SLTSDdtScheme<Type>::SLrDeltaT() const
91 {
92  const surfaceScalarField& phi =
93  mesh().objectRegistry::template
94  lookupObject<surfaceScalarField>(phiName_);
95 
96  const dimensionedScalar& deltaT = mesh().time().deltaT();
97 
98  tmp<volScalarField> trDeltaT
99  (
101  (
102  "rDeltaT",
103  mesh(),
105  extrapolatedCalculatedFvPatchScalarField::typeName
106  )
107  );
108 
109  volScalarField& rDeltaT = trDeltaT.ref();
110 
111  relaxedDiag(rDeltaT, phi);
112 
113  if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
114  {
115  rDeltaT.primitiveFieldRef() = max
116  (
117  rDeltaT.primitiveField()/mesh().V(),
118  scalar(1)/deltaT.value()
119  );
120  }
121  else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
122  {
123  const volScalarField& rho =
124  mesh().objectRegistry::template lookupObject<volScalarField>
125  (
126  rhoName_
127  ).oldTime();
128 
129  rDeltaT.primitiveFieldRef() = max
130  (
131  rDeltaT.primitiveField()/(rho.primitiveField()*mesh().V()),
132  scalar(1)/deltaT.value()
133  );
134  }
135  else
136  {
138  << "Incorrect dimensions of phi: " << phi.dimensions()
139  << abort(FatalError);
140  }
141 
142  rDeltaT.correctBoundaryConditions();
143 
144  return trDeltaT;
145 }
146 
147 
148 template<class Type>
149 tmp<VolField<Type>>
151 (
152  const dimensioned<Type>& dt
153 )
154 {
155  const volScalarField rDeltaT(SLrDeltaT());
156 
157  const word ddtName("ddt("+dt.name()+')');
158 
159  if (mesh().moving())
160  {
161  tmp<VolField<Type>> tdtdt
162  (
164  (
165  ddtName,
166  mesh(),
168  (
169  "0",
170  dt.dimensions()/dimTime,
171  Zero
172  )
173  )
174  );
175 
176  tdtdt.ref().primitiveFieldRef() =
177  rDeltaT.primitiveField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
178 
179  return tdtdt;
180  }
181  else
182  {
183  return tmp<VolField<Type>>
184  (
186  (
187  ddtName,
188  mesh(),
190  (
191  "0",
192  dt.dimensions()/dimTime,
193  Zero
194  ),
196  )
197  );
198  }
199 }
200 
201 
202 template<class Type>
205 (
206  const VolField<Type>& vf
207 )
208 {
209  const volScalarField rDeltaT(SLrDeltaT());
210 
211  const word ddtName("ddt("+vf.name()+')');
212 
213  if (mesh().moving())
214  {
215  return tmp<VolField<Type>>
216  (
218  (
219  ddtName,
220  rDeltaT()*(vf() - vf.oldTime()()*mesh().V0()/mesh().V()),
221  rDeltaT.boundaryField()*
222  (
223  vf.boundaryField() - vf.oldTime().boundaryField()
224  )
225  )
226  );
227  }
228  else
229  {
230  return tmp<VolField<Type>>
231  (
233  (
234  ddtName,
235  rDeltaT*(vf - vf.oldTime())
236  )
237  );
238  }
239 }
240 
241 
242 template<class Type>
245 (
246  const dimensionedScalar& rho,
247  const VolField<Type>& vf
248 )
249 {
250  const volScalarField rDeltaT(SLrDeltaT());
251 
252  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
253 
254  if (mesh().moving())
255  {
256  return tmp<VolField<Type>>
257  (
259  (
260  ddtName,
261  rDeltaT()*rho*(vf() - vf.oldTime()()*mesh().V0()/mesh().V()),
262  rDeltaT.boundaryField()*rho.value()*
263  (
264  vf.boundaryField() - vf.oldTime().boundaryField()
265  )
266  )
267  );
268  }
269  else
270  {
271  return tmp<VolField<Type>>
272  (
274  (
275  ddtName,
276  rDeltaT*rho*(vf - vf.oldTime())
277  )
278  );
279  }
280 }
281 
282 
283 template<class Type>
286 (
287  const volScalarField& rho,
288  const VolField<Type>& vf
289 )
290 {
291  const volScalarField rDeltaT(SLrDeltaT());
292 
293  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
294 
295  if (mesh().moving())
296  {
297  return tmp<VolField<Type>>
298  (
300  (
301  ddtName,
302  rDeltaT()*
303  (
304  rho()*vf()
305  - rho.oldTime()()*vf.oldTime()()*mesh().V0()/mesh().V()
306  ),
307  rDeltaT.boundaryField()*
308  (
309  rho.boundaryField()*vf.boundaryField()
310  - rho.oldTime().boundaryField()
311  *vf.oldTime().boundaryField()
312  )
313  )
314  );
315  }
316  else
317  {
318  return tmp<VolField<Type>>
319  (
321  (
322  ddtName,
323  rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
324  )
325  );
326  }
327 }
328 
329 
330 template<class Type>
333 (
334  const volScalarField& alpha,
335  const volScalarField& rho,
336  const VolField<Type>& vf
337 )
338 {
339  const volScalarField rDeltaT(SLrDeltaT());
340 
341  const word ddtName("ddt("+alpha.name()+','+rho.name()+','+vf.name()+')');
342 
343  if (mesh().moving())
344  {
345  return tmp<VolField<Type>>
346  (
348  (
349  ddtName,
350  rDeltaT()*
351  (
352  alpha()*rho()*vf()
353  - alpha.oldTime()()*rho.oldTime()()*vf.oldTime()()
354  *mesh().Vsc0()/mesh().Vsc()
355  ),
356  rDeltaT.boundaryField()*
357  (
358  alpha.boundaryField()
359  *rho.boundaryField()
360  *vf.boundaryField()
361 
362  - alpha.oldTime().boundaryField()
363  *rho.oldTime().boundaryField()
364  *vf.oldTime().boundaryField()
365  )
366  )
367  );
368  }
369  else
370  {
371  return tmp<VolField<Type>>
372  (
374  (
375  ddtName,
376  rDeltaT
377  *(
378  alpha*rho*vf - alpha.oldTime()*rho.oldTime()*vf.oldTime()
379  )
380  )
381  );
382  }
383 }
384 
385 
386 template<class Type>
389 (
390  const VolField<Type>& vf
391 )
392 {
393  tmp<fvMatrix<Type>> tfvm
394  (
395  new fvMatrix<Type>
396  (
397  vf,
399  )
400  );
401 
402  fvMatrix<Type>& fvm = tfvm.ref();
403 
404  const scalarField rDeltaT(SLrDeltaT()().primitiveField());
405 
406  fvm.diag() = rDeltaT*mesh().V();
407 
408  if (mesh().moving())
409  {
410  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().V0();
411  }
412  else
413  {
414  fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().V();
415  }
416 
417  return tfvm;
418 }
419 
420 
421 template<class Type>
424 (
425  const dimensionedScalar& rho,
426  const VolField<Type>& vf
427 )
428 {
429  tmp<fvMatrix<Type>> tfvm
430  (
431  new fvMatrix<Type>
432  (
433  vf,
434  rho.dimensions()*vf.dimensions()*dimVolume/dimTime
435  )
436  );
437  fvMatrix<Type>& fvm = tfvm.ref();
438 
439  const scalarField rDeltaT(SLrDeltaT()().primitiveField());
440 
441  fvm.diag() = rDeltaT*rho.value()*mesh().V();
442 
443  if (mesh().moving())
444  {
445  fvm.source() = rDeltaT
446  *rho.value()*vf.oldTime().primitiveField()*mesh().V0();
447  }
448  else
449  {
450  fvm.source() = rDeltaT
451  *rho.value()*vf.oldTime().primitiveField()*mesh().V();
452  }
453 
454  return tfvm;
455 }
456 
457 
458 template<class Type>
461 (
462  const volScalarField& rho,
463  const VolField<Type>& vf
464 )
465 {
466  tmp<fvMatrix<Type>> tfvm
467  (
468  new fvMatrix<Type>
469  (
470  vf,
471  rho.dimensions()*vf.dimensions()*dimVolume/dimTime
472  )
473  );
474  fvMatrix<Type>& fvm = tfvm.ref();
475 
476  const scalarField rDeltaT(SLrDeltaT()().primitiveField());
477 
478  fvm.diag() = rDeltaT*rho.primitiveField()*mesh().V();
479 
480  if (mesh().moving())
481  {
482  fvm.source() = rDeltaT
483  *rho.oldTime().primitiveField()
484  *vf.oldTime().primitiveField()*mesh().V0();
485  }
486  else
487  {
488  fvm.source() = rDeltaT
489  *rho.oldTime().primitiveField()
490  *vf.oldTime().primitiveField()*mesh().V();
491  }
492 
493  return tfvm;
494 }
495 
496 
497 template<class Type>
500 (
501  const volScalarField& alpha,
502  const volScalarField& rho,
503  const VolField<Type>& vf
504 )
505 {
506  tmp<fvMatrix<Type>> tfvm
507  (
508  new fvMatrix<Type>
509  (
510  vf,
511  alpha.dimensions()
512  *rho.dimensions()
513  *vf.dimensions()
514  *dimVolume
515  /dimTime
516  )
517  );
518  fvMatrix<Type>& fvm = tfvm.ref();
519 
520  const scalarField rDeltaT(SLrDeltaT()().primitiveField());
521 
522  fvm.diag() =
523  rDeltaT*alpha.primitiveField()*rho.primitiveField()*mesh().Vsc();
524 
525  if (mesh().moving())
526  {
527  fvm.source() = rDeltaT
528  *alpha.oldTime().primitiveField()
529  *rho.oldTime().primitiveField()
530  *vf.oldTime().primitiveField()*mesh().Vsc0();
531  }
532  else
533  {
534  fvm.source() = rDeltaT
535  *alpha.oldTime().primitiveField()
536  *rho.oldTime().primitiveField()
537  *vf.oldTime().primitiveField()*mesh().Vsc();
538  }
539 
540  return tfvm;
541 }
542 
543 
544 template<class Type>
547 (
548  const VolField<Type>& U,
549  const SurfaceField<Type>& Uf
550 )
551 {
552  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
553 
554  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
555  fluxFieldType phiCorr
556  (
557  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
558  );
559 
560  return fluxFieldType::New
561  (
562  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
563  this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
564  *rDeltaT*phiCorr
565  );
566 }
567 
568 
569 template<class Type>
572 (
573  const VolField<Type>& U,
574  const fluxFieldType& phi
575 )
576 {
577  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
578 
579  fluxFieldType phiCorr
580  (
581  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
582  );
583 
584  return fluxFieldType::New
585  (
586  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
587  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
588  *rDeltaT*phiCorr
589  );
590 }
591 
592 
593 template<class Type>
596 (
597  const volScalarField& rho,
598  const VolField<Type>& U,
599  const SurfaceField<Type>& rhoUf
600 )
601 {
602  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
603 
604  if
605  (
606  U.dimensions() == dimVelocity
607  && rhoUf.dimensions() == dimDensity*dimVelocity
608  )
609  {
610  VolField<Type> rhoU0
611  (
612  rho.oldTime()*U.oldTime()
613  );
614 
615  fluxFieldType phiUf0(mesh().Sf() & rhoUf.oldTime());
616  fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
617 
618  return fluxFieldType::New
619  (
620  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + rhoUf.name() + ')',
621  this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime())
622  *rDeltaT*phiCorr
623  );
624  }
625  else if
626  (
627  U.dimensions() == dimDensity*dimVelocity
628  && rhoUf.dimensions() == dimDensity*dimVelocity
629  )
630  {
631  fluxFieldType phiUf0(mesh().Sf() & rhoUf.oldTime());
632  fluxFieldType phiCorr
633  (
634  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
635  );
636 
637  return fluxFieldType::New
638  (
639  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + rhoUf.name() + ')',
640  this->fvcDdtPhiCoeff
641  (
642  U.oldTime(),
643  phiUf0,
644  phiCorr,
645  rho.oldTime()
646  )*rDeltaT*phiCorr
647  );
648  }
649  else
650  {
652  << "dimensions of rhoUf are not correct"
653  << abort(FatalError);
654 
655  return fluxFieldType::null();
656  }
657 }
658 
659 
660 template<class Type>
663 (
664  const volScalarField& rho,
665  const VolField<Type>& U,
666  const fluxFieldType& phi
667 )
668 {
669  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
670 
671  if
672  (
673  U.dimensions() == dimVelocity
674  && phi.dimensions() == rho.dimensions()*dimVolumetricFlux
675  )
676  {
677  VolField<Type> rhoU0
678  (
679  rho.oldTime()*U.oldTime()
680  );
681 
682  fluxFieldType phiCorr
683  (
684  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
685  );
686 
687  return fluxFieldType::New
688  (
689  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
690  this->fvcDdtPhiCoeff
691  (
692  rhoU0,
693  phi.oldTime(),
694  phiCorr,
695  rho.oldTime()
696  )*rDeltaT*phiCorr
697  );
698  }
699  else if
700  (
701  U.dimensions() == rho.dimensions()*dimVelocity
702  && phi.dimensions() == rho.dimensions()*dimVolumetricFlux
703  )
704  {
705  fluxFieldType phiCorr
706  (
707  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
708  );
709 
710  return fluxFieldType::New
711  (
712  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
713  this->fvcDdtPhiCoeff
714  (
715  U.oldTime(),
716  phi.oldTime(),
717  phiCorr,
718  rho.oldTime()
719  )*rDeltaT*phiCorr
720  );
721  }
722  else
723  {
725  << "dimensions of phi are not correct"
726  << abort(FatalError);
727 
728  return fluxFieldType::null();
729  }
730 }
731 
732 
733 template<class Type>
735 (
736  const VolField<Type>&
737 )
738 {
740  (
741  "meshPhi",
742  mesh(),
744  );
745 }
746 
747 
748 template<class Type>
750 (
751  const VolField<Type>&,
752  const label patchi
753 )
754 {
755  return tmp<scalarField>
756  (
757  new scalarField(mesh().boundary()[patchi].size(), 0)
758  );
759 }
760 
761 
762 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
763 
764 } // End namespace fv
765 
766 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
767 
768 } // End namespace Foam
769 
770 // ************************************************************************* //
#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.
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< fluxFieldType > fvcDdtPhiCorr(const VolField< Type > &U, const fluxFieldType &phi)
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:334
A special matrix type and solver, designed for finite volume solutions of scalar equations.
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
const dimensionSet dimless
SurfaceField< scalar > surfaceScalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const dimensionSet dimTime
const dimensionSet dimDensity
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
const dimensionSet dimVolume
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:64
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)