View Issue Details

IDProjectCategoryView StatusLast Update
0003779OpenFOAMBugpublic2022-01-14 15:08
ReportercgoessniAssigned Tohenry 
PrioritynormalSeverityblockReproducibilityrandom
Status closedResolutionunable to reproduce 
Platformamd64OSCentOSOS Version7
Product Version9 
Fixed in Version 
Summary0003779: infinite loop when runnig topoSet inside meshTools/cellClassification/cellClassification.C
DescriptionWhen I run topoSet for a complex geometry, sometimes it would not finish and is stuck inside a while(true) loop:

https://github.com/OpenFOAM/OpenFOAM-9/blob/master/src/meshTools/cellClassification/cellClassification.C#L201

I added some debug prints inside cellClassification.C (attached) and those lines are then repeated in the logfile:

faceTree.findLine = 1 (0.000250031 0.00754 -0.00465) 14085958
pHit = 1 (0.000250031 0.00754 -0.00465) 14085958
pt = (0.000250031 0.00754 -0.00465)
(pt-start) & edgeNormal = 0.000131969
edgeMag = 0.000764
faceTree.findLine = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pHit = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pt = (8.12987e-10 0.00754 -0.00465)
(pt-start) & edgeNormal = 0.000381999
edgeMag = 0.000764
faceTree.findLine = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pHit = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pt = (8.12987e-10 0.00754 -0.00465)
(pt-start) & edgeNormal = 0.000381999
edgeMag = 0.000764
faceTree.findLine = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pHit = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pt = (8.12987e-10 0.00754 -0.00465)
(pt-start) & edgeNormal = 0.000381999
edgeMag = 0.000764
faceTree.findLine = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pHit = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pt = (8.12987e-10 0.00754 -0.00465)
(pt-start) & edgeNormal = 0.000381999
edgeMag = 0.000764
faceTree.findLine = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pHit = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pt = (8.12987e-10 0.00754 -0.00465)
(pt-start) & edgeNormal = 0.000381999
edgeMag = 0.000764
faceTree.findLine = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pHit = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pt = (8.12987e-10 0.00754 -0.00465)
(pt-start) & edgeNormal = 0.000381999
edgeMag = 0.000764
faceTree.findLine = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pHit = 1 (8.13751e-10 0.00754 -0.00465) 14085958
pt = (8.12987e-10 0.00754 -0.00465)
(pt-start) & edgeNormal = 0.000381999
edgeMag = 0.000764

So it seems like the algorithm is stuck somewhere.

This can be avoided by specifying a maximum number of iterations that the while(true) loop must fullfil.
Steps To ReproduceHappens randomly for few case setups.
TagsNo tags attached.

Activities

cgoessni

2022-01-12 08:53

reporter  

cellClassification_workaround.diff.txt (370 bytes)
--- cellClassification.C	2022-01-12 09:12:31.978682898 +0100
+++ cellClassification.C.orig	2022-01-12 09:39:33.503254278 +0100
@@ -198,8 +198,7 @@
         // Current start of pierce test
         point pt = start;
 
-	label iter = 0;
-        while (true && ++iter < 1000)
+        while (true)
         {
             pointIndexHit pHit(faceTree.findLine(pt, end));
 
cellClassification_with_workarouund_and_prints.C (23,227 bytes)

cgoessni

2022-01-12 08:53

reporter   ~0012372

Ops, the diff is the wrong way around, sorry for that.

henry

2022-01-12 09:39

manager   ~0012373

> This can be avoided by specifying a maximum number of iterations that the while(true) loop must fullfil.

This sounds like a hack rather than a solution to the problem. Have you investigated the problem in detail? Can you provide a solution which avoids the need for the hack you propose?

cgoessni

2022-01-12 12:54

reporter   ~0012374

Yes, this is a hack and not intended to be committed.

I did some further analysis:
https://github.com/OpenFOAM/OpenFOAM-9/blob/master/src/meshTools/cellClassification/cellClassification.C#L221
Here, pt is updated using smallVec which is 1e-9 times the edge normal. For my case, the calculation in this line would then update pt to its original value, i.e.
221 pt = pHit.hitPoint() + smallVec;
with
pHit.hitPoint() = pt - smallVec for this specific point. This would then yield an infinite loop.

If I set smallVec to be 1e-10 of the edgeNormal, instead of 1e-9, the problem is not present anymore.
https://github.com/OpenFOAM/OpenFOAM-9/blob/master/src/meshTools/cellClassification/cellClassification.C#L194

henry

2022-01-12 13:05

manager   ~0012375

Unfortunately we don't have any cases which reproduce the problem so we are not in a position to analyse it. My guess is that the 1e-9 was set to that value for some particular case and changing it to 1e-10 may be OK for other cases but perhaps not all.

cgoessni

2022-01-12 13:12

reporter   ~0012376

To avoid such an infinite loop, one could multipe smallVec with a random number x element [0.5,1]

221 pt = pHit.hitPoint() + x * smallVec;

This should ensure for all cases, that the loop will terminate in finite time.

Because as you said, setting the factor from 1e-9 to 1e-10 resolves my specific case, but might not be suited for other cases.

henry

2022-01-12 13:37

manager   ~0012377

Is multiplying smallVec by 0.5 sufficient or do you need to use 1e-10? Using a random number in an already arbitrary tolerance test is "random", I can't see this being a good solution. As I say we don't have any cases which reproduce this problem and hence unable to analyse it.

cgoessni

2022-01-12 14:47

reporter   ~0012378

Using a random underrelaxation factor for smallVec would be a good solution IMHO.

The code already uses some arbitrary, but fixed increment vector (I guess that is what smallVec actually is). However, I have a case where is small increment vector is just the right size to play infinite ping pong in the while(true) loop:

in pseudo code

incr = 1e-9 * normalVector

while (true)
   hitPoint = findLine(pt, end) /* is pt - incr in my case */
   if not hit:
      pt = hitPoint + incr

Now, if we change the pre-factor for incr, there could still be cases which loop forever inside this loop. However, using a random underrelaxation factor for incr in each while loop iteration should make it impossible for this loop to go on forever, while the result should still be the same. Because if we apply incr now one time, or two times (if the underrelaxation factor is 0.5), should not matter for the outcome. Because the outcome of the whole loop is not a vector or a double, but a bool (setting cutFace).

Alternatively, another apporach could be to include a break statement inside the if (!cutFace[facei]) statement? Is the purpose of the while loop to find additional intersected faces? Could a break statement be included there somewhere?

Another approach would be to simply check the special condition that (hitPoint == pt - incr) and break then?

cgoessni

2022-01-12 14:48

reporter   ~0012379

sorry, pseudo-code should read:

while (true)
   hitPoint = findLine(pt, end) /* is pt - incr in my case */
   if hit:
      pt = hitPoint + incr

henry

2022-01-12 15:04

manager   ~0012380

Ideally we would want to replace the arbitrary (random or otherwise) geometric tolerance check with a robust topological test. This would be at least 2 days work (for which we would need to secure funding) and we would need at least one case which reproduces the problem. If you would be interested in funding work on this and can provide a test case please contact the OpenFOAM Foundation and consider the maintenance funding options: https://openfoam.org/news/funding-2022/

Issue History

Date Modified Username Field Change
2022-01-12 08:53 cgoessni New Issue
2022-01-12 08:53 cgoessni File Added: cellClassification_workaround.diff.txt
2022-01-12 08:53 cgoessni File Added: cellClassification_with_workarouund_and_prints.C
2022-01-12 08:53 cgoessni Note Added: 0012372
2022-01-12 09:39 henry Note Added: 0012373
2022-01-12 12:54 cgoessni Note Added: 0012374
2022-01-12 13:05 henry Note Added: 0012375
2022-01-12 13:12 cgoessni Note Added: 0012376
2022-01-12 13:37 henry Note Added: 0012377
2022-01-12 14:47 cgoessni Note Added: 0012378
2022-01-12 14:48 cgoessni Note Added: 0012379
2022-01-12 15:04 henry Note Added: 0012380
2022-01-14 15:08 henry Assigned To => henry
2022-01-14 15:08 henry Status new => closed
2022-01-14 15:08 henry Resolution open => unable to reproduce