Skip to content
Snippets Groups Projects
Commit 7e6faf04 authored by mattijs's avatar mattijs
Browse files

Squashed commit of the following:

commit 510b9353f8cb932a87f1588b17a4dea230c29d3c
Author: mattijs <mattijs>
Date:   Wed Jun 14 16:36:35 2017 +0100

    ENH: timeControl: propagate time-change logic from Time.C

commit 6dc57a8f1e0e7605ea819deb6f11dd4b7874ff30
Author: mattijs <mattijs>
Date:   Thu Jun 1 11:28:56 2017 +0100

    ENH: timeControl: cleanup; avoid division by zero

commit 5ac4bc2dd8c1f4676eef9d7a03215caba23a2e19
Author: mattijs <mattijs>
Date:   Wed May 31 12:04:52 2017 +0100

    ENH: timeControl: allow ramping down as well as up.

commit a6b2db9e791f29258f04f3a9cbd6354aa468977d
Author: mattijs <mattijs>
Date:   Wed May 17 15:29:22 2017 +0100

    ENH: timeControl: limit any timestep change if deltaTCoeff enabled.

commit 004115ee03a6637ae0d23cce303a30d1b3af046f
Author: mattijs <mattijs>
Date:   Wed May 17 11:40:26 2017 +0100

    ENH: setTimeStep: have timeStart, timeEnd controls on time step adjustment.

    Also added smoothly varying time step change (through optional deltaTCoeff)
parent 8a8e5e88
No related merge requests found
...@@ -51,6 +51,7 @@ void Foam::functionObjects::timeControl::readControls() ...@@ -51,6 +51,7 @@ void Foam::functionObjects::timeControl::readControls()
{ {
timeEnd_ = time_.userTimeToTime(timeEnd_); timeEnd_ = time_.userTimeToTime(timeEnd_);
} }
deltaTCoeff_ = dict_.lookupOrDefault("deltaTCoeff", GREAT);
dict_.readIfPresent("nStepsToStartTimeChange", nStepsToStartTimeChange_); dict_.readIfPresent("nStepsToStartTimeChange", nStepsToStartTimeChange_);
} }
...@@ -64,6 +65,292 @@ bool Foam::functionObjects::timeControl::active() const ...@@ -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 * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::functionObjects::timeControl::timeControl Foam::functionObjects::timeControl::timeControl
...@@ -85,7 +372,9 @@ Foam::functionObjects::timeControl::timeControl ...@@ -85,7 +372,9 @@ Foam::functionObjects::timeControl::timeControl
executeControl_(t, dict, "execute"), executeControl_(t, dict, "execute"),
writeControl_(t, dict, "write"), writeControl_(t, dict, "write"),
foPtr_(functionObject::New(name, t, dict_)), foPtr_(functionObject::New(name, t, dict_)),
executeTimeIndex_(-1) executeTimeIndex_(-1),
deltaT0_(0),
seriesDTCoeff_(GREAT)
{ {
readControls(); readControls();
} }
...@@ -113,6 +402,9 @@ bool Foam::functionObjects::timeControl::entriesPresent(const dictionary& dict) ...@@ -113,6 +402,9 @@ bool Foam::functionObjects::timeControl::entriesPresent(const dictionary& dict)
bool Foam::functionObjects::timeControl::execute() bool Foam::functionObjects::timeControl::execute()
{ {
// Store old deltaT for reference
deltaT0_ = time_.deltaTValue();
if (active() && (postProcess || executeControl_.execute())) if (active() && (postProcess || executeControl_.execute()))
{ {
executeTimeIndex_ = time_.timeIndex(); executeTimeIndex_ = time_.timeIndex();
...@@ -152,44 +444,225 @@ bool Foam::functionObjects::timeControl::end() ...@@ -152,44 +444,225 @@ bool Foam::functionObjects::timeControl::end()
} }
bool Foam::functionObjects::timeControl::adjustTimeStep() bool Foam::functionObjects::timeControl::adjustTimeStep()
{ {
if if (active())
(
active()
&& writeControl_.control()
== Foam::timeControl::ocAdjustableRunTime
)
{ {
const label writeTimeIndex = writeControl_.executionIndex(); if
const scalar writeInterval = writeControl_.interval();
scalar timeToNextWrite = max
( (
0.0, writeControl_.control()
(writeTimeIndex + 1)*writeInterval == Foam::timeControl::ocAdjustableRunTime
- (time_.value() - time_.startTime().value()) )
); {
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 //Pout<< "adjustTimeStep by " << foPtr_->name()
// NOTE: Potential problems arise if two function objects dump within // << " from " << oldDeltaT << " to " << wantedDT << endl;
// the same interval
if (nSteps < nStepsToStartTimeChange_)
{
label nStepsToNextWrite = label(nSteps) + 1;
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 (deltaTCoeff_ != GREAT)
if (newDeltaT < deltaT)
{ {
deltaT = max(newDeltaT, 0.2*deltaT); // Clip time step change to deltaTCoeff
const_cast<Time&>(time_).setDeltaT(deltaT, false);
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() ...@@ -198,10 +671,7 @@ bool Foam::functionObjects::timeControl::adjustTimeStep()
} }
bool Foam::functionObjects::timeControl::read bool Foam::functionObjects::timeControl::read(const dictionary& dict)
(
const dictionary& dict
)
{ {
if (dict != dict_) if (dict != dict_)
{ {
...@@ -231,10 +701,7 @@ bool Foam::functionObjects::timeControl::filesModified() const ...@@ -231,10 +701,7 @@ bool Foam::functionObjects::timeControl::filesModified() const
} }
void Foam::functionObjects::timeControl::updateMesh void Foam::functionObjects::timeControl::updateMesh(const mapPolyMesh& mpm)
(
const mapPolyMesh& mpm
)
{ {
if (active()) if (active())
{ {
...@@ -243,10 +710,7 @@ void Foam::functionObjects::timeControl::updateMesh ...@@ -243,10 +710,7 @@ void Foam::functionObjects::timeControl::updateMesh
} }
void Foam::functionObjects::timeControl::movePoints void Foam::functionObjects::timeControl::movePoints(const polyMesh& mesh)
(
const polyMesh& mesh
)
{ {
if (active()) if (active())
{ {
......
...@@ -25,6 +25,13 @@ Class ...@@ -25,6 +25,13 @@ Class
Foam::functionObjects::timeControl Foam::functionObjects::timeControl
Description 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 Note
Since the timeIndex is used directly from Foam::Time, it is unaffected Since the timeIndex is used directly from Foam::Time, it is unaffected
...@@ -85,6 +92,9 @@ class timeControl ...@@ -85,6 +92,9 @@ class timeControl
//- De-activation time - defaults to VGREAT //- De-activation time - defaults to VGREAT
scalar timeEnd_; scalar timeEnd_;
//- Max time step change
scalar deltaTCoeff_;
//- Number of steps before the dump-time during which deltaT //- Number of steps before the dump-time during which deltaT
// may be changed (valid for adjustableRunTime) // may be changed (valid for adjustableRunTime)
label nStepsToStartTimeChange_; label nStepsToStartTimeChange_;
...@@ -103,6 +113,16 @@ class timeControl ...@@ -103,6 +113,16 @@ class timeControl
label executeTimeIndex_; 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 // Private Member Functions
//- Read relevant dictionary entries //- Read relevant dictionary entries
...@@ -111,6 +131,26 @@ class timeControl ...@@ -111,6 +131,26 @@ class timeControl
//- Returns true if within time bounds //- Returns true if within time bounds
bool active() const; 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 //- Disallow default bitwise copy construct
timeControl(const timeControl&); timeControl(const timeControl&);
......
...@@ -89,7 +89,7 @@ bool Foam::functionObjects::setTimeStepFunctionObject::adjustTimeStep() ...@@ -89,7 +89,7 @@ bool Foam::functionObjects::setTimeStepFunctionObject::adjustTimeStep()
index = time_.timeIndex(); index = time_.timeIndex();
// Set time, allow deltaT to be adjusted for writeInterval purposes // 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; return true;
......
...@@ -43,6 +43,9 @@ Usage ...@@ -43,6 +43,9 @@ Usage
{ {
type setTimeStep; type setTimeStep;
libs ("libutilityFunctionObjects.so"); libs ("libutilityFunctionObjects.so");
deltaT table ((0 5e-4)(0.01 1e-3));
... ...
} }
\endverbatim \endverbatim
...@@ -52,6 +55,10 @@ Usage ...@@ -52,6 +55,10 @@ Usage
Property | Description | Required | Default value Property | Description | Required | Default value
type | Type name: setTimeStep | yes | type | Type name: setTimeStep | yes |
enabled | On/off switch | no | 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 \endtable
SourceFiles SourceFiles
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment