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,1,
"getNextTimeStep");
60 Teuchos::as<int>(Teuchos::VERB_NONE);
62 auto changeDT = [] (Scalar dt_old, Scalar dt_new, std::string reason) {
63 std::stringstream message;
65 " (dt = "<<std::scientific<<std::setw(9)<<std::setprecision(3)<<dt_old
66 <<
", new = "<<std::scientific<<std::setw(9)<<std::setprecision(3)<<dt_new
67 <<
") " << reason << std::endl;
71 RCP<SolutionState<Scalar> > workingState=
solutionHistory->getWorkingState();
72 RCP<SolutionStateMetaData<Scalar> > metaData = workingState->getMetaData();
74 const int iStep = metaData->getIStep();
75 int order = metaData->getOrder();
76 Scalar dt = metaData->getDt();
77 bool output = metaData->getOutput();
79 RCP<StepperState<Scalar> > stepperState = workingState->getStepperState();
83 if (getStepType() ==
"Variable") {
85 if (outputAdjustedDt_ ==
true) {
86 if (printChanges) *out << changeDT(dt, dtAfterOutput_,
87 "Reset dt after output.");
89 outputAdjustedDt_ =
false;
94 if (printChanges) *out << changeDT(dt, getInitTimeStep(),
95 "Reset dt to initial dt.");
96 dt = getInitTimeStep();
99 if (dt < getMinTimeStep()) {
100 if (printChanges) *out << changeDT(dt, getMinTimeStep(),
101 "Reset dt to minimum dt.");
102 dt = getMinTimeStep();
114 order = metaData->getOrder();
115 dt = metaData->getDt();
117 if (getStepType() ==
"Variable") {
118 if (dt < getMinTimeStep()) {
119 if (printChanges) *out << changeDT(dt, getMinTimeStep(),
120 "dt is too small. Resetting to minimum dt.");
121 dt = getMinTimeStep();
123 if (dt > getMaxTimeStep()) {
124 if (printChanges) *out << changeDT(dt, getMaxTimeStep(),
125 "dt is too large. Resetting to maximum dt.");
126 dt = getMaxTimeStep();
132 std::vector<int>::const_iterator it =
133 std::find(outputIndices_.begin(), outputIndices_.end(), iStep);
134 if (it != outputIndices_.end()) output =
true;
137 Scalar reltol = 1.0e-6;
138 for (
size_t i=0; i < outputTimes_.size(); ++i) {
139 const Scalar oTime = outputTimes_[i];
140 if (lastTime < oTime && oTime <= lastTime+dt+getMinTimeStep()) {
141 if (std::abs((lastTime+dt-oTime)/(lastTime+dt)) < reltol) {
143 if (getStepType() ==
"Variable") {
144 if (printChanges) *out << changeDT(dt, oTime - lastTime,
145 "Adjusting dt for numerical roundoff to hit the next output time.");
148 outputAdjustedDt_ =
true;
150 dt = oTime - lastTime;
152 }
else if (lastTime*(1.0+reltol) < oTime &&
153 oTime < (lastTime+dt-getMinTimeStep())*(1.0+reltol)) {
155 if (getStepType() ==
"Variable") {
156 if (printChanges) *out << changeDT(dt, oTime - lastTime,
157 "Adjusting dt to hit the next output time.");
161 outputAdjustedDt_ =
true;
163 dt = oTime - lastTime;
166 if (getStepType() ==
"Variable") {
167 if (printChanges) *out << changeDT(dt, (oTime - lastTime)/2.0,
168 "The next output time is within the minimum dt of the next time. "
169 "Adjusting dt to take two steps.");
173 dt = (oTime - lastTime)/2.0;
182 if ((lastTime + dt > getFinalTime() ) ||
183 (std::abs((lastTime+dt-getFinalTime())/(lastTime+dt)) < reltol)) {
184 if (printChanges) *out << changeDT(dt, getFinalTime() - lastTime,
185 "Adjusting dt to hit the final time.");
186 dt = getFinalTime() - lastTime;
190 TEUCHOS_TEST_FOR_EXCEPTION(
191 (lastTime + dt < getInitTime()), std::out_of_range,
192 "Error - Time step does not move time INTO time range.\n"
193 " [timeMin, timeMax] = [" << getInitTime() <<
", "
194 << getFinalTime() <<
"]\n"
195 " T + dt = " << lastTime <<
" + "<< dt <<
" = " << lastTime + dt <<
"\n");
197 TEUCHOS_TEST_FOR_EXCEPTION(
198 (lastTime + dt > getFinalTime()), std::out_of_range,
199 "Error - Time step move time OUT OF time range.\n"
200 " [timeMin, timeMax] = [" << getInitTime() <<
", "
201 << getFinalTime() <<
"]\n"
202 " T + dt = " << lastTime <<
" + "<< dt <<
" = " << lastTime + dt <<
"\n");
204 metaData->setOrder(order);
206 metaData->setTime(lastTime + dt);
207 metaData->setOutput(output);
214 template<
class Scalar>
216 const Scalar relTol = 1.0e-14;
217 bool tir = (getInitTime()*(1.0-relTol) <= time and
218 time < getFinalTime()*(1.0-relTol));
223 template<
class Scalar>
225 bool iir = (getInitIndex() <= iStep and iStep < getFinalIndex());
230 template<
class Scalar>
233 if (numTimeSteps >= 0) {
234 tscPL_->set<
int> (
"Number of Time Steps", numTimeSteps);
235 const int initIndex = getInitIndex();
236 tscPL_->set<
int> (
"Final Time Index", initIndex + numTimeSteps);
237 const double initTime = tscPL_->get<
double>(
"Initial Time");
238 const double finalTime = tscPL_->get<
double>(
"Final Time");
239 double initTimeStep = (finalTime - initTime)/numTimeSteps;
240 if (numTimeSteps == 0) initTimeStep = Scalar(0.0);
241 tscPL_->set<
double> (
"Initial Time Step", initTimeStep);
242 tscPL_->set<
double> (
"Minimum Time Step", initTimeStep);
243 tscPL_->set<
double> (
"Maximum Time Step", initTimeStep);
244 tscPL_->set<std::string>(
"Integrator Step Type",
"Constant");
246 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
247 Teuchos::OSTab ostab(out,1,
"setNumTimeSteps");
248 *out <<
"Warning - Found 'Number of Time Steps' = " << getNumTimeSteps()
249 <<
" Set the following parameters: \n"
250 <<
" 'Final Time Index' = " << getFinalIndex() <<
"\n"
251 <<
" 'Initial Time Step' = " << getInitTimeStep() <<
"\n"
252 <<
" 'Integrator Step Type' = " << getStepType() << std::endl;
257 template<
class Scalar>
260 std::string name =
"Tempus::TimeStepControl";
265 template<
class Scalar>
267 Teuchos::FancyOStream &out,
268 const Teuchos::EVerbosityLevel verbLevel)
const
270 if (verbLevel == Teuchos::VERB_EXTREME) {
271 out << description() <<
"::describe:" << std::endl
272 <<
"pList = " << tscPL_ << std::endl;
277 template <
class Scalar>
279 Teuchos::RCP<Teuchos::ParameterList>
const& pList)
281 if (pList == Teuchos::null) {
283 if (tscPL_ == Teuchos::null) {
284 tscPL_ = Teuchos::parameterList(
"TimeStepControl");
285 *tscPL_ = *(this->getValidParameters());
290 tscPL_->validateParametersAndSetDefaults(*this->getValidParameters(), 0);
293 if (getStepType() ==
"Constant") {
294 const double initTimeStep = tscPL_->get<
double>(
"Initial Time Step");
295 tscPL_->set<
double> (
"Minimum Time Step", initTimeStep);
296 tscPL_->set<
double> (
"Maximum Time Step", initTimeStep);
298 setNumTimeSteps(getNumTimeSteps());
301 setTimeStepControlStrategy();
303 TEUCHOS_TEST_FOR_EXCEPTION(
304 (getInitTime() > getFinalTime() ), std::logic_error,
305 "Error - Inconsistent time range.\n"
306 " (timeMin = "<<getInitTime()<<
") > (timeMax = "<<getFinalTime()<<
")\n");
308 TEUCHOS_TEST_FOR_EXCEPTION(
309 (getMinTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
311 "Error - Negative minimum time step. dtMin = "<<getMinTimeStep()<<
")\n");
312 TEUCHOS_TEST_FOR_EXCEPTION(
313 (getMaxTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
315 "Error - Negative maximum time step. dtMax = "<<getMaxTimeStep()<<
")\n");
316 TEUCHOS_TEST_FOR_EXCEPTION(
317 (getMinTimeStep() > getMaxTimeStep() ), std::logic_error,
318 "Error - Inconsistent time step range.\n"
319 " (dtMin = "<<getMinTimeStep()<<
") > (dtMax = "<<getMaxTimeStep()<<
")\n");
320 TEUCHOS_TEST_FOR_EXCEPTION(
321 (getInitTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
323 "Error - Negative initial time step. dtInit = "<<getInitTimeStep()<<
")\n");
324 TEUCHOS_TEST_FOR_EXCEPTION(
325 (getInitTimeStep() < getMinTimeStep() ||
326 getInitTimeStep() > getMaxTimeStep() ),
328 "Error - Initial time step is out of range.\n"
329 <<
" [dtMin, dtMax] = [" << getMinTimeStep() <<
", "
330 << getMaxTimeStep() <<
"]\n"
331 <<
" dtInit = " << getInitTimeStep() <<
"\n");
333 TEUCHOS_TEST_FOR_EXCEPTION(
334 (getInitIndex() > getFinalIndex() ), std::logic_error,
335 "Error - Inconsistent time index range.\n"
336 " (iStepMin = "<<getInitIndex()<<
") > (iStepMax = "
337 <<getFinalIndex()<<
")\n");
339 TEUCHOS_TEST_FOR_EXCEPTION(
340 (getMaxAbsError() < Teuchos::ScalarTraits<Scalar>::zero() ),
342 "Error - Negative maximum time step. errorMaxAbs = "
343 <<getMaxAbsError()<<
")\n");
344 TEUCHOS_TEST_FOR_EXCEPTION(
345 (getMaxRelError() < Teuchos::ScalarTraits<Scalar>::zero() ),
347 "Error - Negative maximum time step. errorMaxRel = "
348 <<getMaxRelError()<<
")\n");
350 TEUCHOS_TEST_FOR_EXCEPTION(
351 (getMinOrder() < Teuchos::ScalarTraits<Scalar>::zero() ),
353 "Error - Negative minimum order. orderMin = "<<getMinOrder()<<
")\n");
354 TEUCHOS_TEST_FOR_EXCEPTION(
355 (getMaxOrder() < Teuchos::ScalarTraits<Scalar>::zero() ), std::logic_error,
356 "Error - Negative maximum order. orderMax = "<<getMaxOrder()<<
")\n");
357 TEUCHOS_TEST_FOR_EXCEPTION(
358 (getMinOrder() > getMaxOrder() ), std::logic_error,
359 "Error - Inconsistent order range.\n"
360 " (orderMin = "<<getMinOrder()<<
") > (orderMax = "
361 <<getMaxOrder()<<
")\n");
362 TEUCHOS_TEST_FOR_EXCEPTION(
363 (getInitOrder() < getMinOrder() || getInitOrder() > getMaxOrder()),
365 "Error - Initial order is out of range.\n"
366 <<
" [orderMin, orderMax] = [" << getMinOrder() <<
", "
367 << getMaxOrder() <<
"]\n"
368 <<
" order = " << getInitOrder() <<
"\n");
370 TEUCHOS_TEST_FOR_EXCEPTION(
371 (getStepType() !=
"Constant" and getStepType() !=
"Variable"),
373 "Error - 'Integrator Step Type' does not equal none of these:\n"
374 <<
" 'Constant' - Integrator will take constant time step sizes.\n"
375 <<
" 'Variable' - Integrator will allow changes to the time step size.\n"
376 <<
" stepType = " << getStepType() <<
"\n");
381 outputTimes_.clear();
382 std::string str = tscPL_->get<std::string>(
"Output Time List");
383 std::string delimiters(
",");
385 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
387 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
388 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
390 std::string token = str.substr(lastPos,pos-lastPos);
391 outputTimes_.push_back(Scalar(std::stod(token)));
392 if(pos==std::string::npos)
break;
394 lastPos = str.find_first_not_of(delimiters, pos);
395 pos = str.find_first_of(delimiters, lastPos);
398 Scalar outputTimeInterval = tscPL_->get<
double>(
"Output Time Interval");
399 Scalar output_t = getInitTime();
400 while (output_t <= getFinalTime()) {
401 outputTimes_.push_back(output_t);
402 output_t += outputTimeInterval;
406 std::sort(outputTimes_.begin(),outputTimes_.end());
411 outputIndices_.clear();
412 std::string str = tscPL_->get<std::string>(
"Output Index List");
413 std::string delimiters(
",");
414 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
415 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
416 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
417 std::string token = str.substr(lastPos,pos-lastPos);
418 outputIndices_.push_back(
int(std::stoi(token)));
419 if(pos==std::string::npos)
break;
421 lastPos = str.find_first_not_of(delimiters, pos);
422 pos = str.find_first_of(delimiters, lastPos);
425 Scalar outputIndexInterval = tscPL_->get<
int>(
"Output Index Interval");
426 Scalar output_i = getInitIndex();
427 while (output_i <= getFinalIndex()) {
428 outputIndices_.push_back(output_i);
429 output_i += outputIndexInterval;
433 std::sort(outputIndices_.begin(),outputIndices_.end());
439 template<
class Scalar>
444 using Teuchos::ParameterList;
446 if (stepControlStrategy_ == Teuchos::null){
447 stepControlStrategy_ =
451 if (tscs == Teuchos::null) {
454 if (getStepType() ==
"Constant"){
455 stepControlStrategy_->addStrategy(
457 }
else if (getStepType() ==
"Variable") {
460 RCP<ParameterList> tscsPL =
461 Teuchos::sublist(tscPL_,
"Time Step Control Strategy",
true);
463 std::vector<std::string> tscsLists;
467 std::string str = tscsPL->get<std::string>(
"Time Step Control Strategy List");
468 std::string delimiters(
",");
470 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
472 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
473 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
475 std::string token = str.substr(lastPos,pos-lastPos);
476 tscsLists.push_back(token);
477 if(pos==std::string::npos)
break;
479 lastPos = str.find_first_not_of(delimiters, pos);
480 pos = str.find_first_of(delimiters, lastPos);
484 for(
auto el: tscsLists){
486 RCP<Teuchos::ParameterList> pl =
487 Teuchos::rcp(
new ParameterList(tscsPL->sublist(el)));
489 RCP<TimeStepControlStrategy<Scalar>> ts;
492 if(pl->get<std::string>(
"Name") ==
"Integral Controller")
494 else if(pl->get<std::string>(
"Name") ==
"Basic VS")
497 stepControlStrategy_->addStrategy(ts);
503 stepControlStrategy_->addStrategy(tscs);
509 template<
class Scalar>
510 Teuchos::RCP<const Teuchos::ParameterList>
513 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
515 const double stdMax = std::numeric_limits<double>::max();
516 pl->set<
double>(
"Initial Time" , 0.0 ,
"Initial time");
517 pl->set<
double>(
"Final Time" , stdMax ,
"Final time");
518 pl->set<
int> (
"Initial Time Index" , 0 ,
"Initial time index");
519 pl->set<
int> (
"Final Time Index" , 1000000,
"Final time index");
520 pl->set<
double>(
"Minimum Time Step" , 0.0 ,
"Minimum time step size");
521 pl->set<
double>(
"Initial Time Step" , 1.0 ,
"Initial time step size");
522 pl->set<
double>(
"Maximum Time Step" , stdMax ,
"Maximum time step size");
523 pl->set<
int> (
"Minimum Order", 0,
524 "Minimum time-integration order. If set to zero (default), the\n"
525 "Stepper minimum order is used.");
526 pl->set<
int> (
"Initial Order", 0,
527 "Initial time-integration order. If set to zero (default), the\n"
528 "Stepper minimum order is used.");
529 pl->set<
int> (
"Maximum Order", 0,
530 "Maximum time-integration order. If set to zero (default), the\n"
531 "Stepper maximum order is used.");
532 pl->set<
double>(
"Maximum Absolute Error", 1.0e-08,
"Maximum absolute error");
533 pl->set<
double>(
"Maximum Relative Error", 1.0e-08,
"Maximum relative error");
535 pl->set<std::string>(
"Integrator Step Type",
"Variable",
536 "'Integrator Step Type' indicates whether the Integrator will allow "
537 "the time step to be modified.\n"
538 " 'Constant' - Integrator will take constant time step sizes.\n"
539 " 'Variable' - Integrator will allow changes to the time step size.\n");
541 pl->set<std::string>(
"Output Time List",
"",
542 "Comma deliminated list of output times");
543 pl->set<std::string>(
"Output Index List",
"",
544 "Comma deliminated list of output indices");
545 pl->set<
double>(
"Output Time Interval", stdMax,
"Output time interval");
546 pl->set<
int> (
"Output Index Interval", 1000000,
"Output index interval");
548 pl->set<
int> (
"Maximum Number of Stepper Failures", 10,
549 "Maximum number of Stepper failures");
550 pl->set<
int> (
"Maximum Number of Consecutive Stepper Failures", 5,
551 "Maximum number of consecutive Stepper failures");
552 pl->set<
int> (
"Number of Time Steps", -1,
553 "The number of constant time steps. The actual step size gets computed\n"
554 "on the fly given the size of the time domain. Overides and resets\n"
555 " 'Final Time Index' = 'Initial Time Index' + 'Number of Time Steps'\n"
556 " 'Initial Time Step' = "
557 "('Final Time' - 'Initial Time')/'Number of Time Steps'\n"
558 " 'Integrator Step Type' = 'Constant'\n");
560 Teuchos::RCP<Teuchos::ParameterList> tscsPL = Teuchos::parameterList(
"Time Step Control Strategy");
561 tscsPL->set<std::string>(
"Time Step Control Strategy List",
"");
562 pl->set(
"Time Step Control Strategy", *tscsPL);
567 template <
class Scalar>
568 Teuchos::RCP<Teuchos::ParameterList>
575 template <
class Scalar>
576 Teuchos::RCP<Teuchos::ParameterList>
579 Teuchos::RCP<Teuchos::ParameterList> temp_plist = tscPL_;
580 tscPL_ = Teuchos::null;
586 #endif // Tempus_TimeStepControl_impl_hpp