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)
33 outputAdjustedDt_(false),
35 stepControlStrategy_(Teuchos::null),
42 template<
class Scalar>
49 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::TimeStepControl::getNextTimeStep()");
51 RCP<Teuchos::FancyOStream> out = this->getOStream();
52 Teuchos::OSTab ostab(out,0,
"getNextTimeStep");
55 auto changeDT = [] (
int istep, Scalar dt_old, Scalar dt_new,
58 std::stringstream message;
59 message << std::scientific
60 <<std::setw(6)<<std::setprecision(3)<<istep
61 <<
" * (dt = "<<std::setw(9)<<std::setprecision(3)<<dt_old
62 <<
", new = "<<std::setw(9)<<std::setprecision(3)<<dt_new
63 <<
") " << reason << std::endl;
67 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
68 const Scalar lastTime = solutionHistory->getCurrentState()->getTime();
69 const int iStep = workingState->getIndex();
70 int order = workingState->getOrder();
71 Scalar dt = workingState->getTimeStep();
74 RCP<StepperState<Scalar> > stepperState = workingState->getStepperState();
76 if (getStepType() ==
"Variable") {
78 if (outputAdjustedDt_ ==
true) {
79 if (printDtChanges_) *out << changeDT(iStep, dt, dtAfterOutput_,
80 "Reset dt after output.");
82 outputAdjustedDt_ =
false;
87 if (printDtChanges_) *out << changeDT(iStep, dt, getInitTimeStep(),
88 "Reset dt to initial dt.");
89 dt = getInitTimeStep();
92 if (dt < getMinTimeStep()) {
93 if (printDtChanges_) *out << changeDT(iStep, dt, getMinTimeStep(),
94 "Reset dt to minimum dt.");
95 dt = getMinTimeStep();
100 workingState->setTimeStep(dt);
103 stepControlStrategy_->getNextTimeStep(*
this, solutionHistory,
107 order = workingState->getOrder();
108 dt = workingState->getTimeStep();
110 if (getStepType() ==
"Variable") {
111 if (dt < getMinTimeStep()) {
112 if (printDtChanges_) *out << changeDT(iStep, dt, getMinTimeStep(),
113 "dt is too small. Resetting to minimum dt.");
114 dt = getMinTimeStep();
116 if (dt > getMaxTimeStep()) {
117 if (printDtChanges_) *out << changeDT(iStep, dt, getMaxTimeStep(),
118 "dt is too large. Resetting to maximum dt.");
119 dt = getMaxTimeStep();
125 std::vector<int>::const_iterator it =
126 std::find(outputIndices_.begin(), outputIndices_.end(), iStep);
127 if (it != outputIndices_.end()) output =
true;
129 const int iInterval = tscPL_->get<
int>(
"Output Index Interval");
130 if ( (iStep - getInitIndex()) % iInterval == 0) output =
true;
134 Scalar reltol = 1.0e-6;
135 Scalar endTime = lastTime+dt+getMinTimeStep();
138 if (getStepType() ==
"Constant") endTime = lastTime+dt;
139 bool checkOutput =
false;
140 Scalar oTime = getInitTime();
141 for (
size_t i=0; i < outputTimes_.size(); ++i) {
142 oTime = outputTimes_[i];
143 if (lastTime < oTime && oTime <= endTime) {
148 const Scalar tInterval = tscPL_->get<
double>(
"Output Time Interval");
149 Scalar oTime2 = ceil((lastTime-getInitTime())/tInterval)*tInterval
151 if (lastTime < oTime2 && oTime2 <= endTime) {
152 if (checkOutput ==
true) {
153 if (oTime2 < oTime) oTime = oTime2;
160 if (checkOutput ==
true) {
161 const bool outputExactly =
162 tscPL_->get<
bool>(
"Output Exactly On Output Times");
163 if (getStepType() ==
"Variable" && outputExactly ==
true) {
165 if (std::abs((lastTime+dt-oTime)/(lastTime+dt)) < reltol) {
167 if (printDtChanges_) *out << changeDT(iStep, dt, oTime - lastTime,
168 "Adjusting dt for numerical roundoff to hit the next output time.");
171 outputAdjustedDt_ =
true;
173 dt = oTime - lastTime;
174 }
else if (lastTime*(1.0+reltol) < oTime &&
175 oTime < (lastTime+dt-getMinTimeStep())*(1.0+reltol)) {
177 if (printDtChanges_) *out << changeDT(iStep, dt, oTime - lastTime,
178 "Adjusting dt to hit the next output time.");
182 outputAdjustedDt_ =
true;
184 dt = oTime - lastTime;
186 if (printDtChanges_) *out << changeDT(iStep, dt, (oTime - lastTime)/2.0,
187 "The next output time is within the minimum dt of the next time. "
188 "Adjusting dt to take two steps.");
192 dt = (oTime - lastTime)/2.0;
204 if ((lastTime + dt > getFinalTime() ) ||
205 (std::abs((lastTime+dt-getFinalTime())/(lastTime+dt)) < reltol)) {
206 if (printDtChanges_) *out << changeDT(iStep, dt, getFinalTime() - lastTime,
207 "Adjusting dt to hit final time.");
208 dt = getFinalTime() - lastTime;
212 TEUCHOS_TEST_FOR_EXCEPTION( dt <= Scalar(0.0), std::out_of_range,
213 "Error - Time step is not positive. dt = " << dt <<
"\n");
216 TEUCHOS_TEST_FOR_EXCEPTION(
217 (lastTime + dt < getInitTime()), std::out_of_range,
218 "Error - Time step does not move time INTO time range.\n"
219 " [timeMin, timeMax] = [" << getInitTime() <<
", "
220 << getFinalTime() <<
"]\n"
221 " T + dt = " << lastTime <<
" + "<< dt <<
" = " << lastTime + dt <<
"\n");
223 TEUCHOS_TEST_FOR_EXCEPTION(
224 (lastTime + dt > getFinalTime()), std::out_of_range,
225 "Error - Time step move time OUT OF time range.\n"
226 " [timeMin, timeMax] = [" << getInitTime() <<
", "
227 << getFinalTime() <<
"]\n"
228 " T + dt = " << lastTime <<
" + "<< dt <<
" = " << lastTime + dt <<
"\n");
230 workingState->setOrder(order);
231 workingState->setTimeStep(dt);
232 workingState->setTime(lastTime + dt);
233 workingState->setOutput(output);
240 template<
class Scalar>
242 const Scalar relTol = 1.0e-14;
243 bool tir = (getInitTime()*(1.0-relTol) <= time and
244 time < getFinalTime()*(1.0-relTol));
249 template<
class Scalar>
251 bool iir = (getInitIndex() <= iStep and iStep < getFinalIndex());
256 template<
class Scalar>
259 if (numTimeSteps >= 0) {
260 tscPL_->set<
int> (
"Number of Time Steps", numTimeSteps);
261 const int initIndex = getInitIndex();
262 tscPL_->set<
int> (
"Final Time Index", initIndex + numTimeSteps);
263 const double initTime = tscPL_->get<
double>(
"Initial Time");
264 const double finalTime = tscPL_->get<
double>(
"Final Time");
265 double initTimeStep = (finalTime - initTime)/numTimeSteps;
266 if (numTimeSteps == 0) initTimeStep = Scalar(0.0);
267 tscPL_->set<
double> (
"Initial Time Step", initTimeStep);
268 tscPL_->set<
double> (
"Minimum Time Step", initTimeStep);
269 tscPL_->set<
double> (
"Maximum Time Step", initTimeStep);
270 tscPL_->set<std::string>(
"Integrator Step Type",
"Constant");
272 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
273 Teuchos::OSTab ostab(out,1,
"setNumTimeSteps");
274 *out <<
"Warning - Found 'Number of Time Steps' = " << getNumTimeSteps()
275 <<
" Set the following parameters: \n"
276 <<
" 'Final Time Index' = " << getFinalIndex() <<
"\n"
277 <<
" 'Initial Time Step' = " << getInitTimeStep() <<
"\n"
278 <<
" 'Integrator Step Type' = " << getStepType() << std::endl;
283 template<
class Scalar>
286 std::string name =
"Tempus::TimeStepControl";
291 template<
class Scalar>
293 Teuchos::FancyOStream &out,
294 const Teuchos::EVerbosityLevel verbLevel)
const
296 if (verbLevel == Teuchos::VERB_EXTREME) {
297 out << description() <<
"::describe:" << std::endl
298 <<
"pList = " << tscPL_ << std::endl;
303 template <
class Scalar>
305 Teuchos::RCP<Teuchos::ParameterList>
const& pList)
307 if (pList == Teuchos::null) {
309 if (tscPL_ == Teuchos::null) {
310 tscPL_ = Teuchos::parameterList(
"TimeStepControl");
311 *tscPL_ = *(this->getValidParameters());
316 tscPL_->validateParametersAndSetDefaults(*this->getValidParameters(), 0);
319 if (getStepType() ==
"Constant") {
320 const double initTimeStep = tscPL_->get<
double>(
"Initial Time Step");
321 tscPL_->set<
double> (
"Minimum Time Step", initTimeStep);
322 tscPL_->set<
double> (
"Maximum Time Step", initTimeStep);
324 setNumTimeSteps(getNumTimeSteps());
327 setTimeStepControlStrategy();
329 TEUCHOS_TEST_FOR_EXCEPTION(
330 (getInitTime() > getFinalTime() ), std::logic_error,
331 "Error - Inconsistent time range.\n"
332 " (timeMin = "<<getInitTime()<<
") > (timeMax = "<<getFinalTime()<<
")\n");
334 TEUCHOS_TEST_FOR_EXCEPTION(
335 (getMinTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
337 "Error - Negative minimum time step. dtMin = "<<getMinTimeStep()<<
")\n");
338 TEUCHOS_TEST_FOR_EXCEPTION(
339 (getMaxTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
341 "Error - Negative maximum time step. dtMax = "<<getMaxTimeStep()<<
")\n");
342 TEUCHOS_TEST_FOR_EXCEPTION(
343 (getMinTimeStep() > getMaxTimeStep() ), std::logic_error,
344 "Error - Inconsistent time step range.\n"
345 " (dtMin = "<<getMinTimeStep()<<
") > (dtMax = "<<getMaxTimeStep()<<
")\n");
346 TEUCHOS_TEST_FOR_EXCEPTION(
347 (getInitTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
349 "Error - Negative initial time step. dtInit = "<<getInitTimeStep()<<
")\n");
350 TEUCHOS_TEST_FOR_EXCEPTION(
351 (getInitTimeStep() < getMinTimeStep() ||
352 getInitTimeStep() > getMaxTimeStep() ),
354 "Error - Initial time step is out of range.\n"
355 <<
" [dtMin, dtMax] = [" << getMinTimeStep() <<
", "
356 << getMaxTimeStep() <<
"]\n"
357 <<
" dtInit = " << getInitTimeStep() <<
"\n");
359 TEUCHOS_TEST_FOR_EXCEPTION(
360 (getInitIndex() > getFinalIndex() ), std::logic_error,
361 "Error - Inconsistent time index range.\n"
362 " (iStepMin = "<<getInitIndex()<<
") > (iStepMax = "
363 <<getFinalIndex()<<
")\n");
365 TEUCHOS_TEST_FOR_EXCEPTION(
366 (getMaxAbsError() < Teuchos::ScalarTraits<Scalar>::zero() ),
368 "Error - Negative maximum time step. errorMaxAbs = "
369 <<getMaxAbsError()<<
")\n");
370 TEUCHOS_TEST_FOR_EXCEPTION(
371 (getMaxRelError() < Teuchos::ScalarTraits<Scalar>::zero() ),
373 "Error - Negative maximum time step. errorMaxRel = "
374 <<getMaxRelError()<<
")\n");
376 TEUCHOS_TEST_FOR_EXCEPTION(
377 (getMinOrder() < Teuchos::ScalarTraits<Scalar>::zero() ),
379 "Error - Negative minimum order. orderMin = "<<getMinOrder()<<
")\n");
380 TEUCHOS_TEST_FOR_EXCEPTION(
381 (getMaxOrder() < Teuchos::ScalarTraits<Scalar>::zero() ), std::logic_error,
382 "Error - Negative maximum order. orderMax = "<<getMaxOrder()<<
")\n");
383 TEUCHOS_TEST_FOR_EXCEPTION(
384 (getMinOrder() > getMaxOrder() ), std::logic_error,
385 "Error - Inconsistent order range.\n"
386 " (orderMin = "<<getMinOrder()<<
") > (orderMax = "
387 <<getMaxOrder()<<
")\n");
388 TEUCHOS_TEST_FOR_EXCEPTION(
389 (getInitOrder() < getMinOrder() || getInitOrder() > getMaxOrder()),
391 "Error - Initial order is out of range.\n"
392 <<
" [orderMin, orderMax] = [" << getMinOrder() <<
", "
393 << getMaxOrder() <<
"]\n"
394 <<
" order = " << getInitOrder() <<
"\n");
396 TEUCHOS_TEST_FOR_EXCEPTION(
397 (getStepType() !=
"Constant" and getStepType() !=
"Variable"),
399 "Error - 'Integrator Step Type' does not equal none of these:\n"
400 <<
" 'Constant' - Integrator will take constant time step sizes.\n"
401 <<
" 'Variable' - Integrator will allow changes to the time step size.\n"
402 <<
" stepType = " << getStepType() <<
"\n");
407 outputTimes_.clear();
408 std::string str = tscPL_->get<std::string>(
"Output Time List");
409 std::string delimiters(
",");
411 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
413 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
414 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
416 std::string token = str.substr(lastPos,pos-lastPos);
417 outputTimes_.push_back(Scalar(std::stod(token)));
418 if(pos==std::string::npos)
break;
420 lastPos = str.find_first_not_of(delimiters, pos);
421 pos = str.find_first_of(delimiters, lastPos);
425 std::sort(outputTimes_.begin(),outputTimes_.end());
426 outputTimes_.erase(std::unique(outputTimes_.begin(),
427 outputTimes_.end() ),
428 outputTimes_.end() );
433 outputIndices_.clear();
434 std::string str = tscPL_->get<std::string>(
"Output Index List");
435 std::string delimiters(
",");
436 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
437 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
438 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
439 std::string token = str.substr(lastPos,pos-lastPos);
440 outputIndices_.push_back(
int(std::stoi(token)));
441 if(pos==std::string::npos)
break;
443 lastPos = str.find_first_not_of(delimiters, pos);
444 pos = str.find_first_of(delimiters, lastPos);
447 Scalar outputIndexInterval = tscPL_->get<
int>(
"Output Index Interval");
448 Scalar output_i = getInitIndex();
449 while (output_i <= getFinalIndex()) {
450 outputIndices_.push_back(output_i);
451 output_i += outputIndexInterval;
455 std::sort(outputIndices_.begin(),outputIndices_.end());
461 template<
class Scalar>
466 using Teuchos::ParameterList;
468 if (stepControlStrategy_ == Teuchos::null){
469 stepControlStrategy_ =
473 if (tscs == Teuchos::null) {
476 if (getStepType() ==
"Constant"){
477 stepControlStrategy_->addStrategy(
479 }
else if (getStepType() ==
"Variable") {
482 RCP<ParameterList> tscsPL =
483 Teuchos::sublist(tscPL_,
"Time Step Control Strategy",
true);
485 std::vector<std::string> tscsLists;
489 std::string str = tscsPL->get<std::string>(
"Time Step Control Strategy List");
490 std::string delimiters(
",");
492 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
494 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
495 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
497 std::string token = str.substr(lastPos,pos-lastPos);
498 tscsLists.push_back(token);
499 if(pos==std::string::npos)
break;
501 lastPos = str.find_first_not_of(delimiters, pos);
502 pos = str.find_first_of(delimiters, lastPos);
506 for(
auto el: tscsLists){
508 RCP<Teuchos::ParameterList> pl =
509 Teuchos::rcp(
new ParameterList(tscsPL->sublist(el)));
511 RCP<TimeStepControlStrategy<Scalar>> ts;
514 if(pl->get<std::string>(
"Name") ==
"Integral Controller")
516 else if(pl->get<std::string>(
"Name") ==
"Basic VS")
519 stepControlStrategy_->addStrategy(ts);
525 stepControlStrategy_->addStrategy(tscs);
531 template<
class Scalar>
532 Teuchos::RCP<const Teuchos::ParameterList>
535 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
537 const double stdMax = double(1.0e+99);
538 pl->set<
double>(
"Initial Time" , 0.0 ,
"Initial time");
539 pl->set<
double>(
"Final Time" , stdMax ,
"Final time");
540 pl->set<
int> (
"Initial Time Index" , 0 ,
"Initial time index");
541 pl->set<
int> (
"Final Time Index" , 1000000,
"Final time index");
542 pl->set<
double>(
"Minimum Time Step" , 0.0 ,
"Minimum time step size");
543 pl->set<
double>(
"Initial Time Step" , 1.0 ,
"Initial time step size");
544 pl->set<
double>(
"Maximum Time Step" , stdMax ,
"Maximum time step size");
545 pl->set<
int> (
"Minimum Order", 0,
546 "Minimum time-integration order. If set to zero (default), the\n"
547 "Stepper minimum order is used.");
548 pl->set<
int> (
"Initial Order", 0,
549 "Initial time-integration order. If set to zero (default), the\n"
550 "Stepper minimum order is used.");
551 pl->set<
int> (
"Maximum Order", 0,
552 "Maximum time-integration order. If set to zero (default), the\n"
553 "Stepper maximum order is used.");
554 pl->set<
double>(
"Maximum Absolute Error", 1.0e-08,
"Maximum absolute error");
555 pl->set<
double>(
"Maximum Relative Error", 1.0e-08,
"Maximum relative error");
557 pl->set<std::string>(
"Integrator Step Type",
"Variable",
558 "'Integrator Step Type' indicates whether the Integrator will allow "
559 "the time step to be modified.\n"
560 " 'Constant' - Integrator will take constant time step sizes.\n"
561 " 'Variable' - Integrator will allow changes to the time step size.\n");
563 pl->set<
bool>(
"Output Exactly On Output Times",
true,
564 "This determines if the timestep size will be adjusted to exactly land\n"
565 "on the output times for 'Variable' timestepping (default=true).\n"
566 "When set to 'false' or for 'Constant' time stepping, the timestep\n"
567 "following the output time will be flagged for output.\n");
569 pl->set<std::string>(
"Output Time List",
"",
570 "Comma deliminated list of output times");
571 pl->set<std::string>(
"Output Index List",
"",
572 "Comma deliminated list of output indices");
573 pl->set<
double>(
"Output Time Interval", stdMax,
"Output time interval");
574 pl->set<
int> (
"Output Index Interval", 1000000,
"Output index interval");
576 pl->set<
int> (
"Maximum Number of Stepper Failures", 10,
577 "Maximum number of Stepper failures");
578 pl->set<
int> (
"Maximum Number of Consecutive Stepper Failures", 5,
579 "Maximum number of consecutive Stepper failures");
580 pl->set<
int> (
"Number of Time Steps", -1,
581 "The number of constant time steps. The actual step size gets computed\n"
582 "on the fly given the size of the time domain. Overides and resets\n"
583 " 'Final Time Index' = 'Initial Time Index' + 'Number of Time Steps'\n"
584 " 'Initial Time Step' = "
585 "('Final Time' - 'Initial Time')/'Number of Time Steps'\n"
586 " 'Integrator Step Type' = 'Constant'\n");
588 Teuchos::RCP<Teuchos::ParameterList> tscsPL = Teuchos::parameterList(
"Time Step Control Strategy");
589 tscsPL->set<std::string>(
"Time Step Control Strategy List",
"");
590 pl->set(
"Time Step Control Strategy", *tscsPL);
595 template <
class Scalar>
596 Teuchos::RCP<Teuchos::ParameterList>
603 template <
class Scalar>
604 Teuchos::RCP<Teuchos::ParameterList>
607 Teuchos::RCP<Teuchos::ParameterList> temp_plist = tscPL_;
608 tscPL_ = Teuchos::null;
614 #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.
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.