From 9b57ef06cf9cce6ecb3e80a2df71675fb144f43a Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Wed, 9 Aug 2017 15:52:33 +0100 Subject: [PATCH] Lagrangian: Enabled tracking through ACMI patches and minor code improvements Particle collisions with ACMI patches are now handled. The hit detects whether the location is within the overlap or the coupled region and recurses, calling the hit routine appropriate for the region. The low level tracking methods are now more consistently named. There is now a distinction between tracking to a face and hitting it. Function object side effects have been moved out of the base layer and into the parcels on which they are meaningful. --- .../field/nearWallFields/findCellParticle.C | 2 +- .../field/streamLine/streamLineParticle.C | 2 +- .../parcels/Templates/DSMCParcel/DSMCParcel.C | 2 +- src/lagrangian/basic/particle/particle.C | 44 ++++++---- src/lagrangian/basic/particle/particle.H | 33 +++++-- .../basic/particle/particleTemplates.C | 88 +++++++++++++++---- .../KinematicParcel/KinematicParcel.C | 26 ++---- .../KinematicParcel/KinematicParcel.H | 8 -- .../molecule/molecule/molecule.C | 2 +- src/lagrangian/solidParticle/solidParticle.C | 2 +- .../trackedParticle/trackedParticle.C | 3 +- src/sampling/sampledSet/face/faceOnlySet.C | 2 +- 12 files changed, 137 insertions(+), 77 deletions(-) diff --git a/src/functionObjects/field/nearWallFields/findCellParticle.C b/src/functionObjects/field/nearWallFields/findCellParticle.C index e8efde1166..d53a904f3b 100644 --- a/src/functionObjects/field/nearWallFields/findCellParticle.C +++ b/src/functionObjects/field/nearWallFields/findCellParticle.C @@ -105,7 +105,7 @@ bool Foam::findCellParticle::move while (td.keepParticle && !td.switchProcessor && stepFraction() < 1) { const scalar f = 1 - stepFraction(); - trackToFace(f*(end_ - start_), f, td); + trackToAndHitFace(f*(end_ - start_), f, td); } if (stepFraction() == 1 || !td.keepParticle) diff --git a/src/functionObjects/field/streamLine/streamLineParticle.C b/src/functionObjects/field/streamLine/streamLineParticle.C index c587a7279c..58dd0b4030 100644 --- a/src/functionObjects/field/streamLine/streamLineParticle.C +++ b/src/functionObjects/field/streamLine/streamLineParticle.C @@ -198,7 +198,7 @@ bool Foam::streamLineParticle::move dt = maxDt; } - trackToFace(dt*U, 0, td); + trackToAndHitFace(dt*U, 0, td); if ( diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C index 2265ca5f84..7df06d3b98 100644 --- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C +++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C @@ -59,7 +59,7 @@ bool Foam::DSMCParcel::move(TrackData& td, const scalar trackTime) meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking); const scalar f = 1 - p.stepFraction(); - p.trackToFace(f*trackTime*Utracking, f, td); + p.trackToAndHitFace(f*trackTime*Utracking, f, td); if (p.onBoundaryFace() && td.keepParticle) { diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C index 2cef6b8f8b..e0d2d8df9a 100644 --- a/src/lagrangian/basic/particle/particle.C +++ b/src/lagrangian/basic/particle/particle.C @@ -344,6 +344,7 @@ void Foam::particle::changeFace(const label tetTriI) { const label newFaceI = mesh_.cells()[celli_][cellFaceI]; const class face& newFace = mesh_.faces()[newFaceI]; + const label newOwner = mesh_.faceOwner()[newFaceI]; // Exclude the current face if (tetFacei_ == newFaceI) @@ -351,21 +352,21 @@ void Foam::particle::changeFace(const label tetTriI) continue; } - // Loop over the edges, looking for the shared one - label edgeComp = 0; + // Loop over the edges, looking for the shared one. Note that we have to + // match the direction of the edge as well as the end points in order to + // avoid false positives when dealing with coincident ACMI faces. + const label edgeComp = newOwner == celli_ ? -1 : +1; label edgeI = 0; - for (; edgeI < newFace.size(); ++ edgeI) - { - edgeComp = edge::compare(sharedEdge, newFace.faceEdge(edgeI)); - - if (edgeComp != 0) - { - break; - } - } + for + ( + ; + edgeI < newFace.size() + && edge::compare(sharedEdge, newFace.faceEdge(edgeI)) != edgeComp; + ++ edgeI + ); // If the face does not contain the edge, then move on to the next face - if (edgeComp == 0) + if (edgeI >= newFace.size()) { continue; } @@ -643,6 +644,17 @@ Foam::scalar Foam::particle::trackToFace { // The track has hit a face, so set the current face and return facei_ = tetFacei_; + + /* + if (!mesh_.isInternalFace(facei_)) + { + const label patchi = mesh_.boundaryMesh().whichPatch(facei_); + const polyPatch& patch = mesh_.boundaryMesh()[patchi]; + Info<< "Tet face is on " << patch.type() << " patch #" << patchi + << endl << endl; + } + */ + return f; } else @@ -667,8 +679,8 @@ Foam::scalar Foam::particle::trackToStationaryTri if (debug) { - Info<< "Particle " << origId() << endl - << "Tracking from " << x0 << " to " << x0 + x1 << endl; + Info<< "Particle " << origId() << endl << "Tracking from " << x0 + << " along " << x1 << " to " << x0 + x1 << endl; } // Get the tet geometry @@ -778,8 +790,8 @@ Foam::scalar Foam::particle::trackToMovingTri if (debug) { - Info<< "Particle " << origId() << endl - << "Tracking from " << x0 << " to " << x0 + x1 << endl; + Info<< "Particle " << origId() << endl << "Tracking from " << x0 + << " along " << x1 << " to " << x0 + x1 << endl; } // Get the tet geometry diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H index 0c21308758..3a6f8213ef 100644 --- a/src/lagrangian/basic/particle/particle.H +++ b/src/lagrangian/basic/particle/particle.H @@ -57,6 +57,8 @@ class particle; class polyPatch; class cyclicPolyPatch; +class cyclicAMIPolyPatch; +class cyclicACMIPolyPatch; class processorPolyPatch; class symmetryPlanePolyPatch; class symmetryPolyPatch; @@ -294,10 +296,6 @@ protected: // Patch interactions - //- Overridable function to handle the particle hitting a face - template - void hitFace(TrackData& td); - //- Overridable function to handle the particle hitting a // patch. Executed before other patch-hitting functions. // trackFraction is passed in to allow mesh motion to @@ -343,6 +341,16 @@ protected: const vector& direction ); + //- Overridable function to handle the particle hitting a + // cyclicACMIPatch + template + void hitCyclicACMIPatch + ( + const cyclicACMIPolyPatch&, + TrackData& td, + const vector& direction + ); + //- Overridable function to handle the particle hitting a // processorPatch template @@ -585,12 +593,21 @@ public: label& tetTriI ); - //- As non-templated particle::trackToFace, but with additional - // boundary handling. + //- Hit the current face. If the current face is internal than this + // crosses into the next cell. If it is a boundary face then this will + // interact the particle with the relevant patch. template - void trackToFace + void hitFace ( - const vector& displacement, + const vector& direction, + TrackData& td + ); + + //- Convenience function. Cobines trackToFace and hitFace + template + void trackToAndHitFace + ( + const vector& direction, const scalar fraction, TrackData& td ); diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C index 27291a03e8..307dcafff4 100644 --- a/src/lagrangian/basic/particle/particleTemplates.C +++ b/src/lagrangian/basic/particle/particleTemplates.C @@ -27,6 +27,7 @@ License #include "cyclicPolyPatch.H" #include "cyclicAMIPolyPatch.H" +#include "cyclicACMIPolyPatch.H" #include "processorPolyPatch.H" #include "symmetryPlanePolyPatch.H" #include "symmetryPolyPatch.H" @@ -186,38 +187,30 @@ void Foam::particle::writeObjects(const CloudType& c, objectRegistry& obr) template -void Foam::particle::trackToFace +void Foam::particle::hitFace ( - const vector& displacement, - const scalar fraction, + const vector& direction, TrackData& td ) { - // Track - trackToFace(displacement, fraction); - - // If the track is complete, return if (!onFace()) { return; } - - // Hit face/patch processing - typedef typename TrackData::cloudType::particleType particleType; - particleType& p = static_cast(*this); - p.hitFace(td); - if (onInternalFace()) + else if (onInternalFace()) { changeCell(); } - else + else if (onBoundaryFace()) { - label origFacei = facei_; - label patchi = mesh_.boundaryMesh().whichPatch(facei_); + typename TrackData::cloudType::particleType& p = + static_cast(*this); // No action is taken for tetPti_ for tetFacei_ here. These are handled // by the patch interaction call or later during processor transfer. + const label origFacei = facei_; + label patchi = mesh_.boundaryMesh().whichPatch(facei_); const tetIndices faceHitTetIs(celli_, tetFacei_, tetPti_); if @@ -268,13 +261,22 @@ void Foam::particle::trackToFace static_cast(patch), td ); } + else if (isA(patch)) + { + p.hitCyclicACMIPatch + ( + static_cast(patch), + td, + direction + ); + } else if (isA(patch)) { p.hitCyclicAMIPatch ( static_cast(patch), td, - displacement + direction ); } else if (isA(patch)) @@ -301,8 +303,16 @@ void Foam::particle::trackToFace template -void Foam::particle::hitFace(TrackData&) -{} +void Foam::particle::trackToAndHitFace +( + const vector& direction, + const scalar fraction, + TrackData& td +) +{ + trackToFace(direction, fraction); + hitFace(direction, td); +} template @@ -478,6 +488,46 @@ void Foam::particle::hitCyclicAMIPatch } +template +void Foam::particle::hitCyclicACMIPatch +( + const cyclicACMIPolyPatch& cpp, + TrackData& td, + const vector& direction +) +{ + const label localFacei = cpp.whichFace(facei_); + + // If the mask is within the patch tolerance at either end, then we can + // assume an interaction with the appropriate part of the ACMI pair. + const scalar mask = cpp.mask()[localFacei]; + bool couple = mask >= 1 - cpp.tolerance(); + bool nonOverlap = mask <= cpp.tolerance(); + + // If the mask is an intermediate value, then we search for a location on + // the other side of the AMI. If we can't find a location, then we assume + // that we have hit the non-overlap patch. + if (!couple && !nonOverlap) + { + vector pos = position(); + couple = cpp.pointFace(localFacei, direction, pos) >= 0; + nonOverlap = !couple; + } + + if (couple) + { + hitCyclicAMIPatch(cpp, td, direction); + } + else + { + // Move to the face associated with the non-overlap patch and redo the + // face interaction. + tetFacei_ = facei_ = cpp.nonOverlapPatch().start() + localFacei; + hitFace(direction, td); + } +} + + template void Foam::particle::hitProcessorPatch(const processorPolyPatch&, TrackData&) {} diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C index 35f2021187..037ada186c 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -292,8 +292,8 @@ bool Foam::KinematicParcel::move f = min(f, maxCo*l/max(SMALL*l, mag(s))); if (p.active()) { - // Track to the next face - p.trackToFace(f*s, f, td); + // Track to and hit the next face + p.trackToAndHitFace(f*s, f, td); } else { @@ -305,7 +305,6 @@ bool Foam::KinematicParcel::move p.stepFraction() += f; } - const scalar dt = (p.stepFraction() - sfrac)*trackTime; // Avoid problems with extremely small timesteps @@ -332,6 +331,11 @@ bool Foam::KinematicParcel::move p.age() += dt; + if (p.onFace()) + { + td.cloud().functions().postFace(p, p.face(), td.keepParticle); + } + td.cloud().functions().postMove(p, celli, dt, start, td.keepParticle); } @@ -339,22 +343,6 @@ bool Foam::KinematicParcel::move } -template -template -void Foam::KinematicParcel::hitFace(TrackData& td) -{ - typename TrackData::cloudType::parcelType& p = - static_cast(*this); - - td.cloud().functions().postFace(p, p.face(), td.keepParticle); -} - - -template -void Foam::KinematicParcel::hitFace(int& td) -{} - - template template bool Foam::KinematicParcel::hitPatch diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H index a862c4f033..fbdd578e27 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H @@ -597,14 +597,6 @@ public: // Patch interactions - //- Overridable function to handle the particle hitting a face - // without trackData - void hitFace(int& td); - - //- Overridable function to handle the particle hitting a face - template - void hitFace(TrackData& td); - //- Overridable function to handle the particle hitting a patch // Executed before other patch-hitting functions template diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C index f565adbfac..ef11eaddaa 100644 --- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C +++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C @@ -88,7 +88,7 @@ bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime) while (td.keepParticle && !td.switchProcessor && stepFraction() < 1) { const scalar f = 1 - stepFraction(); - trackToFace(f*trackTime*v_, f, td); + trackToAndHitFace(f*trackTime*v_, f, td); } } else if (td.part() == 2) diff --git a/src/lagrangian/solidParticle/solidParticle.C b/src/lagrangian/solidParticle/solidParticle.C index 4b4b8034b9..80737165ec 100644 --- a/src/lagrangian/solidParticle/solidParticle.C +++ b/src/lagrangian/solidParticle/solidParticle.C @@ -58,7 +58,7 @@ bool Foam::solidParticle::move const scalar sfrac = stepFraction(); const scalar f = 1 - stepFraction(); - trackToFace(f*trackTime*U_, f, td); + trackToAndHitFace(f*trackTime*U_, f, td); const scalar dt = (stepFraction() - sfrac)*trackTime; diff --git a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C index 71f44569be..46c13d83a1 100644 --- a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C +++ b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C @@ -139,7 +139,8 @@ bool Foam::trackedParticle::move td.maxLevel_[cell()] = max(td.maxLevel_[cell()], level_); const scalar f = 1 - stepFraction(); - trackToFace(f*(end_ - start_), f, td); + const vector s = end_ - start_; + trackToAndHitFace(f*s, f, td); } } diff --git a/src/sampling/sampledSet/face/faceOnlySet.C b/src/sampling/sampledSet/face/faceOnlySet.C index 429d1c084b..5ce5eaf85c 100644 --- a/src/sampling/sampledSet/face/faceOnlySet.C +++ b/src/sampling/sampledSet/face/faceOnlySet.C @@ -62,7 +62,7 @@ bool Foam::faceOnlySet::trackToBoundary { point oldPoint = trackPt; - singleParticle.trackToFace(end_ - start_, 0, trackData); + singleParticle.trackToAndHitFace(end_ - start_, 0, trackData); trackPt = singleParticle.position();