From cce456fcaa61e73dd053833385e11745e50d0fa4 Mon Sep 17 00:00:00 2001 From: Will Bainbridge 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::move // Initialise the stepFraction moved for the particles forAllIter(typename Cloud, *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::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::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