streamLines: Write out the streamline age

The streamLines function object now writes out "age" by default. This is
calculated as the total integration time from the starting point of the
streamline. This value can be negative if tracking is performed in the
reverse direction. Writing "age" can be deactivated by means of a
"writeAge no;" entry in the function object specification.

src/functionObjects/field/streamLine/streamLine.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -244,6 +244,8 @@ void Foam::functionObjects::streamLine::track()
     // Size to maximum expected sizes.
     allTracks_.clear();
     allTracks_.setCapacity(nSeeds);
+    allAges_.clear();
+    allAges_.setCapacity(nSeeds);
     allScalars_.setSize(vsInterp.size());
     forAll(allScalars_, i)
     {
@@ -272,6 +274,7 @@ void Foam::functionObjects::streamLine::track()
         trackLength_,   // fixed track length
 
         allTracks_,
+        allAges_,
         allScalars_,
         allVectors_
     );
@@ -342,6 +345,8 @@ bool Foam::functionObjects::streamLine::read(const dictionary& dict)
             << exit(FatalIOError);
     }
 
+    writeAge_ = dict.lookupOrDefault<Switch>("writeAge", true);
+
     // The trackForward entry is maintained here for backwards compatibility
     if (!dict.found("direction") && dict.found("trackForward"))
     {
@@ -498,9 +503,9 @@ bool Foam::functionObjects::streamLine::write()
         );
 
         // Distribute the scalars
-        forAll(allScalars_, scalarI)
+        forAll(allScalars_, scalari)
         {
-            allScalars_[scalarI].shrink();
+            allScalars_[scalari].shrink();
             mapDistributeBase::distribute
             (
                 Pstream::commsTypes::scheduled,
@@ -510,15 +515,15 @@ bool Foam::functionObjects::streamLine::write()
                 false,
                 distMap.constructMap(),
                 false,
-                allScalars_[scalarI],
+                allScalars_[scalari],
                 flipOp()
             );
-            allScalars_[scalarI].setCapacity(allScalars_[scalarI].size());
+            allScalars_[scalari].setCapacity(allScalars_[scalari].size());
         }
         // Distribute the vectors
-        forAll(allVectors_, vectorI)
+        forAll(allVectors_, vectori)
         {
-            allVectors_[vectorI].shrink();
+            allVectors_[vectori].shrink();
             mapDistributeBase::distribute
             (
                 Pstream::commsTypes::scheduled,
@@ -528,10 +533,10 @@ bool Foam::functionObjects::streamLine::write()
                 false,
                 distMap.constructMap(),
                 false,
-                allVectors_[vectorI],
+                allVectors_[vectori],
                 flipOp()
             );
-            allVectors_[vectorI].setCapacity(allVectors_[vectorI].size());
+            allVectors_[vectori].setCapacity(allVectors_[vectori].size());
         }
     }
 
@@ -587,21 +592,46 @@ bool Foam::functionObjects::streamLine::write()
 
         // Convert scalar values
 
-        if (allScalars_.size() > 0)
+        if (allScalars_.size() > 0 || writeAge_)
         {
-            List<List<scalarField>> scalarValues(allScalars_.size());
+            List<List<scalarField>> ageAndScalarValues
+            (
+                allScalars_.size() + writeAge_
+            );
+            wordList ageAndScalarNames(allScalars_.size() + writeAge_);
 
-            forAll(allScalars_, scalarI)
+            if (writeAge_)
             {
-                DynamicList<scalarList>& allTrackVals =
+                DynamicList<scalarList>& allTrackAges = allAges_;
-                    allScalars_[scalarI];
-                scalarValues[scalarI].setSize(allTrackVals.size());
 
+                ageAndScalarValues[0].setSize(allTrackAges.size());
+                forAll(allTrackAges, trackI)
+                {
+                    ageAndScalarValues[0][trackI].transfer
+                    (
+                        allTrackAges[trackI]
+                    );
+                }
+
+                ageAndScalarNames[0] = "age";
+            }
+
+            forAll(allScalars_, scalari)
+            {
+                const label ageAndScalari = scalari + writeAge_;
+
+                DynamicList<scalarList>& allTrackVals = allScalars_[scalari];
+
+                ageAndScalarValues[ageAndScalari].setSize(allTrackVals.size());
                 forAll(allTrackVals, trackI)
                 {
-                    scalarList& trackVals = allTrackVals[trackI];
+                    ageAndScalarValues[ageAndScalari][trackI].transfer
-                    scalarValues[scalarI][trackI].transfer(trackVals);
+                    (
+                        allTrackVals[trackI]
+                    );
                 }
+
+                ageAndScalarNames[ageAndScalari] = scalarNames_[scalari];
             }
 
             fileName vtkFile
@@ -610,7 +640,7 @@ bool Foam::functionObjects::streamLine::write()
               / scalarFormatterPtr_().getFileName
                 (
                     tracks[0],
-                    scalarNames_
+                    ageAndScalarNames
                 )
             );
 
@@ -620,8 +650,8 @@ bool Foam::functionObjects::streamLine::write()
             (
                 true,           // writeTracks
                 tracks,
-                scalarNames_,
+                ageAndScalarNames,
-                scalarValues,
+                ageAndScalarValues,
                 OFstream(vtkFile)()
             );
         }
@@ -632,16 +662,17 @@ bool Foam::functionObjects::streamLine::write()
         {
             List<List<vectorField>> vectorValues(allVectors_.size());
 
-            forAll(allVectors_, vectorI)
+            forAll(allVectors_, vectori)
             {
-                DynamicList<vectorList>& allTrackVals =
+                DynamicList<vectorList>& allTrackVals = allVectors_[vectori];
-                    allVectors_[vectorI];
-                vectorValues[vectorI].setSize(allTrackVals.size());
 
+                vectorValues[vectori].setSize(allTrackVals.size());
                 forAll(allTrackVals, trackI)
                 {
-                    vectorList& trackVals = allTrackVals[trackI];
+                    vectorValues[vectori][trackI].transfer
-                    vectorValues[vectorI][trackI].transfer(trackVals);
+                    (
+                        allTrackVals[trackI]
+                    );
                 }
             }
 
@@ -655,6 +686,8 @@ bool Foam::functionObjects::streamLine::write()
                 )
             );
 
+            Info<< "    Writing data to " << vtkFile.path() << endl;
+
             vectorFormatterPtr_().write
             (
                 true,           // writeTracks

src/functionObjects/field/streamLine/streamLine.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,7 +39,7 @@ Description
 
         setFormat       vtk;
         U               U;
-        trackForward    yes;
+        direction       both;
 
         fields
         (
@@ -69,7 +69,9 @@ Usage
         type         | Type name: streamLine   | yes         |
         setFormat    | Output data type        | yes         |
         U            | Tracking velocity field name | no     | U
+        direction    | Direction in which to track | yes     |
         fields       | Fields to sample        | yes         |
+        writeTime    | Write the flow time along the streamline | no | no
         lifetime     | Maximum number of particle tracking steps | yes |
         trackLength  | Tracking segment length | no          |
         nSubCycle    | Number of tracking steps per cell | no|
@@ -187,6 +189,9 @@ private:
         //- Names of vector fields
         wordList vectorNames_;
 
+        //- Write the streamline ages
+        Switch writeAge_;
+
 
         // Demand driven
 
@@ -211,6 +216,9 @@ private:
             //- All tracks. Per particle the points it passed through
             DynamicList<List<point>> allTracks_;
 
+            //- All ages. Per particle the age when it passed through the points
+            DynamicList<List<scalar>> allAges_;
+
             //- Per scalarField, per particle, the sampled value.
             List<DynamicList<scalarList>> allScalars_;
 

src/functionObjects/field/streamLine/streamLineParticle.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 #include "streamLineParticle.H"
 #include "streamLineParticleCloud.H"
 #include "vectorFieldIOField.H"
+#include "scalarFieldIOField.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -88,7 +89,8 @@ Foam::streamLineParticle::streamLineParticle
 )
 :
     particle(mesh, position, celli),
-    lifeTime_(lifeTime)
+    lifeTime_(lifeTime),
+    age_(0)
 {}
 
 
@@ -106,8 +108,8 @@ Foam::streamLineParticle::streamLineParticle
         List<scalarList> sampledScalars;
         List<vectorList> sampledVectors;
 
-        is  >> lifeTime_ >> sampledPositions_ >> sampledScalars
+        is  >> lifeTime_ >> age_ >> sampledPositions_ >> sampledTimes_
-            >> sampledVectors;
+            >> sampledScalars >> sampledVectors;
 
         sampledScalars_.setSize(sampledScalars.size());
         forAll(sampledScalars, i)
@@ -137,8 +139,11 @@ Foam::streamLineParticle::streamLineParticle
 :
     particle(p),
     lifeTime_(p.lifeTime_),
+    age_(p.age_),
     sampledPositions_(p.sampledPositions_),
-    sampledScalars_(p.sampledScalars_)
+    sampledTimes_(p.sampledTimes_),
+    sampledScalars_(p.sampledScalars_),
+    sampledVectors_(p.sampledVectors_)
 {}
 
 
@@ -169,6 +174,7 @@ bool Foam::streamLineParticle::move
 
             // Store current position and sampled velocity.
             sampledPositions_.append(position());
+            sampledTimes_.append(age_);
             vector U = interpolateFields(td, position(), cell(), face());
 
             if (!td.trackForward_)
@@ -205,7 +211,10 @@ bool Foam::streamLineParticle::move
                 dt = maxDt;
             }
 
-            trackToAndHitFace(dt*U, 0, cloud, td);
+            age_ +=
+                (td.trackForward_ ? +1 : -1)
+               *dt
+               *(1 - trackToAndHitFace(dt*U, 0, cloud, td));
 
             if
             (
@@ -237,6 +246,7 @@ bool Foam::streamLineParticle::move
         {
             // Normal exit. Store last position and fields
             sampledPositions_.append(position());
+            sampledTimes_.append(age_);
             interpolateFields(td, position(), cell(), face());
 
             if (debug)
@@ -248,10 +258,16 @@ bool Foam::streamLineParticle::move
         }
 
         // Transfer particle data into trackingData.
-        td.allPositions_.append(vectorList());
+        {
-        vectorList& top = td.allPositions_.last();
+            td.allPositions_.append(vectorList());
-        top.transfer(sampledPositions_);
+            vectorList& top = td.allPositions_.last();
-
+            top.transfer(sampledPositions_);
+        }
+        {
+            td.allTimes_.append(scalarList());
+            scalarList& top = td.allTimes_.last();
+            top.transfer(sampledTimes_);
+        }
         forAll(sampledScalars_, i)
         {
             td.allScalars_[i].append(scalarList());
@@ -406,11 +422,19 @@ void Foam::streamLineParticle::readFields(Cloud<streamLineParticle>& c)
     );
     c.checkFieldIOobject(c, sampledPositions);
 
+    scalarFieldIOField sampledTimes
+    (
+        c.fieldIOobject("sampledTimes", IOobject::MUST_READ),
+        valid
+    );
+    c.checkFieldIOobject(c, sampledTimes);
+
     label i = 0;
     forAllIter(Cloud<streamLineParticle>, c, iter)
     {
         iter().lifeTime_ = lifeTime[i];
         iter().sampledPositions_.transfer(sampledPositions[i]);
+        iter().sampledTimes_.transfer(sampledTimes[i]);
         i++;
     }
 }
@@ -432,17 +456,24 @@ void Foam::streamLineParticle::writeFields(const Cloud<streamLineParticle>& c)
         c.fieldIOobject("sampledPositions", IOobject::NO_READ),
         np
     );
+    scalarFieldIOField sampledTimes
+    (
+        c.fieldIOobject("sampledTimes", IOobject::NO_READ),
+        np
+    );
 
     label i = 0;
     forAllConstIter(Cloud<streamLineParticle>, c, iter)
     {
         lifeTime[i] = iter().lifeTime_;
         sampledPositions[i] = iter().sampledPositions_;
+        sampledTimes[i] = iter().sampledTimes_;
         i++;
     }
 
     lifeTime.write(np > 0);
     sampledPositions.write(np > 0);
+    sampledTimes.write(np > 0);
 }
 
 
@@ -452,7 +483,9 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const streamLineParticle& p)
 {
     os  << static_cast<const particle&>(p)
         << token::SPACE << p.lifeTime_
+        << token::SPACE << p.age_
         << token::SPACE << p.sampledPositions_
+        << token::SPACE << p.sampledTimes_
         << token::SPACE << p.sampledScalars_
         << token::SPACE << p.sampledVectors_;
 

src/functionObjects/field/streamLine/streamLineParticle.H

@@ -86,6 +86,8 @@ public:
 
             DynamicList<vectorList>& allPositions_;
 
+            DynamicList<scalarList>& allTimes_;
+
             List<DynamicList<scalarList>>& allScalars_;
 
             List<DynamicList<vectorList>>& allVectors_;
@@ -104,6 +106,7 @@ public:
                 const label nSubCycle,
                 const scalar trackLength,
                 DynamicList<List<point>>& allPositions,
+                DynamicList<List<scalar>>& allTimes,
                 List<DynamicList<scalarList>>& allScalars,
                 List<DynamicList<vectorList>>& allVectors
             )
@@ -116,6 +119,7 @@ public:
                 nSubCycle_(nSubCycle),
                 trackLength_(trackLength),
                 allPositions_(allPositions),
+                allTimes_(allTimes),
                 allScalars_(allScalars),
                 allVectors_(allVectors)
             {}
@@ -129,9 +133,15 @@ private:
         //- Lifetime of particle. Particle dies when reaches 0.
         label lifeTime_;
 
+        //- Age of the particle
+        scalar age_;
+
         //- Sampled positions
         DynamicList<point> sampledPositions_;
 
+        //- Sampled times
+        DynamicList<scalar> sampledTimes_;
+
         //- Sampled scalars
         List<DynamicList<scalar>> sampledScalars_;