Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

This commit is contained in:
sergio
2017-06-14 14:00:39 -07:00
4 changed files with 555 additions and 44 deletions

View File

@ -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())
);
scalar deltaT = time_.deltaTValue();
scalar nSteps = timeToNextWrite/deltaT - SMALL;
// functionObjects modify deltaT within nStepsToStartTimeChange
// NOTE: Potential problems arise if two function objects dump within
// the same interval
if (nSteps < nStepsToStartTimeChange_)
writeControl_.control()
== Foam::timeControl::ocAdjustableRunTime
)
{
label nStepsToNextWrite = label(nSteps) + 1;
const scalar presentTime = time_.value();
scalar newDeltaT = timeToNextWrite/nStepsToNextWrite;
//const scalar oldDeltaT = time_.deltaTValue();
// Adjust time step
if (newDeltaT < deltaT)
// Call underlying fo to do time step adjusting
foPtr_->adjustTimeStep();
const scalar wantedDT = time_.deltaTValue();
//Pout<< "adjustTimeStep by " << foPtr_->name()
// << " from " << oldDeltaT << " to " << wantedDT << endl;
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)
{
deltaT = max(newDeltaT, 0.2*deltaT);
const_cast<Time&>(time_).setDeltaT(deltaT, false);
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();
if (deltaTCoeff_ != GREAT)
{
// 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())
{

View File

@ -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&);

View File

@ -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;

View File

@ -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