From b76178bd16875f7b183eb131d4ff24691539905e Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Tue, 18 Sep 2012 16:38:59 +0100
Subject: [PATCH] ENH: Corrected cloud penetration for parallel cases
---
.../KinematicCloud/KinematicCloudI.H | 165 ++++++++++--------
1 file changed, 95 insertions(+), 70 deletions(-)
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index a6bdf47cc46..87b9136a4ec 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "fvmSup.H"
+#include "SortableList.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@@ -343,99 +344,123 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
template<class CloudType>
inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
(
- const scalar& prc
+ const scalar& fraction
) const
{
- scalar distance = 0.0;
- scalar mTot = 0.0;
+ if ((fraction < 0) || (fraction > 1))
+ {
+ FatalErrorIn
+ (
+ "inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration"
+ "("
+ "const scalar&"
+ ") const"
+ )
+ << "fraction should be in the range 0 < fraction < 1"
+ << exit(FatalError);
+ }
- label np = this->size();
+ scalar distance = 0.0;
- // arrays containing the parcels mass and
- // distance from injector in ascending order
- scalarField mass(np);
- scalarField dist(np);
+ const label nParcel = this->size();
+ globalIndex globalParcels(nParcel);
+ const label nParcelSum = globalParcels.size();
- if (np > 0)
+ if (nParcelSum == 0)
{
- label n = 0;
+ return distance;
+ }
- // first arrange the parcels in ascending order
- // the first parcel is closest to its injection position
- // and the last one is most far away.
- forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
- {
- const parcelType& p = iter();
- scalar mi = p.nParticle()*p.mass();
- scalar di = mag(p.position() - p.position0());
- mTot += mi;
+ // lists of parcels mass and distance from initial injection point
+ List<scalar> mass(nParcel, 0.0);
+ List<scalar> dist(nParcel, 0.0);
- // insert at the last place
- mass[n] = mi;
- dist[n] = di;
+ 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;
- label i = 0;
- bool found = false;
+ mass[i] = m;
+ dist[i] = d;
- // insert the parcel in the correct place
- // and move the others
- while ((i < n) && (!found))
- {
- if (di < dist[i])
- {
- found = true;
- for (label j=n; j>i; j--)
- {
- mass[j] = mass[j-1];
- dist[j] = dist[j-1];
- }
- mass[i] = mi;
- dist[i] = di;
- }
- i++;
- }
- n++;
- }
+ i++;
}
- reduce(mTot, sumOp<scalar>());
+ // calculate total mass across all processors
+ reduce(mSum, sumOp<scalar>());
+
+ // 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 (np > 0)
+ if (nParcelSum > 1)
{
- scalar mLimit = prc*mTot;
- scalar mOff = (1.0 - prc)*mTot;
+ const scalar mLimit = fraction*mSum;
+ const labelList& indices = allDist.indices();
- if (np > 1)
+ if (mLimit > (mSum - allMass[indices.last()]))
{
- // 'prc' is large enough that the parcel most far
- // away will be used, no need to loop...
- if (mLimit > mTot - mass[np-1])
- {
- distance = dist[np-1];
- }
- else
+ 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)
{
- scalar mOffSum = 0.0;
- label i = np;
+ label indI = indices[i];
+
+ mCurrent += allMass[indI];
- while ((mOffSum < mOff) && (i>0))
+ if (mCurrent > mThreshold)
{
- i--;
- mOffSum += mass[i];
+ i0 = i;
+ break;
}
- distance =
- dist[i+1]
- + (dist[i] - dist[i+1])*(mOffSum - mOff)
- /mass[i+1] ;
}
- }
- else
- {
- distance = dist[0];
+
+ 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]);
+ }
}
}
-
- reduce(distance, maxOp<scalar>());
+ else
+ {
+ distance = allDist.first();
+ }
return distance;
}
--
GitLab