2018-07-16 17:30 BST

View Issue Details Jump to Notes ]
IDProjectCategoryView StatusLast Update
0001147OpenFOAM[All Projects] Bugpublic2018-07-10 11:16
Reporteruser672 
Assigned To 
PrioritynormalSeverityminorReproducibilityalways
StatusnewResolutionopen 
PlatformLinuxOSScientific LinuxOS Version6.4
Product Version 
Target VersionFixed in Version 
Summary0001147: Sample utility: line passing through internal cyclic boundary using midPoint set sampling definition isn't possible
DescriptionNo error happens if "face" set sampling definition is chosen. Same error for either a sampling functionObject, or using the postProcessing sample utility.
Steps To ReproduceIn the attached minimal case I've set up a square baffle in the centre of a channel of square cross-section, and a sampling functionObject to output the velocity along two axial lines; one going through the baffle and one running past it. The line running through the baffle fails if midPoint is used.
Additional InformationSame behaviour observed for both GCC-4.6 and Clang-4.3 builds.
Tagsbaffle, sample
Attached Files
  • tgz file icon sampleCyclic.tgz (3,373 bytes) 2014-02-01 19:47
  • patch file icon bug1147.patch (12,741 bytes) 2015-03-29 21:29 -
    diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
    index c71559d..861b7ec 100644
    --- a/src/lagrangian/basic/particle/particle.C
    +++ b/src/lagrangian/basic/particle/particle.C
    @@ -61,7 +61,8 @@ Foam::particle::particle
         tetFaceI_(tetFaceI),
         tetPtI_(tetPtI),
         origProc_(Pstream::myProcNo()),
    -    origId_(getNewParticleID())
    +    origId_(getNewParticleID()),
    +    markedToTeleport_(false)
     {}
     
     
    @@ -81,7 +82,8 @@ Foam::particle::particle
         tetFaceI_(-1),
         tetPtI_(-1),
         origProc_(Pstream::myProcNo()),
    -    origId_(getNewParticleID())
    +    origId_(getNewParticleID()),
    +    markedToTeleport_(false)
     {
         if (doCellFacePt)
         {
    @@ -100,7 +102,8 @@ Foam::particle::particle(const particle& p)
         tetFaceI_(p.tetFaceI_),
         tetPtI_(p.tetPtI_),
         origProc_(p.origProc_),
    -    origId_(p.origId_)
    +    origId_(p.origId_),
    +    markedToTeleport_(false)
     {}
     
     
    @@ -114,7 +117,8 @@ Foam::particle::particle(const particle& p, const polyMesh& mesh)
         tetFaceI_(p.tetFaceI_),
         tetPtI_(p.tetPtI_),
         origProc_(p.origProc_),
    -    origId_(p.origId_)
    +    origId_(p.origId_),
    +    markedToTeleport_(false)
     {}
     
     
    diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
    index 6cb1838..9d9b099 100644
    --- a/src/lagrangian/basic/particle/particle.H
    +++ b/src/lagrangian/basic/particle/particle.H
    @@ -154,6 +154,8 @@ protected:
     
             //- Local particle id on originating processor
             label origId_;
    +        
    +        bool markedToTeleport_;
     
     
         // Private Member Functions
    @@ -457,6 +459,9 @@ public:
                 //- Is the particle on the boundary/(or outside the domain)?
                 inline bool onBoundary() const;
     
    +            inline bool needsTeleportHandled() const;
    +            inline void resetTeleportHandler();
    +            
                 //- Is this global face an internal face?
                 inline bool internalFace(const label faceI) const;
     
    @@ -511,7 +516,7 @@ public:
                 //  on entry 'stepFraction()' should be set to the fraction of the
                 //  time-step at which the tracking starts.
                 template<class TrackData>
    -            scalar trackToFace(const vector& endPosition, TrackData& td);
    +            scalar trackToFace(const vector& endPosition, TrackData& td, bool handleTeleport=false);
     
                 //- Return the index of the face to be used in the interpolation
                 //  routine
    diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H
    index 803a200..38a0d0b 100644
    --- a/src/lagrangian/basic/particle/particleI.H
    +++ b/src/lagrangian/basic/particle/particleI.H
    @@ -850,7 +850,19 @@ inline void Foam::particle::initCellFacePt()
     
     inline bool Foam::particle::onBoundary() const
     {
    -    return faceI_ != -1 && faceI_ >= mesh_.nInternalFaces();
    +    return faceI_ != -1 && faceI_ >= mesh_.nInternalFaces() && !markedToTeleport_;
    +}
    +
    +
    +inline bool Foam::particle::needsTeleportHandled() const
    +{
    +  return markedToTeleport_;
    +}
    +
    +
    +inline void Foam::particle::resetTeleportHandler()
    +{
    +  markedToTeleport_ = false;
     }
     
     
    diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C
    index 280165d..7e44e58 100644
    --- a/src/lagrangian/basic/particle/particleTemplates.C
    +++ b/src/lagrangian/basic/particle/particleTemplates.C
    @@ -203,7 +203,8 @@ template<class TrackData>
     Foam::scalar Foam::particle::trackToFace
     (
         const vector& endPosition,
    -    TrackData& td
    +    TrackData& td,
    +    bool handleTeleport
     )
     {
         typedef typename TrackData::cloudType cloudType;
    @@ -448,6 +449,9 @@ Foam::scalar Foam::particle::trackToFace
                 // This must be a cell face crossing
                 faceI_ = tetFaceI_;
     
    +            Info << "tetFaceI_: " << tetFaceI_ << endl;
    +            Info << "tri, lambdaMin: " << triI << " " << lambdaMin << endl;
    +
                 // Set the faceHitTetIs to those for the current tet in case a
                 // wall interaction is required with the cell face
                 faceHitTetIs = tetIndices
    @@ -576,14 +580,24 @@ Foam::scalar Foam::particle::trackToFace
                     patchI = patch(faceI_);
                 }
     
    +            Info << "origFaceI: " << origFaceI << endl;
    +            Info << "faceI_: " << faceI_ << endl;
    +            
                 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
     
                 if (isA<wedgePolyPatch>(patch))
                 {
    -                p.hitWedgePatch
    -                (
    -                    static_cast<const wedgePolyPatch&>(patch), td
    -                );
    +                if(!handleTeleport || markedToTeleport_)
    +                {
    +                    p.hitWedgePatch
    +                    (
    +                        static_cast<const wedgePolyPatch&>(patch), td
    +                    );
    +                }
    +                else
    +                {
    +                    markedToTeleport_ = true;
    +                }
                 }
                 else if (isA<symmetryPlanePolyPatch>(patch))
                 {
    @@ -601,26 +615,47 @@ Foam::scalar Foam::particle::trackToFace
                 }
                 else if (isA<cyclicPolyPatch>(patch))
                 {
    -                p.hitCyclicPatch
    -                (
    -                    static_cast<const cyclicPolyPatch&>(patch), td
    -                );
    +                if(!handleTeleport || markedToTeleport_)
    +                {
    +                  p.hitCyclicPatch
    +                  (
    +                      static_cast<const cyclicPolyPatch&>(patch), td
    +                  );
    +                }
    +                else
    +                {
    +                    markedToTeleport_ = true;
    +                }
                 }
                 else if (isA<cyclicAMIPolyPatch>(patch))
                 {
    -                p.hitCyclicAMIPatch
    -                (
    -                    static_cast<const cyclicAMIPolyPatch&>(patch),
    -                    td,
    -                    endPosition - position_
    -                );
    +                if(!handleTeleport || markedToTeleport_)
    +                {
    +                    p.hitCyclicAMIPatch
    +                    (
    +                        static_cast<const cyclicAMIPolyPatch&>(patch),
    +                        td,
    +                        endPosition - position_
    +                    );
    +                }
    +                else
    +                {
    +                    markedToTeleport_ = true;
    +                }
                 }
                 else if (isA<processorPolyPatch>(patch))
                 {
    -                p.hitProcessorPatch
    -                (
    -                    static_cast<const processorPolyPatch&>(patch), td
    -                );
    +                if(!handleTeleport || markedToTeleport_)
    +                {
    +                    p.hitProcessorPatch
    +                    (
    +                        static_cast<const processorPolyPatch&>(patch), td
    +                    );
    +                }
    +                else
    +                {
    +                    markedToTeleport_ = true;
    +                }
                 }
                 else if (isA<wallPolyPatch>(patch))
                 {
    @@ -633,9 +668,14 @@ Foam::scalar Foam::particle::trackToFace
                 {
                     p.hitPatch(patch, td);
                 }
    +
    +            Info << "faceI_ After: " << faceI_ << endl;
    +          
             }
         }
     
    +    Info << "A faceI_: " << faceI_ << endl;
    +    
         if (lambdaMin < SMALL)
         {
             // Apply tracking correction towards tet centre.
    @@ -708,6 +748,8 @@ Foam::scalar Foam::particle::trackToFace
             cloud.trackingRescue();
         }
     
    +    Info << "V faceI_: " << faceI_ << endl;
    +
         return trackFraction;
     }
     
    diff --git a/src/sampling/sampledSet/face/faceOnlySet.C b/src/sampling/sampledSet/face/faceOnlySet.C
    index 85fb7bc..8d0d19b 100644
    --- a/src/sampling/sampledSet/face/faceOnlySet.C
    +++ b/src/sampling/sampledSet/face/faceOnlySet.C
    @@ -60,14 +60,19 @@ bool Foam::faceOnlySet::trackToBoundary
     
         // Alias
         const point& trackPt = singleParticle.position();
    +    bool mustResetTeleport = false;
     
         while(true)
         {
             point oldPoint = trackPt;
     
    -        singleParticle.trackToFace(end_, trackData);
    -
    -        if (singleParticle.face() != -1 && mag(oldPoint - trackPt) > smallDist)
    +        singleParticle.trackToFace(end_, trackData, true);
    +        
    +        if (
    +            singleParticle.face() != -1 && 
    +            (mag(oldPoint - trackPt) > smallDist
    +            || singleParticle.needsTeleportHandled())
    +           )
             {
                 // Reached face. Sample.
                 samplingPts.append(trackPt);
    @@ -86,6 +91,18 @@ bool Foam::faceOnlySet::trackToBoundary
                 // Boundary reached.
                 return true;
             }
    +        if(singleParticle.needsTeleportHandled())
    +        {
    +            if(mustResetTeleport)
    +            {
    +                singleParticle.resetTeleportHandler();
    +                mustResetTeleport = false;
    +            }
    +            else
    +            {
    +                mustResetTeleport = true;
    +            }
    +        }
         }
     }
     
    diff --git a/src/sampling/sampledSet/midPoint/midPointSet.C b/src/sampling/sampledSet/midPoint/midPointSet.C
    index 168d0bb..7e76bfe 100644
    --- a/src/sampling/sampledSet/midPoint/midPointSet.C
    +++ b/src/sampling/sampledSet/midPoint/midPointSet.C
    @@ -62,11 +62,18 @@ void Foam::midPointSet::genSamples()
             {
                 midPoints[midI] =
                     0.5*(operator[](sampleI) + operator[](sampleI+1));
    +                
    +            Info << "faces_[sampleI]: " << faces_[sampleI] << endl;
    +            Info << "faces_[sampleI+1]: " << faces_[sampleI+1] << endl;
     
                 label cell1 = getCell(faces_[sampleI], midPoints[midI]);
                 label cell2 = getCell(faces_[sampleI+1], midPoints[midI]);
     
    -            if (cell1 != cell2)
    +            if (cell1 == -1 || cell2 == -1)
    +            {
    +               Info << "Ignoring intermediate jump" << endl;
    +            }
    +            else if (cell1 != cell2)
                 {
                     FatalErrorIn("midPointSet::genSamples()")
                         << "  sampleI:" << sampleI
    @@ -80,12 +87,14 @@ void Foam::midPointSet::genSamples()
                         << "  cell2:" << cell2
                         << abort(FatalError);
                 }
    +            else
    +            {
    +                midCells[midI] = cell1;
    +                midSegments[midI] = segments_[sampleI];
    +                midCurveDist[midI] = mag(midPoints[midI] - start());
    +                midI++;
    +            }
     
    -            midCells[midI] = cell1;
    -            midSegments[midI] = segments_[sampleI];
    -            midCurveDist[midI] = mag(midPoints[midI] - start());
    -
    -            midI++;
                 sampleI++;
             }
     
    diff --git a/src/sampling/sampledSet/sampledSet/sampledSet.C b/src/sampling/sampledSet/sampledSet/sampledSet.C
    index 1442ca3..2283610 100644
    --- a/src/sampling/sampledSet/sampledSet/sampledSet.C
    +++ b/src/sampling/sampledSet/sampledSet/sampledSet.C
    @@ -30,6 +30,11 @@ License
     #include "writer.H"
     #include "particle.H"
     
    +#include "cyclicPolyPatch.H"
    +#include "cyclicAMIPolyPatch.H"
    +#include "processorPolyPatch.H"
    +#include "wedgePolyPatch.H"
    +
     // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
     
     namespace Foam
    @@ -64,18 +69,33 @@ Foam::label Foam::sampledSet::getCell
                 << abort(FatalError);
         }
     
    -    if (faceI >= mesh().nInternalFaces())
    +    if (!mesh().isInternalFace(faceI))
         {
             label cellI = getBoundaryCell(faceI);
     
             if (!mesh().pointInCell(sample, cellI, searchEngine_.decompMode()))
             {
    -            FatalErrorIn
    -            (
    -                "sampledSet::getCell(const label, const point&)"
    -            )   << "Found cell " << cellI << " using face " << faceI
    -                << ". But cell does not contain point " << sample
    -                << abort(FatalError);
    +            label patchI = mesh_.boundaryMesh().whichPatch(faceI);
    +            const polyPatch& patch = mesh_.boundaryMesh()[patchI];
    +            if(
    +               isA<wedgePolyPatch>(patch) || 
    +               isA<cyclicPolyPatch>(patch) || 
    +               isA<cyclicAMIPolyPatch>(patch) ||
    +               isA<processorPolyPatch>(patch)
    +              )
    +            {
    +                return -1;
    +            }
    +            else
    +            {
    +
    +                FatalErrorIn
    +                (
    +                    "sampledSet::getCell(const label, const point&)"
    +                )   << "Found cell " << cellI << " using face " << faceI
    +                    << ". But cell does not contain point " << sample
    +                    << abort(FatalError);
    +            }
             }
             return cellI;
         }
    
    patch file icon bug1147.patch (12,741 bytes) 2015-03-29 21:29 +

-Relationships
+Relationships

-Notes

~0004549

wyldckat (updater)

@wallace: Many thanks for the simple and quick test case! It's really just a matter of following your instructions!


@Henry: Attached is the file "bug1147.patch", which one crazy hack of a patch. Seriously consider it only a proof of concept, because I couldn't code anything barely more elegant than that in a short time.
And sorry, but I left in a few stray debug lines, to assess if things were working as intended.

The explanation for the reported bug is simple: the method "particle::trackToFace()" always teleports the tracking point to the other side of the cyclic baffle, which after going through the all of the correct steps, it results in "midPointSet::genSamples()" crashing midway of "sampledSet::getCell()", because the second face is always on the other side of the cyclic patch pair, due to the teleport.


Hopefully the following step-by-step makes it a bit clearer, with the details on how to the provided patch tries to overcome this problem:

  1- The story begins in "faceOnlySet::trackToBoundary()", which is where the collection of faces crossed by the tracking line is gathered. The hack here was to use the new option to properly handle teleporting, where:

    a- In the first hit, it stores the first face before teleport, along with a flag indicating that a teleport needs to be performed.

    b- In the second hit, it executes the teleport, stores the face from the other side and resets the flags related to the need for handling a teleport.

    c- The tricky part is how "faceOnlySet::trackToBoundary()" and "particle::trackToFace()" handle the teleportation. This is because in "particle::trackToFace()" it has to postpone the jump to the second pass, which is why it has to rely so heavily on the flag "markedToTeleport_".

  2- Later on, in "midPointSet::genSamples()" we have to handle the new scenario of "cell1" or "cell2" being -1, which occurs when it's a situation of having the two faces be related to the same special pair of patches.

  3- In "sampledSet::getCell()", when it's not an internal face and if it fails to find the cell the point belongs to, then double-check if the "faceI" is part of one of the special patches that do particle teleportation... which if true, then return "-1".

~0004550

henry (reporter)

This is a lot of change for this particular purpose. I am also not happy about the increase in storage for ALL particle/parcel types even though it is only needed for this particular post-processing use. Actually this is a general problem with the Lagrangian library that all data which might be possibly used for any application/model is stored even if it is not needed for the particular case. My feeling is that we need more levels in the class hierachy to handle all the options, i.e. in this case we need a derived variant of particle which stores and handles the "teleportation" data appropriately. I think this will need to be delayed until we start writing the new robust tracking algorithm which handles bad cell which decompose with negative volume pyramids.
+Notes

-Issue History
Date Modified Username Field Change
2014-02-01 19:47 user672 New Issue
2014-02-01 19:47 user672 File Added: sampleCyclic.tgz
2014-02-05 00:03 user672 Tag Attached: Added the used Allrun script
2014-02-05 00:03 user672 Tag Attached: baffle
2014-02-05 00:03 user672 Tag Attached: functionObject
2014-02-05 00:03 user672 Tag Attached: sample
2014-02-05 00:03 user672 Tag Detached: functionObject
2015-03-29 21:29 wyldckat File Added: bug1147.patch
2015-03-29 21:53 wyldckat Note Added: 0004549
2015-03-29 22:03 henry Note Added: 0004550
2018-07-10 11:16 administrator Tag Detached: Added the used Allrun script
+Issue History