Merge branch 'master' of ssh://hunt//home/noisy3/OpenFOAM/OpenFOAM-dev

This commit is contained in:
graham
2011-03-03 11:19:57 +00:00
116 changed files with 278 additions and 221 deletions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -86,18 +86,25 @@ void Foam::Time::adjustDeltaT()
(outputTimeIndex_ + 1)*writeInterval_ - (value() - startTime_)
);
label nStepsToNextWrite = label(timeToNextWrite/deltaT_ - SMALL) + 1;
scalar newDeltaT = timeToNextWrite/nStepsToNextWrite;
scalar nSteps = timeToNextWrite/deltaT_ - SMALL;
// Control the increase of the time step to within a factor of 2
// and the decrease within a factor of 5.
if (newDeltaT >= deltaT_)
// For tiny deltaT the label can overflow!
if (nSteps < labelMax)
{
deltaT_ = min(newDeltaT, 2.0*deltaT_);
}
else
{
deltaT_ = max(newDeltaT, 0.2*deltaT_);
label nStepsToNextWrite = label(nSteps) + 1;
scalar newDeltaT = timeToNextWrite/nStepsToNextWrite;
// Control the increase of the time step to within a factor of 2
// and the decrease within a factor of 5.
if (newDeltaT >= deltaT_)
{
deltaT_ = min(newDeltaT, 2.0*deltaT_);
}
else
{
deltaT_ = max(newDeltaT, 0.2*deltaT_);
}
}
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -30,9 +30,24 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::coupledPolyPatch, 0);
namespace Foam
{
defineTypeNameAndDebug(coupledPolyPatch, 0);
Foam::scalar Foam::coupledPolyPatch::matchTol = 1E-3;
scalar coupledPolyPatch::matchTol = 1E-3;
template<>
const char* NamedEnum<coupledPolyPatch::transformType, 4>::names[] =
{
"unknown",
"rotational",
"translational",
"noOrdering"
};
const NamedEnum<coupledPolyPatch::transformType, 4>
coupledPolyPatch::transformTypeNames;
}
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
@ -204,12 +219,14 @@ void Foam::coupledPolyPatch::calcTransformTensors
const vectorField& nf,
const vectorField& nr,
const scalarField& smallDist,
const scalar absTol
const scalar absTol,
const transformType transform
) const
{
if (debug)
{
Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl
<< " transform:" << transformTypeNames[transform] << nl
<< " (half)size:" << Cf.size() << nl
<< " absTol:" << absTol << nl
<< " smallDist min:" << min(smallDist) << nl
@ -242,9 +259,16 @@ void Foam::coupledPolyPatch::calcTransformTensors
Pout<< " error:" << error << endl;
}
if (sum(mag(nf & nr)) < Cf.size()-error)
if
(
transform == ROTATIONAL
|| (
transform != TRANSLATIONAL
&& (sum(mag(nf & nr)) < Cf.size()-error)
)
)
{
// Rotation, no separation
// Type is rotation or unknown and normals not aligned
// Assume per-face differing transformation, correct later
@ -284,7 +308,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
}
else
{
// No rotation, possible separation
// Translational or (unknown and normals aligned)
forwardT_.setSize(0);
reverseT_.setSize(0);

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -53,6 +53,20 @@ class coupledPolyPatch
:
public polyPatch
{
public:
enum transformType
{
UNKNOWN, // unspecified; automatic ordering
ROTATIONAL, // rotation along coordinate axis
TRANSLATIONAL, // translation
NOORDERING // unspecified, no automatic ordering
};
static const NamedEnum<transformType, 4> transformTypeNames;
private:
// Private data
//- offset (distance) vector from one side of the couple to the other
@ -82,6 +96,8 @@ protected:
//- Calculate the transformation tensors
// smallDist : matching distance per face
// absTol : absolute error in normal
// if transformType = unknown it first tries rotational, then
// translational transform
void calcTransformTensors
(
const vectorField& Cf,
@ -89,7 +105,8 @@ protected:
const vectorField& nf,
const vectorField& nr,
const scalarField& smallDist,
const scalar absTol = matchTol
const scalar absTol = matchTol,
const transformType = UNKNOWN
) const;
//- Initialise the calculation of the patch geometry

View File

@ -44,24 +44,8 @@ namespace Foam
addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, word);
addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, dictionary);
template<>
const char* Foam::NamedEnum
<
Foam::cyclicPolyPatch::transformType,
4
>::names[] =
{
"unknown",
"rotational",
"translational",
"noOrdering"
};
}
const Foam::NamedEnum<Foam::cyclicPolyPatch::transformType, 4>
Foam::cyclicPolyPatch::transformTypeNames;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -271,7 +255,9 @@ void Foam::cyclicPolyPatch::calcTransforms
static_cast<const pointField&>(half1Ctrs),
half0Normals,
half1Normals,
half0Tols
half0Tols,
matchTol,
transform_
);
if (transform_ == ROTATIONAL && !parallel() && forwardT().size() > 1)

View File

@ -64,20 +64,6 @@ class cyclicPolyPatch
:
public coupledPolyPatch
{
public:
enum transformType
{
UNKNOWN, // unspecified; automatic ordering
ROTATIONAL, // rotation along coordinate axis
TRANSLATIONAL, // translation
NOORDERING // unspecified, no automatic ordering
};
static const NamedEnum<transformType, 4> transformTypeNames;
private:
// Private data
//- Name of other half

View File

@ -24,7 +24,6 @@ License
\*---------------------------------------------------------------------------*/
#include "Analytical.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -59,9 +58,9 @@ template<class Type>
typename Foam::IntegrationScheme<Type>::integrationResult
Foam::Analytical<Type>::integrate
(
const Type phi,
const Type& phi,
const scalar dt,
const Type alpha,
const Type& alphaBeta,
const scalar beta
) const
{
@ -69,9 +68,18 @@ Foam::Analytical<Type>::integrate
const scalar expTerm = exp(min(50, -beta*dt));
retValue.average() =
alpha + (phi - alpha)*(1 - expTerm)/(beta*dt + ROOTVSMALL);
retValue.value() = alpha + (phi - alpha)*expTerm;
if (beta > ROOTVSMALL)
{
const Type alpha = alphaBeta/beta;
retValue.average() = alpha + (phi - alpha)*(1 - expTerm)/(beta*dt);
retValue.value() = alpha + (phi - alpha)*expTerm;
}
else
{
retValue.average() = phi;
retValue.value() = phi;
}
return retValue;
}

View File

@ -81,9 +81,9 @@ public:
//- Perform the integration
virtual typename IntegrationScheme<Type>::integrationResult integrate
(
const Type phi,
const Type& phi,
const scalar dt,
const Type alpha,
const Type& alphaBeta,
const scalar beta
) const;
};

View File

@ -24,7 +24,6 @@ License
\*---------------------------------------------------------------------------*/
#include "Euler.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -59,14 +58,14 @@ template<class Type>
typename Foam::IntegrationScheme<Type>::integrationResult
Foam::Euler<Type>::integrate
(
const Type phi,
const Type& phi,
const scalar dt,
const Type alpha,
const Type& alphaBeta,
const scalar beta
) const
{
typename IntegrationScheme<Type>::integrationResult retValue;
retValue.value() = (phi + beta*dt*alpha)/(1.0 + beta*dt);
retValue.value() = (phi + alphaBeta*dt)/(1.0 + beta*dt);
retValue.average() = 0.5*(phi + retValue.value());
return retValue;

View File

@ -78,9 +78,9 @@ public:
//- Perform the integration
virtual typename IntegrationScheme<Type>::integrationResult integrate
(
const Type phi,
const Type& phi,
const scalar dt,
const Type alpha,
const Type& alphaBeta,
const scalar beta
) const;
};

View File

@ -60,9 +60,9 @@ template<class Type>
typename Foam::IntegrationScheme<Type>::integrationResult
Foam::IntegrationScheme<Type>::integrate
(
const Type phi,
const Type& phi,
const scalar dt,
const Type alpha,
const Type& alphaBeta,
const scalar beta
) const
{
@ -71,9 +71,9 @@ Foam::IntegrationScheme<Type>::integrate
"Foam::IntegrationScheme<Type>::integrationResult"
"Foam::IntegrationScheme<Type>::integrate"
"("
"const Type, "
"const Type&, "
"const scalar, "
"const Type, "
"const Type&, "
"const scalar"
") const"
);

View File

@ -183,9 +183,9 @@ public:
//- Perform the Integration
virtual integrationResult integrate
(
const Type phi,
const Type& phi,
const scalar dt,
const Type alpha,
const Type& alphaBeta,
const scalar beta
) const;
};

View File

@ -114,11 +114,11 @@ Foam::CollidingCloud<CloudType>::CollidingCloud
if (this->solution().active())
{
setModels();
}
if (readFields)
{
parcelType::readFields(*this);
if (readFields)
{
parcelType::readFields(*this);
}
}
}

View File

@ -360,11 +360,11 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
if (solution_.active())
{
setModels();
}
if (readFields)
{
parcelType::readFields(*this);
if (readFields)
{
parcelType::readFields(*this);
}
}
if (solution_.resetSourcesOnStartup())
@ -504,7 +504,7 @@ void Foam::KinematicCloud<CloudType>::checkParcelProperties
parcel.rho() = constProps_.rho0();
}
const scalar carrierDt = this->db().time().deltaTValue();
const scalar carrierDt = solution_.deltaT();
parcel.stepFraction() = (carrierDt - lagrangianDt)/carrierDt;
parcel.typeId() = constProps_.parcelTypeId();
}

View File

@ -374,16 +374,38 @@ Foam::KinematicCloud<CloudType>::theta() const
)
);
scalarField& theta = ttheta().internalField();
volScalarField& theta = ttheta();
theta.boundaryField() == 0;
forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
{
const parcelType& p = iter();
const label cellI = p.cell();
if ((p.face() != -1))
{
const label patchI = p.patch(p.face());
if (patchI != -1)
{
scalarField& thetap = theta.boundaryField()[patchI];
const label faceI = p.patchFace(patchI, p.face());
thetap[faceI] += p.nParticle()*p.areaP();
}
}
theta[cellI] += p.nParticle()*p.volume();
}
theta /= mesh().V();
theta.internalField() /= mesh_.V();
forAll(theta.boundaryField(), patchI)
{
scalarField& thetap = theta.boundaryField()[patchI];
if (thetap.size() > 0)
{
thetap /= mesh_.magSf().boundaryField()[patchI];
}
}
return ttheta;
}
@ -420,7 +442,7 @@ Foam::KinematicCloud<CloudType>::alpha() const
alpha[cellI] += p.nParticle()*p.mass();
}
alpha /= (mesh().V()*rho_);
alpha /= (mesh_.V()*rho_);
return talpha;
}
@ -457,7 +479,7 @@ Foam::KinematicCloud<CloudType>::rhoEff() const
rhoEff[cellI] += p.nParticle()*p.mass();
}
rhoEff /= mesh().V();
rhoEff /= mesh_.V();
return trhoEff;
}

View File

@ -117,6 +117,11 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
if (this->solution().active())
{
setModels();
if (readFields)
{
parcelType::readFields(*this, this->composition());
}
}
// Set storage for mass source fields and initialise to zero
@ -142,11 +147,6 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
);
}
if (readFields)
{
parcelType::readFields(*this, this->composition());
}
if (this->solution().resetSourcesOnStartup())
{
resetSourceTerms();

View File

@ -94,11 +94,11 @@ Foam::ReactingMultiphaseCloud<CloudType>::ReactingMultiphaseCloud
if (this->solution().active())
{
setModels();
}
if (readFields)
{
parcelType::readFields(*this, this->composition());
if (readFields)
{
parcelType::readFields(*this, this->composition());
}
}
if (this->solution().resetSourcesOnStartup())

View File

@ -135,11 +135,11 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
if (this->solution().active())
{
setModels();
}
if (readFields)
{
parcelType::readFields(*this);
if (readFields)
{
parcelType::readFields(*this);
}
}
if (this->solution().resetSourcesOnStartup())

View File

@ -192,21 +192,20 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
const parcelType& p = static_cast<const parcelType&>(*this);
const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
forceSuSp Feff = Fcp + Fncp;
Feff.Sp() += ROOTVSMALL;
const forceSuSp Feff = Fcp + Fncp;
// New particle velocity
//~~~~~~~~~~~~~~~~~~~~~~
// Update velocity - treat as 3-D
const vector ap = Uc_ + (Feff.Su() + Su)/Feff.Sp();
const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/mass;
const scalar bp = Feff.Sp()/mass;
Spu = Feff.Sp();
Spu = Feff.Sp()*dt/td.cloud().solution().deltaT();
IntegrationScheme<vector>::integrationResult Ures =
td.cloud().UIntegrator().integrate(U, dt, ap, bp);
td.cloud().UIntegrator().integrate(U, dt, abp, bp);
vector Unew = Ures.value();

View File

@ -333,13 +333,13 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
// Integrate to find the new parcel temperature
IntegrationScheme<scalar>::integrationResult Tres =
td.cloud().TIntegrator().integrate(T, dt, ap, bp);
td.cloud().TIntegrator().integrate(T, dt, ap*bp, bp);
scalar Tnew = max(Tres.value(), td.cloud().constProps().TMin());
dhsTrans += dt*htc*As*(0.5*(T + Tnew) - Tc_);
Cuh = bp;
Cuh = bp*dt/td.cloud().solution().deltaT();
return Tnew;
}

View File

@ -520,7 +520,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
}
const scalar time = this->owner().db().time().value();
const scalar carrierDt = this->owner().db().time().deltaTValue();
const scalar carrierDt = this->owner().solution().deltaT();
const polyMesh& mesh = this->owner().mesh();
// Prepare for next time step

View File

@ -290,7 +290,7 @@ void Foam::PatchInteractionModel<CloudType>::patchData
}
else
{
Up = (Cf - Cf00)/mesh.time().deltaTValue();
Up = (Cf - Cf00)/this->owner().solution().deltaT();
}
if (mag(dn) > SMALL)
@ -312,7 +312,9 @@ void Foam::PatchInteractionModel<CloudType>::patchData
// magOmega = sin(angle between unit normals)
// Normalise omega vector by magOmega, then multiply by
// angle/dt to give the correct angular velocity vector.
omega *= Foam::asin(magOmega)/(magOmega*mesh.time().deltaTValue());
omega *=
Foam::asin(magOmega)
/(magOmega*this->owner().solution().deltaT());
// Project position onto face and calculate this position
// relative to the face centre.

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -709,16 +709,17 @@ Foam::directMappedPatchBase::directMappedPatchBase
}
else
{
FatalErrorIn
FatalIOErrorIn
(
"directMappedPatchBase::directMappedPatchBase\n"
"(\n"
" const polyPatch& pp,\n"
" const dictionary& dict\n"
")\n"
")\n",
dict
) << "Please supply the offsetMode as one of "
<< NamedEnum<offsetMode, 3>::words()
<< exit(FatalError);
<< exit(FatalIOError);
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -53,23 +53,15 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
{
const vector& sample = operator[](probeI);
label faceI = meshSearchEngine.findNearestBoundaryFace(sample);
if (faceI == -1)
{
nearest[probeI].second().first() = Foam::sqr(GREAT);
nearest[probeI].second().second() = Pstream::myProcNo();
}
else
{
const point& fc = mesh.faceCentres()[faceI];
nearest[probeI].first() = pointIndexHit
(
true,
fc,
faceI
);
nearest[probeI].second().first() = magSqr(fc-sample);
nearest[probeI].second().second() = Pstream::myProcNo();
}
const point& fc = mesh.faceCentres()[faceI];
nearest[probeI].first() = pointIndexHit
(
true,
fc,
faceI
);
nearest[probeI].second().first() = magSqr(fc-sample);
nearest[probeI].second().second() = Pstream::myProcNo();
}
@ -92,27 +84,16 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
}
}
// Check if all patchProbes have been found.
forAll(nearest, sampleI)
{
label localI = nearest[sampleI].first().index();
label localI = -1;
if (nearest[sampleI].second().second() == Pstream::myProcNo())
{
localI = nearest[sampleI].first().index();
}
if (localI == -1)
{
if (Pstream::master())
{
WarningIn("patchProbes::findElements()")
<< "Did not find location "
<< nearest[sampleI].second().first()
<< " in any cell. Skipping location." << endl;
}
}
else
{
elementList_[sampleI] = localI;
}
elementList_[sampleI] = localI;
}
}

View File

@ -81,17 +81,18 @@ Foam::radiationCoupledBase::radiationCoupledBase
{
if (!isA<directMappedPatchBase>(patch_.patch()))
{
FatalErrorIn
FatalIOErrorIn
(
"radiationCoupledBase::radiationCoupledBase\n"
"(\n"
" const fvPatch& p,\n"
" const dictionary& dict\n"
")\n"
")\n",
dict
) << "\n patch type '" << patch_.type()
<< "' not type '" << directMappedPatchBase::typeName << "'"
<< "\n for patch " << patch_.name()
<< exit(FatalError);
<< exit(FatalIOError);
}
emissivity_ = scalarField(patch_.size(), 0.0);
@ -102,16 +103,17 @@ Foam::radiationCoupledBase::radiationCoupledBase
{
if(!dict.found("emissivity"))
{
FatalErrorIn
FatalIOErrorIn
(
"radiationCoupledBase::radiationCoupledBase\n"
"(\n"
" const fvPatch& p,\n"
" const dictionary& dict\n"
")\n"
")\n",
dict
) << "\n emissivity key does not exist for patch "
<< patch_.name()
<< exit(FatalError);
<< exit(FatalIOError);
}
else
{