View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0002595OpenFOAMBugpublic2017-06-28 18:092017-07-10 10:24
Reporterkarlvirgil Assigned Towill  
PriorityhighSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
PlatformCentOS 7OSCentOSOS Version7
Product Versiondev 
Fixed in Versiondev 
Summary0002595: ParticleCollector sampling error
DescriptionThe Lagrangian functionObject ParticleCollector is returning incorrect sampling for spray mass crossing through a plane. This behavior started showing after the switch to using barycentric coordinates in this commit:

commit 371762757dbe3cd38a3841a547a9bc8c1aff0b85
Author: Will Bainbridge <http://cfd.direct>
Date: Fri Apr 28 08:03:44 2017 +0100

Steps To ReproduceRun the attached case using a commit after 2017-04-28 and compare to the same case run with a commit before 2017-04-28. The spray injection for the two runs will be the same, but the sampled mass for the later will be higher and actually exceed the amount of spray injected, which is physically impossible

In the case, execute
./run.sh

Then, compare
postProcessing/lagrangian/reactingCloud1/particleCollector/0/particleCollector.dat

The sum of the particleCollector mass after this commit will be 0.1447777, when only 0.11 kg has been injected (per log.fireFoam output).
Additional InformationStudying the ParticleCollector.C file, it appears that the information fed into the functions is the same. The difference seems to come from the "classify" routine, which for the previous commits would correctly return the result, whereas for the later commits it will incorrectly classify the parcel coordinates inaccurately. This seems to be related to the "triangle" class and the switch to barycentric coordinates.
TagsNo tags attached.

Activities

karlvirgil

2017-06-28 18:09

reporter  

case.tgz (45,989 bytes)

will

2017-06-29 09:37

manager   ~0008277

In ParticleCollector::postMove, there exists the following bit of code:

    // Slightly extend end position to avoid falling within tracking tolerances
    const point position1 = position0 + 1.0001*(p.position() - position0);

This extends the tracked path a bit for unspecified tolerance reasons. If the particle ends a track *on* the collector polygon then it gets added both at the end of that track and at the start of the next one, hence the higher total mass.

If we remove this extension (change 1.0001 to 1) then I get a lower total mass of 0.107. This means that some particles which land on the polygon aren't being counted at all. Or maybe both issues are happening and cancelling out.

I don't know how to fix this. There will always be an ambiguity when a track ends (to within floating point accuracy) *on* the polygon. I am very surprised that the old tracking gave the correct total. I don't think that the particle collector can be relied upon to generate totals for either tracking method. That's probably why the totals are never explicitly output; you have to sum them up from the instantaneous values.

The face-set sums, in contrast, are always correct because figuring out whether a particle crosses a face or not is an issue of topology, so it doesn't suffer from the same accuracy issues. If you want a total mass flow over many timesteps then you should use this.

What do we think? Should we remove the magic 1.0001 factor? It would undoubtedly be more accurate, but when it does fail a particle might pass through entirely uncollected. This might be more confusing that it being collected/reported twice.

karlvirgil

2017-07-02 00:43

reporter   ~0008317

Would one possible solution be simply to temporarily store the parcel ID if the mass is collected, and the next time step don't record the mass of the parcel if it has already been collected?

karlvirgil

2017-07-03 08:10

reporter   ~0008321

I've added a patch here that addresses the problem. The patch checks to see if the parcel has been previously collected on the face. If so, it doesn't get collected again. This patch results in the correct summation of mass. One thing the patch doesn't do is clean up the values in collectedParticles_. Ideally this could be done after one time step.

Regarding the face-set sums, I am aware of this feature. However, the ParticleCollector is so much more flexible and convenient for us that we use it almost exclusively. In practice it is very difficult to get face-sets oriented in the correct direction, covering the exact area, etc. Also, it is tedious to define multiple collection areas with face-sets, and the output is not tabulated all together like the ParticleCollector class. Often we have upwards of 100 collection areas in a given simulation, which is trivial for the ParticleCollector class.
patch (2,987 bytes)    
diff --git src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
index 3f278e8..a0880e7 100644
--- src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
+++ src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
@@ -543,7 +543,8 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
     log_(this->coeffDict().lookup("log")),
     outputFilePtr_(),
     timeOld_(owner.mesh().time().value()),
-    hitFaceIDs_()
+    hitFaceIDs_(),
+    collectedParticles_()
 {
     normal_ /= mag(normal_);
 
@@ -594,6 +595,7 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
     mass_.setSize(faces_.size(), 0.0);
     massTotal_.setSize(faces_.size(), 0.0);
     massFlowRate_.setSize(faces_.size(), 0.0);
+    collectedParticles_.setSize(faces_.size());
 
     makeLogFile(faces_, points_, area_);
 }
@@ -626,7 +628,8 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
     log_(pc.log_),
     outputFilePtr_(),
     timeOld_(0.0),
-    hitFaceIDs_()
+    hitFaceIDs_(),
+    collectedParticles_(pc.collectedParticles_)
 {}
 
 
@@ -709,9 +712,21 @@ void Foam::ParticleCollector<CloudType>::postMove
             }
         }
 
-        // Add mass contribution
-        mass_[facei] += m;
+        // If not previously collected, add mass contribution
+        Switch previouslyCollected = false;
+        forAll(collectedParticles_[facei],i)
+        {
+            if (&p == collectedParticles_[facei][i] )
+            {
+                previouslyCollected = true;
+                break;
+            }
+        }
+        if(previouslyCollected)
+        {
+            collectedParticles_[facei].append(&p);
 
+            mass_[facei] += m;
             if (nSector_ == 1)
             {
                 mass_[facei + 1] += m;
@@ -725,6 +740,7 @@ void Foam::ParticleCollector<CloudType>::postMove
             }
         }
     }
+}
 
 
 // ************************************************************************* //
diff --git src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H
index d14195a..621c6e8 100644
--- src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H
+++ src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H
@@ -200,6 +200,9 @@ private:
         //- Work list to store which faces are hit
         mutable DynamicList<label> hitFaceIDs_;
 
+        //- Work list to store which particles are collected
+        List<DynamicList<typename CloudType::parcelType*>> collectedParticles_;
+
 
     // Private Member Functions
patch (2,987 bytes)   

will

2017-07-03 09:07

manager   ~0008322

Comparing particle pointers won't work in parallel. You need to be testing the particle's origProc and origId members. That, presumably, requires two lists or a list of labelPairs, and there would have to be some level of duplication between processors so that a particle isn't counted multiple times on different processes. I haven't fully figured this out, but the easiest way would be to duplicate everything on all cores.

I think the cost of this is prohibitive. Creating a dynamic list that fills up with elements referring to potentially every particle in the simulation isn't a good idea, especially if they cannot be parallelised. Searching this list for every particle that crosses the collection polygon could also be very expensive if a lot of particles are hitting the polygon.

will

2017-07-03 09:07

manager   ~0008323

I don't see why you would want to clear collectedParticles_ at the end of a time-step. Particles might land on the collection polygon at the end of the time-step and get counted both at the end of that step and the start of the next. Surely the point of the collectedParticles_ list is that it persists across time-steps to prevent this issue?

On the other hand, if a particle takes a circular route, it might validly cross the collection polygon twice. How do we discern between this valid case of counting twice from the invalid case created by the tolerance?

will

2017-07-03 09:10

manager   ~0008324

If you wanted to use face-set-sums for this purpose, then I guess you would probably need to explicitly mesh the surfaces that you want to collect through. This would be simple enough with snappy. I suppose the downside is that the mesh in this region gets more distorted as a result of all the snapping.

karlvirgil

2017-07-03 16:35

reporter   ~0008325

I've updated the patch to make sure that collectedParticles_ are only stored for one previous time step, since the problem we are trying to solve just the double counting from time step to time step. This will limit the amount of memory consumption.

As for the parallel issue, I don't really see this as a huge concern. We don't compare position0 and position1 across the processors anyway.

Since this addresses the issue identified by the original bug, could you accept this?
patch-2 (3,856 bytes)    
diff --git src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
index 3f278e8..59304e7 100644
--- src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
+++ src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
@@ -505,7 +505,10 @@ void Foam::ParticleCollector<CloudType>::write()
         mass_[facei] = 0.0;
         massTotal_[facei] = 0.0;
         massFlowRate_[facei] = 0.0;
+        collectedParticles0_[facei] = collectedParticles_[facei];
+        collectedParticles_[facei].clear();
     }
+
 }
 
 
@@ -543,7 +546,9 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
     log_(this->coeffDict().lookup("log")),
     outputFilePtr_(),
     timeOld_(owner.mesh().time().value()),
-    hitFaceIDs_()
+    hitFaceIDs_(),
+    collectedParticles_(),
+    collectedParticles0_()
 {
     normal_ /= mag(normal_);
 
@@ -594,6 +599,8 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
     mass_.setSize(faces_.size(), 0.0);
     massTotal_.setSize(faces_.size(), 0.0);
     massFlowRate_.setSize(faces_.size(), 0.0);
+    collectedParticles_.setSize(faces_.size());
+    collectedParticles0_.setSize(faces_.size());
 
     makeLogFile(faces_, points_, area_);
 }
@@ -626,7 +633,9 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
     log_(pc.log_),
     outputFilePtr_(),
     timeOld_(0.0),
-    hitFaceIDs_()
+    hitFaceIDs_(),
+    collectedParticles_(pc.collectedParticles_),
+    collectedParticles0_(pc.collectedParticles0_)
 {}
 
 
@@ -709,9 +718,29 @@ void Foam::ParticleCollector<CloudType>::postMove
             }
         }
 
-        // Add mass contribution
-        mass_[facei] += m;
+        // If not previously collected, add mass contribution
+        Switch previouslyCollected = false;
+        forAll(collectedParticles0_[facei],i)
+        {
+            if (&p == collectedParticles0_[facei][i] )
+            {
+                previouslyCollected = true;
+                break;
+            }
+        }
+        forAll(collectedParticles_[facei],i)
+        {
+            if (&p == collectedParticles_[facei][i] )
+            {
+                previouslyCollected = true;
+                break;
+            }
+        }
+        if(!previouslyCollected)
+        {
+            collectedParticles_[facei].append(&p);
 
+            mass_[facei] += m;
             if (nSector_ == 1)
             {
                 mass_[facei + 1] += m;
@@ -725,6 +754,7 @@ void Foam::ParticleCollector<CloudType>::postMove
             }
         }
     }
+}
 
 
 // ************************************************************************* //
diff --git src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H
index d14195a..ebe1512 100644
--- src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H
+++ src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H
@@ -200,6 +200,12 @@ private:
         //- Work list to store which faces are hit
         mutable DynamicList<label> hitFaceIDs_;
 
+        //- Work list to store which particles are collected at current time step
+        List<DynamicList<typename CloudType::parcelType*>> collectedParticles_;
+
+        //- Work list to store which particles are collected at previous time step
+        List<DynamicList<typename CloudType::parcelType*>> collectedParticles0_;
+
 
     // Private Member Functions
patch-2 (3,856 bytes)   

will

2017-07-04 20:35

manager   ~0008334

Last edited: 2017-07-04 20:35

 I've tried your patch and I'm now getting 0.107, which still isn't the same as the injected total (0.11).

The problem isn't just the 0.0001 overlap between tracks. There is also an issue in the collectParcelPolygon method. It is using a polygon's triangle decomposition to determine whether there is an intersection, and when the particle hits a diagonal, it can slip through without registering a hit.

Try the patch I've uploaded. I have removed the tolerance, which seems to get rid of the duplicates (though they might still be possible in other cases), and I've also changed the polygon intersection test so that it isn't affected by the diagonals. This gives the right result for me on the supplied case.

patch-1.will (3,104 bytes)

karlvirgil

2017-07-05 17:43

reporter   ~0008344

This patch looks good for this case. Let me test this out on a few other of my cases to make sure I'm still getting the right behavior. I'll get back to you.

karlvirgil

2017-07-05 21:06

reporter   ~0008345

I like your patch. However, here's another case, using your patch, that doesn't collect all the particles. It misses one (the last one to be injected, origId=29) because it lands on the collection plane (y=2m) and doesn't get collected. If I add the 1.0001 factor back in to ParticleCollector.C then it works again. I think what is needed is your fix (which addresses landing on the diagonal) and the 1.0001 factor with my fix to make sure that duplicates aren't counted. The combination of these two fixes works in my cases.
forWill.tgz (15,084 bytes)

karlvirgil

2017-07-05 21:10

reporter  

patch-2.karl (5,915 bytes)

karlvirgil

2017-07-05 21:11

reporter   ~0008346

Here are the combined patches (patch-2.karl)
patch-2-2.karl (5,915 bytes)

will

2017-07-06 13:28

manager   ~0008349

The case attached (forWill.tgz) works for me with just my fix. I'm getting sum(total mass) = 0.1, which is the same as the amount injected.

karlvirgil

2017-07-06 20:28

reporter   ~0008359

For some reason, I get different results from you. I get sum(total mass) = 0.0936, and two parcels are being missed. I've attached the log file and ParticleCollector.dat and ParticleCollector.C (files.tar). I'm using commit 1ff57870007ac317d95a1756b41fd76c5a1e8f26.
files.tar (347,648 bytes)

karlvirgil

2017-07-06 20:49

reporter   ~0008360

I'm wondering if the difference between yours and mine is simply different architectures/compilers/etc. It seems that finding the intersection is very sensitive when the particle endpoint is very close to the plane, and this was likely the reason for the 1.0001 factor to begin with. I think that with the two fixes combined the problem should be solved.

will

2017-07-10 09:34

manager   ~0008373

Maybe. I've tried gcc-4.8.5, though, which is the compiler on CentOS 7, and I still don't see the issue.

I'll push my change, as that's demonstrably better than what we have at present. I'm not going to add the tolerance and collected-particle arrays; these are just going to generate further issues down the line (i.e., parallelism and computational expense).

The method is fundamentally flawed. There is always likely to be some sort of tolerance issue with this sort of geometric test. We should think about how this should be done properly, rather than creating a fragile stack of ad-hoc fixes.

My advice would be that if particles are being missed or counted twice, then the collection surface should be moved slightly. Chances are that it lies in exactly the same plane as a number of faces, and moving it ever so slightly in the direction of it's normal will resolve the problem.

will

2017-07-10 10:24

manager   ~0008374

OK, I've pushed f3ef79b8295d5f2033287598b98fdf530dbcab2b, which is my fix. I appreciate that you may not consider this resolved, but on my machine this makes all the uploaded cases function correctly. Moreover, we now also have a usage recommendation which should resolve the issue; i.e., move the collection surface a small amount. Please reopen if you come across any new cases that aren't resolved by the fix or the recommendation.

Issue History

Date Modified Username Field Change
2017-06-28 18:09 karlvirgil New Issue
2017-06-28 18:09 karlvirgil File Added: case.tgz
2017-06-29 09:37 will Note Added: 0008277
2017-07-02 00:43 karlvirgil Note Added: 0008317
2017-07-03 08:10 karlvirgil File Added: patch
2017-07-03 08:10 karlvirgil Note Added: 0008321
2017-07-03 09:07 will Note Added: 0008322
2017-07-03 09:07 will Note Added: 0008323
2017-07-03 09:10 will Note Added: 0008324
2017-07-03 16:35 karlvirgil File Added: patch-2
2017-07-03 16:35 karlvirgil Note Added: 0008325
2017-07-04 20:35 will File Added: patch-1.will
2017-07-04 20:35 will Note Added: 0008334
2017-07-04 20:35 will Note Edited: 0008334
2017-07-05 17:43 karlvirgil Note Added: 0008344
2017-07-05 21:06 karlvirgil File Added: forWill.tgz
2017-07-05 21:06 karlvirgil Note Added: 0008345
2017-07-05 21:10 karlvirgil File Added: patch-2.karl
2017-07-05 21:11 karlvirgil File Added: patch-2-2.karl
2017-07-05 21:11 karlvirgil Note Added: 0008346
2017-07-06 13:28 will Note Added: 0008349
2017-07-06 20:28 karlvirgil File Added: files.tar
2017-07-06 20:28 karlvirgil Note Added: 0008359
2017-07-06 20:49 karlvirgil Note Added: 0008360
2017-07-10 09:34 will Note Added: 0008373
2017-07-10 10:24 will Assigned To => will
2017-07-10 10:24 will Status new => resolved
2017-07-10 10:24 will Resolution open => fixed
2017-07-10 10:24 will Fixed in Version => dev
2017-07-10 10:24 will Note Added: 0008374