Lagrangian: Removed sub-cycling and improved injection behaviour

Support for mesh sub-cycling has been removed from Lagrangian. Various Lagrangian processes already support sub-dividing the time-step. It is easier and more efficient to extend that support all the way to the high-level cloud objects, rather than to sub-cycle the mesh. This has additional benefits. The cloud now no longer needs to reset the step fraction at the start of every sub-motion. Injection can take advantage of this by modifying particles' step fractions in order to distribute them uniformly through the time-step. This is a simple and efficient method. Previously, injection would track the particles some distance after injection. This was more expensive to do, it resulted in spatial artefacts in the injected Lagrangian field, and it did not correctly handle interactions with patches or parallel transfers. The lack of injection tracking also means that particles injected through patches now start their simulation topologically connected to the face on which they are created. This means that they correctly trigger cloud functions associated with that face and/or patch. The injection models now also return barycentric coordinates directly, rather than the global position of injected particles. For methods that initially generate locations naturally in terms of barycentric coordinates, this prevents unnecessary and potentially fragile conversion from barycentric to Cartesian and back again. These are typically the methods that "fill" some sort of space; e.g., patch or cell-zone injections.
OpenFOAM · Dec 13, 2022 · 06e29c4 · 06e29c4
1 parent d020df4
commit 06e29c4
Show file tree

Hide file tree

Showing 108 changed files with 815 additions and 1,107 deletions.
diff --git a/src/functionObjects/field/nearWallFields/findCellParticle.C b/src/functionObjects/field/nearWallFields/findCellParticle.C
@@ -27,23 +27,6 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::findCellParticle::findCellParticle
-(
-    const polyMesh& mesh,
-    const barycentric& coordinates,
-    const label celli,
-    const label tetFacei,
-    const label tetPtI,
-    const vector& displacement,
-    const label data
-)
-:
-    particle(mesh, coordinates, celli, tetFacei, tetPtI),
-    displacement_(displacement),
-    data_(data)
-{}
-
-
 Foam::findCellParticle::findCellParticle
 (
     const polyMesh& mesh,
@@ -94,8 +77,7 @@ Foam::findCellParticle::findCellParticle(Istream& is, bool readFields)
 bool Foam::findCellParticle::move
 (
     Cloud<findCellParticle>& cloud,
-    trackingData& td,
-    const scalar maxTrackLen
+    trackingData& td
 )
 {
     td.keepParticle = true;

diff --git a/src/functionObjects/field/nearWallFields/findCellParticle.H b/src/functionObjects/field/nearWallFields/findCellParticle.H
@@ -116,18 +116,6 @@ public:
 
     // Constructors
 
-        //- Construct from components
-        findCellParticle
-        (
-            const polyMesh& mesh,
-            const barycentric& coordinates,
-            const label celli,
-            const label tetFacei,
-            const label tetPtI,
-            const vector& displacement,
-            const label data
-        );
-
         //- Construct from a position and a cell, searching for the rest of the
         //  required topology
         findCellParticle
@@ -185,7 +173,7 @@ public:
         // Tracking
 
             //- Track all particles to their end point
-            bool move(Cloud<findCellParticle>&, trackingData&, const scalar);
+            bool move(Cloud<findCellParticle>&, trackingData&);
 
             //- Overridable function to handle the particle hitting a wedge
             void hitWedgePatch(Cloud<findCellParticle>&, trackingData&);

diff --git a/src/functionObjects/field/nearWallFields/nearWallFields.C b/src/functionObjects/field/nearWallFields/nearWallFields.C
@@ -86,7 +86,7 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
                     mesh_,
                     patch.Cf()[patchFacei],
                     patch.faceCells()[patchFacei],
-                    - distance_*nf[patchFacei],
+                  - distance_*nf[patchFacei],
                     globalWalls.toGlobal(nPatchFaces) // passive data
                 )
             );
@@ -123,10 +123,6 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
     // Database to pass into findCellParticle::move
     findCellParticle::trackingData td(cloud, cellToWalls_, cellToSamples_);
 
-    // Track all particles to their end position.
-    scalar maxTrackLen = 2.0*mesh_.bounds().mag();
-
-
     // Debug: collect start points
     pointField start;
     if (debug)
@@ -141,7 +137,8 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
     }
 
 
-    cloud.move(cloud, td, maxTrackLen);
+    // Track
+    cloud.move(cloud, td);
 
 
     // Rework cell-to-globalpatchface into a map

diff --git a/src/functionObjects/field/streamlines/streamlines.C b/src/functionObjects/field/streamlines/streamlines.C
@@ -334,13 +334,13 @@ bool Foam::functionObjects::streamlines::write()
             initialParticles = particles;
         }
 
-        particles.move(particles, td, rootGreat);
+        particles.move(particles, td);
 
         if (trackDirection_ == trackDirection::both)
         {
             particles.IDLList<streamlinesParticle>::operator=(initialParticles);
             td.trackForward_ = !td.trackForward_;
-            particles.move(particles, td, rootGreat);
+            particles.move(particles, td);
         }
     }
 

diff --git a/src/functionObjects/field/streamlines/streamlinesParticle.C b/src/functionObjects/field/streamlines/streamlinesParticle.C
@@ -200,8 +200,7 @@ Foam::streamlinesParticle::streamlinesParticle
 bool Foam::streamlinesParticle::move
 (
     streamlinesCloud& cloud,
-    trackingData& td,
-    const scalar
+    trackingData& td
 )
 {
     td.keepParticle = true;

diff --git a/src/functionObjects/field/streamlines/streamlinesParticle.H b/src/functionObjects/field/streamlines/streamlinesParticle.H
@@ -236,7 +236,7 @@ public:
         // Tracking
 
             //- Track all particles to their end point
-            bool move(streamlinesCloud&, trackingData&, const scalar);
+            bool move(streamlinesCloud&, trackingData&);
 
             //- Overridable function to handle the particle hitting a wedge
             void hitWedgePatch(streamlinesCloud&, trackingData&);

diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
@@ -953,7 +953,7 @@ void Foam::DSMCCloud<ParcelType>::evolve()
     this->inflowBoundary().inflow();
 
     // Move the particles ballistically with their current velocities
-    Cloud<ParcelType>::move(*this, td, mesh_.time().deltaTValue());
+    Cloud<ParcelType>::move(*this, td);
 
     // Update cell occupancy
     buildCellOccupancy();

diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C
@@ -33,8 +33,7 @@ template<class TrackCloudType>
 bool Foam::DSMCParcel<ParcelType>::move
 (
     TrackCloudType& cloud,
-    trackingData& td,
-    const scalar trackTime
+    trackingData& td
 )
 {
     typename TrackCloudType::parcelType& p =
@@ -45,6 +44,8 @@ bool Foam::DSMCParcel<ParcelType>::move
 
     const polyMesh& mesh = cloud.pMesh();
 
+    const scalar trackTime = mesh.time().deltaTValue();
+
     // For reduced-D cases, the velocity used to track needs to be
     // constrained, but the actual U_ of the parcel must not be
     // altered or used, as it is altered by patch interactions an

diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
@@ -158,19 +158,6 @@ public:
 
     // Constructors
 
-        //- Construct from components
-        inline DSMCParcel
-        (
-            const polyMesh& mesh,
-            const barycentric& coordinates,
-            const label celli,
-            const label tetFacei,
-            const label tetPti,
-            const vector& U,
-            const scalar Ei,
-            const label typeId
-        );
-
         //- Construct from a position and a cell, searching for the rest of the
         //  required topology
         inline DSMCParcel
@@ -228,12 +215,7 @@ public:
 
                 //- Move the parcel
                 template<class TrackCloudType>
-                bool move
-                (
-                    TrackCloudType& cloud,
-                    trackingData& td,
-                    const scalar trackTime
-                );
+                bool move(TrackCloudType& cloud, trackingData& td);
 
 
         // Patch interactions

diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H
@@ -52,26 +52,6 @@ inline Foam::DSMCParcel<ParcelType>::constantProperties::constantProperties
 {}
 
 
-template<class ParcelType>
-inline Foam::DSMCParcel<ParcelType>::DSMCParcel
-(
-    const polyMesh& mesh,
-    const barycentric& coordinates,
-    const label celli,
-    const label tetFacei,
-    const label tetPti,
-    const vector& U,
-    const scalar Ei,
-    const label typeId
-)
-:
-    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
-    U_(U),
-    Ei_(Ei),
-    typeId_(typeId)
-{}
-
-
 template<class ParcelType>
 inline Foam::DSMCParcel<ParcelType>::DSMCParcel
 (

diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
@@ -196,7 +196,8 @@ Foam::Cloud<ParticleType>::Cloud
     patchNbrProc_(patchNbrProc(pMesh)),
     patchNbrProcPatch_(patchNbrProcPatch(pMesh)),
     patchNonConformalCyclicPatches_(patchNonConformalCyclicPatches(pMesh)),
-    globalPositionsPtr_()
+    globalPositionsPtr_(),
+    timeIndex_(-1)
 {
     // Ask for the tetBasePtIs and oldCellCentres to trigger all processors to
     // build them, otherwise, if some processors have no particles then there
@@ -256,27 +257,38 @@ void Foam::Cloud<ParticleType>::cloudReset(const Cloud<ParticleType>& c)
 }
 
 
+template<class ParticleType>
+void Foam::Cloud<ParticleType>::changeTimeStep()
+{
+    forAllIter(typename Cloud<ParticleType>, *this, pIter)
+    {
+        pIter().reset(0);
+    }
+
+    timeIndex_ = pMesh_.time().timeIndex();
+}
+
+
 template<class ParticleType>
 template<class TrackCloudType>
 void Foam::Cloud<ParticleType>::move
 (
     TrackCloudType& cloud,
-    typename ParticleType::trackingData& td,
-    const scalar trackTime
+    typename ParticleType::trackingData& td
 )
 {
+    // If the time has changed, modify the particles accordingly
+    if (timeIndex_ != pMesh_.time().timeIndex())
+    {
+        changeTimeStep();
+    }
+
     // Clear the global positions as these are about to change
     globalPositionsPtr_.clear();
 
     // Ensure rays are available for non conformal transfers
     storeRays();
 
-    // Initialise the stepFraction moved for the particles
-    forAllIter(typename Cloud<ParticleType>, *this, pIter)
-    {
-        pIter().reset(0);
-    }
-
     // Create transfer buffers
     PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
 
@@ -300,7 +312,7 @@ void Foam::Cloud<ParticleType>::move
             ParticleType& p = pIter();
 
             // Move the particle
-            bool keepParticle = p.move(cloud, td, trackTime);
+            const bool keepParticle = p.move(cloud, td);
 
             // If the particle is to be kept
             if (keepParticle)

diff --git a/src/lagrangian/basic/Cloud/Cloud.H b/src/lagrangian/basic/Cloud/Cloud.H
@@ -91,6 +91,9 @@ class Cloud
         //- Temporary storage for the global particle positions
         mutable autoPtr<vectorField> globalPositionsPtr_;
 
+        //- Time index
+        mutable label timeIndex_;
+
 
     // Private Member Functions
 
@@ -243,13 +246,16 @@ public:
             //- Reset the particles
             void cloudReset(const Cloud<ParticleType>& c);
 
+            //- Change the particles' state from the end of the previous time
+            //  step to the start of the next time step
+            void changeTimeStep();
+
             //- Move the particles
             template<class TrackCloudType>
             void move
             (
                 TrackCloudType& cloud,
-                typename ParticleType::trackingData& td,
-                const scalar trackTime
+                typename ParticleType::trackingData& td
             );
 
 

diff --git a/src/lagrangian/basic/indexedParticle/indexedParticle.H b/src/lagrangian/basic/indexedParticle/indexedParticle.H
@@ -61,21 +61,6 @@ public:
 
     // Constructors
 
-        //- Construct from components
-        indexedParticle
-        (
-            const polyMesh& mesh,
-            const barycentric& coordinates,
-            const label celli,
-            const label tetFacei,
-            const label tetPti,
-            const label index = 0
-        )
-        :
-            particle(mesh, coordinates, celli, tetFacei, tetPti),
-            index_(index)
-        {}
-
         //- Construct from Istream
         indexedParticle(Istream& is, bool readFields = true)
         :

diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
@@ -480,14 +480,15 @@ Foam::particle::particle
     const barycentric& coordinates,
     const label celli,
     const label tetFacei,
-    const label tetPti
+    const label tetPti,
+    const label facei
 )
 :
     coordinates_(coordinates),
     celli_(celli),
     tetFacei_(tetFacei),
     tetPti_(tetPti),
-    facei_(-1),
+    facei_(facei),
     stepFraction_(1),
     stepFractionBehind_(0),
     nTracksBehind_(0),

diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
@@ -339,7 +339,8 @@ public:
             const barycentric& coordinates,
             const label celli,
             const label tetFacei,
-            const label tetPti
+            const label tetPti,
+            const label facei
         );
 
         //- Construct from a position and a cell, searching for the rest of the
@@ -397,6 +398,9 @@ public:
             //- Return current face particle is on otherwise -1
             inline label face() const;
 
+            //- Return current face particle is on otherwise -1
+            inline label& face();
+
             //- Return the fraction of time-step completed
             inline scalar stepFraction() const;
 
@@ -418,17 +422,6 @@ public:
 
         // Check
 
-            //- Return the step fraction change within the overall time-step.
-            //  Returns the start value and the change as a scalar pair. Always
-            //  return Pair<scalar>(0, 1), unless sub-cycling is in effect, in
-            //  which case the values will reflect the span of the sub-cycle
-            //  within the time-step.
-            inline Pair<scalar> stepFractionSpan(const polyMesh& mesh) const;
-
-            //- Return the current fraction within the timestep. This differs
-            //  from the stored step fraction due to sub-cycling.
-            inline scalar currentTimeFraction(const polyMesh& mesh) const;
-
             //- Return the indices of the current tet that the
             //  particle occupies.
             inline tetIndices currentTetIndices(const polyMesh& mesh) const;