EulerFaD2dt2Scheme.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2017 Volkswagen AG
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "EulerFaD2dt2Scheme.H"
29 #include "faMatrices.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 namespace fa
39 {
40 
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 
43 template<class Type>
44 scalar EulerFaD2dt2Scheme<Type>::deltaT_() const
45 {
46  return mesh().time().deltaTValue();
47 }
48 
49 
50 template<class Type>
51 scalar EulerFaD2dt2Scheme<Type>::deltaT0_() const
52 {
53  return mesh().time().deltaT0Value();
54 }
55 
56 
57 template<class Type>
58 tmp<GeometricField<Type, faPatchField, areaMesh>>
59 EulerFaD2dt2Scheme<Type>::facD2dt2
60 (
61  const dimensioned<Type> dt
62 )
63 {
64 
65  scalar deltaT = deltaT_();
66  scalar deltaT0 = deltaT0_();
67 
68  dimensionedScalar rDeltaT2 =
69  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
70 
71  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
72  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
73 
74  const IOobject d2dt2IOobject
75  (
76  mesh().thisDb().newIOobject
77  (
78  "d2dt2("+dt.name()+')',
79  { IOobject::REGISTER }
80  )
81  );
82 
83  if (mesh().moving())
84  {
85  scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
86 
87  scalarField SS0(mesh().S() + mesh().S0());
88  scalarField S0S00(mesh().S0() + mesh().S00());
89 
90  tmp<GeometricField<Type, faPatchField, areaMesh>> tdt2dt2
91  (
92  new GeometricField<Type, faPatchField, areaMesh>
93  (
94  d2dt2IOobject,
95  mesh(),
96  dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero)
97  )
98  );
99 
100  tdt2dt2.ref().primitiveFieldRef() =
101  halfRdeltaT2*dt.value()
102  *(coefft*SS0 - (coefft*SS0 + coefft00*S0S00) + coefft00*S0S00)
103  /mesh().S();
104 
105  return tdt2dt2;
106  }
107  else
108  {
109  return tmp<GeometricField<Type, faPatchField, areaMesh>>
110  (
111  new GeometricField<Type, faPatchField, areaMesh>
112  (
113  d2dt2IOobject,
114  mesh(),
115  dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero),
116  faPatchFieldBase::calculatedType()
117  )
118  );
119  }
120 }
121 
122 
123 template<class Type>
124 tmp<GeometricField<Type, faPatchField, areaMesh>>
125 EulerFaD2dt2Scheme<Type>::facD2dt2
126 (
127  const GeometricField<Type, faPatchField, areaMesh>& vf
128 )
129 {
130  dimensionedScalar rDeltaT2 =
131  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
132 
133  const IOobject d2dt2IOobject
134  (
135  mesh().thisDb().newIOobject
136  (
137  "d2dt2("+vf.name()+')',
138  { IOobject::REGISTER }
139  )
140  );
141 
142  scalar deltaT = deltaT_();
143  scalar deltaT0 = deltaT0_();
144 
145  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
146  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
147  scalar coefft0 = coefft + coefft00;
148 
149  if (mesh().moving())
150  {
151  scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
152 
153  scalarField SS0(mesh().S() + mesh().S0());
154  scalarField S0S00(mesh().S0() + mesh().S00());
155 
156  return tmp<GeometricField<Type, faPatchField, areaMesh>>
157  (
158  new GeometricField<Type, faPatchField, areaMesh>
159  (
160  d2dt2IOobject,
161  mesh(),
162  rDeltaT2.dimensions()*vf.dimensions(),
163  halfRdeltaT2*
164  (
165  coefft*SS0*vf.primitiveField()
166 
167  - (coefft*SS0 + coefft00*S0S00)
168  *vf.oldTime().primitiveField()
169 
170  + (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
171  )/mesh().S(),
172  rDeltaT2.value()*
173  (
174  coefft*vf.boundaryField()
175  - coefft0*vf.oldTime().boundaryField()
176  + coefft00*vf.oldTime().oldTime().boundaryField()
177  )
178  )
179  );
180  }
181  else
182  {
183  return tmp<GeometricField<Type, faPatchField, areaMesh>>
184  (
185  new GeometricField<Type, faPatchField, areaMesh>
186  (
187  d2dt2IOobject,
188  rDeltaT2*
189  (
190  coefft*vf
191  - coefft0*vf.oldTime()
192  + coefft00*vf.oldTime().oldTime()
193  )
194  )
195  );
196  }
197 }
198 
199 
200 template<class Type>
201 tmp<GeometricField<Type, faPatchField, areaMesh>>
202 EulerFaD2dt2Scheme<Type>::facD2dt2
203 (
204  const dimensionedScalar& rho,
205  const GeometricField<Type, faPatchField, areaMesh>& vf
206 )
207 {
208  dimensionedScalar rDeltaT2 =
209  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
210 
211  const IOobject d2dt2IOobject
212  (
213  mesh().thisDb().newIOobject
214  (
215  "d2dt2("+rho.name()+','+vf.name()+')',
216  { IOobject::REGISTER }
217  )
218  );
219 
220  scalar deltaT = deltaT_();
221  scalar deltaT0 = deltaT0_();
222 
223  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
224  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
225  scalar coefft0 = coefft + coefft00;
226 
227  if (mesh().moving())
228  {
229  scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
230 
231  scalarField SS0rhoRho0
232  (
233  (mesh().S() + mesh().S0())*rho.value()
234  );
235 
236  scalarField S0S00rho0Rho00
237  (
238  (mesh().S0() + mesh().S00())*rho.value()
239  );
240 
241  return tmp<GeometricField<Type, faPatchField, areaMesh>>
242  (
243  new GeometricField<Type, faPatchField, areaMesh>
244  (
245  d2dt2IOobject,
246  mesh(),
247  rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
248  halfRdeltaT2*
249  (
250  coefft*SS0rhoRho0*vf.primitiveField()
251 
252  - (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
253  *vf.oldTime().primitiveField()
254 
255  + (coefft00*S0S00rho0Rho00)
256  *vf.oldTime().oldTime().primitiveField()
257  )/mesh().S(),
258  rDeltaT2.value()*rho.value()*
259  (
260  coefft*vf.boundaryField()
261  - coefft0*vf.oldTime().boundaryField()
262  + coefft00*vf.oldTime().oldTime().boundaryField()
263  )
264  )
265  );
266  }
267  else
268  {
269 
270  return tmp<GeometricField<Type, faPatchField, areaMesh>>
271  (
272  new GeometricField<Type, faPatchField, areaMesh>
273  (
274  d2dt2IOobject,
275  rDeltaT2 * rho.value() *
276  (
277  coefft*vf
278  - coefft0*vf.oldTime()
279  + coefft00*vf.oldTime().oldTime()
280  )
281  )
282  );
283  }
284 }
285 
286 
287 template<class Type>
288 tmp<GeometricField<Type, faPatchField, areaMesh>>
289 EulerFaD2dt2Scheme<Type>::facD2dt2
290 (
291  const areaScalarField& rho,
292  const GeometricField<Type, faPatchField, areaMesh>& vf
293 )
294 {
295  dimensionedScalar rDeltaT2 =
296  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
297 
298  const IOobject d2dt2IOobject
299  (
300  mesh().thisDb().newIOobject
301  (
302  "d2dt2("+rho.name()+','+vf.name()+')',
303  { IOobject::REGISTER }
304  )
305  );
306 
307  scalar deltaT = deltaT_();
308  scalar deltaT0 = deltaT0_();
309 
310  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
311  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
312 
313  if (mesh().moving())
314  {
315  scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
316  scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();
317 
318  scalarField SS0rhoRho0
319  (
320  (mesh().S() + mesh().S0())
321  *(rho.primitiveField() + rho.oldTime().primitiveField())
322  );
323 
324  scalarField S0S00rho0Rho00
325  (
326  (mesh().S0() + mesh().S00())
327  *(
328  rho.oldTime().primitiveField()
329  + rho.oldTime().oldTime().primitiveField()
330  )
331  );
332 
333  return tmp<GeometricField<Type, faPatchField, areaMesh>>
334  (
335  new GeometricField<Type, faPatchField, areaMesh>
336  (
337  d2dt2IOobject,
338  mesh(),
339  rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
340  quarterRdeltaT2*
341  (
342  coefft*SS0rhoRho0*vf.primitiveField()
343 
344  - (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
345  *vf.oldTime().primitiveField()
346 
347  + (coefft00*S0S00rho0Rho00)
348  *vf.oldTime().oldTime().primitiveField()
349  )/mesh().S(),
350  halfRdeltaT2*
351  (
352  coefft
353  *(rho.boundaryField() + rho.oldTime().boundaryField())
354  *vf.boundaryField()
355 
356  - (
357  coefft
358  *(
359  rho.boundaryField()
360  + rho.oldTime().boundaryField()
361  )
362  + coefft00
363  *(
364  rho.oldTime().boundaryField()
365  + rho.oldTime().oldTime().boundaryField()
366  )
367  )*vf.oldTime().boundaryField()
368 
369  + coefft00
370  *(
371  rho.oldTime().boundaryField()
372  + rho.oldTime().oldTime().boundaryField()
373  )*vf.oldTime().oldTime().boundaryField()
374  )
375  )
376  );
377  }
378  else
379  {
380  dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
381 
382  areaScalarField rhoRho0(rho + rho.oldTime());
383  areaScalarField rho0Rho00(rho.oldTime() + rho.oldTime().oldTime());
384 
385  return tmp<GeometricField<Type, faPatchField, areaMesh>>
386  (
387  new GeometricField<Type, faPatchField, areaMesh>
388  (
389  d2dt2IOobject,
390  halfRdeltaT2*
391  (
392  coefft*rhoRho0*vf
393  - (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
394  + coefft00*rho0Rho00*vf.oldTime().oldTime()
395  )
396  )
397  );
398  }
399 }
400 
401 
402 template<class Type>
403 tmp<faMatrix<Type>>
404 EulerFaD2dt2Scheme<Type>::famD2dt2
405 (
406  const GeometricField<Type, faPatchField, areaMesh>& vf
407 )
408 {
409  tmp<faMatrix<Type>> tfam
410  (
411  new faMatrix<Type>
412  (
413  vf,
414  vf.dimensions()*dimArea/dimTime/dimTime
415  )
416  );
417 
418  faMatrix<Type>& fam = tfam.ref();
419 
420  scalar deltaT = deltaT_();
421  scalar deltaT0 = deltaT0_();
422 
423  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
424  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
425  scalar coefft0 = coefft + coefft00;
426 
427  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
428 
429  if (mesh().moving())
430  {
431  scalar halfRdeltaT2 = rDeltaT2/2.0;
432 
433  scalarField SS0(mesh().S() + mesh().S0());
434  scalarField S0S00(mesh().S0() + mesh().S00());
435 
436  fam.diag() = (coefft*halfRdeltaT2)*SS0;
437 
438  fam.source() = halfRdeltaT2*
439  (
440  (coefft*SS0 + coefft00*S0S00)
441  *vf.oldTime().primitiveField()
442 
443  - (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
444  );
445  }
446  else
447  {
448  fam.diag() = (coefft*rDeltaT2)*mesh().S();
449 
450  fam.source() = rDeltaT2*mesh().S()*
451  (
452  coefft0*vf.oldTime().primitiveField()
453  - coefft00*vf.oldTime().oldTime().primitiveField()
454  );
455  }
456 
457  return tfam;
458 }
459 
460 
461 template<class Type>
462 tmp<faMatrix<Type>>
463 EulerFaD2dt2Scheme<Type>::famD2dt2
464 (
465  const dimensionedScalar& rho,
466  const GeometricField<Type, faPatchField, areaMesh>& vf
467 )
468 {
469  tmp<faMatrix<Type>> tfam
470  (
471  new faMatrix<Type>
472  (
473  vf,
474  rho.dimensions()*vf.dimensions()*dimArea
475  /dimTime/dimTime
476  )
477  );
478 
479  faMatrix<Type>& fam = tfam.ref();
480 
481  scalar deltaT = deltaT_();
482  scalar deltaT0 = deltaT0_();
483 
484  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
485  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
486 
487  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
488 
489  if (mesh().moving())
490  {
491  scalar halfRdeltaT2 = 0.5*rDeltaT2;
492 
493  scalarField SS0(mesh().S() + mesh().S0());
494  scalarField S0S00(mesh().S0() + mesh().S00());
495 
496  fam.diag() = rho.value()*(coefft*halfRdeltaT2)*SS0;
497 
498  fam.source() = halfRdeltaT2*rho.value()*
499  (
500  (coefft*SS0 + coefft00*S0S00)
501  *vf.oldTime().primitiveField()
502 
503  - (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
504  );
505  }
506  else
507  {
508  fam.diag() = (coefft*rDeltaT2)*mesh().S()*rho.value();
509 
510  fam.source() = rDeltaT2*mesh().S()*rho.value()*
511  (
512  (coefft + coefft00)*vf.oldTime().primitiveField()
513  - coefft00*vf.oldTime().oldTime().primitiveField()
514  );
515  }
516 
517  return tfam;
518 }
519 
520 
521 template<class Type>
522 tmp<faMatrix<Type>>
523 EulerFaD2dt2Scheme<Type>::famD2dt2
524 (
525  const areaScalarField& rho,
526  const GeometricField<Type, faPatchField, areaMesh>& vf
527 )
528 {
529  tmp<faMatrix<Type>> tfam
530  (
531  new faMatrix<Type>
532  (
533  vf,
534  rho.dimensions()*vf.dimensions()*dimArea/dimTime/dimTime
535  )
536  );
537  faMatrix<Type>& fam = tfam.ref();
538 
539  scalar deltaT =deltaT_();
540  scalar deltaT0 = deltaT0_();
541 
542  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
543  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
544 
545  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
546 
547  if (mesh().moving())
548  {
549  scalar quarterRdeltaT2 = 0.25*rDeltaT2;
550 
551  scalarField SS0rhoRho0
552  (
553  (mesh().S() + mesh().S0())
554  *(rho.primitiveField() + rho.oldTime().primitiveField())
555  );
556 
557  scalarField S0S00rho0Rho00
558  (
559  (mesh().S0() + mesh().S00())
560  *(
561  rho.oldTime().primitiveField()
562  + rho.oldTime().oldTime().primitiveField()
563  )
564  );
565 
566  fam.diag() = (coefft*quarterRdeltaT2)*SS0rhoRho0;
567 
568  fam.source() = quarterRdeltaT2*
569  (
570  (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
571  *vf.oldTime().primitiveField()
572 
573  - (coefft00*S0S00rho0Rho00)
574  *vf.oldTime().oldTime().primitiveField()
575  );
576  }
577  else
578  {
579  scalar halfRdeltaT2 = 0.5*rDeltaT2;
580 
581  scalarField rhoRho0
582  (
583  rho.primitiveField() + rho.oldTime().primitiveField()
584  );
585 
586  scalarField rho0Rho00
587  (
588  rho.oldTime().primitiveField()
589  + rho.oldTime().oldTime().primitiveField()
590  );
591 
592  fam.diag() = (coefft*halfRdeltaT2)*mesh().S()*rhoRho0;
593 
594  fam.source() = halfRdeltaT2*mesh().S()*
595  (
596  (coefft*rhoRho0 + coefft00*rho0Rho00)
597  *vf.oldTime().primitiveField()
598 
599  - (coefft00*rho0Rho00)
600  *vf.oldTime().oldTime().primitiveField()
601  );
602  }
603 
604  return tfam;
605 }
606 
607 
608 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
609 
610 } // End namespace fa
611 
612 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
613 
614 } // End namespace Foam
615 
616 // ************************************************************************* //
Foam::dimensioned::value
const Type & value() const noexcept
Return const reference to value.
Definition: dimensionedType.H:366
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:938
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const noexcept
Return a const-reference to the internal field values.
Definition: GeometricFieldI.H:36
Foam::IOobject::name
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:195
Foam::tmp::ref
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
Definition: tmpI.H:235
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:44
Foam::GeometricField
Generic GeometricField class.
Definition: pointMVCWeight.H:58
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedScalarFwd.H:33
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:46
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const noexcept
Return const reference to dimensions.
Definition: dimensionedType.H:356
Foam::faPatchFieldBase::calculatedType
static const word & calculatedType() noexcept
The type name for calculated patch fields.
Definition: faPatchField.H:174
Foam::fa::EulerFaD2dt2Scheme::famD2dt2
tmp< faMatrix< Type > > famD2dt2(const GeometricField< Type, faPatchField, areaMesh > &)
Definition: EulerFaD2dt2Scheme.C:398
Foam::dimensioned::name
const word & name() const noexcept
Return const reference to name.
Definition: dimensionedType.H:346
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:33
Foam::fa::EulerFaD2dt2Scheme::facD2dt2
tmp< GeometricField< Type, faPatchField, areaMesh > > facD2dt2(const dimensioned< Type >)
Definition: EulerFaD2dt2Scheme.C:53
Foam::faMatrix
A special matrix type and solver, designed for finite area solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: faMatricesFwd.H:37
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
Foam::tmp
A class for managing temporary objects.
Definition: HashPtrTable.H:50
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported...
Definition: IOobject.H:180
fam
Calculate the matrix for the second temporal derivative.
Foam::IOobjectOption::REGISTER
Request registration (bool: true)
Definition: IOobjectOption.H:102
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const noexcept
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:45
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:26
Foam::DimensionedField::dimensions
const dimensionSet & dimensions() const noexcept
Return dimensions.
Definition: DimensionedFieldI.H:34
Foam::areaScalarField
GeometricField< scalar, faPatchField, areaMesh > areaScalarField
Definition: areaFieldsFwd.H:72
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:127