Skip to content
Permalink
Browse files

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 distance 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
  • Loading branch information...
Will Bainbridge
Will Bainbridge committed Sep 20, 2018
1 parent a098cdb commit 1c26ed2d31faf6858a274fae8d6838a568cbefaa
@@ -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()
@@ -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
@@ -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();
@@ -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().reset();
}

// List of lists of particles to be transferred for all of the
@@ -30,7 +30,7 @@ License

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

const Foam::scalar Foam::particle::negativeSpaceDisplacementFactor = 1.01;
const Foam::label Foam::particle::maxNBehind_ = 10;

Foam::label Foam::particle::particleCount_ = 0;

@@ -523,6 +523,8 @@ Foam::particle::particle
tetPti_(tetPti),
facei_(-1),
stepFraction_(0.0),
behind_(0.0),
nBehind_(0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{}
@@ -542,6 +544,8 @@ Foam::particle::particle
tetPti_(-1),
facei_(-1),
stepFraction_(0.0),
behind_(0.0),
nBehind_(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_),
behind_(p.behind_),
nBehind_(p.nBehind_),
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_),
behind_(p.behind_),
nBehind_(p.nBehind_),
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 (nBehind_ < maxNBehind_)
{
f *= trackToTri(f*displacement, f*fraction, tetTriI);

@@ -669,6 +678,25 @@ 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_;

stepFraction_ += f*fraction;

behind_ = 0;
nBehind_ = 0;

return 0;
}


@@ -705,11 +733,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)
{
@@ -779,6 +804,22 @@ Foam::scalar Foam::particle::trackToStationaryTri
// Set the proportion of the track that has been completed
stepFraction_ += fraction*muH*detA;

// Accumulate displacement behind
if (detA <= 0 || nBehind_ > 0)
{
behind_ += muH*detA*mag(displacement);

if (behind_ > 0)
{
behind_ = 0;
nBehind_ = 0;
}
else
{
++ nBehind_;
}
}

return iH != -1 ? 1 - muH*detA : 0;
}

@@ -816,12 +857,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);
@@ -929,6 +967,22 @@ Foam::scalar Foam::particle::trackToMovingTri
// Set the proportion of the track that has been completed
stepFraction_ += fraction*muH*detA[0];

// Accumulate displacement behind
if (detA[0] <= 0 || nBehind_ > 0)
{
behind_ += muH*detA[0]*mag(displacement);

if (behind_ > 0)
{
behind_ = 0;
nBehind_ = 0;
}
else
{
++ nBehind_;
}
}

if (debug)
{
if (iH != -1)
@@ -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 nBehind_ at which tracking is abandoned. See the
// description of nBehind_.
static const label maxNBehind_;


public:
@@ -155,6 +149,18 @@ private:
//- Fraction of time-step completed
scalar stepFraction_;

//- The distance behind the maximum distance reached so far
scalar behind_;

//- The number of tracks carried out that ended in a distance behind the
// maximum distance reached so far. Once this reaches maxNBehind_,
// 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 nBehind_;

//- Originating processor id
label origProc_;

@@ -358,7 +364,8 @@ public:
DefinePropertyList
(
"(coordinatesa coordinatesb coordinatesc coordinatesd) "
"celli tetFacei tetPti facei stepFraction origProc origId"
"celli tetFacei tetPti facei stepFraction "
"behind nBehind origProc origId"
);

//- Cumulative particle counter - used to provide unique ID
@@ -512,6 +519,10 @@ public:

// Track

//- Set step fraction and behind data to zero in preparation for a new
// track
inline void reset();

//- 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
@@ -286,6 +286,14 @@ inline Foam::vector Foam::particle::position() const
}


inline void Foam::particle::reset()
{
stepFraction_ = 0;
nBehind_ = 0;
behind_ = 0;
}


void Foam::particle::patchData(vector& normal, vector& displacement) const
{
if (!onBoundaryFace())
@@ -52,6 +52,8 @@ Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields)
tetPti_(-1),
facei_(-1),
stepFraction_(0.0),
behind_(0.0),
nBehind_(0),
origProc_(Pstream::myProcNo()),
origId_(-1)
{
@@ -61,7 +63,8 @@ Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields)

if (readFields)
{
is >> facei_ >> stepFraction_ >> origProc_ >> origId_;
is >> facei_ >> stepFraction_ >> behind_ >> nBehind_
>> origProc_ >> origId_;
}
}
else
@@ -110,6 +113,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.behind_
<< token::SPACE << p.nBehind_
<< token::SPACE << p.origProc_
<< token::SPACE << p.origId_;
}

0 comments on commit 1c26ed2

Please sign in to comment.
You can’t perform that action at this time.