advectiveFvPatchField.C

Go to the documentation of this file.

00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "advectiveFvPatchField.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 #include <finiteVolume/fvPatchFieldMapper.H>
00029 #include <finiteVolume/volFields.H>
00030 #include <finiteVolume/EulerDdtScheme.H>
00031 #include <finiteVolume/CrankNicholsonDdtScheme.H>
00032 #include <finiteVolume/backwardDdtScheme.H>
00033 
00034 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00035 
00036 namespace Foam
00037 {
00038 
00039 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00040 
00041 template<class Type>
00042 advectiveFvPatchField<Type>::advectiveFvPatchField
00043 (
00044     const fvPatch& p,
00045     const DimensionedField<Type, volMesh>& iF
00046 )
00047 :
00048     mixedFvPatchField<Type>(p, iF),
00049     phiName_("phi"),
00050     rhoName_("rho"),
00051     fieldInf_(pTraits<Type>::zero),
00052     lInf_(0.0)
00053 {
00054     this->refValue() = pTraits<Type>::zero;
00055     this->refGrad() = pTraits<Type>::zero;
00056     this->valueFraction() = 0.0;
00057 }
00058 
00059 
00060 template<class Type>
00061 advectiveFvPatchField<Type>::advectiveFvPatchField
00062 (
00063     const advectiveFvPatchField& ptf,
00064     const fvPatch& p,
00065     const DimensionedField<Type, volMesh>& iF,
00066     const fvPatchFieldMapper& mapper
00067 )
00068 :
00069     mixedFvPatchField<Type>(ptf, p, iF, mapper),
00070     phiName_(ptf.phiName_),
00071     rhoName_(ptf.rhoName_),
00072     fieldInf_(ptf.fieldInf_),
00073     lInf_(ptf.lInf_)
00074 {}
00075 
00076 
00077 template<class Type>
00078 advectiveFvPatchField<Type>::advectiveFvPatchField
00079 (
00080     const fvPatch& p,
00081     const DimensionedField<Type, volMesh>& iF,
00082     const dictionary& dict
00083 )
00084 :
00085     mixedFvPatchField<Type>(p, iF),
00086     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
00087     rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
00088     fieldInf_(pTraits<Type>::zero),
00089     lInf_(0.0)
00090 {
00091     if (dict.found("value"))
00092     {
00093         fvPatchField<Type>::operator=
00094         (
00095             Field<Type>("value", dict, p.size())
00096         );
00097     }
00098     else
00099     {
00100         fvPatchField<Type>::operator=(this->patchInternalField());
00101     }
00102 
00103     this->refValue() = *this;
00104     this->refGrad() = pTraits<Type>::zero;
00105     this->valueFraction() = 0.0;
00106 
00107     if (dict.readIfPresent("lInf", lInf_))
00108     {
00109         dict.lookup("fieldInf") >> fieldInf_;
00110 
00111         if (lInf_ < 0.0)
00112         {
00113             FatalIOErrorIn
00114             (
00115                 "advectiveFvPatchField<Type>::"
00116                 "advectiveFvPatchField"
00117                 "(const fvPatch&, const Field<Type>&, const dictionary&)",
00118                 dict
00119             )   << "unphysical lInf specified (lInf < 0)\n"
00120                 << "    on patch " << this->patch().name()
00121                 << " of field " << this->dimensionedInternalField().name()
00122                 << " in file " << this->dimensionedInternalField().objectPath()
00123                 << exit(FatalIOError);
00124         }
00125     }
00126 }
00127 
00128 
00129 template<class Type>
00130 advectiveFvPatchField<Type>::advectiveFvPatchField
00131 (
00132     const advectiveFvPatchField& ptpsf
00133 )
00134 :
00135     mixedFvPatchField<Type>(ptpsf),
00136     phiName_(ptpsf.phiName_),
00137     rhoName_(ptpsf.rhoName_),
00138     fieldInf_(ptpsf.fieldInf_),
00139     lInf_(ptpsf.lInf_)
00140 {}
00141 
00142 
00143 template<class Type>
00144 advectiveFvPatchField<Type>::advectiveFvPatchField
00145 (
00146     const advectiveFvPatchField& ptpsf,
00147     const DimensionedField<Type, volMesh>& iF
00148 )
00149 :
00150     mixedFvPatchField<Type>(ptpsf, iF),
00151     phiName_(ptpsf.phiName_),
00152     rhoName_(ptpsf.rhoName_),
00153     fieldInf_(ptpsf.fieldInf_),
00154     lInf_(ptpsf.lInf_)
00155 {}
00156 
00157 
00158 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00159 
00160 template<class Type>
00161 tmp<scalarField> advectiveFvPatchField<Type>::advectionSpeed() const
00162 {
00163     const surfaceScalarField& phi =
00164         this->db().objectRegistry::lookupObject<surfaceScalarField>(phiName_);
00165 
00166     fvsPatchField<scalar> phip = this->patch().lookupPatchField
00167     (
00168         phiName_,
00169         reinterpret_cast<const surfaceScalarField*>(0),
00170         reinterpret_cast<const scalar*>(0)
00171     );
00172 
00173     if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
00174     {
00175         const fvPatchScalarField& rhop = this->patch().lookupPatchField
00176         (
00177             rhoName_,
00178             reinterpret_cast<const volScalarField*>(0),
00179             reinterpret_cast<const scalar*>(0)
00180         );
00181 
00182         return phip/(rhop*this->patch().magSf());
00183     }
00184     else
00185     {
00186         return phip/this->patch().magSf();
00187     }
00188 }
00189 
00190 
00191 template<class Type>
00192 void advectiveFvPatchField<Type>::updateCoeffs()
00193 {
00194     if (this->updated())
00195     {
00196         return;
00197     }
00198 
00199     word ddtScheme
00200     (
00201         this->dimensionedInternalField().mesh()
00202        .ddtScheme(this->dimensionedInternalField().name())
00203     );
00204     scalar deltaT = this->db().time().deltaT().value();
00205 
00206     const GeometricField<Type, fvPatchField, volMesh>& field =
00207         this->db().objectRegistry::
        lookupObject<GeometricField<Type, fvPatchField, volMesh> >
00208         (
00209             this->dimensionedInternalField().name()
00210         );
00211 
00212     // Calculate the advection speed of the field wave
00213     // If the wave is incoming set the speed to 0.
00214     scalarField w = Foam::max(advectionSpeed(), scalar(0));
00215 
00216     // Calculate the field wave coefficient alpha (See notes)
00217     scalarField alpha = w*deltaT*this->patch().deltaCoeffs();
00218 
00219     label patchi = this->patch().index();
00220 
00221     // Non-reflecting outflow boundary
00222     // If lInf_ defined setup relaxation to the value fieldInf_.
00223     if (lInf_ > SMALL)
00224     {
00225         // Calculate the field relaxation coefficient k (See notes)
00226         scalarField k = w*deltaT/lInf_;
00227 
00228         if
00229         (
00230             ddtScheme == fv::EulerDdtScheme<scalar>::typeName
00231          || ddtScheme == fv::CrankNicholsonDdtScheme<scalar>::typeName
00232         )
00233         {
00234             this->refValue() =
00235             (
00236                 field.oldTime().boundaryField()[patchi] + k*fieldInf_
00237             )/(1.0 + k);
00238 
00239             this->valueFraction() = (1.0 + k)/(1.0 + alpha + k);
00240         }
00241         else if (ddtScheme == fv::backwardDdtScheme<scalar>::typeName)
00242         {
00243             this->refValue() =
00244             (
00245                 2.0*field.oldTime().boundaryField()[patchi]
00246               - 0.5*field.oldTime().oldTime().boundaryField()[patchi]
00247               + k*fieldInf_
00248             )/(1.5 + k);
00249 
00250             this->valueFraction() = (1.5 + k)/(1.5 + alpha + k);
00251         }
00252         else
00253         {
00254             FatalErrorIn
00255             (
00256                 "advectiveFvPatchField<Type>::updateCoeffs()"
00257             )   << "    Unsupported temporal differencing scheme : "
00258                 << ddtScheme
00259                 << "\n    on patch " << this->patch().name()
00260                 << " of field " << this->dimensionedInternalField().name()
00261                 << " in file " << this->dimensionedInternalField().objectPath()
00262                 << exit(FatalError);
00263         }
00264     }
00265     else
00266     {
00267         if
00268         (
00269             ddtScheme == fv::EulerDdtScheme<scalar>::typeName
00270          || ddtScheme == fv::CrankNicholsonDdtScheme<scalar>::typeName
00271         )
00272         {
00273             this->refValue() = field.oldTime().boundaryField()[patchi];
00274 
00275             this->valueFraction() = 1.0/(1.0 + alpha);
00276         }
00277         else if (ddtScheme == fv::backwardDdtScheme<scalar>::typeName)
00278         {
00279             this->refValue() =
00280             (
00281                 2.0*field.oldTime().boundaryField()[patchi]
00282               - 0.5*field.oldTime().oldTime().boundaryField()[patchi]
00283             )/1.5;
00284 
00285             this->valueFraction() = 1.5/(1.5 + alpha);
00286         }
00287         else
00288         {
00289             FatalErrorIn
00290             (
00291                 "advectiveFvPatchField<Type>::updateCoeffs()"
00292             )   << "    Unsupported temporal differencing scheme : "
00293                 << ddtScheme
00294                 << "\n    on patch " << this->patch().name()
00295                 << " of field " << this->dimensionedInternalField().name()
00296                 << " in file " << this->dimensionedInternalField().objectPath()
00297                 << exit(FatalError);
00298         }
00299     }
00300 
00301     mixedFvPatchField<Type>::updateCoeffs();
00302 }
00303 
00304 
00305 template<class Type>
00306 void advectiveFvPatchField<Type>::write(Ostream& os) const
00307 {
00308     fvPatchField<Type>::write(os);
00309 
00310     if (phiName_ != "phi")
00311     {
00312         os.writeKeyword("phi") << phiName_ << token::END_STATEMENT << nl;
00313     }
00314     if (rhoName_ != "rho")
00315     {
00316         os.writeKeyword("rho") << rhoName_ << token::END_STATEMENT << nl;
00317     }
00318 
00319     if (lInf_ > SMALL)
00320     {
00321         os.writeKeyword("fieldInf") << fieldInf_
00322             << token::END_STATEMENT << nl;
00323         os.writeKeyword("lInf") << lInf_
00324             << token::END_STATEMENT << nl;
00325     }
00326 
00327     this->writeEntry("value", os);
00328 }
00329 
00330 
00331 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00332 
00333 } // End namespace Foam
00334 
00335 // ************************ vim: set sw=4 sts=4 et: ************************ //
00336