Commit 88d644f5 authored by Henry Weller's avatar Henry Weller
Browse files

lagrangian: Move penetration function from KinematicCloud to SprayCloud

Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=992
parent 62eb0e15
......@@ -400,7 +400,6 @@ public:
//- Optional particle forces
// inline const typename parcelType::forceType& forces() const;
inline const forceType& forces() const;
//- Return the optional particle forces
......@@ -499,9 +498,6 @@ public:
//- Total rotational kinetic energy in the system
inline scalar rotationalKineticEnergyOfSystem() const;
//- Penetration for fraction [0-1] of the current total mass
inline scalar penetration(const scalar fraction) const;
//- Mean diameter Dij
inline scalar Dij(const label i, const label j) const;
......
......@@ -355,142 +355,6 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
}
template<class CloudType>
inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
(
const scalar fraction
) const
{
if ((fraction < 0) || (fraction > 1))
{
FatalErrorInFunction
<< "fraction should be in the range 0 < fraction < 1"
<< exit(FatalError);
}
scalar distance = 0.0;
const label nParcel = this->size();
globalIndex globalParcels(nParcel);
const label nParcelSum = globalParcels.size();
if (nParcelSum == 0)
{
return distance;
}
// lists of parcels mass and distance from initial injection point
List<List<scalar>> procMass(Pstream::nProcs());
List<List<scalar>> procDist(Pstream::nProcs());
List<scalar>& mass = procMass[Pstream::myProcNo()];
List<scalar>& dist = procDist[Pstream::myProcNo()];
mass.setSize(nParcel);
dist.setSize(nParcel);
label i = 0;
scalar mSum = 0.0;
forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
{
const parcelType& p = iter();
scalar m = p.nParticle()*p.mass();
scalar d = mag(p.position() - p.position0());
mSum += m;
mass[i] = m;
dist[i] = d;
i++;
}
// calculate total mass across all processors
reduce(mSum, sumOp<scalar>());
Pstream::gatherList(procMass);
Pstream::gatherList(procDist);
if (Pstream::master())
{
// flatten the mass lists
List<scalar> allMass(nParcelSum, 0.0);
SortableList<scalar> allDist(nParcelSum, 0.0);
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
SubList<scalar>
(
allMass,
globalParcels.localSize(procI),
globalParcels.offset(procI)
).assign(procMass[procI]);
// flatten the distance list
SubList<scalar>
(
allDist,
globalParcels.localSize(procI),
globalParcels.offset(procI)
).assign(procDist[procI]);
}
// sort allDist distances into ascending order
// note: allMass masses are left unsorted
allDist.sort();
if (nParcelSum > 1)
{
const scalar mLimit = fraction*mSum;
const labelList& indices = allDist.indices();
if (mLimit > (mSum - allMass[indices.last()]))
{
distance = allDist.last();
}
else
{
// assuming that 'fraction' is generally closer to 1 than 0,
// loop through in reverse distance order
const scalar mThreshold = (1.0 - fraction)*mSum;
scalar mCurrent = 0.0;
label i0 = 0;
forAllReverse(indices, i)
{
label indI = indices[i];
mCurrent += allMass[indI];
if (mCurrent > mThreshold)
{
i0 = i;
break;
}
}
if (i0 == indices.size() - 1)
{
distance = allDist.last();
}
else
{
// linearly interpolate to determine distance
scalar alpha = (mCurrent - mThreshold)/allMass[indices[i0]];
distance =
allDist[i0] + alpha*(allDist[i0+1] - allDist[i0]);
}
}
}
else
{
distance = allDist.first();
}
}
Pstream::scatter(distance);
return distance;
}
template<class CloudType>
inline Foam::cachedRandom& Foam::KinematicCloud<CloudType>::rndGen()
{
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -85,12 +85,6 @@ public:
//- Total linear kinetic energy in the system
virtual scalar linearKineticEnergyOfSystem() const = 0;
//- Total rotational kinetic energy in the system
// virtual scalar rotationalKineticEnergyOfSystem() const = 0;
//- Penetration for percentage of the current total mass
// virtual scalar penetration(const scalar& fraction) const = 0;
//- Mean diameter Dij
virtual scalar Dij(const label i, const label j) const = 0;
......
......@@ -183,6 +183,12 @@ public:
inline scalar averageParcelMass() const;
// Check
//- Penetration for fraction [0-1] of the current total mass
inline scalar penetration(const scalar fraction) const;
// Sub-models
//- Return const-access to the atomization model
......
......@@ -72,4 +72,140 @@ inline Foam::scalar Foam::SprayCloud<CloudType>::averageParcelMass() const
}
template<class CloudType>
inline Foam::scalar Foam::SprayCloud<CloudType>::penetration
(
const scalar fraction
) const
{
if ((fraction < 0) || (fraction > 1))
{
FatalErrorInFunction
<< "fraction should be in the range 0 < fraction < 1"
<< exit(FatalError);
}
scalar distance = 0.0;
const label nParcel = this->size();
globalIndex globalParcels(nParcel);
const label nParcelSum = globalParcels.size();
if (nParcelSum == 0)
{
return distance;
}
// lists of parcels mass and distance from initial injection point
List<List<scalar>> procMass(Pstream::nProcs());
List<List<scalar>> procDist(Pstream::nProcs());
List<scalar>& mass = procMass[Pstream::myProcNo()];
List<scalar>& dist = procDist[Pstream::myProcNo()];
mass.setSize(nParcel);
dist.setSize(nParcel);
label i = 0;
scalar mSum = 0.0;
forAllConstIter(typename SprayCloud<CloudType>, *this, iter)
{
const parcelType& p = iter();
scalar m = p.nParticle()*p.mass();
scalar d = mag(p.position() - p.position0());
mSum += m;
mass[i] = m;
dist[i] = d;
i++;
}
// calculate total mass across all processors
reduce(mSum, sumOp<scalar>());
Pstream::gatherList(procMass);
Pstream::gatherList(procDist);
if (Pstream::master())
{
// flatten the mass lists
List<scalar> allMass(nParcelSum, 0.0);
SortableList<scalar> allDist(nParcelSum, 0.0);
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
SubList<scalar>
(
allMass,
globalParcels.localSize(procI),
globalParcels.offset(procI)
).assign(procMass[procI]);
// flatten the distance list
SubList<scalar>
(
allDist,
globalParcels.localSize(procI),
globalParcels.offset(procI)
).assign(procDist[procI]);
}
// sort allDist distances into ascending order
// note: allMass masses are left unsorted
allDist.sort();
if (nParcelSum > 1)
{
const scalar mLimit = fraction*mSum;
const labelList& indices = allDist.indices();
if (mLimit > (mSum - allMass[indices.last()]))
{
distance = allDist.last();
}
else
{
// assuming that 'fraction' is generally closer to 1 than 0,
// loop through in reverse distance order
const scalar mThreshold = (1.0 - fraction)*mSum;
scalar mCurrent = 0.0;
label i0 = 0;
forAllReverse(indices, i)
{
label indI = indices[i];
mCurrent += allMass[indI];
if (mCurrent > mThreshold)
{
i0 = i;
break;
}
}
if (i0 == indices.size() - 1)
{
distance = allDist.last();
}
else
{
// linearly interpolate to determine distance
scalar alpha = (mCurrent - mThreshold)/allMass[indices[i0]];
distance =
allDist[i0] + alpha*(allDist[i0+1] - allDist[i0]);
}
}
}
else
{
distance = allDist.first();
}
}
Pstream::scatter(distance);
return distance;
}
// ************************************************************************* //
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment