diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLine.C b/src/postProcessing/functionObjects/field/streamLine/streamLine.C
index 8336fd8f9a9144b8f3465e144421547b81b63a1a..1a7c86c1a49153b45f15cbcde112a2fae655281e 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLine.C
+++ b/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;
//
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLine.H b/src/postProcessing/functionObjects/field/streamLine/streamLine.H
index 55edc7efe091e92e0681fe66f51cba081af01d2e..ed39204d2230ece5b940b9881cb714940b514d4b 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLine.H
+++ b/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_;
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
index 56b05c88caf9b1f39c3ebf21f79258629363ce7f..7230a21345accc9cacd715176a2fe9cc8d940b72 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
+++ b/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)
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
index d556e4791decf15ab076f885143fe1adb7b3f1d1..bd16038d83cbf29bec5f763bf7154b005563912d 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
+++ b/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
(
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
index c81f8fd5885ea5463ed2cb12e3032fe154880fb1..eed783b8a95f21d2c362180b881e9f8ef612f72f 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
+++ b/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;
diff --git a/tutorials/incompressible/simpleFoam/pitzDaily/system/controlDict b/tutorials/incompressible/simpleFoam/pitzDaily/system/controlDict
index 86cf5ca88a82ecc234f0b72118b4ba84f5590c39..e608ae6e77483c7f15b6135d93f7648c7d6b170d 100644
--- a/tutorials/incompressible/simpleFoam/pitzDaily/system/controlDict
+++ b/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;