1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "CoEulerDdtScheme.H"
28 #include "surfaceInterpolate.H"
30 #include "fvMatrices.H"
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
45 tmp<volScalarField> CoEulerDdtScheme<Type>::CorDeltaT() const
47 surfaceScalarField cofrDeltaT = CofrDeltaT();
49 tmp<volScalarField> tcorDeltaT
56 cofrDeltaT.instance(),
60 dimensionedScalar("CorDeltaT", cofrDeltaT.dimensions(), 0.0),
61 zeroGradientFvPatchScalarField::typeName
65 volScalarField& corDeltaT = tcorDeltaT();
67 const unallocLabelList& owner = mesh().owner();
68 const unallocLabelList& neighbour = mesh().neighbour();
72 corDeltaT[owner[faceI]] =
73 max(corDeltaT[owner[faceI]], cofrDeltaT[faceI]);
75 corDeltaT[neighbour[faceI]] =
76 max(corDeltaT[neighbour[faceI]], cofrDeltaT[faceI]);
79 forAll(corDeltaT.boundaryField(), patchi)
81 const fvsPatchScalarField& pcofrDeltaT =
82 cofrDeltaT.boundaryField()[patchi];
84 const fvPatch& p = pcofrDeltaT.patch();
85 const unallocLabelList& faceCells = p.patch().faceCells();
87 forAll(pcofrDeltaT, patchFacei)
89 corDeltaT[faceCells[patchFacei]] = max
91 corDeltaT[faceCells[patchFacei]],
92 pcofrDeltaT[patchFacei]
97 corDeltaT.correctBoundaryConditions();
99 //corDeltaT = max(corDeltaT, max(corDeltaT)/100.0);
106 tmp<surfaceScalarField> CoEulerDdtScheme<Type>::CofrDeltaT() const
108 const dimensionedScalar& deltaT = mesh().time().deltaT();
110 const surfaceScalarField& phi =
111 static_cast<const objectRegistry&>(mesh())
112 .lookupObject<surfaceScalarField>(phiName_);
114 if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
116 surfaceScalarField Co
118 mesh().surfaceInterpolation::deltaCoeffs()
119 *(mag(phi)/mesh().magSf())
123 return max(Co/maxCo_, scalar(1))/deltaT;
125 else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
127 const volScalarField& rho =
128 static_cast<const objectRegistry&>(mesh())
129 .lookupObject<volScalarField>(rhoName_).oldTime();
131 surfaceScalarField Co
133 mesh().surfaceInterpolation::deltaCoeffs()
134 *(mag(phi)/(fvc::interpolate(rho)*mesh().magSf()))
138 return max(Co/maxCo_, scalar(1))/deltaT;
142 FatalErrorIn("CoEulerDdtScheme<Type>::CofrDeltaT() const")
143 << "Incorrect dimensions of phi: " << phi.dimensions()
144 << abort(FatalError);
146 return tmp<surfaceScalarField>(NULL);
152 tmp<GeometricField<Type, fvPatchField, volMesh> >
153 CoEulerDdtScheme<Type>::fvcDdt
155 const dimensioned<Type>& dt
158 volScalarField rDeltaT = CorDeltaT();
162 "ddt("+dt.name()+')',
163 mesh().time().timeName(),
169 tmp<GeometricField<Type, fvPatchField, volMesh> > tdtdt
171 new GeometricField<Type, fvPatchField, volMesh>
178 dt.dimensions()/dimTime,
184 tdtdt().internalField() =
185 rDeltaT.internalField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
191 return tmp<GeometricField<Type, fvPatchField, volMesh> >
193 new GeometricField<Type, fvPatchField, volMesh>
200 dt.dimensions()/dimTime,
203 calculatedFvPatchField<Type>::typeName
211 tmp<GeometricField<Type, fvPatchField, volMesh> >
212 CoEulerDdtScheme<Type>::fvcDdt
214 const GeometricField<Type, fvPatchField, volMesh>& vf
217 volScalarField rDeltaT = CorDeltaT();
221 "ddt("+vf.name()+')',
222 mesh().time().timeName(),
228 return tmp<GeometricField<Type, fvPatchField, volMesh> >
230 new GeometricField<Type, fvPatchField, volMesh>
234 rDeltaT.dimensions()*vf.dimensions(),
235 rDeltaT.internalField()*
238 - vf.oldTime().internalField()*mesh().V0()/mesh().V()
240 rDeltaT.boundaryField()*
242 vf.boundaryField() - vf.oldTime().boundaryField()
249 return tmp<GeometricField<Type, fvPatchField, volMesh> >
251 new GeometricField<Type, fvPatchField, volMesh>
254 rDeltaT*(vf - vf.oldTime())
262 tmp<GeometricField<Type, fvPatchField, volMesh> >
263 CoEulerDdtScheme<Type>::fvcDdt
265 const dimensionedScalar& rho,
266 const GeometricField<Type, fvPatchField, volMesh>& vf
269 volScalarField rDeltaT = CorDeltaT();
273 "ddt("+rho.name()+','+vf.name()+')',
274 mesh().time().timeName(),
280 return tmp<GeometricField<Type, fvPatchField, volMesh> >
282 new GeometricField<Type, fvPatchField, volMesh>
286 rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
287 rDeltaT.internalField()*rho.value()*
290 - vf.oldTime().internalField()*mesh().V0()/mesh().V()
292 rDeltaT.boundaryField()*rho.value()*
294 vf.boundaryField() - vf.oldTime().boundaryField()
301 return tmp<GeometricField<Type, fvPatchField, volMesh> >
303 new GeometricField<Type, fvPatchField, volMesh>
306 rDeltaT*rho*(vf - vf.oldTime())
314 tmp<GeometricField<Type, fvPatchField, volMesh> >
315 CoEulerDdtScheme<Type>::fvcDdt
317 const volScalarField& rho,
318 const GeometricField<Type, fvPatchField, volMesh>& vf
321 volScalarField rDeltaT = CorDeltaT();
325 "ddt("+rho.name()+','+vf.name()+')',
326 mesh().time().timeName(),
332 return tmp<GeometricField<Type, fvPatchField, volMesh> >
334 new GeometricField<Type, fvPatchField, volMesh>
338 rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
339 rDeltaT.internalField()*
341 rho.internalField()*vf.internalField()
342 - rho.oldTime().internalField()
343 *vf.oldTime().internalField()*mesh().V0()/mesh().V()
345 rDeltaT.boundaryField()*
347 rho.boundaryField()*vf.boundaryField()
348 - rho.oldTime().boundaryField()
349 *vf.oldTime().boundaryField()
356 return tmp<GeometricField<Type, fvPatchField, volMesh> >
358 new GeometricField<Type, fvPatchField, volMesh>
361 rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
370 CoEulerDdtScheme<Type>::fvmDdt
372 GeometricField<Type, fvPatchField, volMesh>& vf
375 tmp<fvMatrix<Type> > tfvm
380 vf.dimensions()*dimVol/dimTime
384 fvMatrix<Type>& fvm = tfvm();
386 scalarField rDeltaT = CorDeltaT()().internalField();
388 fvm.diag() = rDeltaT*mesh().V();
392 fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
396 fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V();
405 CoEulerDdtScheme<Type>::fvmDdt
407 const dimensionedScalar& rho,
408 GeometricField<Type, fvPatchField, volMesh>& vf
411 tmp<fvMatrix<Type> > tfvm
416 rho.dimensions()*vf.dimensions()*dimVol/dimTime
419 fvMatrix<Type>& fvm = tfvm();
421 scalarField rDeltaT = CorDeltaT()().internalField();
423 fvm.diag() = rDeltaT*rho.value()*mesh().V();
427 fvm.source() = rDeltaT
428 *rho.value()*vf.oldTime().internalField()*mesh().V0();
432 fvm.source() = rDeltaT
433 *rho.value()*vf.oldTime().internalField()*mesh().V();
442 CoEulerDdtScheme<Type>::fvmDdt
444 const volScalarField& rho,
445 GeometricField<Type, fvPatchField, volMesh>& vf
448 tmp<fvMatrix<Type> > tfvm
453 rho.dimensions()*vf.dimensions()*dimVol/dimTime
456 fvMatrix<Type>& fvm = tfvm();
458 scalarField rDeltaT = CorDeltaT()().internalField();
460 fvm.diag() = rDeltaT*rho.internalField()*mesh().V();
464 fvm.source() = rDeltaT
465 *rho.oldTime().internalField()
466 *vf.oldTime().internalField()*mesh().V0();
470 fvm.source() = rDeltaT
471 *rho.oldTime().internalField()
472 *vf.oldTime().internalField()*mesh().V();
480 tmp<typename CoEulerDdtScheme<Type>::fluxFieldType>
481 CoEulerDdtScheme<Type>::fvcDdtPhiCorr
483 const volScalarField& rA,
484 const GeometricField<Type, fvPatchField, volMesh>& U,
485 const fluxFieldType& phi
490 "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
491 mesh().time().timeName(),
497 return tmp<fluxFieldType>
503 dimensioned<typename flux<Type>::type>
506 rA.dimensions()*phi.dimensions()/dimTime,
507 pTraits<typename flux<Type>::type>::zero
514 volScalarField rDeltaT = CorDeltaT();
516 return tmp<fluxFieldType>
521 fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
523 fvc::interpolate(rDeltaT*rA)*phi.oldTime()
524 - (fvc::interpolate(rDeltaT*rA*U.oldTime()) & mesh().Sf())
533 tmp<typename CoEulerDdtScheme<Type>::fluxFieldType>
534 CoEulerDdtScheme<Type>::fvcDdtPhiCorr
536 const volScalarField& rA,
537 const volScalarField& rho,
538 const GeometricField<Type, fvPatchField, volMesh>& U,
539 const fluxFieldType& phi
545 + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
546 mesh().time().timeName(),
552 return tmp<fluxFieldType>
558 dimensioned<typename flux<Type>::type>
561 rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
562 pTraits<typename flux<Type>::type>::zero
569 volScalarField rDeltaT = CorDeltaT();
573 U.dimensions() == dimVelocity
574 && phi.dimensions() == dimVelocity*dimArea
577 return tmp<fluxFieldType>
582 fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
584 fvc::interpolate(rDeltaT*rA*rho.oldTime())*phi.oldTime()
585 - (fvc::interpolate(rDeltaT*rA*rho.oldTime()*U.oldTime())
593 U.dimensions() == dimVelocity
594 && phi.dimensions() == dimDensity*dimVelocity*dimArea
597 return tmp<fluxFieldType>
605 phi.oldTime()/fvc::interpolate(rho.oldTime())
608 fvc::interpolate(rDeltaT*rA*rho.oldTime())
609 *phi.oldTime()/fvc::interpolate(rho.oldTime())
613 rDeltaT*rA*rho.oldTime()*U.oldTime()
622 U.dimensions() == dimDensity*dimVelocity
623 && phi.dimensions() == dimDensity*dimVelocity*dimArea
626 return tmp<fluxFieldType>
631 fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())
633 fvc::interpolate(rDeltaT*rA)*phi.oldTime()
635 fvc::interpolate(rDeltaT*rA*U.oldTime())&mesh().Sf()
645 "CoEulerDdtScheme<Type>::fvcDdtPhiCorr"
646 ) << "dimensions of phi are not correct"
647 << abort(FatalError);
649 return fluxFieldType::null();
656 tmp<surfaceScalarField> CoEulerDdtScheme<Type>::meshPhi
658 const GeometricField<Type, fvPatchField, volMesh>&
661 return tmp<surfaceScalarField>
663 new surfaceScalarField
668 mesh().time().timeName(),
672 dimensionedScalar("0", dimVolume/dimTime, 0.0)
678 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
680 } // End namespace fv
682 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
684 } // End namespace Foam
686 // ************************************************************************* //