ENH: streamline : added subcycling

src/postProcessing/functionObjects/field/streamLine/streamLine.C

@@ -53,16 +53,10 @@ void Foam::streamLine::track()
         initialParticles
     );
 
-    //Pout<< "Seeding particles." << endl;
-
     const sampledSet& seedPoints = sampledSetPtr_();
 
     forAll(seedPoints, i)
     {
-        //Pout<< "Seeded particle at " << seedPoints[i]
-        //    << " at cell:" << seedPoints.cells()[i]
-        //    << endl;
-
         particles.addParticle
         (
             new streamLineParticle
@@ -228,15 +222,14 @@ void Foam::streamLine::track()
         vvInterp,
         UIndex,         // index of U in vvInterp
         trackForward_,  // track in +u direction
+        nSubCycle_,     // step through cells in steps?
         allTracks_,
         allScalars_,
         allVectors_
     );
 
     // Track
-    //Pout<< "Tracking particles." << endl;
     particles.move(td);
-    //Pout<< "Finished tracking particles." << endl;
 }
 
 
@@ -296,6 +289,19 @@ void Foam::streamLine::read(const dictionary& dict)
         UName_ = dict.lookupOrDefault<word>("U", "U");
         dict.lookup("trackForward") >> trackForward_;
         dict.lookup("lifeTime") >> lifeTime_;
+        if (lifeTime_ < 1)
+        {
+            FatalErrorIn(":streamLine::read(const dictionary&)")
+                << "Illegal value " << lifeTime_ << " for lifeTime"
+                << exit(FatalError);
+        }
+
+        dict.lookup("nSubCycle") >> nSubCycle_;
+        if (nSubCycle_ < 1)
+        {
+            nSubCycle_ = 1;
+        }
+
         cloudName_ = dict.lookupOrDefault<word>("cloudName", "streamLine");
         dict.lookup("seedSampleSet") >> seedSet_;
 
@@ -323,7 +329,6 @@ void Foam::streamLine::execute()
 {
 //    const Time& runTime = const_cast<Time&>(obr_.time());
 //    Pout<< "**streamLine::execute : time:" << runTime.timeName() << endl;
-//    Pout<< "**streamLine::execute : time:" << runTime.timeIndex() << endl;
 //
 //    bool isOutputTime = false;
 //

src/postProcessing/functionObjects/field/streamLine/streamLine.H

@@ -92,6 +92,9 @@ class streamLine
         //- Maximum lifetime (= number of cells) of particle
         label lifeTime_;
 
+        //- Number of subcycling steps
+        label nSubCycle_;
+
         //- Optional specified name of particles
         word cloudName_;
 

src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C

@@ -36,6 +36,24 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+// Estimate dt to cross cell in a few steps
+Foam::scalar Foam::streamLineParticle::calcSubCycleDeltaT
+(
+    streamLineParticle::trackData& td,
+    const scalar dt,
+    const vector& U
+) const
+{
+    streamLineParticle testParticle(*this);
+    bool oldKeepParticle = td.keepParticle;
+    bool oldSwitchProcessor = td.switchProcessor;
+    scalar fraction = testParticle.trackToFace(position()+dt*U, td);
+    td.keepParticle = oldKeepParticle;
+    td.switchProcessor = oldSwitchProcessor;
+    // Adapt the dt to subdivide the trajectory into 4 substeps.
+    return dt *fraction/td.nSubCycle_;
+}
+
 Foam::vector Foam::streamLineParticle::interpolateFields
 (
     const streamLineParticle::trackData& td,
@@ -170,34 +188,56 @@ bool Foam::streamLineParticle::move(streamLineParticle::trackData& td)
     && lifeTime_ > 0
     )
     {
-        // TBD: implement subcycling so step through cells in more than
-        //      one step.
-
         // set the lagrangian time-step
         scalar dt = min(dtMax, tEnd);
 
-        // Store current position and sampled velocity.
-        --lifeTime_;
-        sampledPositions_.append(position());
-        vector U = interpolateFields(td, position(), cell());
-
-        if (!td.trackForward_)
+        // Cross cell in steps
+        for (label subIter = 0; subIter < td.nSubCycle_; subIter++)
         {
-            U = -U;
-        }
+            --lifeTime_;
 
-        dt *= trackToFace(position()+dt*U, td);
+            // Store current position and sampled velocity.
+            sampledPositions_.append(position());
+            vector U = interpolateFields(td, position(), cell());
 
-        tEnd -= dt;
-        stepFraction() = 1.0 - tEnd/deltaT;
+            if (!td.trackForward_)
+            {
+                U = -U;
+            }
 
-        if (tEnd <= ROOTVSMALL)
-        {
-            // Force removal
-            lifeTime_ = 0;
+
+            if (subIter == 0 && td.nSubCycle_ > 1)
+            {
+                // Adapt dt to cross cell in a few steps
+                dt = calcSubCycleDeltaT(td, dt, U);
+            }
+            else if (subIter == td.nSubCycle_-1)
+            {
+                // Do full step on last subcycle
+                dt = min(dtMax, tEnd);
+            }
+
+
+            scalar fraction = trackToFace(position()+dt*U, td);
+            dt *= fraction;
+
+            tEnd -= dt;
+            stepFraction() = 1.0 - tEnd/deltaT;
+
+            if (tEnd <= ROOTVSMALL)
+            {
+                // Force removal
+                lifeTime_ = 0;
+            }
+
+            if (!td.keepParticle || td.switchProcessor || lifeTime_ == 0)
+            {
+                break;
+            }
         }
     }
 
+
     if (!td.keepParticle || lifeTime_ == 0)
     {
         if (lifeTime_ == 0)

src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H

@@ -72,6 +72,7 @@ public:
         const PtrList<interpolationCellPoint<vector> >& vvInterp_;
         const label UIndex_;
         const bool trackForward_;
+        const label nSubCycle_;
 
         DynamicList<vectorList>& allPositions_;
         List<DynamicList<scalarList> >& allScalars_;
@@ -87,6 +88,7 @@ public:
                 const PtrList<interpolationCellPoint<vector> >& vvInterp,
                 const label UIndex,
                 const bool trackForward,
+                const label nSubCycle,
                 DynamicList<List<point> >& allPositions,
                 List<DynamicList<scalarList> >& allScalars,
                 List<DynamicList<vectorList> >& allVectors
@@ -97,6 +99,7 @@ public:
                 vvInterp_(vvInterp),
                 UIndex_(UIndex),
                 trackForward_(trackForward),
+                nSubCycle_(nSubCycle),
                 allPositions_(allPositions),
                 allScalars_(allScalars),
                 allVectors_(allVectors)
@@ -123,6 +126,15 @@ private:
 
     // Private Member Functions
 
+        //- Estimate dt to cross from current face to next one in nSubCycle
+        //  steps.
+        scalar calcSubCycleDeltaT
+        (
+            streamLineParticle::trackData& td,
+            const scalar dt,
+            const vector& U
+        ) const;
+
         //- Interpolate all quantities; return interpolated velocity.
         vector interpolateFields
         (

tutorials/incompressible/simpleFoam/motorBike/system/controlDict

@@ -83,7 +83,10 @@ functions
         fields (p U k);
 
         // Cells particles can travel before being removed
-        lifeTime        1000;
+        lifeTime        10000;
+
+        // Number of steps per cell (estimate). Set to 1 to disable subcycling.
+        nSubCycle 5;
 
         // Cloud name to use
         cloudName       particleTracks;

tutorials/incompressible/simpleFoam/pitzDaily/system/controlDict

@@ -85,7 +85,10 @@ functions
         fields (p k U);
 
         // Cells particles can travel before being removed
-        lifeTime        1000;
+        lifeTime        10000;
+
+        // Number of steps per cell (estimate). Set to 1 to disable subcycling.
+        nSubCycle 5;
 
         // Cloud name to use
         cloudName       particleTracks;