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.

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,
+    trackingData& td
-    const scalar maxTrackLen
 )
 {
     td.keepParticle = true;

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&);

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

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);
         }
     }
 

src/functionObjects/field/streamlines/streamlinesParticle.C

@@ -200,8 +200,7 @@ Foam::streamlinesParticle::streamlinesParticle
 bool Foam::streamlinesParticle::move
 (
     streamlinesCloud& cloud,
-    trackingData& td,
+    trackingData& td
-    const scalar
 )
 {
     td.keepParticle = true;

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&);

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();

src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C

@@ -33,8 +33,7 @@ template<class TrackCloudType>
 bool Foam::DSMCParcel<ParcelType>::move
 (
     TrackCloudType& cloud,
-    trackingData& td,
+    trackingData& td
-    const scalar trackTime
 )
 {
     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

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
+                bool move(TrackCloudType& cloud, trackingData& td);
-                (
-                    TrackCloudType& cloud,
-                    trackingData& td,
-                    const scalar trackTime
-                );
 
 
         // Patch interactions

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
 (

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,
+    typename ParticleType::trackingData& td
-    const scalar trackTime
 )
 {
+    // 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)

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,
+                typename ParticleType::trackingData& td
-                const scalar trackTime
             );
 
 

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)
         :

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),

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;

src/lagrangian/basic/particle/particleI.H

@@ -77,11 +77,7 @@ inline void Foam::particle::movingTetGeometry
     const vector ccOld = mesh.oldCellCentres()[celli_];
     const vector ccNew = mesh.cellCentres()[celli_];
 
-    // Old and new points and cell centres are not sub-cycled. If we are sub-
+    const scalar f0 = stepFraction_, f1 = fraction;
-    // cycling, then we have to account for the timestep change here by
-    // modifying the fractions that we take of the old and new geometry.
-    const Pair<scalar> s = stepFractionSpan(mesh);
-    const scalar f0 = s[0] + stepFraction_*s[1], f1 = fraction*s[1];
 
     centre[0] = ccOld + f0*(ccNew - ccOld);
     base[0] = ptsOld[triIs[0]] + f0*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
@@ -159,6 +155,12 @@ inline Foam::label Foam::particle::face() const
 }
 
 
+inline Foam::label& Foam::particle::face()
+{
+    return facei_;
+}
+
+
 inline Foam::scalar Foam::particle::stepFraction() const
 {
     return stepFraction_;
@@ -195,44 +197,6 @@ inline Foam::label& Foam::particle::origId()
 }
 
 
-inline Foam::Pair<Foam::scalar> Foam::particle::stepFractionSpan
-(
-    const polyMesh& mesh
-) const
-{
-    if (mesh.time().subCycling())
-    {
-        const TimeState& tsNew = mesh.time();
-        const TimeState& tsOld = mesh.time().prevTimeState();
-
-        const scalar tFrac =
-        (
-            (tsNew.value() - tsNew.deltaTValue())
-          - (tsOld.value() - tsOld.deltaTValue())
-        )/tsOld.deltaTValue();
-
-        const scalar dtFrac = tsNew.deltaTValue()/tsOld.deltaTValue();
-
-        return Pair<scalar>(tFrac, dtFrac);
-    }
-    else
-    {
-        return Pair<scalar>(0, 1);
-    }
-}
-
-
-inline Foam::scalar Foam::particle::currentTimeFraction
-(
-    const polyMesh& mesh
-) const
-{
-    const Pair<scalar> s = stepFractionSpan(mesh);
-
-    return s[0] + stepFraction_*s[1];
-}
-
-
 inline Foam::tetIndices Foam::particle::currentTetIndices
 (
     const polyMesh& mesh

src/lagrangian/basic/particle/particleTemplates.C

@@ -347,8 +347,7 @@ bool Foam::particle::hitNonConformalCyclicPatch
     const remote receiveProcFace =
         nccpp.ray
         (
-            stepFractionSpan(td.mesh)[0]
+            stepFraction_,
-          + stepFraction_*stepFractionSpan(td.mesh)[1],
             nccpp.origPatch().whichFace(facei_),
             sendPos,
             displacement - fraction*sendDisplacement,

src/lagrangian/basic/passiveParticle/passiveParticle.H

@@ -67,7 +67,7 @@ public:
             const label tetPti
         )
         :
-            particle(mesh, coordinates, celli, tetFacei, tetPti)
+            particle(mesh, coordinates, celli, tetFacei, tetPti, -1)
         {}
 
 

src/lagrangian/molecularDynamics/molecule/molecule/molecule.C

@@ -68,8 +68,7 @@ Foam::tensor Foam::molecule::rotationTensorZ(scalar phi) const
 bool Foam::molecule::move
 (
     moleculeCloud& cloud,
-    trackingData& td,
+    trackingData& td
-    const scalar trackTime
 )
 {
     td.keepParticle = true;
@@ -77,6 +76,8 @@ bool Foam::molecule::move
 
     const constantProperties& constProps(cloud.constProps(id_));
 
+    const scalar trackTime = td.mesh.time().deltaTValue();
+
     if (td.part() == trackingData::tpVelocityHalfStep0)
     {
         // First leapfrog velocity adjust part, required before tracking+force

src/lagrangian/molecularDynamics/molecule/molecule/molecule.H

@@ -259,25 +259,6 @@ public:
 
     // Constructors
 
-        //- Construct from components
-        inline molecule
-        (
-            const polyMesh& mesh,
-            const barycentric& coordinates,
-            const label celli,
-            const label tetFacei,
-            const label tetPti,
-            const tensor& Q,
-            const vector& v,
-            const vector& a,
-            const vector& pi,
-            const vector& tau,
-            const vector& specialPosition,
-            const constantProperties& constProps,
-            const label special,
-            const label id
-        );
-
         //- Construct from a position and a cell, searching for the rest of the
         //  required topology
         inline molecule
@@ -316,7 +297,7 @@ public:
 
         // Tracking
 
-            bool move(moleculeCloud&, trackingData&, const scalar trackTime);
+            bool move(moleculeCloud&, trackingData&);
 
             virtual void transformProperties(const transformer&);
 

src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H

@@ -219,43 +219,6 @@ inline Foam::molecule::constantProperties::constantProperties
 }
 
 
-inline Foam::molecule::molecule
-(
-    const polyMesh& mesh,
-    const barycentric& coordinates,
-    const label celli,
-    const label tetFacei,
-    const label tetPti,
-    const tensor& Q,
-    const vector& v,
-    const vector& a,
-    const vector& pi,
-    const vector& tau,
-    const vector& specialPosition,
-    const constantProperties& constProps,
-    const label special,
-    const label id
-
-)
-:
-    particle(mesh, coordinates, celli, tetFacei, tetPti),
-    Q_(Q),
-    v_(v),
-    a_(a),
-    pi_(pi),
-    tau_(tau),
-    specialPosition_(specialPosition),
-    potentialEnergy_(0.0),
-    rf_(Zero),
-    special_(special),
-    id_(id),
-    siteForces_(constProps.nSites(), Zero),
-    sitePositions_(constProps.nSites())
-{
-    setSitePositions(mesh, constProps);
-}
-
-
 inline Foam::molecule::molecule
 (
     const polyMesh& mesh,

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C

@@ -1116,18 +1116,18 @@ void Foam::moleculeCloud::evolve()
     molecule::trackingData td(*this);
 
     td.part() = molecule::trackingData::tpVelocityHalfStep0;
-    Cloud<molecule>::move(*this, td, mesh_.time().deltaTValue());
+    Cloud<molecule>::move(*this, td);
 
     td.part() = molecule::trackingData::tpLinearTrack;
-    Cloud<molecule>::move(*this, td, mesh_.time().deltaTValue());
+    Cloud<molecule>::move(*this, td);
 
     td.part() = molecule::trackingData::tpRotationalTrack;
-    Cloud<molecule>::move(*this, td, mesh_.time().deltaTValue());
+    Cloud<molecule>::move(*this, td);
 
     calculateForce();
 
     td.part() = molecule::trackingData::tpVelocityHalfStep1;
-    Cloud<molecule>::move(*this, td, mesh_.time().deltaTValue());
+    Cloud<molecule>::move(*this, td);
 }
 
 

src/lagrangian/parcel/clouds/Templates/CollidingCloud/CollidingCloud.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,34 +45,6 @@ void Foam::CollidingCloud<CloudType>::setModels()
 }
 
 
-template<class CloudType>
-template<class TrackCloudType>
-void Foam::CollidingCloud<CloudType>::moveCollide
-(
-    TrackCloudType& cloud,
-    typename parcelType::trackingData& td,
-    const scalar deltaT
-)
-{
-    td.part() = parcelType::trackingData::tpVelocityHalfStep;
-    CloudType::move(cloud, td, deltaT);
-
-    td.part() = parcelType::trackingData::tpLinearTrack;
-    CloudType::move(cloud, td, deltaT);
-
-    // td.part() = parcelType::trackingData::tpRotationalTrack;
-    // CloudType::move(cloud, td, deltaT);
-
-    this->updateCellOccupancy();
-
-    this->collision().collide();
-
-    td.part() = parcelType::trackingData::tpVelocityHalfStep;
-    CloudType::move(cloud, td, deltaT);
-}
-
-
-
 template<class CloudType>
 void Foam::CollidingCloud<CloudType>::cloudReset(CollidingCloud<CloudType>& c)
 {
@@ -216,29 +188,40 @@ void  Foam::CollidingCloud<CloudType>::motion
     // + calculate forces in new position
     // + apply half deltaV with new force
 
-    label nSubCycles = collision().nSubCycles();
+    const label nSubCycles = collision().nSubCycles();
 
     if (nSubCycles > 1)
     {
         Info<< "    " << nSubCycles << " move-collide subCycles" << endl;
-
-        subCycleTime moveCollideSubCycle
-        (
-            const_cast<Time&>(this->db().time()),
-            nSubCycles
-        );
-
-        while(!(++moveCollideSubCycle).end())
-        {
-            moveCollide(cloud, td, this->db().time().deltaTValue());
-        }
-
-        moveCollideSubCycle.endSubCycle();
     }
-    else
+
+    for (label subCyclei = 0; subCyclei < nSubCycles; ++ subCyclei)
     {
-        moveCollide(cloud, td, this->db().time().deltaTValue());
+        td.stepFractionRange() =
+            Pair<scalar>
+            (
+                scalar(subCyclei)/nSubCycles,
+                scalar(subCyclei + 1)/nSubCycles
+            );
+
+        td.part() = parcelType::trackingData::tpVelocityHalfStep;
+        CloudType::move(cloud, td);
+
+        td.part() = parcelType::trackingData::tpLinearTrack;
+        CloudType::move(cloud, td);
+
+        // td.part() = parcelType::trackingData::tpRotationalTrack;
+        // CloudType::move(cloud, td);
+
+        this->updateCellOccupancy();
+
+        this->collision().collide();
+
+        td.part() = parcelType::trackingData::tpVelocityHalfStep;
+        CloudType::move(cloud, td);
     }
+
+    td.stepFractionRange() = Pair<scalar>(0, 1);
 }
 
 

src/lagrangian/parcel/clouds/Templates/CollidingCloud/CollidingCloud.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -116,15 +116,6 @@ protected:
 
         // Cloud evolution functions
 
-            //- Move-collide particles
-            template<class TrackCloudType>
-            void moveCollide
-            (
-                TrackCloudType& cloud,
-                typename parcelType::trackingData& td,
-                const scalar deltaT
-            );
-
             //- Reset state of cloud
             void cloudReset(CollidingCloud<CloudType>& c);
 

src/lagrangian/parcel/clouds/Templates/MPPICCloud/MPPICCloud.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -196,7 +196,7 @@ void Foam::MPPICCloud<CloudType>::motion
 
     // force calculation and tracking
     td.part() = parcelType::trackingData::tpPredictTrack;
-    CloudType::move(cloud, td, this->db().time().deltaTValue());
+    CloudType::move(cloud, td);
 
 
     // Preliminary
@@ -227,11 +227,11 @@ void Foam::MPPICCloud<CloudType>::motion
 
         // calculate the damping velocity corrections without moving the parcels
         td.part() = parcelType::trackingData::tpDampingNoTrack;
-        CloudType::move(cloud, td, this->db().time().deltaTValue());
+        CloudType::move(cloud, td);
 
         // correct the parcel positions and velocities
         td.part() = parcelType::trackingData::tpCorrectTrack;
-        CloudType::move(cloud, td, this->db().time().deltaTValue());
+        CloudType::move(cloud, td);
 
         // finalise and free memory
         dampingModel_->cacheFields(false);
@@ -254,9 +254,9 @@ void Foam::MPPICCloud<CloudType>::motion
         td.updateAverages(cloud);
         packingModel_->cacheFields(true);
         td.part() = parcelType::trackingData::tpPackingNoTrack;
-        CloudType::move(cloud, td, this->db().time().deltaTValue());
+        CloudType::move(cloud, td);
         td.part() = parcelType::trackingData::tpCorrectTrack;
-        CloudType::move(cloud, td, this->db().time().deltaTValue());
+        CloudType::move(cloud, td);
         packingModel_->cacheFields(false);
     }
 

src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloud.C

@@ -94,29 +94,57 @@ void Foam::MomentumCloud<CloudType>::solve
     typename parcelType::trackingData& td
 )
 {
+    this->changeTimeStep();
+
     if (solution_.steadyState())
     {
         cloud.storeState();
+    }
 
-        cloud.preEvolve();
+    cloud.preEvolve();
 
-        evolveCloud(cloud, td);
+    if (solution_.coupled())
+    {
+        cloud.resetSourceTerms();
+    }
 
-        if (solution_.coupled())
+    if (solution_.transient())
+    {
+        label preInjectionSize = this->size();
+
+        this->surfaceFilm().inject(cloud);
+
+        // Update the cellOccupancy if the size of the cloud has changed
+        // during the injection.
+        if (preInjectionSize != this->size())
         {
-            cloud.relaxSources(cloud.cloudCopy());
+            updateCellOccupancy();
+            preInjectionSize = this->size();
         }
+
+        injectors_.inject(cloud, td);
+
+        // Assume that motion will update the cellOccupancy as necessary
+        // before it is required.
+        cloud.motion(cloud, td);
+
+        stochasticCollision().update(td);
     }
     else
     {
-        cloud.preEvolve();
+        injectors_.injectSteadyState(cloud, td);
 
-        evolveCloud(cloud, td);
+        CloudType::move(cloud, td);
+    }
 
-        if (solution_.coupled())
+    if (solution_.coupled() && solution_.transient())
-        {
+    {
-            cloud.scaleSources();
+        cloud.scaleSources();
-        }
+    }
+
+    if (solution_.coupled() && solution_.steadyState())
+    {
+        cloud.relaxSources(cloud.cloudCopy());
     }
 
     cloud.info();
@@ -175,53 +203,6 @@ void Foam::MomentumCloud<CloudType>::updateCellOccupancy()
 }
 
 
-template<class CloudType>
-template<class TrackCloudType>
-void Foam::MomentumCloud<CloudType>::evolveCloud
-(
-    TrackCloudType& cloud,
-    typename parcelType::trackingData& td
-)
-{
-    if (solution_.coupled())
-    {
-        cloud.resetSourceTerms();
-    }
-
-    if (solution_.transient())
-    {
-        label preInjectionSize = this->size();
-
-        this->surfaceFilm().inject(cloud);
-
-        // Update the cellOccupancy if the size of the cloud has changed
-        // during the injection.
-        if (preInjectionSize != this->size())
-        {
-            updateCellOccupancy();
-            preInjectionSize = this->size();
-        }
-
-        injectors_.inject(cloud, td);
-
-
-        // Assume that motion will update the cellOccupancy as necessary
-        // before it is required.
-        cloud.motion(cloud, td);
-
-        stochasticCollision().update(td, solution_.trackTime());
-    }
-    else
-    {
-//        this->surfaceFilm().injectSteadyState(cloud);
-
-        injectors_.injectSteadyState(cloud, td, solution_.trackTime());
-
-        CloudType::move(cloud, td, solution_.trackTime());
-    }
-}
-
-
 template<class CloudType>
 void Foam::MomentumCloud<CloudType>::postEvolve()
 {
@@ -553,8 +534,7 @@ Foam::MomentumCloud<CloudType>::~MomentumCloud()
 template<class CloudType>
 void Foam::MomentumCloud<CloudType>::setParcelThermoProperties
 (
-    parcelType& parcel,
+    parcelType& parcel
-    const scalar lagrangianDt
 )
 {
     parcel.rho() = constProps_.rho0();
@@ -565,13 +545,9 @@ template<class CloudType>
 void Foam::MomentumCloud<CloudType>::checkParcelProperties
 (
     parcelType& parcel,
-    const scalar lagrangianDt,
     const bool fullyDescribed
 )
 {
-    const scalar carrierDt = this->mesh().time().deltaTValue();
-    parcel.stepFraction() = (carrierDt - lagrangianDt)/carrierDt;
-
     if (parcel.typeId() == -1)
     {
         parcel.typeId() = constProps_.parcelTypeId();
@@ -694,7 +670,7 @@ void Foam::MomentumCloud<CloudType>::motion
     typename parcelType::trackingData& td
 )
 {
-    CloudType::move(cloud, td, solution_.trackTime());
+    CloudType::move(cloud, td);
 
     updateCellOccupancy();
 }
@@ -731,9 +707,7 @@ void Foam::MomentumCloud<CloudType>::patchData
             const vector& Uw0 =
                 U_.oldTime().boundaryField()[patchi][patchFacei];
 
-            const scalar f = p.currentTimeFraction(this->mesh());
+            const vector Uw = Uw0 + p.stepFraction()*(Uw1 - Uw0);
-
-            const vector Uw = Uw0 + f*(Uw1 - Uw0);
 
             const tensor nnw = nw*nw;
 

src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloud.H

@@ -253,14 +253,6 @@ protected:
             //  already been used
             void updateCellOccupancy();
 
-            //- Evolve the cloud
-            template<class TrackCloudType>
-            void evolveCloud
-            (
-                TrackCloudType& cloud,
-                typename parcelType::trackingData& td
-            );
-
             //- Post-evolve
             void postEvolve();
 
@@ -538,17 +530,12 @@ public:
         // Cloud evolution functions
 
             //- Set parcel thermo properties
-            void setParcelThermoProperties
+            void setParcelThermoProperties(parcelType& parcel);
-            (
-                parcelType& parcel,
-                const scalar lagrangianDt
-            );
 
             //- Check parcel properties
             void checkParcelProperties
             (
                 parcelType& parcel,
-                const scalar lagrangianDt,
                 const bool fullyDescribed
             );
 

src/lagrangian/parcel/clouds/Templates/ReactingCloud/ReactingCloud.C

@@ -195,11 +195,10 @@ Foam::ReactingCloud<CloudType>::~ReactingCloud()
 template<class CloudType>
 void Foam::ReactingCloud<CloudType>::setParcelThermoProperties
 (
-    parcelType& parcel,
+    parcelType& parcel
-    const scalar lagrangianDt
 )
 {
-    CloudType::setParcelThermoProperties(parcel, lagrangianDt);
+    CloudType::setParcelThermoProperties(parcel);
 
     parcel.Y() = this->composition().YMixture0();
 }
@@ -209,11 +208,10 @@ template<class CloudType>
 void Foam::ReactingCloud<CloudType>::checkParcelProperties
 (
     parcelType& parcel,
-    const scalar lagrangianDt,
     const bool fullyDescribed
 )
 {
-    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
+    CloudType::checkParcelProperties(parcel, fullyDescribed);
 
     if (fullyDescribed)
     {

src/lagrangian/parcel/clouds/Templates/ReactingCloud/ReactingCloud.H

@@ -259,17 +259,12 @@ public:
         // Cloud evolution functions
 
             //- Set parcel thermo properties
-            void setParcelThermoProperties
+            void setParcelThermoProperties(parcelType& parcel);
-            (
-                parcelType& parcel,
-                const scalar lagrangianDt
-            );
 
             //- Check parcel properties
             void checkParcelProperties
             (
                 parcelType& parcel,
-                const scalar lagrangianDt,
                 const bool fullyDescribed
             );
 

src/lagrangian/parcel/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C

@@ -152,11 +152,10 @@ Foam::ReactingMultiphaseCloud<CloudType>::~ReactingMultiphaseCloud()
 template<class CloudType>
 void Foam::ReactingMultiphaseCloud<CloudType>::setParcelThermoProperties
 (
-    parcelType& parcel,
+    parcelType& parcel
-    const scalar lagrangianDt
 )
 {
-    CloudType::setParcelThermoProperties(parcel, lagrangianDt);
+    CloudType::setParcelThermoProperties(parcel);
 
     label idGas = this->composition().idGas();
     label idLiquid = this->composition().idLiquid();
@@ -172,11 +171,10 @@ template<class CloudType>
 void Foam::ReactingMultiphaseCloud<CloudType>::checkParcelProperties
 (
     parcelType& parcel,
-    const scalar lagrangianDt,
     const bool fullyDescribed
 )
 {
-    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
+    CloudType::checkParcelProperties(parcel, fullyDescribed);
 
     if (fullyDescribed)
     {

src/lagrangian/parcel/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.H

@@ -255,17 +255,12 @@ public:
         // Cloud evolution functions
 
             //- Set parcel thermo properties
-            void setParcelThermoProperties
+            void setParcelThermoProperties(parcelType& parcel);
-            (
-                parcelType& parcel,
-                const scalar lagrangianDt
-            );
 
             //- Check parcel properties
             void checkParcelProperties
             (
                 parcelType& parcel,
-                const scalar lagrangianDt,
                 const bool fullyDescribed
             );
 

src/lagrangian/parcel/clouds/Templates/SprayCloud/SprayCloud.C

@@ -140,11 +140,10 @@ Foam::SprayCloud<CloudType>::~SprayCloud()
 template<class CloudType>
 void Foam::SprayCloud<CloudType>::setParcelThermoProperties
 (
-    parcelType& parcel,
+    parcelType& parcel
-    const scalar lagrangianDt
 )
 {
-    CloudType::setParcelThermoProperties(parcel, lagrangianDt);
+    CloudType::setParcelThermoProperties(parcel);
 
     const liquidMixtureProperties& liqMix = this->composition().liquids();
 
@@ -164,11 +163,10 @@ template<class CloudType>
 void Foam::SprayCloud<CloudType>::checkParcelProperties
 (
     parcelType& parcel,
-    const scalar lagrangianDt,
     const bool fullyDescribed
 )
 {
-    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
+    CloudType::checkParcelProperties(parcel, fullyDescribed);
 
     // store the injection position and initial drop size
     parcel.position0() = parcel.position(this->mesh());

src/lagrangian/parcel/clouds/Templates/SprayCloud/SprayCloud.H

@@ -209,17 +209,12 @@ public:
         // Cloud evolution functions
 
             //- Set parcel thermo properties
-            void setParcelThermoProperties
+            void setParcelThermoProperties(parcelType& parcel);
-            (
-                parcelType& parcel,
-                const scalar lagrangianDt
-            );
 
             //- Check parcel properties
             void checkParcelProperties
             (
                 parcelType& parcel,
-                const scalar lagrangianDt,
                 const bool fullyDescribed
             );
 

src/lagrangian/parcel/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -359,11 +359,10 @@ Foam::ThermoCloud<CloudType>::~ThermoCloud()
 template<class CloudType>
 void Foam::ThermoCloud<CloudType>::setParcelThermoProperties
 (
-    parcelType& parcel,
+    parcelType& parcel
-    const scalar lagrangianDt
 )
 {
-    CloudType::setParcelThermoProperties(parcel, lagrangianDt);
+    CloudType::setParcelThermoProperties(parcel);
 
     parcel.T() = constProps_.T0();
     parcel.Cp() = constProps_.Cp0();
@@ -374,11 +373,10 @@ template<class CloudType>
 void Foam::ThermoCloud<CloudType>::checkParcelProperties
 (
     parcelType& parcel,
-    const scalar lagrangianDt,
     const bool fullyDescribed
 )
 {
-    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
+    CloudType::checkParcelProperties(parcel, fullyDescribed);
 }
 
 

src/lagrangian/parcel/clouds/Templates/ThermoCloud/ThermoCloud.H

@@ -348,17 +348,12 @@ public:
         // Cloud evolution functions
 
             //- Set parcel thermo properties
-            void setParcelThermoProperties
+            void setParcelThermoProperties(parcelType& parcel);
-            (
-                parcelType& parcel,
-                const scalar lagrangianDt
-            );
 
             //- Check parcel properties
             void checkParcelProperties
             (
                 parcelType& parcel,
-                const scalar lagrangianDt,
                 const bool fullyDescribed
             );
 

src/lagrangian/parcel/parcels/Templates/CollidingParcel/CollidingParcel.C

@@ -48,8 +48,7 @@ template<class TrackCloudType>
 bool Foam::CollidingParcel<ParcelType>::move
 (
     TrackCloudType& cloud,
-    trackingData& td,
+    trackingData& td
-    const scalar trackTime
 )
 {
     typename TrackCloudType::parcelType& p =
@@ -62,9 +61,11 @@ bool Foam::CollidingParcel<ParcelType>::move
             // First and last leapfrog velocity adjust part, required
             // before and after tracking and force calculation
 
-            p.U() += 0.5*trackTime*p.f()/p.mass();
+            const Pair<scalar>& sfr = td.stepFractionRange();
+            const scalar dt = (sfr.second() - sfr.first())*td.trackTime()/2;
 
-            p.angularMomentum() += 0.5*trackTime*p.torque();
+            p.U() += dt*p.f()/p.mass();
+            p.angularMomentum() += dt*p.torque();
 
             td.keepParticle = true;
             td.sendToProc = -1;
@@ -74,7 +75,7 @@ bool Foam::CollidingParcel<ParcelType>::move
 
         case trackingData::tpLinearTrack:
         {
-            ParcelType::move(cloud, td, trackTime);
+            ParcelType::move(cloud, td);
 
             break;
         }
@@ -85,13 +86,6 @@ bool Foam::CollidingParcel<ParcelType>::move
 
             break;
         }
-
-        default:
-        {
-            FatalErrorInFunction
-                << td.part() << " is an invalid part of the tracking method."
-                << abort(FatalError);
-        }
     }
 
     return td.keepParticle;

src/lagrangian/parcel/parcels/Templates/CollidingParcel/CollidingParcel.H

@@ -219,7 +219,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
@@ -286,12 +287,7 @@ public:
 
             //- Move the parcel
             template<class TrackCloudType>
-            bool move
+            bool move(TrackCloudType& cloud, trackingData& td);
-            (
-                TrackCloudType& cloud,
-                trackingData& td,
-                const scalar trackTime
-            );
 
             //- Transform the physical properties of the particle
             //  according to the given transformation tensor

src/lagrangian/parcel/parcels/Templates/CollidingParcel/CollidingParcelI.H

@@ -66,10 +66,11 @@ inline Foam::CollidingParcel<ParcelType>::CollidingParcel
     const barycentric& coordinates,
     const label celli,
     const label tetFacei,
-    const label tetPti
+    const label tetPti,
+    const label facei
 )
 :
-    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti, facei),
     f_(Zero),
     angularMomentum_(Zero),
     torque_(Zero),

src/lagrangian/parcel/parcels/Templates/MPPICParcel/MPPICParcel.C

@@ -45,8 +45,7 @@ template<class TrackCloudType>
 bool Foam::MPPICParcel<ParcelType>::move
 (
     TrackCloudType& cloud,
-    trackingData& td,
+    trackingData& td
-    const scalar trackTime
 )
 {
     typename TrackCloudType::parcelType& p =
@@ -56,14 +55,14 @@ bool Foam::MPPICParcel<ParcelType>::move
     {
         case trackingData::tpPredictTrack:
         {
-            ParcelType::move(cloud, td, trackTime);
+            ParcelType::move(cloud, td);
 
             break;
         }
         case trackingData::tpDampingNoTrack:
         {
             p.UCorrect() =
-                cloud.dampingModel().velocityCorrection(p, trackTime);
+                cloud.dampingModel().velocityCorrection(p, td.trackTime());
 
             td.keepParticle = true;
             td.sendToProc = -1;
@@ -73,7 +72,7 @@ bool Foam::MPPICParcel<ParcelType>::move
         case trackingData::tpPackingNoTrack:
         {
             p.UCorrect() =
-                cloud.packingModel().velocityCorrection(p, trackTime);
+                cloud.packingModel().velocityCorrection(p, td.trackTime());
 
             td.keepParticle = true;
             td.sendToProc = -1;
@@ -87,7 +86,7 @@ bool Foam::MPPICParcel<ParcelType>::move
 
             Swap(p.U(), p.UCorrect());
 
-            ParcelType::move(cloud, td, trackTime);
+            ParcelType::move(cloud, td);
 
             Swap(p.U(), p.UCorrect());
 

src/lagrangian/parcel/parcels/Templates/MPPICParcel/MPPICParcel.H

@@ -191,7 +191,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
@@ -237,12 +238,7 @@ public:
 
             //- Move the parcel
             template<class TrackCloudType>
-            bool move
+            bool move(TrackCloudType& cloud, trackingData& td);
-            (
-                TrackCloudType& cloud,
-                trackingData& td,
-                const scalar trackTime
-            );
 
 
     // Friend Functions

src/lagrangian/parcel/parcels/Templates/MPPICParcel/MPPICParcelI.H

@@ -34,10 +34,11 @@ inline Foam::MPPICParcel<ParcelType>::MPPICParcel
     const barycentric& coordinates,
     const label celli,
     const label tetFacei,
-    const label tetPti
+    const label tetPti,
+    const label facei
 )
 :
-    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti, facei),
     UCorrect_(Zero)
 {}
 

src/lagrangian/parcel/parcels/Templates/MomentumParcel/MomentumParcel.C

@@ -257,8 +257,7 @@ template<class TrackCloudType>
 bool Foam::MomentumParcel<ParcelType>::move
 (
     TrackCloudType& cloud,
-    trackingData& td,
+    trackingData& td
-    const scalar trackTime
 )
 {
     typename TrackCloudType::parcelType& p =
@@ -272,7 +271,12 @@ bool Foam::MomentumParcel<ParcelType>::move
     const scalarField& cellLengthScale = cloud.cellLengthScale();
     const scalar maxCo = cloud.solution().maxCo();
 
-    while (ttd.keepParticle && ttd.sendToProc == -1 && p.stepFraction() < 1)
+    while
+    (
+        ttd.keepParticle
+     && ttd.sendToProc == -1
+     && p.stepFraction() < ttd.stepFractionRange().second()
+    )
     {
         if (p.moving() && p.onFace())
         {
@@ -284,7 +288,7 @@ bool Foam::MomentumParcel<ParcelType>::move
         const scalar sfrac = p.stepFraction();
 
         // Total displacement over the time-step
-        const vector s = trackTime*U_;
+        const vector s = ttd.trackTime()*U_;
 
         // Cell length scale
         const scalar l = cellLengthScale[p.cell()];
@@ -295,7 +299,7 @@ bool Foam::MomentumParcel<ParcelType>::move
         // Fraction of the displacement to track in this loop. This is limited
         // to ensure that the both the time and distance tracked is less than
         // maxCo times the total value.
-        scalar f = 1 - p.stepFraction();
+        scalar f = ttd.stepFractionRange().second() - p.stepFraction();
         f = min(f, maxCo);
         f = min(f, maxCo/min(max(mag(s)/l, rootSmall), rootGreat));
         if (p.moving())
@@ -313,7 +317,7 @@ bool Foam::MomentumParcel<ParcelType>::move
             p.stepFraction() += f;
         }
 
-        const scalar dt = (p.stepFraction() - sfrac)*trackTime;
+        const scalar dt = (p.stepFraction() - sfrac)*ttd.trackTime();
 
         // Avoid problems with extremely small timesteps
         if (dt > rootVSmall)

src/lagrangian/parcel/parcels/Templates/MomentumParcel/MomentumParcel.H

@@ -191,6 +191,16 @@ public:
             const vector& g_;
 
 
+            // Tracking Parameters
+
+                //- Track time. Total across the entire step.
+                scalar trackTime_;
+
+                //- Step fraction range to track between. Allows for the
+                //  creation of sub-steps.
+                Pair<scalar> stepFractionRange_;
+
+
     public:
 
         // Constructors
@@ -232,8 +242,20 @@ public:
             //- Access the continuous phase viscosity
             inline scalar& muc();
 
-            // Return const access to the gravitational acceleration vector
+            // Return the gravitational acceleration vector
             inline const vector& g() const;
+
+            //- Return the tracking time
+            inline scalar trackTime() const;
+
+            //- Access the tracking time
+            inline scalar& trackTime();
+
+            //- Return the step fraction range to track between
+            inline Pair<scalar> stepFractionRange() const;
+
+            //- Access the step fraction range to track between
+            inline Pair<scalar>& stepFractionRange();
     };
 
 
@@ -355,7 +377,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
@@ -536,12 +559,7 @@ public:
 
             //- Move the parcel
             template<class TrackCloudType>
-            bool move
+            bool move(TrackCloudType& cloud, trackingData& td);
-            (
-                TrackCloudType& cloud,
-                trackingData& td,
-                const scalar trackTime
-            );
 
 
         // Transformations

src/lagrangian/parcel/parcels/Templates/MomentumParcel/MomentumParcelI.H

@@ -77,10 +77,11 @@ inline Foam::MomentumParcel<ParcelType>::MomentumParcel
     const barycentric& coordinates,
     const label celli,
     const label tetFacei,
-    const label tetPti
+    const label tetPti,
+    const label facei
 )
 :
-    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti, facei),
     moving_(true),
     typeId_(-1),
     nParticle_(0),

src/lagrangian/parcel/parcels/Templates/MomentumParcel/MomentumParcelTrackingDataI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,7 +58,9 @@ inline Foam::MomentumParcel<ParcelType>::trackingData::trackingData
     rhoc_(Zero),
     Uc_(Zero),
     muc_(Zero),
-    g_(cloud.g().value())
+    g_(cloud.g().value()),
+    trackTime_(cloud.solution().trackTime()),
+    stepFractionRange_(0, 1)
 {}
 
 
@@ -86,14 +88,6 @@ Foam::MomentumParcel<ParcelType>::trackingData::muInterp() const
 }
 
 
-template<class ParcelType>
-inline const Foam::vector&
-Foam::MomentumParcel<ParcelType>::trackingData::g() const
-{
-    return g_;
-}
-
-
 template<class ParcelType>
 inline Foam::scalar
 Foam::MomentumParcel<ParcelType>::trackingData::rhoc() const
@@ -138,4 +132,44 @@ inline Foam::scalar& Foam::MomentumParcel<ParcelType>::trackingData::muc()
 }
 
 
+template<class ParcelType>
+inline const Foam::vector&
+Foam::MomentumParcel<ParcelType>::trackingData::g() const
+{
+    return g_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar
+Foam::MomentumParcel<ParcelType>::trackingData::trackTime() const
+{
+    return trackTime_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar&
+Foam::MomentumParcel<ParcelType>::trackingData::trackTime()
+{
+    return trackTime_;
+}
+
+
+template<class ParcelType>
+inline Foam::Pair<Foam::scalar>
+Foam::MomentumParcel<ParcelType>::trackingData::stepFractionRange() const
+{
+    return stepFractionRange_;
+}
+
+
+template<class ParcelType>
+inline Foam::Pair<Foam::scalar>&
+Foam::MomentumParcel<ParcelType>::trackingData::stepFractionRange()
+{
+    return stepFractionRange_;
+}
+
+
 // ************************************************************************* //

src/lagrangian/parcel/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H

@@ -290,7 +290,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

src/lagrangian/parcel/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H

@@ -73,10 +73,11 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
     const barycentric& coordinates,
     const label celli,
     const label tetFacei,
-    const label tetPti
+    const label tetPti,
+    const label facei
 )
 :
-    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti, facei),
     YGas_(0),
     YLiquid_(0),
     YSolid_(0),

src/lagrangian/parcel/parcels/Templates/ReactingParcel/ReactingParcel.H

@@ -197,7 +197,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

src/lagrangian/parcel/parcels/Templates/ReactingParcel/ReactingParcelI.H

@@ -68,10 +68,11 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
     const barycentric& coordinates,
     const label celli,
     const label tetFacei,
-    const label tetPti
+    const label tetPti,
+    const label facei
 )
 :
-    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti, facei),
     mass0_(0.0),
     Y_(0)
 {}

src/lagrangian/parcel/parcels/Templates/SprayParcel/SprayParcel.H

@@ -199,7 +199,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

src/lagrangian/parcel/parcels/Templates/SprayParcel/SprayParcelI.H

@@ -112,10 +112,11 @@ inline Foam::SprayParcel<ParcelType>::SprayParcel
     const barycentric& coordinates,
     const label celli,
     const label tetFacei,
-    const label tetPti
+    const label tetPti,
+    const label facei
 )
 :
-    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti, facei),
     d0_(this->d()),
     position0_(this->position(mesh)),
     sigma_(0.0),

src/lagrangian/parcel/parcels/Templates/ThermoParcel/ThermoParcel.H

@@ -314,7 +314,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

src/lagrangian/parcel/parcels/Templates/ThermoParcel/ThermoParcelI.H

@@ -79,10 +79,11 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
     const barycentric& coordinates,
     const label celli,
     const label tetFacei,
-    const label tetPti
+    const label tetPti,
+    const label facei
 )
 :
-    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti, facei),
     T_(0.0),
     Cp_(0.0)
 {}

src/lagrangian/parcel/submodels/Momentum/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C

@@ -144,22 +144,19 @@ bool Foam::PairSpringSliderDashpot<CloudType>::controlsTimestep() const
 template<class CloudType>
 Foam::label Foam::PairSpringSliderDashpot<CloudType>::nSubCycles() const
 {
-    if (!(this->owner().size()))
+    if (!this->owner().size())
     {
         return 1;
     }
 
-    scalar RMin;
+    scalar RMin, rhoMax, UMagMax;
-    scalar rhoMax;
-    scalar UMagMax;
-
     findMinMaxProperties(RMin, rhoMax, UMagMax);
 
     // Note:  pi^(7/5)*(5/4)^(2/5) = 5.429675
-    scalar minCollisionDeltaT =
+    const scalar minCollisionDeltaT =
         5.429675
        *RMin
-       *pow(rhoMax/(Estar_*sqrt(UMagMax) + vSmall), 0.4)
+       *pow(rhoMax/(Estar_*sqrt(UMagMax) + rootVSmall), 0.4)
        /collisionResolutionSteps_;
 
     return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);

src/lagrangian/parcel/submodels/Momentum/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C

@@ -311,22 +311,19 @@ bool Foam::WallLocalSpringSliderDashpot<CloudType>::controlsTimestep() const
 template<class CloudType>
 Foam::label Foam::WallLocalSpringSliderDashpot<CloudType>::nSubCycles() const
 {
-    if (!(this->owner().size()))
+    if (!this->owner().size())
     {
         return 1;
     }
 
-    scalar rMin;
+    scalar rMin, rhoMax, UMagMax;
-    scalar rhoMax;
-    scalar UMagMax;
-
     findMinMaxProperties(rMin, rhoMax, UMagMax);
 
     // Note:  pi^(7/5)*(5/4)^(2/5) = 5.429675
-    scalar minCollisionDeltaT =
+    const scalar minCollisionDeltaT =
         5.429675
        *rMin
-       *pow(rhoMax/(Estar_[maxEstarIndex_]*sqrt(UMagMax) + vSmall), 0.4)
+       *pow(rhoMax/(Estar_[maxEstarIndex_]*sqrt(UMagMax) + rootVSmall), 0.4)
        /collisionResolutionSteps_;
 
     return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);

src/lagrangian/parcel/submodels/Momentum/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C

@@ -241,22 +241,19 @@ bool Foam::WallSpringSliderDashpot<CloudType>::controlsTimestep() const
 template<class CloudType>
 Foam::label Foam::WallSpringSliderDashpot<CloudType>::nSubCycles() const
 {
-    if (!(this->owner().size()))
+    if (!this->owner().size())
     {
         return 1;
     }
 
-    scalar rMin;
+    scalar rMin, rhoMax, UMagMax;
-    scalar rhoMax;
-    scalar UMagMax;
-
     findMinMaxProperties(rMin, rhoMax, UMagMax);
 
     // Note:  pi^(7/5)*(5/4)^(2/5) = 5.429675
-    scalar minCollisionDeltaT =
+    const scalar minCollisionDeltaT =
         5.429675
        *rMin
-       *pow(rhoMax/(Estar_*sqrt(UMagMax) + vSmall), 0.4)
+       *pow(rhoMax/(Estar_*sqrt(UMagMax) + rootVSmall), 0.4)
        /collisionResolutionSteps_;
 
     return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);

src/lagrangian/parcel/submodels/Momentum/InjectionModel/CellZoneInjection/CellZoneInjection.C

@@ -42,7 +42,7 @@ void Foam::CellZoneInjection<CloudType>::setPositions
     const label nCells = cellZoneCells.size();
     Random& rnd = this->owner().rndGen();
 
-    DynamicList<vector> positions(nCells);          // initial size only
+    DynamicList<barycentric> injectorCoordinates(nCells); // initial size only
     DynamicList<label> injectorCells(nCells);       // initial size only
     DynamicList<label> injectorTetFaces(nCells);    // initial size only
     DynamicList<label> injectorTetPts(nCells);      // initial size only
@@ -95,8 +95,8 @@ void Foam::CellZoneInjection<CloudType>::setPositions
                     break;
                 }
             }
-            positions.append(cellTetIs[tetI].tet(mesh).randomPoint(rnd));
 
+            injectorCoordinates.append(barycentric01(rnd));
             injectorCells.append(celli);
             injectorTetFaces.append(cellTetIs[tetI].face());
             injectorTetPts.append(cellTetIs[tetI].tetPt());
@@ -104,22 +104,26 @@ void Foam::CellZoneInjection<CloudType>::setPositions
     }
 
     // Parallel operation manipulations
-    globalIndex globalPositions(positions.size());
+    globalIndex globalPositions(injectorCoordinates.size());
-    List<vector> allPositions(globalPositions.size(), point::max);
+    List<barycentric> allCoordinates
+    (
+        globalPositions.size(),
+        barycentric::uniform(NaN)
+    );
     List<label> allInjectorCells(globalPositions.size(), -1);
     List<label> allInjectorTetFaces(globalPositions.size(), -1);
     List<label> allInjectorTetPts(globalPositions.size(), -1);
 
     // Gather all positions on to all processors
-    SubList<vector>
+    SubList<barycentric>
     (
-        allPositions,
+        allCoordinates,
         globalPositions.localSize(Pstream::myProcNo()),
         globalPositions.offset(Pstream::myProcNo())
-    ) = positions;
+    ) = injectorCoordinates;
 
-    Pstream::listCombineGather(allPositions, minEqOp<point>());
+    Pstream::listCombineGather(allCoordinates, minEqOp<barycentric>());
-    Pstream::listCombineScatter(allPositions);
+    Pstream::listCombineScatter(allCoordinates);
 
     // Gather local cell tet and tet-point Ids, but leave non-local ids set -1
     SubList<label>
@@ -142,20 +146,10 @@ void Foam::CellZoneInjection<CloudType>::setPositions
     ) = injectorTetPts;
 
     // Transfer data
-    positions_.transfer(allPositions);
+    injectorCoordinates_.transfer(allCoordinates);
     injectorCells_.transfer(allInjectorCells);
     injectorTetFaces_.transfer(allInjectorTetFaces);
     injectorTetPts_.transfer(allInjectorTetPts);
-
-
-    if (debug)
-    {
-        OFstream points("points.obj");
-        forAll(positions_, i)
-        {
-            meshTools::writeOBJ(points, positions_[i]);
-        }
-    }
 }
 
 
@@ -172,7 +166,7 @@ Foam::CellZoneInjection<CloudType>::CellZoneInjection
     InjectionModel<CloudType>(dict, owner, modelName, typeName),
     cellZoneName_(this->coeffDict().lookup("cellZone")),
     numberDensity_(this->coeffDict().template lookup<scalar>("numberDensity")),
-    positions_(),
+    injectorCoordinates_(),
     injectorCells_(),
     injectorTetFaces_(),
     injectorTetPts_(),
@@ -199,7 +193,7 @@ Foam::CellZoneInjection<CloudType>::CellZoneInjection
     InjectionModel<CloudType>(im),
     cellZoneName_(im.cellZoneName_),
     numberDensity_(im.numberDensity_),
-    positions_(im.positions_),
+    injectorCoordinates_(im.injectorCoordinates_),
     injectorCells_(im.injectorCells_),
     injectorTetFaces_(im.injectorTetFaces_),
     injectorTetPts_(im.injectorTetPts_),
@@ -252,10 +246,11 @@ void Foam::CellZoneInjection<CloudType>::topoChange()
         setPositions(cellZoneCells);
 
         Info<< "    number density      = " << numberDensity_ << nl
-            << "    number of particles = " << positions_.size() << endl;
+            << "    number of particles = " << injectorCoordinates_.size()
+            << endl;
 
         // Construct parcel diameters
-        diameters_.setSize(positions_.size());
+        diameters_.setSize(injectorCoordinates_.size());
         forAll(diameters_, i)
         {
             diameters_[i] = sizeDistribution_->sample();
@@ -284,7 +279,7 @@ Foam::label Foam::CellZoneInjection<CloudType>::parcelsToInject
 {
     if ((0.0 >= time0) && (0.0 < time1))
     {
-        return positions_.size();
+        return injectorCoordinates_.size();
     }
     else
     {
@@ -318,14 +313,15 @@ void Foam::CellZoneInjection<CloudType>::setPositionAndCell
     const label parcelI,
     const label nParcels,
     const scalar time,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
-    position = positions_[parcelI];
+    coordinates = injectorCoordinates_[parcelI];
-    cellOwner = injectorCells_[parcelI];
+    celli = injectorCells_[parcelI];
     tetFacei = injectorTetFaces_[parcelI];
     tetPti = injectorTetPts_[parcelI];
 }

src/lagrangian/parcel/submodels/Momentum/InjectionModel/CellZoneInjection/CellZoneInjection.H

@@ -70,16 +70,16 @@ class CellZoneInjection
         //- Number density
         const scalar numberDensity_;
 
-        //- Field of parcel positions
+        //- Field of parcel coordinates
-        List<vector> positions_;
+        List<barycentric> injectorCoordinates_;
 
-        //- List of cell labels corresponding to injector positions
+        //- List of cell labels corresponding to injector coordinates
         labelList injectorCells_;
 
-        //- List of tetFace labels corresponding to injector positions
+        //- List of tetFace labels corresponding to injector coordinates
         labelList injectorTetFaces_;
 
-        //- List of tetPt labels corresponding to injector positions
+        //- List of tetPt labels corresponding to injector coordinates
         labelList injectorTetPts_;
 
         //- Field of parcel diameters
@@ -154,10 +154,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             //- Set the parcel properties

src/lagrangian/parcel/submodels/Momentum/InjectionModel/ConeInjection/ConeInjection.C

@@ -138,6 +138,7 @@ Foam::ConeInjection<CloudType>::ConeInjection
             this->coeffDict()
         )
     ),
+    injectorCoordinates_(barycentric::uniform(NaN)),
     injectorCell_(-1),
     injectorTetFace_(-1),
     injectorTetPt_(-1),
@@ -211,6 +212,7 @@ Foam::ConeInjection<CloudType>::ConeInjection
     position_(im.position_),
     positionIsConstant_(im.positionIsConstant_),
     direction_(im.direction_),
+    injectorCoordinates_(im.injectorCoordinates_),
     injectorCell_(im.injectorCell_),
     injectorTetFace_(im.injectorTetFace_),
     injectorTetPt_(im.injectorTetPt_),
@@ -245,10 +247,11 @@ void Foam::ConeInjection<CloudType>::topoChange()
         vector position = position_.value(0);
         this->findCellAtPosition
         (
+            position,
+            injectorCoordinates_,
             injectorCell_,
             injectorTetFace_,
-            injectorTetPt_,
+            injectorTetPt_
-            position
         );
     }
 }
@@ -307,10 +310,11 @@ void Foam::ConeInjection<CloudType>::setPositionAndCell
     const label parcelI,
     const label,
     const scalar time,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
     Random& rndGen = this->owner().rndGen();
@@ -321,10 +325,11 @@ void Foam::ConeInjection<CloudType>::setPositionAndCell
     {
         case imPoint:
         {
-            position = position_.value(t);
+            const point pos = position_.value(t);
             if (positionIsConstant_)
             {
-                cellOwner = injectorCell_;
+                coordinates = injectorCoordinates_;
+                celli = injectorCell_;
                 tetFacei = injectorTetFace_;
                 tetPti = injectorTetPt_;
             }
@@ -332,10 +337,11 @@ void Foam::ConeInjection<CloudType>::setPositionAndCell
             {
                 this->findCellAtPosition
                 (
-                    cellOwner,
+                    pos,
+                    coordinates,
+                    celli,
                     tetFacei,
                     tetPti,
-                    position,
                     false
                 );
             }
@@ -350,13 +356,14 @@ void Foam::ConeInjection<CloudType>::setPositionAndCell
             const vector t2 = normalised(n ^ t1);
             const vector tanVec = t1*cos(beta) + t2*sin(beta);
             const scalar d = sqrt((1 - frac)*sqr(dInner_) + frac*sqr(dOuter_));
-            position = position_.value(t) + d/2*tanVec;
+            const point pos = position_.value(t) + d/2*tanVec;
             this->findCellAtPosition
             (
-                cellOwner,
+                pos,
+                coordinates,
+                celli,
                 tetFacei,
                 tetPti,
-                position,
                 false
             );
             break;

src/lagrangian/parcel/submodels/Momentum/InjectionModel/ConeInjection/ConeInjection.H

@@ -195,6 +195,9 @@ private:
         //- Centreline direction in which to inject
         const TimeFunction1<vector> direction_;
 
+        //- Coordinates corresponding to the injector position
+        barycentric injectorCoordinates_;
+
         //- Cell label corresponding to the injector position
         label injectorCell_;
 
@@ -310,10 +313,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             //- Set the parcel properties

src/lagrangian/parcel/submodels/Momentum/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.C

@@ -67,6 +67,7 @@ Foam::FieldActivatedInjection<CloudType>::FieldActivatedInjection
             IOobject::NO_WRITE
         )
     ),
+    injectorCoordinates_(positions_.size()),
     injectorCells_(positions_.size()),
     injectorTetFaces_(positions_.size()),
     injectorTetPts_(positions_.size()),
@@ -112,6 +113,7 @@ Foam::FieldActivatedInjection<CloudType>::FieldActivatedInjection
     thresholdField_(im.thresholdField_),
     positionsFile_(im.positionsFile_),
     positions_(im.positions_),
+    injectorCoordinates_(im.injectorCoordinates_),
     injectorCells_(im.injectorCells_),
     injectorTetFaces_(im.injectorTetFaces_),
     injectorTetPts_(im.injectorTetPts_),
@@ -140,10 +142,11 @@ void Foam::FieldActivatedInjection<CloudType>::topoChange()
     {
         this->findCellAtPosition
         (
+            positions_[i],
+            injectorCoordinates_[i],
             injectorCells_[i],
             injectorTetFaces_[i],
-            injectorTetPts_[i],
+            injectorTetPts_[i]
-            positions_[i]
         );
     }
 }
@@ -198,14 +201,15 @@ void Foam::FieldActivatedInjection<CloudType>::setPositionAndCell
     const label parcelI,
     const label,
     const scalar,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
-    position = positions_[parcelI];
+    coordinates = injectorCoordinates_[parcelI];
-    cellOwner = injectorCells_[parcelI];
+    celli = injectorCells_[parcelI];
     tetFacei = injectorTetFaces_[parcelI];
     tetPti = injectorTetPts_[parcelI];
 }

src/lagrangian/parcel/submodels/Momentum/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.H

@@ -87,6 +87,9 @@ class FieldActivatedInjection
             //- Field of injector (x,y,z) positions
             GlobalIOField<vector> positions_;
 
+            //- List of coordinates corresponding to injector positions
+            List<barycentric> injectorCoordinates_;
+
             //- List of cell labels corresponding to injector positions
             labelList injectorCells_;
 
@@ -172,10 +175,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             //- Set the parcel properties

src/lagrangian/parcel/submodels/Momentum/InjectionModel/InflationInjection/InflationInjection.C

@@ -435,20 +435,20 @@ void Foam::InflationInjection<CloudType>::setPositionAndCell
     const label parcelI,
     const label,
     const scalar,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
-    position = newParticles_[parcelI].first().first();
-
     this->findCellAtPosition
     (
-        cellOwner,
+        newParticles_[parcelI].first().first(),
+        coordinates,
+        celli,
         tetFacei,
         tetPti,
-        position,
         false
     );
 }

src/lagrangian/parcel/submodels/Momentum/InjectionModel/InflationInjection/InflationInjection.H

@@ -167,10 +167,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             //- Set the parcel properties

src/lagrangian/parcel/submodels/Momentum/InjectionModel/InjectionModel/InjectionModel.C

@@ -91,96 +91,115 @@ bool Foam::InjectionModel<CloudType>::prepareForNextTimeStep
 template<class CloudType>
 bool Foam::InjectionModel<CloudType>::findCellAtPosition
 (
+    const point& position,
+    barycentric& coordinates,
     label& celli,
     label& tetFacei,
     label& tetPti,
-    vector& position,
     bool errorOnNotFound
 )
 {
-    const volVectorField& cellCentres = this->owner().mesh().C();
+    // Subroutine for finding the cell
-
+    auto findProcAndCell = [this](const point& pos)
-    const vector p0 = position;
-
-    this->owner().mesh().findCellFacePt
-    (
-        position,
-        celli,
-        tetFacei,
-        tetPti
-    );
-
-    label proci = -1;
-
-    if (celli >= 0)
     {
-        proci = Pstream::myProcNo();
+        // Find the containing cell
-    }
+        label celli = this->owner().mesh().findCell(pos);
-
-    reduce(proci, maxOp<label>());
-
-    // Ensure that only one processor attempts to insert this Parcel
-
-    if (proci != Pstream::myProcNo())
-    {
-        celli = -1;
-        tetFacei = -1;
-        tetPti = -1;
-    }
-
-    // Last chance - find nearest cell and try that one - the point is
-    // probably on an edge
-    if (proci == -1)
-    {
-        celli = this->owner().mesh().findNearestCell(position);
-
-        if (celli >= 0)
-        {
-            position += small*(cellCentres[celli] - position);
-
-            this->owner().mesh().findCellFacePt
-            (
-                position,
-                celli,
-                tetFacei,
-                tetPti
-            );
-
-            if (celli > 0)
-            {
-                proci = Pstream::myProcNo();
-            }
-        }
 
+        // Synchronise so only a single processor finds this position
+        label proci = celli >= 0 ? Pstream::myProcNo() : -1;
         reduce(proci, maxOp<label>());
-
         if (proci != Pstream::myProcNo())
         {
             celli = -1;
-            tetFacei = -1;
-            tetPti = -1;
         }
+
+        return labelPair(proci, celli);
+    };
+
+    point pos = position;
+
+    // Try and find the cell at the given position
+    const labelPair procAndCelli = findProcAndCell(pos);
+    label proci = procAndCelli.first();
+    celli = procAndCelli.second();
+
+    // Didn't find it. The point may be awkwardly on an edge or face. Try
+    // again, but move the point into its nearest cell a little bit.
+    if (proci == -1)
+    {
+        pos += small*(this->owner().mesh().C()[celli] - pos);
+        const labelPair procAndCelli = findProcAndCell(pos);
+        proci = procAndCelli.first();
+        celli = procAndCelli.second();
     }
 
+    // Didn't find it. Error or return false.
     if (proci == -1)
     {
         if (errorOnNotFound)
         {
             FatalErrorInFunction
                 << "Cannot find parcel injection cell. "
-                << "Parcel position = " << p0 << nl
+                << "Parcel position = " << position << nl
                 << exit(FatalError);
         }
-        else
+
-        {
+        return false;
-            return false;
+    }
-        }
+
+    // Found it. Construct the barycentric coordinates.
+    if (proci == Pstream::myProcNo())
+    {
+        particle p(this->owner().mesh(), pos, celli);
+        coordinates = p.coordinates();
+        celli = p.cell();
+        tetFacei = p.tetFace();
+        tetPti = p.tetPt();
     }
 
     return true;
 }
 
 
+template<class CloudType>
+void Foam::InjectionModel<CloudType>::constrainPosition
+(
+    typename CloudType::parcelType::trackingData& td,
+    typename CloudType::parcelType& parcel
+)
+{
+    const vector d = parcel.deviationFromMeshCentre(td.mesh);
+
+    if (d == vector::zero)
+    {
+        return;
+    }
+
+    const label facei = parcel.face();
+
+    // If the parcel is not on a face, then just track it to the mesh centre
+    if (facei == -1)
+    {
+        parcel.track(td.mesh, - d, 0);
+    }
+
+    // If the parcel is on a face, then track in two steps, going slightly into
+    // the current cell. This prevents a boundary hit from ending the track
+    // prematurely.
+    if (facei != -1)
+    {
+        const vector pc =
+            td.mesh.cellCentres()[parcel.cell()] - parcel.position(td.mesh);
+
+        parcel.track(td.mesh, - d/2 + rootSmall*pc, 0);
+        parcel.track(td.mesh, - d/2 - rootSmall*pc, 0);
+    }
+
+    // Restore any face-association that got changed during tracking
+    parcel.face() = facei;
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::InjectionModel<CloudType>::setNumberOfParticles
 (
@@ -411,6 +430,8 @@ void Foam::InjectionModel<CloudType>::inject
     typename CloudType::parcelType::trackingData& td
 )
 {
+    const polyMesh& mesh = this->owner().mesh();
+
     const scalar time = this->owner().db().time().value();
 
     // Prepare for next time step
@@ -421,12 +442,9 @@ void Foam::InjectionModel<CloudType>::inject
 
     if (prepareForNextTimeStep(time, newParcels, newVolumeFraction))
     {
-        const scalar trackTime = this->owner().solution().trackTime();
-        const polyMesh& mesh = this->owner().mesh();
-
         // Duration of injection period during this timestep
         const scalar deltaT =
-            max(0.0, min(trackTime, min(time - SOI_, timeEnd() - time0_)));
+            max(0.0, min(td.trackTime(), min(time - SOI_, timeEnd() - time0_)));
 
         // Pad injection time if injection starts during this timestep
         const scalar padTime = max(0.0, SOI_ - time0_);
@@ -439,44 +457,50 @@ void Foam::InjectionModel<CloudType>::inject
                 // Calculate the pseudo time of injection for parcel 'parcelI'
                 scalar timeInj = time0_ + padTime + deltaT*parcelI/newParcels;
 
-                // Determine the injection position and owner cell,
+                // Determine the injection coordinates and owner cell,
                 // tetFace and tetPt
-                label celli = -1;
+                barycentric coordinates = barycentric::uniform(NaN);
-                label tetFacei = -1;
+                label celli = -1, tetFacei = -1, tetPti = -1, facei = -1;
-                label tetPti = -1;
-
-                vector pos = Zero;
-
                 setPositionAndCell
                 (
                     parcelI,
                     newParcels,
                     timeInj,
-                    pos,
+                    coordinates,
                     celli,
                     tetFacei,
-                    tetPti
+                    tetPti,
+                    facei
                 );
 
                 if (celli > -1)
                 {
                     // Lagrangian timestep
-                    const scalar dt = time - timeInj;
+                    const scalar dt = timeInj - time0_;
-
-                    // Apply corrections to position for 2-D cases
-                    meshTools::constrainToMeshCentre(mesh, pos);
 
                     // Create a new parcel
-                    parcelType* pPtr = new parcelType(mesh, pos, celli);
+                    parcelType* pPtr =
+                        new parcelType
+                        (
+                            mesh,
+                            coordinates,
+                            celli,
+                            tetFacei,
+                            tetPti,
+                            facei
+                        );
+
+                    // Correct the position for reduced-dimension cases
+                    constrainPosition(td, *pPtr);
 
                     // Check/set new parcel thermo properties
-                    cloud.setParcelThermoProperties(*pPtr, dt);
+                    cloud.setParcelThermoProperties(*pPtr);
 
                     // Assign new parcel properties in injection model
                     setProperties(parcelI, newParcels, timeInj, *pPtr);
 
                     // Check/set new parcel injection properties
-                    cloud.checkParcelProperties(*pPtr, dt, fullyDescribed());
+                    cloud.checkParcelProperties(*pPtr, fullyDescribed());
 
                     // Apply correction to velocity for 2-D cases
                     meshTools::constrainDirection
@@ -496,18 +520,14 @@ void Foam::InjectionModel<CloudType>::inject
                             pPtr->rho()
                         );
 
+                    // Modify the step fraction so that the particles are
+                    // injected continually through the time-step
+                    pPtr->stepFraction() = dt/td.trackTime();
+
+                    // Add the new parcel
                     parcelsAdded ++;
-
                     massAdded += pPtr->nParticle()*pPtr->mass();
-
+                    cloud.addParticle(pPtr);
-                    if (pPtr->move(cloud, td, dt))
-                    {
-                        cloud.addParticle(pPtr);
-                    }
-                    else
-                    {
-                        delete pPtr;
-                    }
                 }
             }
         }
@@ -522,8 +542,7 @@ template<class TrackCloudType>
 void Foam::InjectionModel<CloudType>::injectSteadyState
 (
     TrackCloudType& cloud,
-    typename CloudType::parcelType::trackingData& td,
+    typename CloudType::parcelType::trackingData& td
-    const scalar trackTime
 )
 {
     const polyMesh& mesh = this->owner().mesh();
@@ -544,41 +563,47 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
         // Volume to inject is split equally amongst all parcel streams
         scalar newVolumeFraction = 1.0/scalar(newParcels);
 
-        // Determine the injection position and owner cell,
+        // Determine the injection coordinates and owner cell,
         // tetFace and tetPt
-        label celli = -1;
+        barycentric coordinates = barycentric::uniform(NaN);
-        label tetFacei = -1;
+        label celli = -1, tetFacei = -1, tetPti = -1, facei = -1;
-        label tetPti = -1;
-
-        vector pos = Zero;
-
         setPositionAndCell
         (
             parcelI,
             newParcels,
-            0.0,
+            0,
-            pos,
+            coordinates,
             celli,
             tetFacei,
-            tetPti
+            tetPti,
+            facei
         );
 
         if (celli > -1)
         {
-            // Apply corrections to position for 2-D cases
-            meshTools::constrainToMeshCentre(mesh, pos);
-
             // Create a new parcel
-            parcelType* pPtr = new parcelType(mesh, pos, celli);
+            parcelType* pPtr =
+                new parcelType
+                (
+                    mesh,
+                    coordinates,
+                    celli,
+                    tetFacei,
+                    tetPti,
+                    facei
+                );
+
+            // Correct the position for reduced-dimension cases
+            constrainPosition(td, *pPtr);
 
             // Check/set new parcel thermo properties
-            cloud.setParcelThermoProperties(*pPtr, 0.0);
+            cloud.setParcelThermoProperties(*pPtr);
 
             // Assign new parcel properties in injection model
             setProperties(parcelI, newParcels, 0.0, *pPtr);
 
             // Check/set new parcel injection properties
-            cloud.checkParcelProperties(*pPtr, 0.0, fullyDescribed());
+            cloud.checkParcelProperties(*pPtr, fullyDescribed());
 
             // Apply correction to velocity for 2-D cases
             meshTools::constrainDirection(mesh, mesh.solutionD(), pPtr->U());
@@ -593,11 +618,13 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
                     pPtr->rho()
                 );
 
-            // Add the new parcel
+            // Initial step fraction
-            cloud.addParticle(pPtr);
+            pPtr->stepFraction() = 0;
 
+            // Add the new parcel
+            parcelsAdded ++;
             massAdded += pPtr->nParticle()*pPtr->mass();
-            parcelsAdded++;
+            cloud.addParticle(pPtr);
         }
     }
 

src/lagrangian/parcel/submodels/Momentum/InjectionModel/InjectionModel/InjectionModel.H

@@ -53,6 +53,7 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "CloudSubModelBase.H"
 #include "vector.H"
+#include "particle.H"
 #include "TimeFunction1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -142,7 +143,7 @@ protected:
         virtual bool validInjection(const label parcelI) = 0;
 
         //- Determine properties for next time step/injection interval
-        virtual bool prepareForNextTimeStep
+        bool prepareForNextTimeStep
         (
             const scalar time,
             label& newParcels,
@@ -152,17 +153,26 @@ protected:
         //- Find the cell that contains the supplied position
         //  Will modify position slightly towards the owner cell centroid to
         //  ensure that it lies in a cell and not edge/face
-        virtual bool findCellAtPosition
+        bool findCellAtPosition
         (
+            const point& position,
+            barycentric& coordinates,
             label& celli,
             label& tetFacei,
             label& tetPti,
-            vector& position,
             bool errorOnNotFound = true
         );
 
+        //- Constrain a parcel's position appropriately to the geometric
+        //  dimensions of the mesh
+        void constrainPosition
+        (
+            typename CloudType::parcelType::trackingData& td,
+            typename CloudType::parcelType& parcel
+        );
+
         //- Set number of particles to inject given parcel properties
-        virtual scalar setNumberOfParticles
+        scalar setNumberOfParticles
         (
             const label parcels,
             const scalar volumeFraction,
@@ -171,7 +181,7 @@ protected:
         );
 
         //- Post injection checks
-        virtual void postInjectCheck
+        void postInjectCheck
         (
             const label parcelsAdded,
             const scalar massAdded
@@ -309,8 +319,7 @@ public:
             void injectSteadyState
             (
                 TrackCloudType& cloud,
-                typename CloudType::parcelType::trackingData& td,
+                typename CloudType::parcelType::trackingData& td
-                const scalar trackTime
             );
 
 
@@ -322,10 +331,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             ) = 0;
 
             //- Set the parcel properties

src/lagrangian/parcel/submodels/Momentum/InjectionModel/InjectionModel/InjectionModelList.C

@@ -198,13 +198,12 @@ template<class TrackCloudType>
 void Foam::InjectionModelList<CloudType>::injectSteadyState
 (
     TrackCloudType& cloud,
-    typename CloudType::parcelType::trackingData& td,
+    typename CloudType::parcelType::trackingData& td
-    const scalar trackTime
 )
 {
     forAll(*this, i)
     {
-        this->operator[](i).injectSteadyState(cloud, td, trackTime);
+        this->operator[](i).injectSteadyState(cloud, td);
     }
 }
 

src/lagrangian/parcel/submodels/Momentum/InjectionModel/InjectionModel/InjectionModelList.H

@@ -118,8 +118,7 @@ public:
             void injectSteadyState
             (
                 TrackCloudType& cloud,
-                typename CloudType::parcelType::trackingData& td,
+                typename CloudType::parcelType::trackingData& td
-                const scalar trackTime
             );
 
 

src/lagrangian/parcel/submodels/Momentum/InjectionModel/ManualInjection/ManualInjection.C

@@ -53,6 +53,7 @@ Foam::ManualInjection<CloudType>::ManualInjection
         )
     ),
     diameters_(positions_.size()),
+    injectorCoordinates_(positions_.size(), barycentric::uniform(NaN)),
     injectorCells_(positions_.size(), -1),
     injectorTetFaces_(positions_.size(), -1),
     injectorTetPts_(positions_.size(), -1),
@@ -93,6 +94,7 @@ Foam::ManualInjection<CloudType>::ManualInjection
     positionsFile_(im.positionsFile_),
     positions_(im.positions_),
     diameters_(im.diameters_),
+    injectorCoordinates_(im.injectorCoordinates_),
     injectorCells_(im.injectorCells_),
     injectorTetFaces_(im.injectorTetFaces_),
     injectorTetPts_(im.injectorTetPts_),
@@ -124,10 +126,11 @@ void Foam::ManualInjection<CloudType>::topoChange()
         (
             !this->findCellAtPosition
             (
+                positions_[pI],
+                injectorCoordinates_[pI],
                 injectorCells_[pI],
                 injectorTetFaces_[pI],
                 injectorTetPts_[pI],
-                positions_[pI],
                 !ignoreOutOfBounds_
             )
         )
@@ -141,6 +144,7 @@ void Foam::ManualInjection<CloudType>::topoChange()
     {
         inplaceSubset(keep, positions_);
         inplaceSubset(keep, diameters_);
+        inplaceSubset(keep, injectorCoordinates_);
         inplaceSubset(keep, injectorCells_);
         inplaceSubset(keep, injectorTetFaces_);
         inplaceSubset(keep, injectorTetPts_);
@@ -203,14 +207,15 @@ void Foam::ManualInjection<CloudType>::setPositionAndCell
     const label parcelI,
     const label,
     const scalar,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
-    position = positions_[parcelI];
+    coordinates = injectorCoordinates_[parcelI];
-    cellOwner = injectorCells_[parcelI];
+    celli = injectorCells_[parcelI];
     tetFacei = injectorTetFaces_[parcelI];
     tetPti = injectorTetPts_[parcelI];
 }

src/lagrangian/parcel/submodels/Momentum/InjectionModel/ManualInjection/ManualInjection.H

@@ -74,6 +74,9 @@ class ManualInjection
         //- Field of parcel diameters
         scalarList diameters_;
 
+        //- List of coordinates of injector positions
+        List<barycentric> injectorCoordinates_;
+
         //- List of cell labels corresponding to injector positions
         labelList injectorCells_;
 
@@ -149,10 +152,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             //- Set the parcel properties

src/lagrangian/parcel/submodels/Momentum/InjectionModel/MomentumLookupTableInjection/MomentumLookupTableInjection.C

@@ -55,6 +55,7 @@ Foam::MomentumLookupTableInjection<CloudType>::MomentumLookupTableInjection
             IOobject::NO_WRITE
         )
     ),
+    injectorCoordinates_(0),
     injectorCells_(0),
     injectorTetFaces_(0),
     injectorTetPts_(0)
@@ -90,6 +91,7 @@ Foam::MomentumLookupTableInjection<CloudType>::MomentumLookupTableInjection
     parcelsPerSecond_(im.parcelsPerSecond_),
     randomise_(im.randomise_),
     injectors_(im.injectors_),
+    injectorCoordinates_(im.injectorCoordinates_),
     injectorCells_(im.injectorCells_),
     injectorTetFaces_(im.injectorTetFaces_),
     injectorTetPts_(im.injectorTetPts_)
@@ -113,10 +115,11 @@ void Foam::MomentumLookupTableInjection<CloudType>::topoChange()
     {
         this->findCellAtPosition
         (
+            injectors_[i].x(),
+            injectorCoordinates_[i],
             injectorCells_[i],
             injectorTetFaces_[i],
-            injectorTetPts_[i],
+            injectorTetPts_[i]
-            injectors_[i].x()
         );
     }
 }
@@ -173,10 +176,11 @@ void Foam::MomentumLookupTableInjection<CloudType>::setPositionAndCell
     const label parcelI,
     const label nParcels,
     const scalar time,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
     label injectorI = 0;
@@ -190,8 +194,8 @@ void Foam::MomentumLookupTableInjection<CloudType>::setPositionAndCell
         injectorI = int64_t(parcelI)*int64_t(injectors_.size())/nParcels;
     }
 
-    position = injectors_[injectorI].x();
+    coordinates = injectorCoordinates_[injectorI];
-    cellOwner = injectorCells_[injectorI];
+    celli = injectorCells_[injectorI];
     tetFacei = injectorTetFaces_[injectorI];
     tetPti = injectorTetPts_[injectorI];
 }

src/lagrangian/parcel/submodels/Momentum/InjectionModel/MomentumLookupTableInjection/MomentumLookupTableInjection.H

@@ -88,6 +88,9 @@ class MomentumLookupTableInjection
         //- List of injectors
         momentumParcelInjectionDataIOList injectors_;
 
+        //- List of coordinates corresponding to injector positions
+        List<barycentric> injectorCoordinates_;
+
         //- List of cell labels corresponding to injector positions
         labelList injectorCells_;
 
@@ -157,10 +160,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             //- Set the parcel properties

src/lagrangian/parcel/submodels/Momentum/InjectionModel/NoInjection/NoInjection.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,7 +91,8 @@ void Foam::NoInjection<CloudType>::setPositionAndCell
     const label,
     const label,
     const scalar,
-    vector&,
+    barycentric&,
+    label&,
     label&,
     label&,
     label&

src/lagrangian/parcel/submodels/Momentum/InjectionModel/NoInjection/NoInjection.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -100,10 +100,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             virtual void setProperties

src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C

@@ -223,20 +223,22 @@ void Foam::PatchFlowRateInjection<CloudType>::setPositionAndCell
     const label,
     const label,
     const scalar,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
     patchInjectionBase::setPositionAndCell
     (
         this->owner().mesh(),
         this->owner().rndGen(),
-        position,
+        coordinates,
-        cellOwner,
+        celli,
         tetFacei,
-        tetPti
+        tetPti,
+        facei
     );
 }
 

src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.H

@@ -157,10 +157,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             virtual void setProperties

src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/PatchInjection.C

@@ -170,20 +170,22 @@ void Foam::PatchInjection<CloudType>::setPositionAndCell
     const label,
     const label,
     const scalar,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
     patchInjectionBase::setPositionAndCell
     (
         this->owner().mesh(),
         this->owner().rndGen(),
-        position,
+        coordinates,
-        cellOwner,
+        celli,
         tetFacei,
-        tetPti
+        tetPti,
+        facei
     );
 }
 

src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/PatchInjection.H

@@ -149,10 +149,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             virtual void setProperties

src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/patchInjectionBase.C

@@ -42,13 +42,9 @@ Foam::patchInjectionBase::patchInjectionBase
 :
     patchName_(patchName),
     patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
-    patchArea_(0.0),
+    sumProcArea_(),
-    patchNormal_(),
+    sumFaceArea_(),
-    cellOwners_(),
+    sumFaceTriArea_()
-    triFace_(),
-    triToFace_(),
-    triCumulativeMagSf_(),
-    sumTriMagSf_(Pstream::nProcs() + 1, 0.0)
 {
     if (patchId_ < 0)
     {
@@ -66,13 +62,9 @@ Foam::patchInjectionBase::patchInjectionBase(const patchInjectionBase& pib)
 :
     patchName_(pib.patchName_),
     patchId_(pib.patchId_),
-    patchArea_(pib.patchArea_),
+    sumProcArea_(pib.sumProcArea_),
-    patchNormal_(pib.patchNormal_),
+    sumFaceArea_(pib.sumFaceArea_),
-    cellOwners_(pib.cellOwners_),
+    sumFaceTriArea_(pib.sumFaceTriArea_)
-    triFace_(pib.triFace_),
-    triToFace_(pib.triToFace_),
-    triCumulativeMagSf_(pib.triCumulativeMagSf_),
-    sumTriMagSf_(pib.sumTriMagSf_)
 {}
 
 
@@ -88,65 +80,49 @@ void Foam::patchInjectionBase::topoChange(const polyMesh& mesh)
 {
     // Set/cache the injector cells
     const polyPatch& patch = mesh.boundaryMesh()[patchId_];
-    const pointField& points = patch.points();
 
-    cellOwners_ = patch.faceCells();
+    // Initialise
-
+    sumProcArea_.resize(Pstream::nProcs(), 0);
-    // Triangulate the patch faces and create addressing
+    sumFaceArea_.resize(patch.size(), 0);
-    DynamicList<label> triToFace(2*patch.size());
+    sumFaceTriArea_.resize(patch.size());
-    DynamicList<scalar> triMagSf(2*patch.size());
+    forAll(patch, patchFacei)
-    DynamicList<triFace> triFace(2*patch.size());
-
-    // Set zero value at the start of the tri area list
-    triMagSf.append(0.0);
-
-    // Create a triangulation engine
-    polygonTriangulate triEngine;
-
-    forAll(patch, facei)
     {
-        const face& f = patch[facei];
+        sumFaceTriArea_[patchFacei].resize(patch[patchFacei].nTriangles(), 0);
+    }
 
-        triEngine.triangulate(UIndirectList<point>(points, f));
+    // Cumulatively sum up the areas through the local patch
+    scalar patchArea = 0;
+    forAll(patch, patchFacei)
+    {
+        const label facei = patchFacei + patch.start();
+        const label celli = patch.faceCells()[patchFacei];
 
-        forAll(triEngine.triPoints(), i)
+        scalar patchFaceArea = 0;
+        for
+        (
+            label patchFaceTrii = 0;
+            patchFaceTrii < patch[patchFacei].nTriangles();
+            ++ patchFaceTrii
+        )
         {
-            triToFace.append(facei);
+            const tetIndices tet(celli, facei, patchFaceTrii + 1);
-            triFace.append(triEngine.triPoints(i, f));
+
-            triMagSf.append(triEngine.triPoints(i, f).mag(points));
+            patchFaceArea += tet.faceTri(mesh).mag();
+            sumFaceTriArea_[patchFacei][patchFaceTrii] = patchFaceArea;
         }
+
+        patchArea += patchFaceArea;
+        sumFaceArea_[patchFacei] = patchArea;
     }
 
-    forAll(sumTriMagSf_, i)
+    // Cumulatively sum the total areas across the processors
+    sumProcArea_[Pstream::myProcNo()] = patchArea;
+    Pstream::listCombineGather(sumProcArea_, maxEqOp<scalar>());
+    Pstream::listCombineScatter(sumProcArea_);
+    for (label proci = 1; proci < Pstream::nProcs(); proci ++)
     {
-        sumTriMagSf_[i] = 0.0;
+        sumProcArea_[proci] += sumProcArea_[proci - 1];
     }
-
-    sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
-
-    Pstream::listCombineGather(sumTriMagSf_, maxEqOp<scalar>());
-    Pstream::listCombineScatter(sumTriMagSf_);
-
-    for (label i = 1; i < triMagSf.size(); i++)
-    {
-        triMagSf[i] += triMagSf[i-1];
-    }
-
-    // Transfer to persistent storage
-    triFace_.transfer(triFace);
-    triToFace_.transfer(triToFace);
-    triCumulativeMagSf_.transfer(triMagSf);
-
-    // Convert sumTriMagSf_ into cumulative sum of areas per proc
-    for (label i = 1; i < sumTriMagSf_.size(); i++)
-    {
-        sumTriMagSf_[i] += sumTriMagSf_[i-1];
-    }
-
-    const scalarField magSf(patch.magFaceAreas());
-    patchArea_ = sum(magSf);
-    patchNormal_ = patch.faceAreas()/magSf;
-    reduce(patchArea_, sumOp<scalar>());
 }
 
 
@@ -154,125 +130,75 @@ void Foam::patchInjectionBase::setPositionAndCell
 (
     const fvMesh& mesh,
     Random& rnd,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
-    scalar areaFraction = rnd.globalScalar01()*patchArea_;
+    // Linear searching for area fractions
-
+    auto findArea = []
-    if (cellOwners_.size() > 0)
+    (
+        const scalarList& areas,
+        scalar& area
+    )
     {
-        // Determine which processor to inject from
+        for (label i = areas.size(); i > 0; i --)
-        label proci = 0;
-        forAllReverse(sumTriMagSf_, i)
         {
-            if (areaFraction >= sumTriMagSf_[i])
+            if (area >= areas[i - 1])
             {
-                proci = i;
+                area -= areas[i - 1];
-                break;
+                return i;
             }
         }
 
+        return 0;
+    };
+
+    const polyPatch& patch = mesh.boundaryMesh()[patchId_];
+
+    scalar area = rnd.globalScalar01()*sumProcArea_.last();
+
+    if (patch.size() > 0)
+    {
+        // Determine which processor to inject from
+        const label proci = findArea(sumProcArea_, area);
+
+        // If this processor...
         if (Pstream::myProcNo() == proci)
         {
-            // Find corresponding decomposed face triangle
+            // Determine the face to inject from
-            label trii = 0;
+            const label patchFacei = findArea(sumFaceArea_, area);
-            scalar offset = sumTriMagSf_[proci];
-            forAllReverse(triCumulativeMagSf_, i)
-            {
-                if (areaFraction > triCumulativeMagSf_[i] + offset)
-                {
-                    trii = i;
-                    break;
-                }
-            }
 
-            // Set cellOwner
+            // Determine the face-tri to inject from
-            label facei = triToFace_[trii];
+            const label patchFaceTrii =
-            cellOwner = cellOwners_[facei];
+                findArea(sumFaceTriArea_[patchFacei], area);
 
-            // Find random point in triangle
+            // Set the topology
-            const polyPatch& patch = mesh.boundaryMesh()[patchId_];
+            const barycentric2D r = barycentric2D01(rnd);
-            const pointField& points = patch.points();
+            coordinates = barycentric(0, r.a(), r.b(), r.c());
-            const face& tf = triFace_[trii];
+            celli = mesh.faceOwner()[patch.start() + patchFacei];
-            const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
+            tetFacei = patch.start() + patchFacei;
-            const point pf(tri.randomPoint(rnd));
+            tetPti = patchFaceTrii + 1;
-
+            facei = patch.start() + patchFacei;
-            // Position perturbed away from face (into domain)
-            const scalar a = rnd.scalarAB(0.1, 0.5);
-            const vector& pc = mesh.cellCentres()[cellOwner];
-            const vector d =
-                mag((pf - pc) & patchNormal_[facei])*patchNormal_[facei];
-
-            position = pf - a*d;
-
-            // Try to find tetFacei and tetPti in the current position
-            mesh.findTetFacePt(cellOwner, position, tetFacei, tetPti);
-
-            // tetFacei and tetPti not found, check if the cell has changed
-            if (tetFacei == -1 ||tetPti == -1)
-            {
-                mesh.findCellFacePt(position, cellOwner, tetFacei, tetPti);
-            }
-
-            // Both searches failed, choose a random position within
-            // the original cell
-            if (tetFacei == -1 ||tetPti == -1)
-            {
-                // Reset cellOwner
-                cellOwner = cellOwners_[facei];
-                const scalarField& V = mesh.V();
-
-                // Construct cell tet indices
-                const List<tetIndices> cellTetIs =
-                    polyMeshTetDecomposition::cellTetIndices(mesh, cellOwner);
-
-                // Construct cell tet volume fractions
-                scalarList cTetVFrac(cellTetIs.size(), 0.0);
-                for (label teti=1; teti<cellTetIs.size()-1; teti++)
-                {
-                    cTetVFrac[teti] =
-                        cTetVFrac[teti-1]
-                      + cellTetIs[teti].tet(mesh).mag()/V[cellOwner];
-                }
-                cTetVFrac.last() = 1;
-
-                // Set new particle position
-                const scalar volFrac = rnd.sample01<scalar>();
-                label teti = 0;
-                forAll(cTetVFrac, vfI)
-                {
-                    if (cTetVFrac[vfI] > volFrac)
-                    {
-                        teti = vfI;
-                        break;
-                    }
-                }
-                position = cellTetIs[teti].tet(mesh).randomPoint(rnd);
-                tetFacei = cellTetIs[teti].face();
-                tetPti = cellTetIs[teti].tetPt();
-            }
         }
         else
         {
-            cellOwner = -1;
+            coordinates = barycentric::uniform(NaN);
+            celli = -1;
             tetFacei = -1;
             tetPti = -1;
-
+            facei = -1;
-            // Dummy position
-            position = pTraits<vector>::max;
         }
     }
     else
     {
-        cellOwner = -1;
+        coordinates = barycentric::uniform(NaN);
+        celli = -1;
         tetFacei = -1;
         tetPti = -1;
-
+        facei = -1;
-        // Dummy position
-        position = pTraits<vector>::max;
     }
 }
 

src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/patchInjectionBase.H

@@ -46,6 +46,7 @@ SourceFiles
 #include "vectorList.H"
 #include "faceList.H"
 #include "triFaceList.H"
+#include "barycentric.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,26 +74,14 @@ protected:
         //- Patch ID
         const label patchId_;
 
-        //- Patch area - total across all processors
+        //- ...
-        scalar patchArea_;
+        scalarList sumProcArea_;
 
-        //- Patch face normal directions
+        //- ...
-        vectorList patchNormal_;
+        scalarList sumFaceArea_;
 
-        //- List of cell labels corresponding to injector positions
+        //- ...
-        labelList cellOwners_;
+        scalarListList sumFaceTriArea_;
-
-        //- Decomposed patch faces as a list of triangles
-        triFaceList triFace_;
-
-        //- Addressing from per triangle to patch face
-        labelList triToFace_;
-
-        //- Cumulative triangle area per triangle face
-        scalarList triCumulativeMagSf_;
-
-        //- Cumulative area fractions per processor
-        scalarList sumTriMagSf_;
 
 
 public:
@@ -120,10 +109,11 @@ public:
         (
             const fvMesh& mesh,
             Random& rnd,
-            vector& position,
+            barycentric& coordinates,
-            label& cellOwner,
+            label& celli,
             label& tetFacei,
-            label& tetPti
+            label& tetPti,
+            label& facei
         );
 };
 

src/lagrangian/parcel/submodels/Momentum/StochasticCollision/NoStochasticCollision/NoStochasticCollision.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,8 +30,7 @@ License
 template<class CloudType>
 void Foam::NoStochasticCollision<CloudType>::collide
 (
-    typename CloudType::parcelType::trackingData&,
+    typename CloudType::parcelType::trackingData&
-    const scalar
 )
 {}
 

src/lagrangian/parcel/submodels/Momentum/StochasticCollision/NoStochasticCollision/NoStochasticCollision.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,11 +52,7 @@ protected:
     // Protected Member Functions
 
         //- Update the model
-        virtual void collide
+        virtual void collide(typename CloudType::parcelType::trackingData& td);
-        (
-            typename CloudType::parcelType::trackingData& td,
-            const scalar dt
-        );
 
 
 public:

src/lagrangian/parcel/submodels/Momentum/StochasticCollision/StochasticCollisionModel/StochasticCollisionModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,11 +71,10 @@ Foam::StochasticCollisionModel<CloudType>::~StochasticCollisionModel()
 template<class CloudType>
 void Foam::StochasticCollisionModel<CloudType>::update
 (
-    typename CloudType::parcelType::trackingData& td,
+    typename CloudType::parcelType::trackingData& td
-    const scalar dt
 )
 {
-    this->collide(td, dt);
+    this->collide(td);
 }
 
 

src/lagrangian/parcel/submodels/Momentum/StochasticCollision/StochasticCollisionModel/StochasticCollisionModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,11 +58,7 @@ class StochasticCollisionModel
 protected:
 
     //- Main collision routine
-    virtual void collide
+    virtual void collide(typename CloudType::parcelType::trackingData& td) = 0;
-    (
-        typename CloudType::parcelType::trackingData& td,
-        const scalar dt
-    ) = 0;
 
 
 public:
@@ -119,11 +115,7 @@ public:
     // Member Functions
 
         //- Update the model
-        void update
+        void update(typename CloudType::parcelType::trackingData& td);
-        (
-            typename CloudType::parcelType::trackingData& td,
-            const scalar dt
-        );
 };
 
 

src/lagrangian/parcel/submodels/Momentum/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C

@@ -201,14 +201,14 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackCloudType& cloud)
                         new parcelType(this->owner().pMesh(), pos, celli);
 
                     // Check/set new parcel thermo properties
-                    cloud.setParcelThermoProperties(*pPtr, 0.0);
+                    cloud.setParcelThermoProperties(*pPtr);
 
                     setParcelProperties(*pPtr, j);
 
                     if (pPtr->nParticle() > 0.001)
                     {
                         // Check new parcel properties
-                        cloud.checkParcelProperties(*pPtr, 0.0, false);
+                        cloud.checkParcelProperties(*pPtr, false);
 
                         // Add the new parcel to the cloud
                         cloud.addParticle(pPtr);

src/lagrangian/parcel/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.C

@@ -54,6 +54,7 @@ Foam::ReactingLookupTableInjection<CloudType>::ReactingLookupTableInjection
             IOobject::NO_WRITE
         )
     ),
+    injectorCoordinates_(0),
     injectorCells_(0),
     injectorTetFaces_(0),
     injectorTetPts_(0)
@@ -61,6 +62,7 @@ Foam::ReactingLookupTableInjection<CloudType>::ReactingLookupTableInjection
     duration_ = owner.db().time().userTimeToTime(duration_);
 
     // Set/cache the injector cells
+    injectorCoordinates_.setSize(injectors_.size());
     injectorCells_.setSize(injectors_.size());
     injectorTetFaces_.setSize(injectors_.size());
     injectorTetPts_.setSize(injectors_.size());
@@ -89,6 +91,7 @@ Foam::ReactingLookupTableInjection<CloudType>::ReactingLookupTableInjection
     parcelsPerSecond_(im.parcelsPerSecond_),
     randomise_(im.randomise_),
     injectors_(im.injectors_),
+    injectorCoordinates_(im.injectorCoordinates_),
     injectorCells_(im.injectorCells_),
     injectorTetFaces_(im.injectorTetFaces_),
     injectorTetPts_(im.injectorTetPts_)
@@ -112,10 +115,11 @@ void Foam::ReactingLookupTableInjection<CloudType>::topoChange()
     {
         this->findCellAtPosition
         (
+            injectors_[i].x(),
+            injectorCoordinates_[i],
             injectorCells_[i],
             injectorTetFaces_[i],
-            injectorTetPts_[i],
+            injectorTetPts_[i]
-            injectors_[i].x()
         );
     }
 }
@@ -172,10 +176,11 @@ void Foam::ReactingLookupTableInjection<CloudType>::setPositionAndCell
     const label parcelI,
     const label nParcels,
     const scalar time,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
     label injectorI = 0;
@@ -189,8 +194,8 @@ void Foam::ReactingLookupTableInjection<CloudType>::setPositionAndCell
         injectorI = parcelI*injectorCells_.size()/nParcels;
     }
 
-    position = injectors_[injectorI].x();
+    coordinates = injectorCoordinates_[injectorI];
-    cellOwner = injectorCells_[injectorI];
+    celli = injectorCells_[injectorI];
     tetFacei = injectorTetFaces_[injectorI];
     tetPti = injectorTetPts_[injectorI];
 }

src/lagrangian/parcel/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.H

@@ -87,6 +87,9 @@ class ReactingLookupTableInjection
         //- List of injectors
         reactingParcelInjectionDataIOList injectors_;
 
+        //- List of coordinates corresponding to injector positions
+        List<barycentric> injectorCoordinates_;
+
         //- List of cell labels corresponding to injector positions
         labelList injectorCells_;
 
@@ -157,10 +160,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             //- Set the parcel properties

src/lagrangian/parcel/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.C

@@ -55,6 +55,7 @@ ReactingMultiphaseLookupTableInjection
             IOobject::NO_WRITE
         )
     ),
+    injectorCoordinates_(0),
     injectorCells_(0),
     injectorTetFaces_(0),
     injectorTetPts_(0)
@@ -62,6 +63,7 @@ ReactingMultiphaseLookupTableInjection
     duration_ = owner.db().time().userTimeToTime(duration_);
 
     // Set/cache the injector cells
+    injectorCoordinates_.setSize(injectors_.size());
     injectorCells_.setSize(injectors_.size());
     injectorTetFaces_.setSize(injectors_.size());
     injectorTetPts_.setSize(injectors_.size());
@@ -91,6 +93,7 @@ ReactingMultiphaseLookupTableInjection
     parcelsPerSecond_(im.parcelsPerSecond_),
     randomise_(im.randomise_),
     injectors_(im.injectors_),
+    injectorCoordinates_(im.injectorCoordinates_),
     injectorCells_(im.injectorCells_),
     injectorTetFaces_(im.injectorTetFaces_),
     injectorTetPts_(im.injectorTetPts_)
@@ -115,10 +118,11 @@ void Foam::ReactingMultiphaseLookupTableInjection<CloudType>::topoChange()
     {
         this->findCellAtPosition
         (
+            injectors_[i].x(),
+            injectorCoordinates_[i],
             injectorCells_[i],
             injectorTetFaces_[i],
-            injectorTetPts_[i],
+            injectorTetPts_[i]
-            injectors_[i].x()
         );
     }
 }
@@ -178,10 +182,11 @@ void Foam::ReactingMultiphaseLookupTableInjection<CloudType>::setPositionAndCell
     const label parcelI,
     const label nParcels,
     const scalar time,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
     label injectorI = 0;
@@ -195,8 +200,8 @@ void Foam::ReactingMultiphaseLookupTableInjection<CloudType>::setPositionAndCell
         injectorI = parcelI*injectorCells_.size()/nParcels;
     }
 
-    position = injectors_[injectorI].x();
+    coordinates = injectorCoordinates_[injectorI];
-    cellOwner = injectorCells_[injectorI];
+    celli = injectorCells_[injectorI];
     tetFacei = injectorTetFaces_[injectorI];
     tetPti = injectorTetPts_[injectorI];
 }

src/lagrangian/parcel/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.H

@@ -90,6 +90,9 @@ class ReactingMultiphaseLookupTableInjection
         //- List of injectors
         reactingMultiphaseParcelInjectionDataIOList injectors_;
 
+        //- List of coordinates corresponding to injector positions
+        List<barycentric> injectorCoordinates_;
+
         //- List of cell labels corresponding to injector positions
         labelList injectorCells_;
 
@@ -159,10 +162,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
-                label& cellOwner,
+                label& celli,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             //- Set the parcel properties

src/lagrangian/parcel/submodels/ReactingMultiphase/StochasticCollision/SuppressionCollision/SuppressionCollision.C

@@ -31,8 +31,7 @@ License
 template<class CloudType>
 void Foam::SuppressionCollision<CloudType>::collide
 (
-    typename CloudType::parcelType::trackingData& td,
+    typename CloudType::parcelType::trackingData& td
-    const scalar dt
 )
 {
     const parcelCloud& sc =
@@ -41,7 +40,7 @@ void Foam::SuppressionCollision<CloudType>::collide
 
     volScalarField vDotSweep(sc.vDotSweep());
 
-    dimensionedScalar Dt("dt", dimTime, dt);
+    dimensionedScalar Dt("dt", dimTime, td.trackTime());
     volScalarField P(typedName("p"), 1.0 - exp(-vDotSweep*Dt));
 
     forAllIter(typename CloudType, this->owner(), iter)

src/lagrangian/parcel/submodels/ReactingMultiphase/StochasticCollision/SuppressionCollision/SuppressionCollision.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -62,11 +62,7 @@ protected:
    // Protected Member Functions
 
         //- Update the model
-        virtual void collide
+        virtual void collide(typename CloudType::parcelType::trackingData& td);
-        (
-            typename CloudType::parcelType::trackingData& td,
-            const scalar dt
-        );
 
 
 public:

src/lagrangian/parcel/submodels/Spray/StochasticCollision/ORourkeCollision/ORourkeCollision.C

@@ -35,8 +35,7 @@ using namespace Foam::constant::mathematical;
 template<class CloudType>
 void Foam::ORourkeCollision<CloudType>::collide
 (
-    typename CloudType::parcelType::trackingData& td,
+    typename CloudType::parcelType::trackingData& td
-    const scalar dt
 )
 {
     const liquidMixtureProperties& liquids =
@@ -78,7 +77,8 @@ void Foam::ORourkeCollision<CloudType>::collide
                     scalar m1 = p1.nParticle()*p1.mass();
                     scalar m2 = p2.nParticle()*p2.mass();
 
-                    bool massChanged = collideParcels(dt, p1, p2, m1, m2);
+                    bool massChanged =
+                        collideParcels(td.trackTime(), p1, p2, m1, m2);
 
                     if (massChanged)
                     {

src/lagrangian/parcel/submodels/Spray/StochasticCollision/ORourkeCollision/ORourkeCollision.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,11 +63,7 @@ protected:
     // Protected Member Functions
 
         //- Main collision routine
-        virtual void collide
+        virtual void collide(typename CloudType::parcelType::trackingData& td);
-        (
-            typename CloudType::parcelType::trackingData& td,
-            const scalar dt
-        );
 
         //- Collide parcels and return true if mass has changed
         virtual bool collideParcels

src/lagrangian/parcel/submodels/Spray/StochasticCollision/TrajectoryCollision/TrajectoryCollision.C

@@ -30,11 +30,10 @@ License
 template<class CloudType>
 void Foam::TrajectoryCollision<CloudType>::collide
 (
-    typename CloudType::parcelType::trackingData& td,
+    typename CloudType::parcelType::trackingData& td
-    const scalar dt
 )
 {
-    ORourkeCollision<CloudType>::collide(td, dt);
+    ORourkeCollision<CloudType>::collide(td);
 }
 
 

src/lagrangian/parcel/submodels/Spray/StochasticCollision/TrajectoryCollision/TrajectoryCollision.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,11 +65,7 @@ protected:
     // Protected Member Functions
 
         //- Main collision routine
-        virtual void collide
+        virtual void collide(typename CloudType::parcelType::trackingData& td);
-        (
-            typename CloudType::parcelType::trackingData& td,
-            const scalar dt
-        );
 
         //- Collide parcels and return true if mass has changed
         virtual bool collideParcels

src/lagrangian/parcel/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.C

@@ -55,6 +55,7 @@ Foam::ThermoLookupTableInjection<CloudType>::ThermoLookupTableInjection
             IOobject::NO_WRITE
         )
     ),
+    injectorCoordinates_(0),
     injectorCells_(0),
     injectorTetFaces_(0),
     injectorTetPts_(0)
@@ -62,6 +63,7 @@ Foam::ThermoLookupTableInjection<CloudType>::ThermoLookupTableInjection
     duration_ = owner.db().time().userTimeToTime(duration_);
 
     // Set/cache the injector cells
+    injectorCoordinates_.setSize(injectors_.size());
     injectorCells_.setSize(injectors_.size());
     injectorTetFaces_.setSize(injectors_.size());
     injectorTetPts_.setSize(injectors_.size());
@@ -90,6 +92,7 @@ Foam::ThermoLookupTableInjection<CloudType>::ThermoLookupTableInjection
     parcelsPerSecond_(im.parcelsPerSecond_),
     randomise_(im.randomise_),
     injectors_(im.injectors_),
+    injectorCoordinates_(im.injectorCoordinates_),
     injectorCells_(im.injectorCells_),
     injectorTetFaces_(im.injectorTetFaces_),
     injectorTetPts_(im.injectorTetPts_)
@@ -113,10 +116,11 @@ void Foam::ThermoLookupTableInjection<CloudType>::topoChange()
     {
         this->findCellAtPosition
         (
+            injectors_[i].x(),
+            injectorCoordinates_[i],
             injectorCells_[i],
             injectorTetFaces_[i],
-            injectorTetPts_[i],
+            injectorTetPts_[i]
-            injectors_[i].x()
         );
     }
 }
@@ -173,10 +177,11 @@ void Foam::ThermoLookupTableInjection<CloudType>::setPositionAndCell
     const label parcelI,
     const label nParcels,
     const scalar time,
-    vector& position,
+    barycentric& coordinates,
-    label& cellOwner,
+    label& celli,
     label& tetFacei,
-    label& tetPti
+    label& tetPti,
+    label& facei
 )
 {
     label injectorI = 0;
@@ -190,8 +195,8 @@ void Foam::ThermoLookupTableInjection<CloudType>::setPositionAndCell
         injectorI = parcelI*injectorCells_.size()/nParcels;
     }
 
-    position = injectors_[injectorI].x();
+    coordinates = injectorCoordinates_[injectorI];
-    cellOwner = injectorCells_[injectorI];
+    celli = injectorCells_[injectorI];
     tetFacei = injectorTetFaces_[injectorI];
     tetPti = injectorTetPts_[injectorI];
 }

src/lagrangian/parcel/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.H

@@ -86,6 +86,9 @@ class ThermoLookupTableInjection
         //- List of injectors
         thermoParcelInjectionDataIOList injectors_;
 
+        //- List of coordinates corresponding to injector positions
+        List<barycentric> injectorCoordinates_;
+
         //- List of cell labels corresponding to injector positions
         labelList injectorCells_;
 
@@ -156,10 +159,11 @@ public:
                 const label parcelI,
                 const label nParcels,
                 const scalar time,
-                vector& position,
+                barycentric& coordinates,
                 label& cellOwner,
                 label& tetFacei,
-                label& tetPti
+                label& tetPti,
+                label& facei
             );
 
             //- Set the parcel properties

src/lagrangian/solidParticle/solidParticle.C

@@ -38,13 +38,14 @@ namespace Foam
 bool Foam::solidParticle::move
 (
     solidParticleCloud& cloud,
-    trackingData& td,
+    trackingData& td
-    const scalar trackTime
 )
 {
     td.keepParticle = true;
     td.sendToProc = -1;
 
+    const scalar trackTime = td.mesh.time().deltaTValue();
+
     while (td.keepParticle && td.sendToProc == -1 && stepFraction() < 1)
     {
         if (debug)
@@ -54,7 +55,6 @@ bool Foam::solidParticle::move
                 << " stepFraction() = " << stepFraction() << endl;
         }
 
-
         const scalar sfrac = stepFraction();
 
         const scalar f = 1 - stepFraction();

src/lagrangian/solidParticle/solidParticle.H

@@ -132,18 +132,6 @@ public:
 
     // Constructors
 
-        //- Construct from components
-        inline solidParticle
-        (
-            const polyMesh& mesh,
-            const barycentric& coordinates,
-            const label celli,
-            const label tetFacei,
-            const label tetPti,
-            const scalar d,
-            const vector& U
-        );
-
         //- Construct from Istream
         solidParticle(Istream& is, bool readFields = true);
 
@@ -174,7 +162,7 @@ public:
         // Tracking
 
             //- Move
-            bool move(solidParticleCloud&, trackingData&, const scalar);
+            bool move(solidParticleCloud&, trackingData&);
 
 
         // Patch interactions

src/lagrangian/solidParticle/solidParticleCloud.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,7 +84,7 @@ void Foam::solidParticleCloud::move(const dimensionedVector& g)
     solidParticle::trackingData
         td(*this, rhoInterp, UInterp, nuInterp, g.value());
 
-    Cloud<solidParticle>::move(*this, td, mesh_.time().deltaTValue());
+    Cloud<solidParticle>::move(*this, td);
 }