View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0003056OpenFOAMBugpublic2018-08-27 19:232018-09-24 10:09
Reporterwyldckat Assigned Towill  
PrioritylowSeverityblockReproducibilitysometimes
Status resolvedResolutionfixed 
Fixed in Versiondev 
Summary0003056: Found an "Achilles heel" in the Barycentric implementation
DescriptionThis is a somewhat very specific corner case, which so far I've only been able to trigger this in the following conditions:

  1. If the mesh has cells that let the particles travel in reverse;
  2. And if a particle starts exactly from the centre of a cell;
  3. And if the particle happens to travel in a critical direction.

Then that particle gets stuck traveling back and forth between the same group of tets, blocking the solver in an infinite loop at 'Foam::particle::trackToFace'.

More specifically regarding condition 1, the provided test case demonstrates a situation where a polyhedral cell decomposes into tetrahedral cells in such a way that one of those tets has a negative determinant, which leads the particle getting stuck in an infinite travel loop between tets... essentially because it was meant to travel through that tet and instead is sent back to where it started from.

I have studied the code to try and fix this, but the only way I've been able to break the infinite loop is to count and limit for how long the particle stays in the same exact place while traveling between several tets in method 'Foam::particle::trackToFace', e.g. if it travels between 10 tets and never leaves the same position, then break the loop. Once it exits the method, it reports back that the tracking had to be canceled and is then handled accordingly by the caller method.

Another possible workaround was to try and change the precision of the position of the particle (because in the tests I've done I start from the centre of the cells), but it does not always work, since it can still get stuck when traveling in another direction.

I did not manage to understand if removing negative determinant tets from the search space would somewhat solve this issue, because when setting the debug flag for 'particle', we can see the tet-particles it travels in-between and doesn't seem to always land on the negative tets...
Steps To ReproduceThe attached package 'particleTrackingTestCase_v1.tar.gz' provides the case 'particleTrackingTestCase', which allows to reproduce the error, as well as demonstrates a workaround that happened to work. The mesh is a biopsy with 90 cells from a much larger mesh on which this was first found.
The large mesh works just fine with the existing flaws, both for fluid flow and conventional particle transport, but breaks the barycentric algorithm on this specific scenario (and with other directions).

Scripts for running the case:

  * To run the case in a way that gets the particle stuck in a travel loop, run the script 'Allrun.infiniteLoop'.
  * To clean up the case, run 'Allclean'.
  * To run the case without problems, run 'Allrun.runs'.

The only difference is that the write precision is set to 16 on the infinite loop and 8 for running without problems.

The work flow of the case is as follows:

  1. The flow field is as good as frozen, with flow velocity set to 0 on the whole field and with no pressure gradients.
  2. The solver 'pisoFoam' is used to run for 2 seconds, with deltaT set to 1, writing every 1 second.
  3. The function object 'icoUncoupledKinematicCloud' is used to transport the particles for just one time step.
  4. The particles are injected with the injector 'manualInjection', with the positions file set to 'C' and with a fixed velocity vector.
  5. The 'C' file is generated with "postProcess -func writeCellCentres" and then manipulated with the 'sed' command (done by the 'Allrun.*' scripts).

You can add the following block to 'controlDict' to get a log of where the particles travel through:

    DebugSwitches
    {
        particle 1;
    }

With or without that flag, the following warning is repeated several times until it silences it:

    No base point for face 122, 4(139 94 163 36), produces a valid tet decomposition.

which is likely related to the error given by 'checkMesh -allTopology -allGeometry':

    Error in face tets: 1 faces with low quality or negative volume decomposition tets.

and related to the problematic tet where the particle gets stuck to.
Additional InformationFirst reproduced this issue with OpenFOAM 5.x, but also occurs with the latest OpenFOAM 6 (b38715c4b19) and OpenFOAM dev (8ed92de98).

This was only possible to trigger because of the runs done with an in-house solver that does a sort of ray-tracing with OpenFOAM's particle transport algorithms. It was this way that we tripped over an old bug when running in parallel and then we provided the bug fix for it in report #2315 for the old transport algorithm.

Based on the information provided at https://cfd.direct/openfoam/free-software/barycentric-tracking/ - I was not expecting to trip over this bug :(
TagsNo tags attached.

Activities

wyldckat

2018-08-27 19:23

updater  

particleTrackingTestCase_v1.tar.gz (13,160 bytes)

will

2018-09-07 09:18

manager   ~0010056

> Then that particle gets stuck traveling back and forth between the same group
> of tets, blocking the solver in an infinite loop at
> 'Foam::particle::trackToFace'.

You have found the fundamental problem with tracking through negative space; that it opens the possibility of closed loops which never terminate.

> I have studied the code to try and fix this, but the only way I've been able
> to break the infinite loop is to count and limit for how long the particle
> stays in the same exact place while travelling between several tets in method
> 'Foam::particle::trackToFace',

I have avoided counters thus far as that sort of logic has masked all sorts of sins in the past, and also because these loops can cross between cyclics/AMI-s/processors/..., so the counter actually needs to be stored by the particle, which increases storage.

The hack that is in there at the moment is a factor which artificially decreases the (negative) amount of time that a particle spends in negative space. This means that as a particle moves around a loop the positive and negative space fractions don't cancel out and the track terminates. This isn't really any better than a counter in that the particle still ends up in the wrong place, but it does at least avoid additional storage.

Your case hangs despite this protection because you injected exactly at the cell centre, so the loop is of zero size, so even with the factor there is no difference in the time-accumulated in positive and negative space (they are both zero).

> I did not manage to understand if removing negative determinant tets from the
> search space would somewhat solve this issue, because when setting the debug
> flag for 'particle', we can see the tet-particles it travels in-between and
> doesn't seem to always land on the negative tets...

This isn't guaranteed to help. A particle can be injected into a positive tet and still subsequently enter one of these loops.

...

I'm happy to try a counter. That will fix your problem, and it should put a stop to these sorts of hangs (the time-hack doesn't really work on moving meshes either). The exact break condition will need a bit of care to write and locate within the various loops, if it is to catch all eventualities.

I would still like to come up with some means of fixing the deficit of time-not-tracked as a result of abandoned tracks, even if that means more storage. That issue can be considered separately from this one, though.

wyldckat

2018-09-07 19:21

updater   ~0010060

I still have to double-check what's going on when I run the case in parallel, because with the original case it locks up at the very first time step.
I'll try to diagnose the parallel case as soon as I can, to assess if the same issue is occurring and resulting in the prolonged/endless loop through a processor patch. Because if this is the situation, then that way we should have the two test cases needed for double-checking if the counter mechanism always works in serial and parallel modes.

will

2018-09-21 10:09

manager   ~0010065

Last edited: 2018-09-21 10:13

 This patch is my attempt at fixing this with a counter. It applies to dev 3a7b7619. I will run most of the Lagrangian tutorials to completion over the weekend to make sure that it doesn't introduce any other failure mechanisms, and I'll push it next week if all goes well. I thought the patch might be appreciated in the meantime, should you want run any more of your own tests.

You might notice that it wasn't enough just to add a counter. If it was a steady tracking case, or something with long time-steps, just cutting off after a fixed number of inverted tetrahedra might stop valid tracks prematurely. The logic I decided on was to cut off after a number of steps behind the maximum currently achieved step fraction. Once the particle achieves a new maximum step fraction the counter resets. This means it is only when we are not progressing do we stop the track.

Unfortunately, this meant adding storage for the maximum currently achieved step fraction. I can't think of another way which avoids this extra storage. I would be very pleased if someone else did.

0001-particle-Added-logic-to-break-closed-loops.patch (16,177 bytes)    
From cce456fcaa61e73dd053833385e11745e50d0fa4 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Thu, 20 Sep 2018 10:56:16 +0100
Subject: [PATCH] particle: Added logic to break closed loops

When a part of the tetrahedral decomposition is inverted, tracking along
a straight line can result in a closed loop which never ends.

This change adds a limit to the number of tracks that are done that end
before or at the maximum time already achieved. This breaks these closed
loops and prevents the simulation from hanging. The particles do,
however, end up in an incorrect position as a result of the tracking
being abandoned at an intermediate point in the step. A warning is
printed to indicate when this is occuring.

This resolves bug report https://bugs.openfoam.org/view.php?id=3056
---
 .../polyMeshTetDecomposition/tetIndices.C     |  7 ---
 .../polyMeshTetDecomposition/tetIndices.H     | 10 -----
 .../polyMeshTetDecomposition/tetIndicesI.H    | 16 +++----
 src/lagrangian/basic/Cloud/Cloud.C            |  2 +-
 src/lagrangian/basic/particle/particle.C      | 43 +++++++++++++------
 src/lagrangian/basic/particle/particle.H      | 41 ++++++++++++------
 src/lagrangian/basic/particle/particleI.H     | 30 ++++++++++---
 src/lagrangian/basic/particle/particleIO.C    |  8 +++-
 .../Templates/KinematicCloud/KinematicCloud.C |  2 +-
 .../KinematicParcel/KinematicParcel.C         |  2 +-
 10 files changed, 100 insertions(+), 61 deletions(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.C
index f2c7603b8..cd3232174 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.C
@@ -25,13 +25,6 @@ License
 
 #include "tetIndices.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-const Foam::label Foam::tetIndices::maxNWarnings = 100;
-
-Foam::label Foam::tetIndices::nWarnings = 0;
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::tetIndices::tetIndices()
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
index 6c7f0ecdf..4563af371 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
@@ -94,16 +94,6 @@ class tetIndices
         label tetPti_;
 
 
-    // Private static data
-
-        //- Maximum number of bad tet warnings
-        static const label maxNWarnings;
-
-        //- Current number of bad tet warnings. Warnings stop when this reaches
-        //  the maximum number.
-        static label nWarnings;
-
-
 public:
 
     // Constructors
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H
index e1a7cffb0..2d8e45841 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H
@@ -69,20 +69,18 @@ inline Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const
 
     if (faceBasePtI < 0)
     {
-        faceBasePtI = 0;
-        if (nWarnings < maxNWarnings)
+        static label badFacei = -1;
+
+        if (badFacei != face())
         {
             WarningInFunction
                 << "No base point for face " << face() << ", " << f
                 << ", produces a valid tet decomposition." << endl;
-            ++ nWarnings;
-        }
-        if (nWarnings == maxNWarnings)
-        {
-            Warning
-                << "Suppressing any further warnings." << endl;
-            ++ nWarnings;
+
+            badFacei = face();
         }
+
+        faceBasePtI = 0;
     }
 
     label facePtI = (tetPt() + faceBasePtI) % f.size();
diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 0cc8ac155..3215e5cc0 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -168,7 +168,7 @@ void Foam::Cloud<ParticleType>::move
     // Initialise the stepFraction moved for the particles
     forAllIter(typename Cloud<ParticleType>, *this, pIter)
     {
-        pIter().stepFraction() = 0;
+        pIter().setStepFraction(0);
     }
 
     // List of lists of particles to be transferred for all of the
diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index 79640c09a..8ff70a334 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.C
@@ -30,7 +30,7 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const Foam::scalar Foam::particle::negativeSpaceDisplacementFactor = 1.01;
+const Foam::label Foam::particle::maxNStepFractionsBelowMax_ = 10;
 
 Foam::label Foam::particle::particleCount_ = 0;
 
@@ -523,6 +523,8 @@ Foam::particle::particle
     tetPti_(tetPti),
     facei_(-1),
     stepFraction_(0.0),
+    maxStepFraction_(0.0),
+    nStepFractionsBelowMax_(0),
     origProc_(Pstream::myProcNo()),
     origId_(getNewParticleID())
 {}
@@ -542,6 +544,8 @@ Foam::particle::particle
     tetPti_(-1),
     facei_(-1),
     stepFraction_(0.0),
+    maxStepFraction_(0.0),
+    nStepFractionsBelowMax_(0),
     origProc_(Pstream::myProcNo()),
     origId_(getNewParticleID())
 {
@@ -564,6 +568,8 @@ Foam::particle::particle(const particle& p)
     tetPti_(p.tetPti_),
     facei_(p.facei_),
     stepFraction_(p.stepFraction_),
+    maxStepFraction_(p.maxStepFraction_),
+    nStepFractionsBelowMax_(p.nStepFractionsBelowMax_),
     origProc_(p.origProc_),
     origId_(p.origId_)
 {}
@@ -578,6 +584,8 @@ Foam::particle::particle(const particle& p, const polyMesh& mesh)
     tetPti_(p.tetPti_),
     facei_(p.facei_),
     stepFraction_(p.stepFraction_),
+    maxStepFraction_(p.maxStepFraction_),
+    nStepFractionsBelowMax_(p.nStepFractionsBelowMax_),
     origProc_(p.origProc_),
     origId_(p.origId_)
 {}
@@ -648,7 +656,8 @@ Foam::scalar Foam::particle::trackToFace
 
     facei_ = -1;
 
-    while (true)
+    // Loop the tets in the current cell
+    while (nStepFractionsBelowMax_ < maxNStepFractionsBelowMax_)
     {
         f *= trackToTri(f*displacement, f*fraction, tetTriI);
 
@@ -669,6 +678,22 @@ Foam::scalar Foam::particle::trackToFace
             changeTet(tetTriI);
         }
     }
+
+    // Warn if stuck, and incorrectly advance the step fraction to completion
+    static label stuckID = -1, stuckProc = -1;
+    if (origId_ != stuckID && origProc_ != stuckProc)
+    {
+        WarningInFunction
+            << "Particle #" << origId_ << " got stuck at " << position()
+            << endl;
+    }
+
+    stuckID = origId_;
+    stuckProc = origProc_;
+
+    addToStepFraction(f*fraction);
+
+    return 0;
 }
 
 
@@ -705,11 +730,8 @@ Foam::scalar Foam::particle::trackToStationaryTri
             << "Start local coordinates = " << y0 << endl;
     }
 
-    // Get the factor by which the displacement is increased
-    const scalar f = detA >= 0 ? 1 : negativeSpaceDisplacementFactor;
-
     // Calculate the local tracking displacement
-    barycentric Tx1(f*x1 & T);
+    barycentric Tx1(x1 & T);
 
     if (debug)
     {
@@ -777,7 +799,7 @@ Foam::scalar Foam::particle::trackToStationaryTri
     }
 
     // Set the proportion of the track that has been completed
-    stepFraction_ += fraction*muH*detA;
+    addToStepFraction(fraction*muH*detA);
 
     return iH != -1 ? 1 - muH*detA : 0;
 }
@@ -816,12 +838,9 @@ Foam::scalar Foam::particle::trackToMovingTri
             << "Start local coordinates = " << y0[0] << endl;
     }
 
-    // Get the factor by which the displacement is increased
-    const scalar f = detA[0] >= 0 ? 1 : negativeSpaceDisplacementFactor;
-
     // Get the relative global position
     const vector x0Rel = x0 - centre[0];
-    const vector x1Rel = f*x1 - centre[1];
+    const vector x1Rel = x1 - centre[1];
 
     // Form the determinant and hit equations
     cubicEqn detAEqn(sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
@@ -927,7 +946,7 @@ Foam::scalar Foam::particle::trackToMovingTri
     tetTriI = iH;
 
     // Set the proportion of the track that has been completed
-    stepFraction_ += fraction*muH*detA[0];
+    addToStepFraction(fraction*muH*detA[0]);
 
     if (debug)
     {
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index 5949f23a6..1b3ba2cc1 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -94,15 +94,9 @@ class particle
         //- Size in bytes of the fields
         static const std::size_t sizeofFields_;
 
-        //- The factor by which the displacement is increased when passing
-        //  through negative space. This should be slightly bigger than one.
-        //  This is needed as a straight trajectory can form a closed loop
-        //  through regions of overlapping positive and negative space, leading
-        //  to a track which never ends. By increasing the rate of displacement
-        //  through negative regions, the change in track fraction over this
-        //  part of the loop no longer exactly cancels the change over the
-        //  positive part, and the track therefore terminates.
-        static const scalar negativeSpaceDisplacementFactor;
+        //- The value of nStepFractionsBelowMax_ at which tracking is
+        //  abandoned. See the description of nStepFractionsBelowMax_.
+        static const label maxNStepFractionsBelowMax_;
 
 
 public:
@@ -155,6 +149,18 @@ private:
         //- Fraction of time-step completed
         scalar stepFraction_;
 
+        //- The maximum fraction of time-step completed so far
+        scalar maxStepFraction_;
+
+        //- The number of tracks carried out that ended in a step fraction
+        //  below the maximum. Once this reaches maxNStepFractionsBelowMax_,
+        //  tracking is abandoned for the current step. This is needed because
+        //  when tetrahedra are inverted a straight trajectory can form a closed
+        //  loop through regions of overlapping positive and negative space.
+        //  Without this break clause, such loops can result in a valid track
+        //  which never ends.
+        label nStepFractionsBelowMax_;
+
         //- Originating processor id
         label origProc_;
 
@@ -358,7 +364,9 @@ public:
         DefinePropertyList
         (
             "(coordinatesa coordinatesb coordinatesc coordinatesd) "
-            "celli tetFacei tetPti facei stepFraction origProc origId"
+            "celli tetFacei tetPti facei "
+            "stepFraction maxStepFraction nStepFractionsBelowMax "
+            "origProc origId"
         );
 
         //- Cumulative particle counter - used to provide unique ID
@@ -454,9 +462,6 @@ public:
             //- Return the fraction of time-step completed
             inline scalar stepFraction() const;
 
-            //- Return the fraction of time-step completed
-            inline scalar& stepFraction();
-
             //- Return the originating processor ID
             inline label origProc() const;
 
@@ -512,6 +517,16 @@ public:
 
     // Track
 
+        //- Set the value of the step fraction. Typically used to initialise at
+        //  the start of a tracking step. This is a set method, rather than an
+        //  accessor, as it also resets the maximum step fraction and the
+        //  number of tracks behind.
+        inline void setStepFraction(const scalar stepFraction);
+
+        //- Add a value to the step fraction. Not an accessor for the same
+        //  reasons as setStepFraction.
+        inline void addToStepFraction(const scalar dStepFraction);
+
         //- Track along the displacement for a given fraction of the overall
         //  step. End when the track is complete, or when a boundary is hit.
         //  On exit, stepFraction_ will have been incremented to the current
diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H
index 6d7878579..ae08c8f08 100644
--- a/src/lagrangian/basic/particle/particleI.H
+++ b/src/lagrangian/basic/particle/particleI.H
@@ -169,12 +169,6 @@ inline Foam::scalar Foam::particle::stepFraction() const
 }
 
 
-inline Foam::scalar& Foam::particle::stepFraction()
-{
-    return stepFraction_;
-}
-
-
 inline Foam::label Foam::particle::origProc() const
 {
     return origProc_;
@@ -286,6 +280,30 @@ inline Foam::vector Foam::particle::position() const
 }
 
 
+inline void Foam::particle::setStepFraction(const scalar stepFraction)
+{
+    stepFraction_ = stepFraction;
+    maxStepFraction_ = stepFraction;
+    nStepFractionsBelowMax_ = 0;
+}
+
+
+inline void Foam::particle::addToStepFraction(const scalar dStepFraction)
+{
+    stepFraction_ += dStepFraction;
+
+    if (stepFraction_ > maxStepFraction_)
+    {
+        maxStepFraction_ = stepFraction_;
+        nStepFractionsBelowMax_ = 0;
+    }
+    else
+    {
+        ++ nStepFractionsBelowMax_;
+    }
+}
+
+
 void Foam::particle::patchData(vector& normal, vector& displacement) const
 {
     if (!onBoundaryFace())
diff --git a/src/lagrangian/basic/particle/particleIO.C b/src/lagrangian/basic/particle/particleIO.C
index 6c99e0bd7..bf89392c0 100644
--- a/src/lagrangian/basic/particle/particleIO.C
+++ b/src/lagrangian/basic/particle/particleIO.C
@@ -52,6 +52,8 @@ Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields)
     tetPti_(-1),
     facei_(-1),
     stepFraction_(0.0),
+    maxStepFraction_(0.0),
+    nStepFractionsBelowMax_(0),
     origProc_(Pstream::myProcNo()),
     origId_(-1)
 {
@@ -61,7 +63,9 @@ Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields)
 
         if (readFields)
         {
-            is  >> facei_ >> stepFraction_ >> origProc_ >> origId_;
+            is  >> facei_
+                >> stepFraction_ >> maxStepFraction_ >> nStepFractionsBelowMax_
+                >> origProc_ >> origId_;
         }
     }
     else
@@ -110,6 +114,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const particle& p)
             << token::SPACE << p.tetPti_
             << token::SPACE << p.facei_
             << token::SPACE << p.stepFraction_
+            << token::SPACE << p.maxStepFraction_
+            << token::SPACE << p.nStepFractionsBelowMax_
             << token::SPACE << p.origProc_
             << token::SPACE << p.origId_;
     }
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 3b65e031e..d02677844 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -565,7 +565,7 @@ void Foam::KinematicCloud<CloudType>::checkParcelProperties
 )
 {
     const scalar carrierDt = mesh_.time().deltaTValue();
-    parcel.stepFraction() = (carrierDt - lagrangianDt)/carrierDt;
+    parcel.setStepFraction((carrierDt - lagrangianDt)/carrierDt);
 
     if (parcel.typeId() == -1)
     {
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 87b2b05bc..dade674c8 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -326,7 +326,7 @@ bool Foam::KinematicParcel<ParcelType>::move
             // the same relative to the face that it is on. The local
             // coordinates therefore do not change. We still advance in time and
             // perform the relevant interactions with the fixed particle.
-            p.stepFraction() += f;
+            p.addToStepFraction(f);
         }
 
         const scalar dt = (p.stepFraction() - sfrac)*trackTime;
-- 
2.18.0
0001-particle-Added-logic-to-break-closed-loops.patch (16,177 bytes)   

will

2018-09-24 10:09

manager   ~0010067

Fix pushed to dev as 1c26ed2d. The particleTrackingTestCase_v1.tar.gz case now runs, and will spit a warning to indicate that a particle got stuck.

The fix is roughly as outlined in my earlier commit, except that it is working with a maximum achieved distance, rather than the step fraction. Often we track without changing the step fraction, so this change was needed to fix those cases, too. All post-processing tracks are done this way, so that the mesh doesn't move. The pros and cons of the change haven't changed; the implementation is similar, and the memory usage is identical.

I also cleaned up the warnings a bit. Rather than printing N times and then going silent, they now print only when the particle or tet ID is different from last time. I think this is a bit more log-friendly.

Issue History

Date Modified Username Field Change
2018-08-27 19:23 wyldckat New Issue
2018-08-27 19:23 wyldckat File Added: particleTrackingTestCase_v1.tar.gz
2018-09-07 09:18 will Note Added: 0010056
2018-09-07 19:21 wyldckat Note Added: 0010060
2018-09-21 10:09 will File Added: 0001-particle-Added-logic-to-break-closed-loops.patch
2018-09-21 10:09 will Note Added: 0010065
2018-09-21 10:13 will Note Edited: 0010065
2018-09-24 10:09 will Assigned To => will
2018-09-24 10:09 will Status new => resolved
2018-09-24 10:09 will Resolution open => fixed
2018-09-24 10:09 will Fixed in Version => dev
2018-09-24 10:09 will Note Added: 0010067