diff --git a/src/OpenFOAM/db/functionObjects/timeControl/timeControlFunctionObject.C b/src/OpenFOAM/db/functionObjects/timeControl/timeControlFunctionObject.C
index 8e99aebe071bd8ba74829165bca0195bd61ae7dc..27f07f281c6c347cfae2daac49dbc036bcf7a299 100644
--- a/src/OpenFOAM/db/functionObjects/timeControl/timeControlFunctionObject.C
+++ b/src/OpenFOAM/db/functionObjects/timeControl/timeControlFunctionObject.C
@@ -51,6 +51,7 @@ void Foam::functionObjects::timeControl::readControls()
{
timeEnd_ = time_.userTimeToTime(timeEnd_);
}
+ deltaTCoeff_ = dict_.lookupOrDefault("deltaTCoeff", GREAT);
dict_.readIfPresent("nStepsToStartTimeChange", nStepsToStartTimeChange_);
}
@@ -64,6 +65,292 @@ bool Foam::functionObjects::timeControl::active() const
}
+Foam::scalar Foam::functionObjects::timeControl::calcExpansion
+(
+ const scalar startRatio,
+ const scalar y,
+ const label n
+)
+{
+ scalar ratio = startRatio;
+
+ // This function is used to calculate the 'expansion ratio' or
+ // 'time step growth factor'. Inputs:
+ // - ratio: wanted ratio
+ // - y: overall time required (divided by deltaT) to get from 1 to ratio
+ // - n: number of steps in which to get to wanted ratio
+
+ // Let a0 = first term in geometric progression (GP), an = last term,
+ // n = number of terms, ratio = ratio of GP, Sn = sum of series for
+ // first n terms.
+ // We have an=a0*ratio^(n-1) Sn = a0(ratio^n -1)/(ratio -1)
+ // -> Sn=an/ratio^(n-1)*(ratio^n -1)/(ratio -1), assume y=Sn/an
+ // -> y*ratio^(n-1)=(ratio^n -1)/(ratio -1)
+ // -> y*ratio^n - y*ratio^(n-1) -ratio^n +1 = 0
+
+ // Solve using Newton-Raphson iteration to determine the expansion
+ // ratio giving the desired overall timestep change.
+ // Note max iteration count to avoid hanging; function generally
+ // converges in 5 iterations or so.
+ for (label iter = 0; iter < 100; iter++)
+ {
+ // Dimensionless equation
+ scalar f = (y-1)*pow(ratio, n)+1-y*pow(ratio, n-1);
+ scalar dfdratio = (y-1)*n*pow(ratio, n-1)-y*(n-1)*pow(ratio, n-2);
+ scalar dratio = f/(dfdratio + SMALL);
+ ratio -= dratio;
+ // Exit if function is satisfied
+ if (mag(f) < 1e-6)
+ {
+ return ratio;
+ }
+ }
+
+ if (debug)
+ {
+ WarningInFunction
+ << "Did not converge to find new timestep growth factor given "
+ << "overall factor " << y << " and " << n << " steps to do it in."
+ << endl << " Returning current best guess " << ratio << endl;
+ }
+ return ratio;
+}
+
+
+void Foam::functionObjects::timeControl::calcDeltaTCoeff
+(
+ scalar& requiredDeltaTCoeff,
+ const scalar wantedDT,
+ const label nSteps,
+ const scalar presentTime,
+ const scalar timeToNextWrite,
+ const bool rampDirectionUp
+)
+{
+ const scalar writeInterval = writeControl_.interval();
+
+
+ // Set new deltaT approximated to last step value
+ // adjust to include geometric series sum of nSteps terms
+ scalar newDeltaT = deltaT0_;
+
+ if (seriesDTCoeff_ != GREAT)
+ {
+ newDeltaT *= seriesDTCoeff_;
+ }
+
+ // Recalculate new deltaTCoeff_ based on rounded steps
+ requiredDeltaTCoeff = Foam::exp(Foam::log(wantedDT/newDeltaT)/nSteps);
+
+ // Calculate total time required with given dT increment
+ // to reach specified dT value
+ scalar requiredTimeInterval = newDeltaT;
+
+ // For nSteps = 1 during we get coeff = 1.0. Adding SMALL in denominator
+ // might lead to SMALL requiredTimeInterval, and further unnecessary
+ // calculations
+ if (requiredDeltaTCoeff != 1.0)
+ {
+ requiredTimeInterval *=
+ (pow(requiredDeltaTCoeff, nSteps) - 1)
+ /(requiredDeltaTCoeff - 1);
+ }
+
+ // Nearest time which is multiple of wantedDT, after
+ // requiredTimeInterval
+ scalar timeToNextMultiple = -presentTime;
+
+ if (rampDirectionUp)
+ {
+ timeToNextMultiple +=
+ ceil
+ (
+ (presentTime + requiredTimeInterval)
+ /wantedDT
+ )
+ *wantedDT;
+ }
+ else
+ {
+ timeToNextMultiple +=
+ floor
+ (
+ (presentTime - requiredTimeInterval)
+ /wantedDT
+ )
+ *wantedDT;
+ }
+
+ // Check nearest time and decide parameters
+ if (timeToNextWrite <= timeToNextMultiple)
+ {
+ // As the ramping can't be achieved, use current, unadapted deltaT
+ // (from e.g. maxCo calculation in solver) and adjust when
+ // writeInterval is closer than this deltaT
+ newDeltaT = deltaT0_;
+ if (timeToNextWrite < newDeltaT)
+ {
+ newDeltaT = timeToNextWrite;
+ }
+
+ // Corner case when required time span is less than writeInterval
+ if (requiredTimeInterval > writeInterval)
+ {
+ WarningInFunction
+ << "With given ratio needed time span "
+ << requiredTimeInterval
+ << " exceeds available writeInterval "
+ << writeInterval << nl
+ << "Disabling all future time step ramping"
+ << endl;
+ deltaTCoeff_ = GREAT;
+ newDeltaT = wantedDT;
+ }
+ else
+ {
+ // Reset the series ratio, as this could change later
+ seriesDTCoeff_ = GREAT;
+ if (debug)
+ {
+ Info<< "Disabling ramping until time "
+ << presentTime+timeToNextWrite << endl;
+ }
+ }
+ requiredDeltaTCoeff = newDeltaT/deltaT0_;
+ }
+ else
+ {
+ // Find the minimum number of time steps (with constant deltaT
+ // variation) to get us to the writeInterval and multiples of wantedDT.
+ // Increases the number of nSteps to do the adjustment in until it
+ // has reached one that is possible. Returns the new expansionratio.
+
+ scalar y = timeToNextMultiple/wantedDT;
+
+ label requiredSteps = nSteps;
+
+ scalar ratioEstimate = deltaTCoeff_;
+ scalar ratioMax = deltaTCoeff_;
+
+ if (seriesDTCoeff_ != GREAT)
+ {
+ ratioEstimate = seriesDTCoeff_;
+ }
+
+ if (!rampDirectionUp)
+ {
+ ratioEstimate = 1/ratioEstimate;
+ requiredSteps *= -1;
+ }
+ // Provision for fall-back value if we can't satisfy requirements
+ bool searchConverged = false;
+ for (label iter = 0; iter < 100; iter++)
+ {
+ // Calculate the new expansion and the overall time step changes
+ const scalar newRatio = calcExpansion
+ (
+ ratioEstimate,
+ y,
+ requiredSteps
+ );
+ const scalar deltaTLoop = wantedDT/pow(newRatio, requiredSteps-1);
+ scalar firstDeltaRatio = deltaTLoop/deltaT0_;
+ if (debug)
+ {
+ // Avoid division by zero for ratio = 1.0
+ scalar Sn =
+ deltaTLoop
+ *(pow(newRatio, requiredSteps)-1)
+ /(newRatio-1+SMALL);
+
+ Info<< " nSteps " << requiredSteps
+ << " ratioEstimate " << ratioEstimate
+ << " a0 " << deltaTLoop
+ << " firstDeltaRatio " << firstDeltaRatio
+ << " Sn " << Sn << " Sn+Time " << Sn+presentTime
+ << " seriesRatio "<< newRatio << endl;
+ }
+
+ // If series ratio becomes unity check next deltaT multiple
+ // Do the same if first ratio is decreasing when ramping up!
+ if
+ (
+ (firstDeltaRatio < 1.0 && rampDirectionUp)
+ || (firstDeltaRatio > 1.0 && !rampDirectionUp)
+ || newRatio == 1.0
+ )
+ {
+ y += 1;
+ requiredSteps = nSteps;
+ if (debug)
+ {
+ Info<< "firstDeltaRatio " << firstDeltaRatio << " rampDir"
+ << rampDirectionUp << " newRatio " << newRatio
+ << " y " << y << " steps " << requiredSteps <<endl;
+ }
+ continue;
+ }
+
+ if (firstDeltaRatio > ratioMax)
+ {
+ requiredSteps++;
+ if (debug)
+ {
+ Info<< "First ratio " << firstDeltaRatio
+ << " exceeds threshold " << ratioMax << nl
+ << "Increasing required steps " << requiredSteps
+ << endl;
+ }
+ }
+ else if (newRatio > ratioMax)
+ {
+ y += 1;
+ requiredSteps = nSteps;
+ if (debug)
+ {
+ Info<< "Series ratio " << newRatio << "exceeds threshold "
+ << ratioMax << nl << "Consider next deltaT multiple "
+ << y*wantedDT + presentTime << endl;
+ }
+ }
+ else
+ {
+ seriesDTCoeff_ = newRatio;
+ // Reset future series expansion at last step
+ if (requiredSteps == 1)
+ {
+ seriesDTCoeff_ = GREAT;
+ }
+ if (debug)
+ {
+ Info<< "All conditions satisfied deltaT0_:" << deltaT0_
+ << " calculated deltaTCoeff:" << firstDeltaRatio
+ << " series ratio set to:" << seriesDTCoeff_ << endl;
+ }
+ searchConverged = true;
+ requiredDeltaTCoeff = firstDeltaRatio;
+ break;
+ }
+ }
+
+ if (!searchConverged)
+ {
+ if (debug)
+ {
+ // Not converged. Return fall-back value
+ WarningInFunction
+ << "Did not converge to find new timestep growth factor"
+ << " given overall factor " << y
+ << " and " << requiredSteps << " steps to do it in."
+ << nl << " Falling back to non-adjusted deltatT "
+ << deltaT0_ << endl;
+ }
+ requiredDeltaTCoeff = 1.0;
+ }
+ }
+}
+
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::functionObjects::timeControl::timeControl
@@ -85,7 +372,9 @@ Foam::functionObjects::timeControl::timeControl
executeControl_(t, dict, "execute"),
writeControl_(t, dict, "write"),
foPtr_(functionObject::New(name, t, dict_)),
- executeTimeIndex_(-1)
+ executeTimeIndex_(-1),
+ deltaT0_(0),
+ seriesDTCoeff_(GREAT)
{
readControls();
}
@@ -113,6 +402,9 @@ bool Foam::functionObjects::timeControl::entriesPresent(const dictionary& dict)
bool Foam::functionObjects::timeControl::execute()
{
+ // Store old deltaT for reference
+ deltaT0_ = time_.deltaTValue();
+
if (active() && (postProcess || executeControl_.execute()))
{
executeTimeIndex_ = time_.timeIndex();
@@ -152,44 +444,225 @@ bool Foam::functionObjects::timeControl::end()
}
-
bool Foam::functionObjects::timeControl::adjustTimeStep()
{
- if
- (
- active()
- && writeControl_.control()
- == Foam::timeControl::ocAdjustableRunTime
- )
+ if (active())
{
- const label writeTimeIndex = writeControl_.executionIndex();
- const scalar writeInterval = writeControl_.interval();
-
- scalar timeToNextWrite = max
+ if
(
- 0.0,
- (writeTimeIndex + 1)*writeInterval
- - (time_.value() - time_.startTime().value())
- );
+ writeControl_.control()
+ == Foam::timeControl::ocAdjustableRunTime
+ )
+ {
+ const scalar presentTime = time_.value();
- scalar deltaT = time_.deltaTValue();
+ //const scalar oldDeltaT = time_.deltaTValue();
- scalar nSteps = timeToNextWrite/deltaT - SMALL;
+ // Call underlying fo to do time step adjusting
+ foPtr_->adjustTimeStep();
+ const scalar wantedDT = time_.deltaTValue();
- // functionObjects modify deltaT within nStepsToStartTimeChange
- // NOTE: Potential problems arise if two function objects dump within
- // the same interval
- if (nSteps < nStepsToStartTimeChange_)
- {
- label nStepsToNextWrite = label(nSteps) + 1;
+ //Pout<< "adjustTimeStep by " << foPtr_->name()
+ // << " from " << oldDeltaT << " to " << wantedDT << endl;
- scalar newDeltaT = timeToNextWrite/nStepsToNextWrite;
+ const label writeTimeIndex = writeControl_.executionIndex();
+ const scalar writeInterval = writeControl_.interval();
+
+ // Get time to next writeInterval
+ scalar timeToNextWrite =
+ (writeTimeIndex+1)*writeInterval
+ - (presentTime - time_.startTime().value());
+ if (timeToNextWrite <= 0.0)
+ {
+ timeToNextWrite = writeInterval;
+ }
+
+ // The target DT (coming from fixed dT or maxCo etc) may not be
+ // an integral multiple of writeInterval. If so ignore any step
+ // calculation and give priority to writeInterval.
+ // Note looser tolerance on the relative intervalError - SMALL is
+ // too strict.
+ scalar intervalError =
+ remainder(writeInterval, wantedDT)/writeInterval;
+ if
+ (
+ mag(intervalError) > ROOTSMALL
+ || deltaTCoeff_ == GREAT
+ )
+ {
+ scalar deltaT = time_.deltaTValue();
+ scalar nSteps = timeToNextWrite/deltaT;
+
+ // For tiny deltaT the label can overflow!
+ if (nSteps < labelMax)
+ {
+ // nSteps can be < 1 so make sure at least 1
+
+ label nStepsToNextWrite = max(1, round(nSteps));
+ scalar newDeltaT = timeToNextWrite/nStepsToNextWrite;
+
+ // Backwards compatibility: clip deltaT to 0.2 .. 2
+ scalar clipThreshold = 2;
+ if (deltaTCoeff_ != GREAT)
+ {
+ clipThreshold = deltaTCoeff_;
+ }
+ // Adjust time step
+ if (newDeltaT >= deltaT)
+ {
+ deltaT = min(newDeltaT, clipThreshold*deltaT);
+ }
+ else
+ {
+ clipThreshold = 1/clipThreshold;
+ deltaT = max(newDeltaT, clipThreshold*deltaT);
+ }
+ const_cast<Time&>(time_).setDeltaT(deltaT, false);
+ }
+ }
+ else
+ {
+ // Set initial approximation to user defined ratio.
+ scalar requiredDeltaTCoeff = deltaTCoeff_;
+
+ // Use if we have series expansion ratio from last attempt.
+ // Starting from user defined coefficient, might get different
+ // steps and convergence time (could be recursive - worst case)
+ // This is more likely to happen when coefficient limit is high
+ if (seriesDTCoeff_ != GREAT)
+ {
+ requiredDeltaTCoeff = seriesDTCoeff_;
+ }
+ // Avoid divide by zero if we need ratio = 1
+ if (1/Foam::log(requiredDeltaTCoeff)> labelMax)
+ {
+ requiredDeltaTCoeff = deltaTCoeff_;
+ }
+
+ // Try and observe the deltaTCoeff and the writeInterval
+ // and the wanted deltaT. Below code calculates the
+ // minimum number of time steps in which to end up on
+ // writeInterval+n*deltaT with the minimum possible deltaTCoeff.
+ // This is non-trivial!
+
+ // Approx steps required from present dT to required dT
+ label nSteps = 0;
+ if (wantedDT < deltaT0_)
+ {
+ nSteps = label
+ (
+ floor
+ (
+ Foam::log(wantedDT/deltaT0_)
+ /Foam::log(requiredDeltaTCoeff + SMALL)
+ )
+ );
+ }
+ else
+ {
+ nSteps = label
+ (
+ ceil
+ (
+ Foam::log(wantedDT/deltaT0_)
+ /Foam::log(requiredDeltaTCoeff + SMALL)
+ )
+ );
+ }
+ // Negative steps -> ramp down needed,
+ // zero steps -> we are at writeInterval-multiple time!
+ if (nSteps < 1)
+ {
+ // Not enough steps to do gradual time step adjustment
+ // if earlier deltaT larger than wantedDT
+ if (wantedDT < deltaT0_)
+ {
+ requiredDeltaTCoeff = 1/requiredDeltaTCoeff;
+ calcDeltaTCoeff
+ (
+ requiredDeltaTCoeff,
+ wantedDT,
+ nSteps,
+ presentTime,
+ timeToNextWrite,
+ false
+ );
+ }
+ else
+ {
+ if (timeToNextWrite > wantedDT)
+ {
+ requiredDeltaTCoeff = wantedDT/deltaT0_;
+ }
+ else
+ {
+ requiredDeltaTCoeff = timeToNextWrite/deltaT0_;
+ }
+ }
+ // When allowed deltaTCoeff can't give required ramp
+ if (requiredDeltaTCoeff > deltaTCoeff_ && debug)
+ {
+ WarningInFunction
+ << "Required deltaTCoeff "<< requiredDeltaTCoeff
+ << " is larger than allowed value "<< deltaTCoeff_
+ << endl;
+ }
+ }
+ else
+ {
+ calcDeltaTCoeff
+ (
+ requiredDeltaTCoeff,
+ wantedDT,
+ nSteps,
+ presentTime,
+ timeToNextWrite,
+ true
+ );
+ }
+
+ // Calculate deltaT to be used based on ramp
+ scalar newDeltaT = deltaT0_*requiredDeltaTCoeff;
+
+ // Set time, deltaT adjustments for writeInterval purposes
+ // are already taken care. Hence disable auto-update
+ const_cast<Time&>( time_).setDeltaT(newDeltaT, false);
+
+ if (seriesDTCoeff_ < 1.0)
+ {
+ requiredDeltaTCoeff = 1/requiredDeltaTCoeff;
+ seriesDTCoeff_ = 1/seriesDTCoeff_;
+ }
+ }
+ }
+ else
+ {
+ foPtr_->adjustTimeStep();
+ const scalar wantedDT = time_.deltaTValue();
- // Adjust time step
- if (newDeltaT < deltaT)
+ if (deltaTCoeff_ != GREAT)
{
- deltaT = max(newDeltaT, 0.2*deltaT);
- const_cast<Time&>(time_).setDeltaT(deltaT, false);
+ // Clip time step change to deltaTCoeff
+
+ scalar requiredDeltaTCoeff =
+ (
+ max
+ (
+ 1.0/deltaTCoeff_,
+ min(deltaTCoeff_, wantedDT/deltaT0_)
+ )
+ );
+
+ scalar newDeltaT = deltaT0_*requiredDeltaTCoeff;
+
+ // Set time, deltaT adjustments for writeInterval purposes
+ // are already taken care. Hence disable auto-update
+ const_cast<Time&>( time_).setDeltaT(newDeltaT, false);
+ }
+ else
+ {
+ // Set the DT without any ramping
+ const_cast<Time&>( time_).setDeltaT(wantedDT, false);
}
}
}
@@ -198,10 +671,7 @@ bool Foam::functionObjects::timeControl::adjustTimeStep()
}
-bool Foam::functionObjects::timeControl::read
-(
- const dictionary& dict
-)
+bool Foam::functionObjects::timeControl::read(const dictionary& dict)
{
if (dict != dict_)
{
@@ -231,10 +701,7 @@ bool Foam::functionObjects::timeControl::filesModified() const
}
-void Foam::functionObjects::timeControl::updateMesh
-(
- const mapPolyMesh& mpm
-)
+void Foam::functionObjects::timeControl::updateMesh(const mapPolyMesh& mpm)
{
if (active())
{
@@ -243,10 +710,7 @@ void Foam::functionObjects::timeControl::updateMesh
}
-void Foam::functionObjects::timeControl::movePoints
-(
- const polyMesh& mesh
-)
+void Foam::functionObjects::timeControl::movePoints(const polyMesh& mesh)
{
if (active())
{
diff --git a/src/OpenFOAM/db/functionObjects/timeControl/timeControlFunctionObject.H b/src/OpenFOAM/db/functionObjects/timeControl/timeControlFunctionObject.H
index aac5343c424ededd4a4f2eec0a1f94dd033d7ec9..4b4bca78e1617f296dbaa7dbbbb6601a1b0eddf0 100644
--- a/src/OpenFOAM/db/functionObjects/timeControl/timeControlFunctionObject.H
+++ b/src/OpenFOAM/db/functionObjects/timeControl/timeControlFunctionObject.H
@@ -25,6 +25,13 @@ Class
Foam::functionObjects::timeControl
Description
+ Wrapper around functionObjects to add time control.
+
+ Adds
+ - timeStart, timeEnd activation
+ - deltaT modification to hit writeInterval (note: if adjustableRunTime
+ also in controlDict this has to be an integer multiple)
+ - smooth deltaT variation through optional deltaTCoeff parameter
Note
Since the timeIndex is used directly from Foam::Time, it is unaffected
@@ -85,6 +92,9 @@ class timeControl
//- De-activation time - defaults to VGREAT
scalar timeEnd_;
+ //- Max time step change
+ scalar deltaTCoeff_;
+
//- Number of steps before the dump-time during which deltaT
// may be changed (valid for adjustableRunTime)
label nStepsToStartTimeChange_;
@@ -103,6 +113,16 @@ class timeControl
label executeTimeIndex_;
+ // For time-step change limiting (deltaTCoeff supplied) only:
+
+ //- Store the old deltaT value
+ scalar deltaT0_;
+
+ //- Store the series expansion coefficient value
+ scalar seriesDTCoeff_;
+
+
+
// Private Member Functions
//- Read relevant dictionary entries
@@ -111,6 +131,26 @@ class timeControl
//- Returns true if within time bounds
bool active() const;
+ //- Return expansion ratio (deltaT change) that gives overall time
+ static scalar calcExpansion
+ (
+ const scalar startRatio,
+ const scalar y,
+ const label n
+ );
+
+ //- Calculate deltaT change such that the next write interval is
+ // obeyed. Updates requiredDeltaTCoeff
+ void calcDeltaTCoeff
+ (
+ scalar& requiredDeltaTCoeff,
+ const scalar wantedDT,
+ const label nSteps,
+ const scalar presentTime,
+ const scalar timeToNextWrite,
+ const bool rampDirectionUp
+ );
+
//- Disallow default bitwise copy construct
timeControl(const timeControl&);
diff --git a/src/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.C b/src/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.C
index 67a8a60a3b0def6e60e383ad3330fb403822db9f..c17fa5ecbdfe9ac9be95fa514742f7d08e0c903c 100644
--- a/src/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.C
+++ b/src/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.C
@@ -89,7 +89,7 @@ bool Foam::functionObjects::setTimeStepFunctionObject::adjustTimeStep()
index = time_.timeIndex();
// Set time, allow deltaT to be adjusted for writeInterval purposes
- const_cast<Time&>(time_).setDeltaT(newDeltaT, true);
+ const_cast<Time&>(time_).setDeltaT(newDeltaT, false);
}
return true;
diff --git a/src/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.H b/src/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.H
index 70925dc0a046abd07bc161dd89f4718ce2985237..f331073336a475397c488240fa3bf90bbd569062 100644
--- a/src/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.H
+++ b/src/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.H
@@ -43,6 +43,9 @@ Usage
{
type setTimeStep;
libs ("libutilityFunctionObjects.so");
+
+ deltaT table ((0 5e-4)(0.01 1e-3));
+
...
}
\endverbatim
@@ -52,6 +55,10 @@ Usage
Property | Description | Required | Default value
type | Type name: setTimeStep | yes |
enabled | On/off switch | no | yes
+ deltaT | Time step value | yes |
+ timeStart | Start time | no | 0
+ timeEnd | End time | no | GREAT
+ deltaTCoeff | Time step change limit | no | GREAT
\endtable
SourceFiles