ENH: Updated cloud penetration calc

This commit is contained in:
andy
2013-02-01 16:29:07 +00:00
parent ed96224c12
commit b3d5349edb

View File

@ -372,8 +372,14 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
}
// lists of parcels mass and distance from initial injection point
List<scalar> mass(nParcel, 0.0);
List<scalar> dist(nParcel, 0.0);
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;
@ -392,75 +398,86 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
// calculate total mass across all processors
reduce(mSum, sumOp<scalar>());
Pstream::gatherList(procMass);
Pstream::gatherList(procDist);
// flatten the mass list
List<scalar> allMass(nParcelSum, 0.0);
SubList<scalar>
(
allMass,
globalParcels.localSize(Pstream::myProcNo()),
globalParcels.offset(Pstream::myProcNo())
).assign(mass);
// flatten the distance list
SortableList<scalar> allDist(nParcelSum, 0.0);
SubList<scalar>
(
allDist,
globalParcels.localSize(Pstream::myProcNo()),
globalParcels.offset(Pstream::myProcNo())
).assign(dist);
// sort allDist distances into ascending order
// note: allMass masses are left unsorted
allDist.sort();
if (nParcelSum > 1)
if (Pstream::master())
{
const scalar mLimit = fraction*mSum;
const labelList& indices = allDist.indices();
if (mLimit > (mSum - allMass[indices.last()]))
// flatten the mass lists
List<scalar> allMass(nParcelSum, 0.0);
SortableList<scalar> allDist(nParcelSum, 0.0);
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
distance = allDist.last();
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]);
}
else
// sort allDist distances into ascending order
// note: allMass masses are left unsorted
allDist.sort();
if (nParcelSum > 1)
{
// 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;
const scalar mLimit = fraction*mSum;
const labelList& indices = allDist.indices();
forAllReverse(indices, i)
{
label indI = indices[i];
mCurrent += allMass[indI];
if (mCurrent > mThreshold)
{
i0 = i;
break;
}
}
if (i0 == indices.size() - 1)
if (mLimit > (mSum - allMass[indices.last()]))
{
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]);
// 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();
}
}
else
{
distance = allDist.first();
}
Pstream::scatter(distance);
return distance;
}