9 #ifndef Tempus_TimeStepControl_impl_hpp
10 #define Tempus_TimeStepControl_impl_hpp
13 #include "Teuchos_ScalarTraits.hpp"
14 #include "Teuchos_StandardParameterEntryValidators.hpp"
15 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
16 #include "Teuchos_TimeMonitor.hpp"
25 #include "Thyra_VectorStdOps.hpp"
29 template<
class Scalar>
31 Teuchos::RCP<Teuchos::ParameterList> pList)
32 : outputAdjustedDt_(false), dtAfterOutput_(0.0)
37 template<
class Scalar>
39 : tscPL_ (tsc_.tscPL_ ),
40 outputIndices_ (tsc_.outputIndices_ ),
41 outputTimes_ (tsc_.outputTimes_ ),
42 outputAdjustedDt_ (tsc_.outputAdjustedDt_ ),
43 dtAfterOutput_ (tsc_.dtAfterOutput_ ),
44 stepControlStrategy_(tsc_.stepControlStrategy_ )
48 template<
class Scalar>
55 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::TimeStepControl::getNextTimeStep()");
57 RCP<Teuchos::FancyOStream> out = this->getOStream();
58 Teuchos::OSTab ostab(out,0,
"getNextTimeStep");
60 Teuchos::as<int>(Teuchos::VERB_NONE);
62 auto changeDT = [] (
int istep, Scalar dt_old, Scalar dt_new,
65 std::stringstream message;
66 message << std::scientific
67 <<std::setw(6)<<std::setprecision(3)<<istep
68 <<
" * (dt = "<<std::setw(9)<<std::setprecision(3)<<dt_old
69 <<
", new = "<<std::setw(9)<<std::setprecision(3)<<dt_new
70 <<
") " << reason << std::endl;
74 RCP<SolutionState<Scalar> > workingState=
solutionHistory->getWorkingState();
75 RCP<SolutionStateMetaData<Scalar> > metaData = workingState->getMetaData();
77 const int iStep = metaData->getIStep();
78 int order = metaData->getOrder();
79 Scalar dt = metaData->getDt();
82 RCP<StepperState<Scalar> > stepperState = workingState->getStepperState();
84 if (getStepType() ==
"Variable") {
86 if (outputAdjustedDt_ ==
true) {
87 if (printChanges) *out << changeDT(iStep, dt, dtAfterOutput_,
88 "Reset dt after output.");
90 outputAdjustedDt_ =
false;
95 if (printChanges) *out << changeDT(iStep, dt, getInitTimeStep(),
96 "Reset dt to initial dt.");
97 dt = getInitTimeStep();
100 if (dt < getMinTimeStep()) {
101 if (printChanges) *out << changeDT(iStep, dt, getMinTimeStep(),
102 "Reset dt to minimum dt.");
103 dt = getMinTimeStep();
115 order = metaData->getOrder();
116 dt = metaData->getDt();
118 if (getStepType() ==
"Variable") {
119 if (dt < getMinTimeStep()) {
120 if (printChanges) *out << changeDT(iStep, dt, getMinTimeStep(),
121 "dt is too small. Resetting to minimum dt.");
122 dt = getMinTimeStep();
124 if (dt > getMaxTimeStep()) {
125 if (printChanges) *out << changeDT(iStep, dt, getMaxTimeStep(),
126 "dt is too large. Resetting to maximum dt.");
127 dt = getMaxTimeStep();
133 std::vector<int>::const_iterator it =
134 std::find(outputIndices_.begin(), outputIndices_.end(), iStep);
135 if (it != outputIndices_.end()) output =
true;
137 const int iInterval = tscPL_->get<
int>(
"Output Index Interval");
138 if ( (iStep - getInitIndex()) % iInterval == 0) output =
true;
142 Scalar reltol = 1.0e-6;
143 const Scalar endTime = lastTime+dt+getMinTimeStep();
144 bool checkOutput =
false;
145 Scalar oTime = getInitTime();
146 for (
size_t i=0; i < outputTimes_.size(); ++i) {
147 oTime = outputTimes_[i];
148 if (lastTime < oTime && oTime <= endTime) {
153 const Scalar tInterval = tscPL_->get<
double>(
"Output Time Interval");
154 Scalar oTime2 = ceil((lastTime-getInitTime())/tInterval)*tInterval
156 if (lastTime < oTime2 && oTime2 <= endTime) {
157 if (checkOutput ==
true) {
158 if (oTime2 < oTime) oTime = oTime2;
165 if (checkOutput ==
true) {
166 const bool outputExactly =
167 tscPL_->get<
bool>(
"Output Exactly On Output Times");
168 if (getStepType() ==
"Variable" && outputExactly ==
true) {
170 if (std::abs((lastTime+dt-oTime)/(lastTime+dt)) < reltol) {
172 if (printChanges) *out << changeDT(iStep, dt, oTime - lastTime,
173 "Adjusting dt for numerical roundoff to hit the next output time.");
176 outputAdjustedDt_ =
true;
178 dt = oTime - lastTime;
179 }
else if (lastTime*(1.0+reltol) < oTime &&
180 oTime < (lastTime+dt-getMinTimeStep())*(1.0+reltol)) {
182 if (printChanges) *out << changeDT(iStep, dt, oTime - lastTime,
183 "Adjusting dt to hit the next output time.");
187 outputAdjustedDt_ =
true;
189 dt = oTime - lastTime;
191 if (printChanges) *out << changeDT(iStep, dt, (oTime - lastTime)/2.0,
192 "The next output time is within the minimum dt of the next time. "
193 "Adjusting dt to take two steps.");
197 dt = (oTime - lastTime)/2.0;
209 if ((lastTime + dt > getFinalTime() ) ||
210 (std::abs((lastTime+dt-getFinalTime())/(lastTime+dt)) < reltol)) {
211 if (printChanges) *out << changeDT(iStep, dt, getFinalTime() - lastTime,
212 "Adjusting dt to hit the final time.");
213 dt = getFinalTime() - lastTime;
217 TEUCHOS_TEST_FOR_EXCEPTION( dt <= Scalar(0.0), std::out_of_range,
218 "Error - Time step is not positive. dt = " << dt <<
"\n");
221 TEUCHOS_TEST_FOR_EXCEPTION(
222 (lastTime + dt < getInitTime()), std::out_of_range,
223 "Error - Time step does not move time INTO time range.\n"
224 " [timeMin, timeMax] = [" << getInitTime() <<
", "
225 << getFinalTime() <<
"]\n"
226 " T + dt = " << lastTime <<
" + "<< dt <<
" = " << lastTime + dt <<
"\n");
228 TEUCHOS_TEST_FOR_EXCEPTION(
229 (lastTime + dt > getFinalTime()), std::out_of_range,
230 "Error - Time step move time OUT OF time range.\n"
231 " [timeMin, timeMax] = [" << getInitTime() <<
", "
232 << getFinalTime() <<
"]\n"
233 " T + dt = " << lastTime <<
" + "<< dt <<
" = " << lastTime + dt <<
"\n");
235 metaData->setOrder(order);
237 metaData->setTime(lastTime + dt);
238 metaData->setOutput(output);
245 template<
class Scalar>
247 const Scalar relTol = 1.0e-14;
248 bool tir = (getInitTime()*(1.0-relTol) <= time and
249 time < getFinalTime()*(1.0-relTol));
254 template<
class Scalar>
256 bool iir = (getInitIndex() <= iStep and iStep < getFinalIndex());
261 template<
class Scalar>
264 if (numTimeSteps >= 0) {
265 tscPL_->set<
int> (
"Number of Time Steps", numTimeSteps);
266 const int initIndex = getInitIndex();
267 tscPL_->set<
int> (
"Final Time Index", initIndex + numTimeSteps);
268 const double initTime = tscPL_->get<
double>(
"Initial Time");
269 const double finalTime = tscPL_->get<
double>(
"Final Time");
270 double initTimeStep = (finalTime - initTime)/numTimeSteps;
271 if (numTimeSteps == 0) initTimeStep = Scalar(0.0);
272 tscPL_->set<
double> (
"Initial Time Step", initTimeStep);
273 tscPL_->set<
double> (
"Minimum Time Step", initTimeStep);
274 tscPL_->set<
double> (
"Maximum Time Step", initTimeStep);
275 tscPL_->set<std::string>(
"Integrator Step Type",
"Constant");
277 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
278 Teuchos::OSTab ostab(out,1,
"setNumTimeSteps");
279 *out <<
"Warning - Found 'Number of Time Steps' = " << getNumTimeSteps()
280 <<
" Set the following parameters: \n"
281 <<
" 'Final Time Index' = " << getFinalIndex() <<
"\n"
282 <<
" 'Initial Time Step' = " << getInitTimeStep() <<
"\n"
283 <<
" 'Integrator Step Type' = " << getStepType() << std::endl;
288 template<
class Scalar>
291 std::string name =
"Tempus::TimeStepControl";
296 template<
class Scalar>
298 Teuchos::FancyOStream &out,
299 const Teuchos::EVerbosityLevel verbLevel)
const
301 if (verbLevel == Teuchos::VERB_EXTREME) {
302 out << description() <<
"::describe:" << std::endl
303 <<
"pList = " << tscPL_ << std::endl;
308 template <
class Scalar>
310 Teuchos::RCP<Teuchos::ParameterList>
const& pList)
312 if (pList == Teuchos::null) {
314 if (tscPL_ == Teuchos::null) {
315 tscPL_ = Teuchos::parameterList(
"TimeStepControl");
316 *tscPL_ = *(this->getValidParameters());
321 tscPL_->validateParametersAndSetDefaults(*this->getValidParameters(), 0);
324 if (getStepType() ==
"Constant") {
325 const double initTimeStep = tscPL_->get<
double>(
"Initial Time Step");
326 tscPL_->set<
double> (
"Minimum Time Step", initTimeStep);
327 tscPL_->set<
double> (
"Maximum Time Step", initTimeStep);
329 setNumTimeSteps(getNumTimeSteps());
332 setTimeStepControlStrategy();
334 TEUCHOS_TEST_FOR_EXCEPTION(
335 (getInitTime() > getFinalTime() ), std::logic_error,
336 "Error - Inconsistent time range.\n"
337 " (timeMin = "<<getInitTime()<<
") > (timeMax = "<<getFinalTime()<<
")\n");
339 TEUCHOS_TEST_FOR_EXCEPTION(
340 (getMinTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
342 "Error - Negative minimum time step. dtMin = "<<getMinTimeStep()<<
")\n");
343 TEUCHOS_TEST_FOR_EXCEPTION(
344 (getMaxTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
346 "Error - Negative maximum time step. dtMax = "<<getMaxTimeStep()<<
")\n");
347 TEUCHOS_TEST_FOR_EXCEPTION(
348 (getMinTimeStep() > getMaxTimeStep() ), std::logic_error,
349 "Error - Inconsistent time step range.\n"
350 " (dtMin = "<<getMinTimeStep()<<
") > (dtMax = "<<getMaxTimeStep()<<
")\n");
351 TEUCHOS_TEST_FOR_EXCEPTION(
352 (getInitTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
354 "Error - Negative initial time step. dtInit = "<<getInitTimeStep()<<
")\n");
355 TEUCHOS_TEST_FOR_EXCEPTION(
356 (getInitTimeStep() < getMinTimeStep() ||
357 getInitTimeStep() > getMaxTimeStep() ),
359 "Error - Initial time step is out of range.\n"
360 <<
" [dtMin, dtMax] = [" << getMinTimeStep() <<
", "
361 << getMaxTimeStep() <<
"]\n"
362 <<
" dtInit = " << getInitTimeStep() <<
"\n");
364 TEUCHOS_TEST_FOR_EXCEPTION(
365 (getInitIndex() > getFinalIndex() ), std::logic_error,
366 "Error - Inconsistent time index range.\n"
367 " (iStepMin = "<<getInitIndex()<<
") > (iStepMax = "
368 <<getFinalIndex()<<
")\n");
370 TEUCHOS_TEST_FOR_EXCEPTION(
371 (getMaxAbsError() < Teuchos::ScalarTraits<Scalar>::zero() ),
373 "Error - Negative maximum time step. errorMaxAbs = "
374 <<getMaxAbsError()<<
")\n");
375 TEUCHOS_TEST_FOR_EXCEPTION(
376 (getMaxRelError() < Teuchos::ScalarTraits<Scalar>::zero() ),
378 "Error - Negative maximum time step. errorMaxRel = "
379 <<getMaxRelError()<<
")\n");
381 TEUCHOS_TEST_FOR_EXCEPTION(
382 (getMinOrder() < Teuchos::ScalarTraits<Scalar>::zero() ),
384 "Error - Negative minimum order. orderMin = "<<getMinOrder()<<
")\n");
385 TEUCHOS_TEST_FOR_EXCEPTION(
386 (getMaxOrder() < Teuchos::ScalarTraits<Scalar>::zero() ), std::logic_error,
387 "Error - Negative maximum order. orderMax = "<<getMaxOrder()<<
")\n");
388 TEUCHOS_TEST_FOR_EXCEPTION(
389 (getMinOrder() > getMaxOrder() ), std::logic_error,
390 "Error - Inconsistent order range.\n"
391 " (orderMin = "<<getMinOrder()<<
") > (orderMax = "
392 <<getMaxOrder()<<
")\n");
393 TEUCHOS_TEST_FOR_EXCEPTION(
394 (getInitOrder() < getMinOrder() || getInitOrder() > getMaxOrder()),
396 "Error - Initial order is out of range.\n"
397 <<
" [orderMin, orderMax] = [" << getMinOrder() <<
", "
398 << getMaxOrder() <<
"]\n"
399 <<
" order = " << getInitOrder() <<
"\n");
401 TEUCHOS_TEST_FOR_EXCEPTION(
402 (getStepType() !=
"Constant" and getStepType() !=
"Variable"),
404 "Error - 'Integrator Step Type' does not equal none of these:\n"
405 <<
" 'Constant' - Integrator will take constant time step sizes.\n"
406 <<
" 'Variable' - Integrator will allow changes to the time step size.\n"
407 <<
" stepType = " << getStepType() <<
"\n");
412 outputTimes_.clear();
413 std::string str = tscPL_->get<std::string>(
"Output Time List");
414 std::string delimiters(
",");
416 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
418 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
419 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
421 std::string token = str.substr(lastPos,pos-lastPos);
422 outputTimes_.push_back(Scalar(std::stod(token)));
423 if(pos==std::string::npos)
break;
425 lastPos = str.find_first_not_of(delimiters, pos);
426 pos = str.find_first_of(delimiters, lastPos);
430 std::sort(outputTimes_.begin(),outputTimes_.end());
431 outputTimes_.erase(std::unique(outputTimes_.begin(),
432 outputTimes_.end() ),
433 outputTimes_.end() );
438 outputIndices_.clear();
439 std::string str = tscPL_->get<std::string>(
"Output Index List");
440 std::string delimiters(
",");
441 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
442 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
443 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
444 std::string token = str.substr(lastPos,pos-lastPos);
445 outputIndices_.push_back(
int(std::stoi(token)));
446 if(pos==std::string::npos)
break;
448 lastPos = str.find_first_not_of(delimiters, pos);
449 pos = str.find_first_of(delimiters, lastPos);
452 Scalar outputIndexInterval = tscPL_->get<
int>(
"Output Index Interval");
453 Scalar output_i = getInitIndex();
454 while (output_i <= getFinalIndex()) {
455 outputIndices_.push_back(output_i);
456 output_i += outputIndexInterval;
460 std::sort(outputIndices_.begin(),outputIndices_.end());
466 template<
class Scalar>
471 using Teuchos::ParameterList;
473 if (stepControlStrategy_ == Teuchos::null){
474 stepControlStrategy_ =
478 if (tscs == Teuchos::null) {
481 if (getStepType() ==
"Constant"){
482 stepControlStrategy_->addStrategy(
484 }
else if (getStepType() ==
"Variable") {
487 RCP<ParameterList> tscsPL =
488 Teuchos::sublist(tscPL_,
"Time Step Control Strategy",
true);
490 std::vector<std::string> tscsLists;
494 std::string str = tscsPL->get<std::string>(
"Time Step Control Strategy List");
495 std::string delimiters(
",");
497 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
499 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
500 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
502 std::string token = str.substr(lastPos,pos-lastPos);
503 tscsLists.push_back(token);
504 if(pos==std::string::npos)
break;
506 lastPos = str.find_first_not_of(delimiters, pos);
507 pos = str.find_first_of(delimiters, lastPos);
511 for(
auto el: tscsLists){
513 RCP<Teuchos::ParameterList> pl =
514 Teuchos::rcp(
new ParameterList(tscsPL->sublist(el)));
516 RCP<TimeStepControlStrategy<Scalar>> ts;
519 if(pl->get<std::string>(
"Name") ==
"Integral Controller")
521 else if(pl->get<std::string>(
"Name") ==
"Basic VS")
524 stepControlStrategy_->addStrategy(ts);
530 stepControlStrategy_->addStrategy(tscs);
536 template<
class Scalar>
537 Teuchos::RCP<const Teuchos::ParameterList>
540 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
542 const double stdMax = double(1.0e+99);
543 pl->set<
double>(
"Initial Time" , 0.0 ,
"Initial time");
544 pl->set<
double>(
"Final Time" , stdMax ,
"Final time");
545 pl->set<
int> (
"Initial Time Index" , 0 ,
"Initial time index");
546 pl->set<
int> (
"Final Time Index" , 1000000,
"Final time index");
547 pl->set<
double>(
"Minimum Time Step" , 0.0 ,
"Minimum time step size");
548 pl->set<
double>(
"Initial Time Step" , 1.0 ,
"Initial time step size");
549 pl->set<
double>(
"Maximum Time Step" , stdMax ,
"Maximum time step size");
550 pl->set<
int> (
"Minimum Order", 0,
551 "Minimum time-integration order. If set to zero (default), the\n"
552 "Stepper minimum order is used.");
553 pl->set<
int> (
"Initial Order", 0,
554 "Initial time-integration order. If set to zero (default), the\n"
555 "Stepper minimum order is used.");
556 pl->set<
int> (
"Maximum Order", 0,
557 "Maximum time-integration order. If set to zero (default), the\n"
558 "Stepper maximum order is used.");
559 pl->set<
double>(
"Maximum Absolute Error", 1.0e-08,
"Maximum absolute error");
560 pl->set<
double>(
"Maximum Relative Error", 1.0e-08,
"Maximum relative error");
562 pl->set<std::string>(
"Integrator Step Type",
"Variable",
563 "'Integrator Step Type' indicates whether the Integrator will allow "
564 "the time step to be modified.\n"
565 " 'Constant' - Integrator will take constant time step sizes.\n"
566 " 'Variable' - Integrator will allow changes to the time step size.\n");
568 pl->set<
bool>(
"Output Exactly On Output Times",
true,
569 "This determines if the timestep size will be adjusted to exactly land\n"
570 "on the output times for 'Variable' timestepping (default=true).\n"
571 "When set to 'false' or for 'Constant' time stepping, the timestep\n"
572 "following the output time will be flagged for output.\n");
574 pl->set<std::string>(
"Output Time List",
"",
575 "Comma deliminated list of output times");
576 pl->set<std::string>(
"Output Index List",
"",
577 "Comma deliminated list of output indices");
578 pl->set<
double>(
"Output Time Interval", stdMax,
"Output time interval");
579 pl->set<
int> (
"Output Index Interval", 1000000,
"Output index interval");
581 pl->set<
int> (
"Maximum Number of Stepper Failures", 10,
582 "Maximum number of Stepper failures");
583 pl->set<
int> (
"Maximum Number of Consecutive Stepper Failures", 5,
584 "Maximum number of consecutive Stepper failures");
585 pl->set<
int> (
"Number of Time Steps", -1,
586 "The number of constant time steps. The actual step size gets computed\n"
587 "on the fly given the size of the time domain. Overides and resets\n"
588 " 'Final Time Index' = 'Initial Time Index' + 'Number of Time Steps'\n"
589 " 'Initial Time Step' = "
590 "('Final Time' - 'Initial Time')/'Number of Time Steps'\n"
591 " 'Integrator Step Type' = 'Constant'\n");
593 Teuchos::RCP<Teuchos::ParameterList> tscsPL = Teuchos::parameterList(
"Time Step Control Strategy");
594 tscsPL->set<std::string>(
"Time Step Control Strategy List",
"");
595 pl->set(
"Time Step Control Strategy", *tscsPL);
600 template <
class Scalar>
601 Teuchos::RCP<Teuchos::ParameterList>
608 template <
class Scalar>
609 Teuchos::RCP<Teuchos::ParameterList>
612 Teuchos::RCP<Teuchos::ParameterList> temp_plist = tscPL_;
613 tscPL_ = Teuchos::null;
619 #endif // Tempus_TimeStepControl_impl_hpp
StepControlStrategy class for TimeStepControl.
virtual void setNumTimeSteps(int numTimeSteps)
virtual void getNextTimeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory, Status &integratorStatus)
Determine the time step size.
std::string description() const
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
virtual void initialize(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
virtual void setTimeStepControlStrategy(Teuchos::RCP< TimeStepControlStrategy< Scalar > > tscs=Teuchos::null)
Set the TimeStepControlStrategy.
Status
Status for the Integrator, the Stepper and the SolutionState.
virtual bool indexInRange(const int iStep) const
Check if time step index is within minimum and maximum index.
TimeStepControl manages the time step size. There several mechanicisms that effect the time step size...
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl)
StepControlStrategy class for TimeStepControl.
StepControlStrategy class for TimeStepControl.
virtual bool timeInRange(const Scalar time) const
Check if time is within minimum and maximum time.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
StepControlStrategy class for TimeStepControl.
StepControlStrategy class for TimeStepControl.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
TimeStepControl(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Constructor.