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-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 "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<GeometricField<Type, fvPatchField, volMesh>>
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  {
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  {
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 (
207 )
208 {
209  const volScalarField rDeltaT(SLrDeltaT());
210 
211  const word ddtName("ddt("+vf.name()+')');
212 
213  if (mesh().moving())
214  {
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  {
231  (
233  (
234  ddtName,
235  rDeltaT*(vf - vf.oldTime())
236  )
237  );
238  }
239 }
240 
241 
242 template<class Type>
245 (
246  const dimensionedScalar& rho,
248 )
249 {
250  const volScalarField rDeltaT(SLrDeltaT());
251 
252  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
253 
254  if (mesh().moving())
255  {
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  {
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,
289 )
290 {
291  const volScalarField rDeltaT(SLrDeltaT());
292 
293  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
294 
295  if (mesh().moving())
296  {
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  {
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,
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  {
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  {
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 (
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,
427 )
428 {
429  tmp<fvMatrix<Type>> tfvm
430  (
431  new fvMatrix<Type>
432  (
433  vf,
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,
464 )
465 {
466  tmp<fvMatrix<Type>> tfvm
467  (
468  new fvMatrix<Type>
469  (
470  vf,
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,
504 )
505 {
506  tmp<fvMatrix<Type>> tfvm
507  (
508  new fvMatrix<Type>
509  (
510  vf,
511  alpha.dimensions()*rho.dimensions()*vf.dimensions()*dimVol/dimTime
512  )
513  );
514  fvMatrix<Type>& fvm = tfvm.ref();
515 
516  const scalarField rDeltaT(SLrDeltaT()().primitiveField());
517 
518  fvm.diag() =
519  rDeltaT*alpha.primitiveField()*rho.primitiveField()*mesh().Vsc();
520 
521  if (mesh().moving())
522  {
523  fvm.source() = rDeltaT
524  *alpha.oldTime().primitiveField()
525  *rho.oldTime().primitiveField()
526  *vf.oldTime().primitiveField()*mesh().Vsc0();
527  }
528  else
529  {
530  fvm.source() = rDeltaT
531  *alpha.oldTime().primitiveField()
532  *rho.oldTime().primitiveField()
533  *vf.oldTime().primitiveField()*mesh().Vsc();
534  }
535 
536  return tfvm;
537 }
538 
539 
540 template<class Type>
543 (
546 )
547 {
548  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
549 
550  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
551  fluxFieldType phiCorr
552  (
553  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
554  );
555 
556  return fluxFieldType::New
557  (
558  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
559  this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)
560  *rDeltaT*phiCorr
561  );
562 }
563 
564 
565 template<class Type>
568 (
570  const fluxFieldType& phi
571 )
572 {
573  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
574 
575  fluxFieldType phiCorr
576  (
577  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
578  );
579 
580  return fluxFieldType::New
581  (
582  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
583  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
584  *rDeltaT*phiCorr
585  );
586 }
587 
588 
589 template<class Type>
592 (
593  const volScalarField& rho,
596 )
597 {
598  const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
599 
600  if
601  (
602  U.dimensions() == dimVelocity
604  )
605  {
607  (
608  rho.oldTime()*U.oldTime()
609  );
610 
611  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
612  fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0));
613 
614  return fluxFieldType::New
615  (
616  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
617  this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime())
618  *rDeltaT*phiCorr
619  );
620  }
621  else if
622  (
625  )
626  {
627  fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime());
628  fluxFieldType phiCorr
629  (
630  phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
631  );
632 
633  return fluxFieldType::New
634  (
635  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
636  this->fvcDdtPhiCoeff
637  (
638  U.oldTime(),
639  phiUf0,
640  phiCorr,
641  rho.oldTime()
642  )*rDeltaT*phiCorr
643  );
644  }
645  else
646  {
648  << "dimensions of Uf are not correct"
649  << abort(FatalError);
650 
651  return fluxFieldType::null();
652  }
653 }
654 
655 
656 template<class Type>
659 (
660  const volScalarField& rho,
662  const fluxFieldType& phi
663 )
664 {
665  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
666 
667  if
668  (
669  U.dimensions() == dimVelocity
670  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
671  )
672  {
674  (
675  rho.oldTime()*U.oldTime()
676  );
677 
678  fluxFieldType phiCorr
679  (
680  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), rhoU0)
681  );
682 
683  return fluxFieldType::New
684  (
685  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
686  this->fvcDdtPhiCoeff
687  (
688  rhoU0,
689  phi.oldTime(),
690  phiCorr,
691  rho.oldTime()
692  )*rDeltaT*phiCorr
693  );
694  }
695  else if
696  (
697  U.dimensions() == rho.dimensions()*dimVelocity
698  && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
699  )
700  {
701  fluxFieldType phiCorr
702  (
703  phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
704  );
705 
706  return fluxFieldType::New
707  (
708  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
709  this->fvcDdtPhiCoeff
710  (
711  U.oldTime(),
712  phi.oldTime(),
713  phiCorr,
714  rho.oldTime()
715  )*rDeltaT*phiCorr
716  );
717  }
718  else
719  {
721  << "dimensions of phi are not correct"
722  << abort(FatalError);
723 
724  return fluxFieldType::null();
725  }
726 }
727 
728 
729 template<class Type>
731 (
733 )
734 {
736  (
737  "meshPhi",
738  mesh(),
740  );
741 }
742 
743 
744 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
745 
746 } // End namespace fv
747 
748 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
749 
750 } // End namespace Foam
751 
752 // ************************************************************************* //
fvsPatchField< scalar > fvsPatchScalarField
rho oldTime()
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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.
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
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)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual tmp< GeometricField< Type, fvPatchField, volMesh > > fvcDdt(const dimensioned< Type > &)
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
virtual tmp< fluxFieldType > fvcDdtUfCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
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
label patchi
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
virtual tmp< surfaceScalarField > meshPhi(const GeometricField< Type, fvPatchField, volMesh > &)
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.
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
virtual tmp< fluxFieldType > fvcDdtPhiCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const fluxFieldType &phi)
tmp< volScalarField > trDeltaT
Definition: createRDeltaT.H:3
virtual tmp< fvMatrix< Type > > fvmDdt(const GeometricField< Type, fvPatchField, volMesh > &)
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Namespace for OpenFOAM.
const dimensionSet dimVelocity