ENH: extractEulerianParticles - refactoring and robustness improvements

src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -264,214 +264,163 @@ void Foam::functionObjects::extractEulerianParticles::setBlockedFaces
 }
 
 
-void Foam::functionObjects::extractEulerianParticles::collectParticles
+void Foam::functionObjects::extractEulerianParticles::collectParticle
 (
     const scalar time,
-    const boolList& collectParticle
+    const label regioni
 )
 {
-    DebugInFunction << "collectParticle: " << collectParticle << endl;
+    DebugInFunction << "collectParticle: " << regioni << endl;
 
-    // Collect particles on local processors that have passed through faceZone
+    const label particlei = regionToParticleMap_[regioni];
-    forAll(collectParticle, regioni)
+    eulerianParticle p = particles_[particlei];
+
+    if (p.faceIHit != -1 && nInjectorLocations_)
     {
-        if (!collectParticle[regioni])
+        // Use coarse face index for tag output
+        label coarseFacei = fineToCoarseAddr_[p.faceIHit];
+        p.faceIHit = globalCoarseFaces_.toGlobal(coarseFacei);
+    }
+
+    reduce(p, sumParticleOp<eulerianParticle>());
+
+    const scalar pDiameter = cbrt(6.0*p.V/constant::mathematical::pi);
+
+    if ((pDiameter > minDiameter_) && (pDiameter < maxDiameter_))
+    {
+        if (Pstream::master())
         {
-            continue;
+            const scalar d = cbrt(6*p.V/constant::mathematical::pi);
-        }
+            const point position = p.VC/(p.V + ROOTVSMALL);
-
+            const vector U = p.VU/(p.V + ROOTVSMALL);
-        Map<label>::const_iterator iter = regionToParticleMap_.find(regioni);
+            label tag = -1;
-        eulerianParticle p = particles_[iter()];
+            if (nInjectorLocations_)
-
-        if (p.faceIHit != -1 && nInjectorLocations_)
-        {
-            // Use coarse face index for tag output
-            label coarseFacei = fineToCoarseAddr_[p.faceIHit];
-            p.faceIHit = globalCoarseFaces_.toGlobal(coarseFacei);
-        }
-
-        reduce(p, sumParticleOp<eulerianParticle>());
-
-        const scalar pDiameter = cbrt(6.0*p.V/constant::mathematical::pi);
-
-        if ((pDiameter > minDiameter_) && (pDiameter < maxDiameter_))
-        {
-            if (Pstream::master())
             {
-                const scalar d = cbrt(6*p.V/constant::mathematical::pi);
+                tag = p.faceIHit;
-                const point position = p.VC/(p.V + ROOTVSMALL);
-                const vector U = p.VU/(p.V + ROOTVSMALL);
-                label tag = -1;
-                if (nInjectorLocations_)
-                {
-                    tag = p.faceIHit;
-                }
-
-                injectedParticle* ip = new injectedParticle
-                (
-                    mesh_,
-                    position,
-                    tag,
-                    time,
-                    d,
-                    U,
-                    false // not looking to set cell owner etc.
-                );
-
-                cloud_.addParticle(ip);
             }
 
-            nCollectedParticles_++;
+            injectedParticle* ip = new injectedParticle
-        }
+            (
-        else
+                mesh_,
-        {
+                position,
-            // Discard particles over/under diameter threshold
+                tag,
-            nDiscardedParticles_++;
+                time,
-            discardedVolume_ += p.V;
+                d,
+                U,
+                false // not looking to set cell owner etc.
+            );
+
+            cloud_.addParticle(ip);
+
+            collectedVolume_ += p.V;
         }
+
+        ++nCollectedParticles_;
+    }
+    else
+    {
+        // Discard particles over/under diameter threshold
+        ++nDiscardedParticles_;
+        discardedVolume_ += p.V;
     }
 }
 
 
 void Foam::functionObjects::extractEulerianParticles::calculateAddressing
 (
-    const label nRegionsOld,
+    const label nNewRegions,
-    const label nRegionsNew,
     const scalar time,
     labelList& regionFaceIDs
 )
 {
     DebugInFunction << endl;
 
-    // New region can only point to one old region
+    // Determine mapping between old and new regions so that we can
-    // Old region can only point to one new region.  If old region intersects
+    // accumulate particle info
-    // multiple new regions, select max of new region indices.
+    labelList oldToNewRegion(particles_.size(), -1);
-
+    labelList newToNewRegion(identity(nNewRegions));
-    labelList oldToNewRegion(nRegionsOld, -1);
-    labelList newToOldRegion(nRegionsNew, -1);
 
     forAll(regionFaceIDs, facei)
     {
         label newRegioni = regionFaceIDs[facei];
         label oldRegioni = regions0_[facei];
 
-        if (newRegioni != -1)
+        if (newRegioni != -1 && oldRegioni != -1)
         {
-            newToOldRegion[newRegioni] = oldRegioni;
+            // If old region has split into multiple regions we need to
-
+            // renumber new regions to maintain connectivity with old regions
-            if (oldRegioni != -1)
+            newToNewRegion[newRegioni] =
-            {
+                max(newRegioni, oldToNewRegion[oldRegioni]);
-                // New region linked to old (so can accumulate particle data)
+            oldToNewRegion[oldRegioni] = newRegioni;
-                // Old region might already be mapped to a new region
-                oldToNewRegion[oldRegioni] =
-                    max(newRegioni, oldToNewRegion[oldRegioni]);
-            }
         }
     }
 
-    // Need to re-number the new region indices based on knowledge of which
+    // Create map from new regions to slots in particles list
-    // old region they intersect.  After operations, there should be a
+    // - filter through new-to-new addressing to identify new particles
-    // one-to-one match between the old and new regions.
+    Pstream::listCombineGather(newToNewRegion, maxEqOp<label>());
+    Pstream::listCombineScatter(newToNewRegion);
 
-    // Ensure all old regions point to the same new regions (if present)
+    label nParticle = -1;
+    labelHashSet newRegions;
+    Map<label> newRegionToParticleMap;
+    forAll(newToNewRegion, newRegioni0)
+    {
+        label newRegioni = newToNewRegion[newRegioni0];
+        if (newRegions.insert(newRegioni))
+        {
+            ++nParticle;
+        }
+
+        // New particle slot
+        newRegionToParticleMap.insert(newRegioni0, nParticle);
+    }
+
+    // Accumulate old region data or create a new particle if there is no
+    // mapping from the old-to-new region
     Pstream::listCombineGather(oldToNewRegion, maxEqOp<label>());
     Pstream::listCombineScatter(oldToNewRegion);
-
+    List<eulerianParticle> newParticles(newRegionToParticleMap.size());
-    // Any new region that points to an old region should be renumbered to the
+    forAll(oldToNewRegion, oldRegioni)
-    // new region specified by the oldToNewRegion index
-
-    if (oldToNewRegion.size())
     {
-        // Create corrected new to new addressing
+        label newRegioni = oldToNewRegion[oldRegioni];
-        labelList newToNewRegionCorr(newToOldRegion.size(), -1);
+        if (newRegioni == -1)
-        forAll(newToOldRegion, newRegioni)
         {
-            label oldRegioni = newToOldRegion[newRegioni];
+            // No mapping from old-to-new - collect new particle
-            if (oldRegioni != -1)
+            DebugInfo
-            {
+                << "Collecting particle from oldRegion:" << oldRegioni
-                label newRegionICorr = oldToNewRegion[oldRegioni];
+                << endl;
-                newToNewRegionCorr[newRegioni] = newRegionICorr;
-                newToOldRegion[newRegioni] = oldRegioni;
-            }
-        }
 
-        // Renumber the new (current) face region IDs
+            collectParticle(time, oldRegioni);
-        forAll(regionFaceIDs, facei)
-        {
-            label newRegioni = regionFaceIDs[facei];
-
-            if (newRegioni != -1)
-            {
-                label newRegionICorr = newToNewRegionCorr[newRegioni];
-
-                // Note: newRegionICorr can be -1 if the region is new, since
-                // the address corrections are based on inverting the
-                // old-to-new addressing
-                if (newRegionICorr != -1)
-                {
-                    regionFaceIDs[facei] = newRegionICorr;
-                }
-            }
-        }
-
-
-        boolList collectParticleFlag(nRegionsOld, true);
-        forAll(oldToNewRegion, oldRegioni)
-        {
-            label newRegioni = oldToNewRegion[oldRegioni];
-            if (newRegioni != -1)
-            {
-                collectParticleFlag[oldRegioni] = false;
-            }
-        }
-
-        // Collect particles whose IDs are no longer active
-        collectParticles(time, collectParticleFlag);
-    }
-
-
-    // Re-order collection bins
-
-    Map<label> newRegionToParticleMap(nRegionsNew);
-    List<eulerianParticle> newParticles(nRegionsNew);
-    label particlei = 0;
-
-    forAll(newToOldRegion, newRegioni)
-    {
-        label oldRegioni = newToOldRegion[newRegioni];
-        if (oldRegioni == -1)
-        {
-            // No mapping from old to new - this is a new particle
-            newRegionToParticleMap.insert(newRegioni, particlei);
-            particlei++;
         }
         else
         {
-            // Update addressing for old to new regions
+            // Combine existing particle into new particle
+            label newParticlei = newRegionToParticleMap[newRegioni];
             label oldParticlei = regionToParticleMap_[oldRegioni];
-            if (newRegionToParticleMap.insert(newRegioni, particlei))
+
-            {
+            DebugInfo
-                // First time this particle has been seen
+                << "Combining newRegioni: " << newRegioni
-                newParticles[particlei] = particles_[oldParticlei];
+                << "(p:" << newParticlei << ") and "
-                particlei++;
+                << "oldRegioni: " << oldRegioni
-            }
+                << "(p:" << oldParticlei << ")"
-            else
+                << endl;
-            {
+
-                // Combine with existing contribution
+            newParticles[newParticlei] =
-                label newParticlei = newRegionToParticleMap[newRegioni];
+                sumParticleOp<eulerianParticle>()
-                newParticles[newParticlei] =
+                (
-                    sumParticleOp<eulerianParticle>()
+                    newParticles[newParticlei],
-                    (
+                    particles_[oldParticlei]
-                        newParticles[newParticlei],
+                );
-                        particles_[oldParticlei]
-                    );
-            }
         }
     }
 
     // Reset the particles list and addressing for latest available info
     particles_.transfer(newParticles);
     regionToParticleMap_ = newRegionToParticleMap;
+
+    // Reset the region IDs for the next integration step
+    // - these become the oldRegioni's
+    regions0_ = regionFaceIDs;
 }
 
 
@@ -494,7 +443,6 @@ void Foam::functionObjects::extractEulerianParticles::accumulateParticleInfo
     forAll(regionFaceIDs, localFacei)
     {
         const label newRegioni = regionFaceIDs[localFacei];
-
         if (newRegioni != -1)
         {
             const label particlei = regionToParticleMap_[newRegioni];
@@ -547,14 +495,14 @@ Foam::functionObjects::extractEulerianParticles::extractEulerianParticles
     fineToCoarseAddr_(),
     globalCoarseFaces_(),
     regions0_(),
-    nRegions0_(0),
     particles_(),
     regionToParticleMap_(),
     minDiameter_(ROOTVSMALL),
     maxDiameter_(GREAT),
-    nCollectedParticles_(0),
+    nCollectedParticles_(getProperty<label>("nCollectedParticles", 0)),
-    nDiscardedParticles_(0),
+    collectedVolume_(getProperty<scalar>("collectedVolume", 0)),
-    discardedVolume_(0)
+    nDiscardedParticles_(getProperty<label>("nDiscardedParticles", 0)),
+    discardedVolume_(getProperty<scalar>("discardedVolume", 0))
 {
     if (mesh_.nSolutionD() != 3)
     {
@@ -567,12 +515,6 @@ Foam::functionObjects::extractEulerianParticles::extractEulerianParticles
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::functionObjects::extractEulerianParticles::~extractEulerianParticles()
-{}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 bool Foam::functionObjects::extractEulerianParticles::read
@@ -644,9 +586,9 @@ bool Foam::functionObjects::extractEulerianParticles::execute()
 
     // Calculate the addressing between the old and new region information
     // Also collects particles that have traversed the faceZone
+    // - Note: may also update regionFaceIDs
     calculateAddressing
     (
-        nRegions0_,
         nRegionsNew,
         mesh_.time().value(),
         regionFaceIDs
@@ -656,13 +598,11 @@ bool Foam::functionObjects::extractEulerianParticles::execute()
     tmp<surfaceScalarField> tphi = phiU();
     accumulateParticleInfo(alphaf, tphi(), regionFaceIDs, fz);
 
-    // Reset the blocked faces for the next integration step
+    Log << "    Collected particles   : " << nCollectedParticles_ << nl
-    nRegions0_ = nRegionsNew;
+        << "    Collected volume      : " << collectedVolume_ << nl
-    regions0_ = regionFaceIDs;
+        << "    Discarded particles   : " << nDiscardedParticles_ << nl
-
+        << "    Discarded volume      : " << discardedVolume_ << nl
-    Log << "    Collected particles: " << nCollectedParticles_ << nl
+        << "    Particles in progress : " << particles_.size() << nl
-        << "    Discarded particles: " << nDiscardedParticles_ << nl
-        << "    Discarded volume: " << discardedVolume_ << nl
         << endl;
 
     return true;
@@ -675,6 +615,11 @@ bool Foam::functionObjects::extractEulerianParticles::write()
 
     cloud_.write();
 
+    setProperty("nCollectedParticles", nCollectedParticles_);
+    setProperty("collectedVolume", collectedVolume_);
+    setProperty("nDiscardedParticles", nDiscardedParticles_);
+    setProperty("discardedVolume", discardedVolume_);
+
     return true;
 }
 

src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,6 +23,7 @@ License
 
 Class
     Foam::functionObjects::extractEulerianParticles
+
 Group
     grpFieldFunctionObjects
 
@@ -73,10 +74,8 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "polyMesh.H"
 #include "surfaceFieldsFwd.H"
-#include "vectorList.H"
 #include "globalIndex.H"
 #include "eulerianParticle.H"
-#include "IOdictionary.H"
 #include "injectedParticleCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -153,9 +152,6 @@ protected:
             //- Region indices in faceZone faces from last iteration
             labelList regions0_;
 
-            //- Number of regions from last iteration
-            label nRegions0_;
-
             //- Particle properties (partial, being accumulated)
             List<eulerianParticle> particles_;
 
@@ -176,6 +172,9 @@ protected:
             //- Total number of collected particles
             label nCollectedParticles_;
 
+            //- Total collected volume [m3]
+            scalar collectedVolume_;
+
             //- Total number of discarded particles
             label nDiscardedParticles_;
 
@@ -206,19 +205,19 @@ protected:
         );
 
         //- Calculate the addressing between regions between iterations
+        //- Returns the number of active regions (particles)
         virtual void calculateAddressing
         (
-            const label nRegionsOld,
             const label nRegionsNew,
             const scalar time,
             labelList& regionFaceIDs
         );
 
         //- Collect particles that have passed through the faceZone
-        virtual void collectParticles
+        virtual void collectParticle
         (
             const scalar time,
-            const boolList& collectParticle
+            const label regioni
         );
 
         //- Process latest region information
@@ -265,7 +264,7 @@ public:
 
 
     //- Destructor
-    virtual ~extractEulerianParticles();
+    virtual ~extractEulerianParticles() = default;
 
 
     // Member Functions