From b3d5349edb2747c888a895e4bbea77baa3930240 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Fri, 1 Feb 2013 16:29:07 +0000
Subject: [PATCH] ENH: Updated cloud penetration calc

---
 .../KinematicCloud/KinematicCloudI.H          | 129 ++++++++++--------
 1 file changed, 73 insertions(+), 56 deletions(-)

diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index be67f182926..68acada6cad 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -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();
-        }
-        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;
+            SubList<scalar>
+            (
+                allMass,
+                globalParcels.localSize(procI),
+                globalParcels.offset(procI)
+            ).assign(procMass[procI]);
 
-            forAllReverse(indices, i)
-            {
-                label indI = indices[i];
+            // flatten the distance list
+            SubList<scalar>
+            (
+                allDist,
+                globalParcels.localSize(procI),
+                globalParcels.offset(procI)
+            ).assign(procDist[procI]);
+        }
 
-                mCurrent += allMass[indI];
+        // sort allDist distances into ascending order
+        // note: allMass masses are left unsorted
+        allDist.sort();
 
-                if (mCurrent > mThreshold)
-                {
-                    i0 = i;
-                    break;
-                }
-            }
+        if (nParcelSum > 1)
+        {
+            const scalar mLimit = fraction*mSum;
+            const labelList& indices = allDist.indices();
 
-            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;
 }
-- 
GitLab