Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_TimeStepControl_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ****************************************************************************
3 // Tempus: Copyright (2017) Sandia Corporation
4 //
5 // Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6 // ****************************************************************************
7 // @HEADER
8 
9 #ifndef Tempus_TimeStepControl_impl_hpp
10 #define Tempus_TimeStepControl_impl_hpp
11 
12 #include "Teuchos_TimeMonitor.hpp"
13 
18 #include "Tempus_TimeEventRange.hpp"
19 #include "Tempus_TimeEventRangeIndex.hpp"
20 #include "Tempus_TimeEventList.hpp"
21 #include "Tempus_TimeEventListIndex.hpp"
22 
23 
24 namespace Tempus {
25 
26 template<class Scalar>
28  : isInitialized_ (false),
29  initTime_ (0.0),
30  finalTime_ (1.0e+99),
31  minTimeStep_ (0.0),
32  initTimeStep_ (1.0e+99),
33  maxTimeStep_ (1.0e+99),
34  initIndex_ (0),
35  finalIndex_ (1000000),
36  maxAbsError_ (1.0e-08),
37  maxRelError_ (1.0e-08),
38  maxFailures_ (10),
39  maxConsecFailures_ (5),
40  numTimeSteps_ (-1),
41  printDtChanges_ (true),
42  teAdjustedDt_ (false),
43  dtAfterTimeEvent_ (0.0)
44 {
46  setTimeEvents();
47  this->initialize();
48 }
49 
50 
51 template<class Scalar>
53  Scalar initTime,
54  Scalar finalTime,
55  Scalar minTimeStep,
56  Scalar initTimeStep,
57  Scalar maxTimeStep,
58  int initIndex,
59  int finalIndex,
60  Scalar maxAbsError,
61  Scalar maxRelError,
62  int maxFailures,
63  int maxConsecFailures,
64  int numTimeSteps,
65  bool printDtChanges,
66  bool outputExactly,
67  std::vector<int> outputIndices,
68  std::vector<Scalar> outputTimes,
69  int outputIndexInterval,
70  Scalar outputTimeInterval,
72  Teuchos::RCP<TimeStepControlStrategy<Scalar>> stepControlStrategy)
73  : isInitialized_ (false ),
74  initTime_ (initTime ),
75  finalTime_ (finalTime ),
76  minTimeStep_ (minTimeStep ),
77  initTimeStep_ (initTimeStep ),
78  maxTimeStep_ (maxTimeStep ),
79  initIndex_ (initIndex ),
80  finalIndex_ (finalIndex ),
81  maxAbsError_ (maxAbsError ),
82  maxRelError_ (maxRelError ),
83  maxFailures_ (maxFailures ),
84  maxConsecFailures_ (maxConsecFailures ),
85  printDtChanges_ (printDtChanges ),
86  teAdjustedDt_ (false ),
87  dtAfterTimeEvent_ (0.0 )
88 {
89  using Teuchos::rcp;
90 
91  setTimeStepControlStrategy(stepControlStrategy);
92  setNumTimeSteps(numTimeSteps);
93 
94  auto tec = rcp(new TimeEventComposite<Scalar>());
95  tec->setName("Time Step Control Events");
96  auto tes = timeEvent->getTimeEvents();
97  for(auto& e : tes) tec->add(e);
98 
99  // Add a range of times.
100  auto teRange = rcp(new TimeEventRange<Scalar>(
101  initTime, finalTime, outputTimeInterval, "Output Time Interval", outputExactly));
102  tec->add(teRange);
103 
104  // Add a list of times.
105  if (!outputTimes.empty())
106  tec->add(rcp(new TimeEventList<Scalar>(
107  outputTimes, "Output Time List", outputExactly)));
108 
109  // Add a range of indices.
110  tec->add(rcp(new TimeEventRangeIndex<Scalar>(
111  initIndex, finalIndex, outputIndexInterval, "Output Index Interval")));
112 
113  // Add a list of indices.
114  if (!outputTimes.empty())
115  tec->add(rcp(new TimeEventListIndex<Scalar>(
116  outputIndices, "Output Index List")));
117 
118  setTimeEvents(tec);
119 
120  this->initialize();
121 }
122 
123 
124 template<class Scalar>
126 {
128  (getInitTime() > getFinalTime() ), std::logic_error,
129  "Error - Inconsistent time range.\n"
130  " (timeMin = "<<getInitTime()<<") > (timeMax = "<<getFinalTime()<<")\n");
131 
133  (getMinTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
134  std::logic_error,
135  "Error - Negative minimum time step. dtMin = "<<getMinTimeStep()<<")\n");
137  (getMaxTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
138  std::logic_error,
139  "Error - Negative maximum time step. dtMax = "<<getMaxTimeStep()<<")\n");
141  (getMinTimeStep() > getMaxTimeStep() ), std::logic_error,
142  "Error - Inconsistent time step range.\n"
143  " (dtMin = "<<getMinTimeStep()<<") > (dtMax = "<<getMaxTimeStep()<<")\n");
145  (getInitTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
146  std::logic_error,
147  "Error - Negative initial time step. dtInit = "<<getInitTimeStep()<<")\n");
149  (getInitTimeStep() < getMinTimeStep() ||
150  getInitTimeStep() > getMaxTimeStep() ),
151  std::out_of_range,
152  "Error - Initial time step is out of range.\n"
153  << " [dtMin, dtMax] = [" << getMinTimeStep() << ", "
154  << getMaxTimeStep() << "]\n"
155  << " dtInit = " << getInitTimeStep() << "\n");
156 
158  (getInitIndex() > getFinalIndex() ), std::logic_error,
159  "Error - Inconsistent time index range.\n"
160  " (iStepMin = "<<getInitIndex()<<") > (iStepMax = "
161  <<getFinalIndex()<<")\n");
162 
164  (getMaxAbsError() < Teuchos::ScalarTraits<Scalar>::zero() ),
165  std::logic_error,
166  "Error - Negative maximum time step. errorMaxAbs = "
167  <<getMaxAbsError()<<")\n");
169  (getMaxRelError() < Teuchos::ScalarTraits<Scalar>::zero() ),
170  std::logic_error,
171  "Error - Negative maximum time step. errorMaxRel = "
172  <<getMaxRelError()<<")\n");
173 
175  (stepControlStrategy_ == Teuchos::null), std::logic_error,
176  "Error - Strategy is unset!\n");
177 
178  stepControlStrategy_->initialize();
179 
181  (getStepType() != "Constant" && getStepType() != "Variable"),
182  std::out_of_range,
183  "Error - 'Step Type' does not equal one of these:\n"
184  << " 'Constant' - Integrator will take constant time step sizes.\n"
185  << " 'Variable' - Integrator will allow changes to the time step size.\n"
186  << " stepType = " << getStepType() << "\n");
187 
188  isInitialized_ = true; // Only place where this is set to true!
189 }
190 
191 
192 template<class Scalar>
194 printDtChanges(int istep, Scalar dt_old, Scalar dt_new, std::string reason) const
195 {
196  if (!getPrintDtChanges()) return;
197 
198  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
199  out->setOutputToRootOnly(0);
200  Teuchos::OSTab ostab(out,0,"printDtChanges");
201 
202  std::stringstream message;
203  message << std::scientific
204  <<std::setw(6)<<std::setprecision(3)<<istep
205  << " * (dt = "<<std::setw(9)<<std::setprecision(3)<<dt_old
206  << ", new = "<<std::setw(9)<<std::setprecision(3)<<dt_new
207  << ") " << reason << std::endl;
208  *out << message.str();
209 }
210 
211 
212 template<class Scalar>
214 {
215  if ( !isInitialized_ ) {
216  this->describe( *(this->getOStream()), Teuchos::VERB_MEDIUM);
217  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
218  "Error - " << this->description() << " is not initialized!");
219  }
220 }
221 
222 
223 template<class Scalar>
225  const Teuchos::RCP<SolutionHistory<Scalar> > & solutionHistory,
226  Status & integratorStatus)
227 {
228  using Teuchos::RCP;
229 
230  checkInitialized();
231 
232  TEMPUS_FUNC_TIME_MONITOR("Tempus::TimeStepControl::setNextTimeStep()");
233  {
234  RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
235  const Scalar lastTime = solutionHistory->getCurrentState()->getTime();
236  const int iStep = workingState->getIndex();
237  Scalar dt = workingState->getTimeStep();
238  Scalar time = workingState->getTime();
239 
240  RCP<StepperState<Scalar> > stepperState = workingState->getStepperState();
241 
242  // If last time step was adjusted for time event, reinstate previous dt.
243  if (getStepType() == "Variable") {
244  if (teAdjustedDt_ == true) {
245  printDtChanges(iStep, dt, dtAfterTimeEvent_, "Reset dt after time event.");
246  dt = dtAfterTimeEvent_;
247  time = lastTime + dt;
248  teAdjustedDt_ = false;
249  dtAfterTimeEvent_ = 0.0;
250  }
251 
252  if (dt <= 0.0) {
253  printDtChanges(iStep, dt, getInitTimeStep(), "Reset dt to initial dt.");
254  dt = getInitTimeStep();
255  time = lastTime + dt;
256  }
257 
258  if (dt < getMinTimeStep()) {
259  printDtChanges(iStep, dt, getMinTimeStep(), "Reset dt to minimum dt.");
260  dt = getMinTimeStep();
261  time = lastTime + dt;
262  }
263  }
264 
265  // Update dt for the step control strategy to be informed
266  workingState->setTimeStep(dt);
267  workingState->setTime(time);
268 
269  // Call the step control strategy (to update dt if needed)
270  stepControlStrategy_->setNextTimeStep(*this, solutionHistory,
271  integratorStatus);
272 
273  // Get the dt (It was probably changed by stepControlStrategy_.)
274  dt = workingState->getTimeStep();
275  time = workingState->getTime();
276 
277  if (getStepType() == "Variable") {
278  if (dt < getMinTimeStep()) { // decreased below minimum dt
279  printDtChanges(iStep, dt, getMinTimeStep(),
280  "dt is too small. Resetting to minimum dt.");
281  dt = getMinTimeStep();
282  time = lastTime + dt;
283  }
284  if (dt > getMaxTimeStep()) { // increased above maximum dt
285  printDtChanges(iStep, dt, getMaxTimeStep(),
286  "dt is too large. Resetting to maximum dt.");
287  dt = getMaxTimeStep();
288  time = lastTime + dt;
289  }
290  }
291 
292 
293  // Check if this index is a TimeEvent and whether it is an output event.
294  bool outputDueToIndex = false;
295  std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > constrainingTEs;
296  if (timeEvent_->isIndex(iStep, constrainingTEs)) {
297  for(auto& e : constrainingTEs) {
298  if (e->getName() == "Output Index Interval" ||
299  e->getName() == "Output Index List" )
300  { outputDueToIndex = true; }
301  }
302  }
303 
304  // Check if during this time step is there a TimeEvent.
305  Scalar endTime = lastTime+dt;
306  // For "Variable", need to check t+dt+t_min in case we need to split
307  // the time step into two steps (e.g., time step lands within dt_min).
308  if (getStepType() == "Variable") endTime = lastTime+dt+getMinTimeStep();
309 
310  bool teThisStep = timeEvent_->eventInRange(lastTime, endTime, constrainingTEs);
311 
312  bool outputDueToTime = false;
313  Scalar tone = endTime;
314  bool landOnExactly = false;
315  if (teThisStep) {
316  for (auto& e : constrainingTEs) {
317  if (e->getName() == "Output Time Interval" ||
318  e->getName() == "Output Time List" )
319  { outputDueToTime = true; }
320 
321  if (e->getLandOnExactly() == true) {
322  landOnExactly = true;
323  tone = e->timeOfNextEvent(lastTime);
324  break;
325  }
326  }
327  }
328 
329  Scalar reltol = 1.0e-6;
330  if (teThisStep && getStepType() == "Variable") {
331  if (landOnExactly == true) {
332  // Adjust time step to hit TimeEvent.
333  if ( time > tone ) {
334  // Next TimeEvent is not near next time. It is more than
335  // getMinTimeStep() away from it.
336  printDtChanges(iStep, dt, tone - lastTime,
337  "Adjusting dt to hit the next TimeEvent.");
338  teAdjustedDt_ = true;
339  dtAfterTimeEvent_ = dt;
340  dt = tone - lastTime;
341  time = lastTime + dt;
342  } else if (std::fabs(time-tone) < reltol*std::fabs(time)) {
343  // Next TimeEvent IS VERY near next time. It is less than
344  // reltol away from it, e.g., adjust for numerical roundoff.
345  printDtChanges(iStep, dt, tone - lastTime,
346  "Adjusting dt for numerical roundoff to hit the next TimeEvent.");
347  teAdjustedDt_ = true;
348  dtAfterTimeEvent_ = dt;
349  dt = tone - lastTime;
350  time = lastTime + dt;
351  } else if (std::fabs(time + getMinTimeStep() - tone)
352  < reltol*std::fabs(tone)) {
353  // Next TimeEvent IS near next time. It is less than
354  // getMinTimeStep() away from it. Take two time steps
355  // to get to next TimeEvent.
356  printDtChanges(iStep, dt, (tone - lastTime)/2.0,
357  "The next TimeEvent is within the minimum dt of the next time. "
358  "Adjusting dt to take two steps.");
359  dt = (tone - lastTime)/2.0;
360  time = lastTime + dt;
361  // This time step is no longer an output step due the time constraint.
362  outputDueToTime = false;
363  }
364  }
365  }
366 
367  // Adjust time step to hit final time or correct for small
368  // numerical differences.
369  if (getStepType() == "Variable") {
370  if ( (lastTime+dt > getFinalTime()) ||
371  (std::fabs((lastTime+dt-getFinalTime()))
372  < reltol*std::fabs(lastTime+dt)) ) {
373  printDtChanges(iStep, dt, getFinalTime() - lastTime,
374  "Adjusting dt to hit final time.");
375  dt = getFinalTime() - lastTime;
376  time = lastTime + dt;
377  }
378  }
379 
380  // Check for negative time step.
381  TEUCHOS_TEST_FOR_EXCEPTION( dt <= Scalar(0.0), std::out_of_range,
382  "Error - Time step is not positive. dt = " << dt <<"\n");
383 
384  // Time step always needs to keep time within range.
386  (lastTime + dt < getInitTime()), std::out_of_range,
387  "Error - Time step does not move time INTO time range.\n"
388  " [timeMin, timeMax] = [" << getInitTime() << ", "
389  << getFinalTime() << "]\n"
390  " T + dt = " << lastTime <<" + "<< dt <<" = " << lastTime + dt <<"\n");
391 
392  if (getStepType() == "Variable") {
394  (lastTime + dt > getFinalTime()), std::out_of_range,
395  "Error - Time step move time OUT OF time range.\n"
396  " [timeMin, timeMax] = [" << getInitTime() << ", "
397  << getFinalTime() << "]\n"
398  " T + dt = " << lastTime <<" + "<< dt <<" = " << lastTime + dt <<"\n");
399  }
400 
401  workingState->setTimeStep(dt);
402  workingState->setTime(time);
403  workingState->setOutput(outputDueToIndex || outputDueToTime);
404  }
405  return;
406 }
407 
408 
410 template<class Scalar>
411 bool TimeStepControl<Scalar>::timeInRange(const Scalar time) const
412 {
413  // Get absolute tolerance 1.0e-(i+14), i.e., 14 digits of accuracy.
414  const int relTol = 14;
415  const int i =
416  (getInitTime() == 0) ? 0 : 1 + (int)std::floor(std::log10(std::fabs(getInitTime()) ) );
417  const Scalar absTolInit = std::pow(10, i-relTol);
418  const int j =
419  (getFinalTime() == 0) ? 0 : 1 + (int)std::floor(std::log10(std::fabs(getFinalTime()) ) );
420  const Scalar absTolFinal = std::pow(10, j-relTol);
421 
422  const bool test1 = getInitTime() - absTolInit <= time;
423  const bool test2 = time < getFinalTime() - absTolFinal;
424 
425  return (test1 && test2);
426 }
427 
428 
430 template<class Scalar>
431 bool TimeStepControl<Scalar>::indexInRange(const int iStep) const
432 {
433  bool iir = (getInitIndex() <= iStep && iStep < getFinalIndex());
434  return iir;
435 }
436 
437 
438 template<class Scalar>
440 {
442  (stepControlStrategy_ == Teuchos::null), std::logic_error,
443  "Error - getStepType() - Strategy is unset!\n");
444  return stepControlStrategy_->getStepType();
445 }
446 
447 
448 template<class Scalar>
450 {
451  auto event = timeEvent_->find("Output Time Interval");
452  if (event != Teuchos::null) return event->getLandOnExactly();
453 
454  event = timeEvent_->find("Output Time List");
455  if (event != Teuchos::null) return event->getLandOnExactly();
456 
457  return true;
458 }
459 
460 
461 template<class Scalar>
463 {
464  std::vector<int> indices;
465  if ( timeEvent_->find("Output Index List") == Teuchos::null)
466  return indices;
467  auto teli = Teuchos::rcp_dynamic_cast<TimeEventListIndex<Scalar> >(
468  timeEvent_->find("Output Index List"));
469  return teli->getIndexList();
470 }
471 
472 
473 template<class Scalar>
474 std::vector<Scalar> TimeStepControl<Scalar>::getOutputTimes() const
475 {
476  std::vector<Scalar> times;
477  if ( timeEvent_->find("Output Time List") == Teuchos::null)
478  return times;
479  auto tel = Teuchos::rcp_dynamic_cast<TimeEventList<Scalar> >(
480  timeEvent_->find("Output Time List"));
481  return tel->getTimeList();
482 }
483 
484 
485 template<class Scalar>
487 {
488  if ( timeEvent_->find("Output Index Interval") == Teuchos::null)
489  return 1000000;
490  auto teri = Teuchos::rcp_dynamic_cast<TimeEventRangeIndex<Scalar> >(
491  timeEvent_->find("Output Index Interval"));
492  return teri->getIndexStride();
493 }
494 
495 
496 template<class Scalar>
498 {
499  if ( timeEvent_->find("Output Time Interval") == Teuchos::null)
500  return 1.0e+99;
501  auto ter = Teuchos::rcp_dynamic_cast<TimeEventRange<Scalar> >(
502  timeEvent_->find("Output Time Interval"));
503  return ter->getTimeStride();
504 }
505 
506 
507 template<class Scalar>
509 {
510  TEUCHOS_TEST_FOR_EXCEPTION( (getStepType() != "Constant"), std::out_of_range,
511  "Error - Can only use setNumTimeSteps() when 'Step Type' == 'Constant'.\n");
512 
513  numTimeSteps_ = numTimeSteps;
514  if (numTimeSteps_ >= 0) {
515  setFinalIndex(getInitIndex() + numTimeSteps_);
516  Scalar initTimeStep;
517  if (numTimeSteps_ == 0)
518  initTimeStep = Scalar(0.0);
519  else
520  initTimeStep = (getFinalTime() - getInitTime())/numTimeSteps_;
521  setInitTimeStep(initTimeStep);
522  setMinTimeStep (initTimeStep);
523  setMaxTimeStep (initTimeStep);
524 
525  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
526  out->setOutputToRootOnly(0);
527  Teuchos::OSTab ostab(out,1,"setNumTimeSteps");
528  *out << "Warning - setNumTimeSteps() Setting 'Number of Time Steps' = " << getNumTimeSteps()
529  << " Set the following parameters: \n"
530  << " 'Final Time Index' = " << getFinalIndex() << "\n"
531  << " 'Initial Time Step' = " << getInitTimeStep() << "\n"
532  << " 'Step Type' = " << getStepType() << std::endl;
533 
534  isInitialized_ = false;
535  }
536 }
537 
538 
539 template<class Scalar>
541 {
542  auto event = timeEvent_->find("Output Time Interval");
543  if (event != Teuchos::null) {
544  auto ter = Teuchos::rcp_dynamic_cast<TimeEventRange<Scalar> >(event);
545  ter->setLandOnExactly(b);
546  }
547  event = timeEvent_->find("Output Time List");
548  if (event != Teuchos::null) {
549  auto tel = Teuchos::rcp_dynamic_cast<TimeEventList<Scalar> >(event);
550  tel->setLandOnExactly(b);
551  }
552 }
553 
554 
555 template<class Scalar>
556 void TimeStepControl<Scalar>::setOutputIndices(std::vector<int> outputIndices)
557 {
559  outputIndices, "Output Index List")));
560  isInitialized_ = false;
561 }
562 
563 
564 template<class Scalar>
565 void TimeStepControl<Scalar>::setOutputTimes(std::vector<Scalar> outputTimes)
566 {
567  timeEvent_->add(Teuchos::rcp(new Tempus::TimeEventList<Scalar>(
568  outputTimes, "Output Time List", getOutputExactly())));
569  isInitialized_ = false;
570 }
571 
572 
573 template<class Scalar>
575 {
576  timeEvent_->add(Teuchos::rcp(new TimeEventRangeIndex<Scalar>(
577  getInitIndex(), getFinalIndex(), i, "Output Index Interval")));
578  isInitialized_ = false;
579 }
580 
581 
582 template<class Scalar>
584 {
585  timeEvent_->add(Teuchos::rcp(new TimeEventRange<Scalar>(
586  getInitTime(), getFinalTime(), t, "Output Time Interval", getOutputExactly())));
587  isInitialized_ = false;
588 }
589 
590 
591 template<class Scalar>
594 {
595  using Teuchos::rcp;
596 
597  if ( timeEvent != Teuchos::null ) {
598  timeEvent_ = timeEvent;
599  } else {
600  auto tec = rcp(new TimeEventComposite<Scalar>());
601  tec->setName("Time Step Control Events");
602  tec->add(rcp(new TimeEventRangeIndex<Scalar>(
603  getInitIndex(), getFinalIndex(), 1000000, "Output Index Interval")));
604  tec->add(rcp(new TimeEventRange<Scalar>(
605  getInitTime(), getFinalTime(), 1.0e+99, "Output Time Interval", true)));
606  timeEvent_ = tec;
607  }
608  isInitialized_ = false;
609 }
610 
611 
612 template<class Scalar>
615 {
616  using Teuchos::rcp;
617 
618  if ( tscs != Teuchos::null ) {
619  stepControlStrategy_ = tscs;
620  } else {
621  stepControlStrategy_ =
622  rcp(new TimeStepControlStrategyConstant<Scalar>(getInitTimeStep()));
623  }
624  isInitialized_ = false;
625 }
626 
627 
628 template<class Scalar>
630 {
631  std::string name = "Tempus::TimeStepControl";
632  return(name);
633 }
634 
635 
636 template<class Scalar>
639  const Teuchos::EVerbosityLevel verbLevel) const
640 {
641  auto l_out = Teuchos::fancyOStream( out.getOStream() );
642  Teuchos::OSTab ostab(*l_out, 2, this->description());
643  l_out->setOutputToRootOnly(0);
644 
645  *l_out << "\n--- " << this->description() << " ---" <<std::endl;
646 
647  if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_MEDIUM)) {
648  std::vector<int> idx = getOutputIndices();
649  std::ostringstream listIdx;
650  if (!idx.empty()) {
651  for(std::size_t i = 0; i < idx.size()-1; ++i) listIdx << idx[i] << ", ";
652  listIdx << idx[idx.size()-1];
653  }
654 
655  std::vector<Scalar> times = getOutputTimes();
656  std::ostringstream listTimes;
657  if (!times.empty()) {
658  for(std::size_t i = 0; i < times.size()-1; ++i)
659  listTimes << times[i] << ", ";
660  listTimes << times[times.size()-1];
661  }
662 
663  *l_out << " stepType = " << getStepType() << std::endl
664  << " initTime = " << getInitTime() << std::endl
665  << " finalTime = " << getFinalTime() << std::endl
666  << " minTimeStep = " << getMinTimeStep() << std::endl
667  << " initTimeStep = " << getInitTimeStep() << std::endl
668  << " maxTimeStep = " << getMaxTimeStep() << std::endl
669  << " initIndex = " << getInitIndex() << std::endl
670  << " finalIndex = " << getFinalIndex() << std::endl
671  << " maxAbsError = " << getMaxAbsError() << std::endl
672  << " maxRelError = " << getMaxRelError() << std::endl
673  << " maxFailures = " << getMaxFailures() << std::endl
674  << " maxConsecFailures = " << getMaxConsecFailures() << std::endl
675  << " numTimeSteps = " << getNumTimeSteps() << std::endl
676  << " printDtChanges = " << getPrintDtChanges() << std::endl
677  << " outputExactly = " << getOutputExactly() << std::endl
678  << " outputIndices = " << listIdx.str() << std::endl
679  << " outputTimes = " << listTimes.str() << std::endl
680  << " outputIndexInterval= " << getOutputIndexInterval() << std::endl
681  << " outputTimeInterval = " << getOutputTimeInterval() << std::endl
682  << " outputAdjustedDt = " << teAdjustedDt_ << std::endl
683  << " dtAfterOutput = " << dtAfterTimeEvent_ << std::endl;
684  stepControlStrategy_->describe(*l_out, verbLevel);
685  timeEvent_->describe(*l_out, verbLevel);
686  }
687  *l_out << std::string(this->description().length()+8, '-') <<std::endl;
688 
689 }
690 
691 
692 template<class Scalar>
695 {
696  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList("Time Step Control");
697 
698  pl->set<double>("Initial Time" , getInitTime() , "Initial time");
699  pl->set<double>("Final Time" , getFinalTime() , "Final time");
700  pl->set<double>("Minimum Time Step" , getMinTimeStep() , "Minimum time step size");
701  pl->set<double>("Initial Time Step" , getInitTimeStep(), "Initial time step size");
702  pl->set<double>("Maximum Time Step" , getMaxTimeStep() , "Maximum time step size");
703  pl->set<int> ("Initial Time Index" , getInitIndex() , "Initial time index");
704  pl->set<int> ("Final Time Index" , getFinalIndex() , "Final time index");
705  pl->set<int> ("Number of Time Steps", getNumTimeSteps(),
706  "The number of constant time steps. The actual step size gets computed\n"
707  "on the fly given the size of the time domain. Overides and resets\n"
708  " 'Final Time Index' = 'Initial Time Index' + 'Number of Time Steps'\n"
709  " 'Initial Time Step' = "
710  "('Final Time' - 'Initial Time')/'Number of Time Steps'\n");
711  pl->set<double>("Maximum Absolute Error", getMaxAbsError() , "Maximum absolute error");
712  pl->set<double>("Maximum Relative Error", getMaxRelError() , "Maximum relative error");
713 
714  pl->set<bool> ("Print Time Step Changes", getPrintDtChanges(),
715  "Print timestep size when it changes");
716 
717  auto event = timeEvent_->find("Output Index Interval");
718  if (event != Teuchos::null) {
719  auto teri = Teuchos::rcp_dynamic_cast<TimeEventRangeIndex<Scalar> >(event);
720  pl->set<int> ("Output Index Interval", teri->getIndexStride(), "Output Index Interval");
721  }
722 
723  event = timeEvent_->find("Output Time Interval");
724  if (event != Teuchos::null) {
725  auto ter = Teuchos::rcp_dynamic_cast<TimeEventRange<Scalar> >(event);
726  pl->set<double>("Output Time Interval", ter->getTimeStride(), "Output time interval");
727  }
728 
729  pl->set<bool>("Output Exactly On Output Times", getOutputExactly(),
730  "This determines if the timestep size will be adjusted to exactly land\n"
731  "on the output times for 'Variable' timestepping (default=true).\n"
732  "When set to 'false' or for 'Constant' time stepping, the timestep\n"
733  "following the output time will be flagged for output.\n");
734 
735  {
736  event = timeEvent_->find("Output Index List");
737  std::ostringstream list;
738  if (event != Teuchos::null) {
739  auto teli = Teuchos::rcp_dynamic_cast<TimeEventListIndex<Scalar> >(event);
740  std::vector<int> idx = teli->getIndexList();
741  if (!idx.empty()) {
742  for(std::size_t i = 0; i < idx.size()-1; ++i) list << idx[i] << ", ";
743  list << idx[idx.size()-1];
744  }
745  }
746  pl->set<std::string>("Output Index List", list.str(),
747  "Comma deliminated list of output indices");
748  }
749 
750  {
751  event = timeEvent_->find("Output Time List");
752  std::ostringstream list;
753  if (event != Teuchos::null) {
754  auto teli = Teuchos::rcp_dynamic_cast<TimeEventList<Scalar> >(event);
755  std::vector<Scalar> times = teli->getTimeList();
756  if (!times.empty()) {
757  for(std::size_t i = 0; i < times.size()-1; ++i) list << times[i] << ", ";
758  list << times[times.size()-1];
759  }
760  }
761  pl->set<std::string>("Output Time List", list.str(),
762  "Comma deliminated list of output times");
763  }
764 
765  { // Do not duplicate the above "Output" events in "Time Step Control Events".
766  auto tecTmp = Teuchos::rcp(new TimeEventComposite<Scalar>());
767  tecTmp->setTimeEvents(timeEvent_->getTimeEvents());
768  tecTmp->remove("Output Index Interval");
769  tecTmp->remove("Output Time Interval");
770  tecTmp->remove("Output Index List");
771  tecTmp->remove("Output Time List");
772  if (tecTmp->getSize() > 0)
773  pl->set("Time Step Control Events", *tecTmp->getValidParameters());
774  }
775 
776  pl->set<int> ("Maximum Number of Stepper Failures", getMaxFailures(),
777  "Maximum number of Stepper failures");
778  pl->set<int> ("Maximum Number of Consecutive Stepper Failures", getMaxConsecFailures(),
779  "Maximum number of consecutive Stepper failures");
780 
781  pl->set("Time Step Control Strategy", *stepControlStrategy_->getValidParameters());
782 
783  return pl;
784 }
785 
786 
787 // Nonmember constructor - ParameterList
788 // ------------------------------------------------------------------------
789 template <class Scalar>
791  Teuchos::RCP<Teuchos::ParameterList> const& pList, bool runInitialize)
792 {
793  using Teuchos::RCP;
794  using Teuchos::rcp;
796 
797  auto tsc = Teuchos::rcp(new TimeStepControl<Scalar>());
798  if (pList == Teuchos::null || pList->numParams() == 0) return tsc;
799 
800  auto tscValidPL = Teuchos::rcp_const_cast<ParameterList>(tsc->getValidParameters());
801  // Handle optional "Time Step Control Events" sublist.
802  if (pList->isSublist("Time Step Control Events")) {
803  RCP<ParameterList> tsctePL =
804  Teuchos::sublist(pList, "Time Step Control Events", true);
805  tscValidPL->set("Time Step Control Events", *tsctePL);
806  }
807 
808  pList->validateParametersAndSetDefaults(*tscValidPL, 0);
809 
810  tsc->setInitTime( pList->get<double>("Initial Time"));
811  tsc->setFinalTime( pList->get<double>("Final Time"));
812  tsc->setMinTimeStep( pList->get<double>("Minimum Time Step"));
813  tsc->setInitTimeStep( pList->get<double>("Initial Time Step"));
814  tsc->setMaxTimeStep( pList->get<double>("Maximum Time Step"));
815  tsc->setInitIndex( pList->get<int> ("Initial Time Index"));
816  tsc->setFinalIndex( pList->get<int> ("Final Time Index"));
817  tsc->setMaxAbsError( pList->get<double>("Maximum Absolute Error"));
818  tsc->setMaxRelError( pList->get<double>("Maximum Relative Error"));
819  tsc->setMaxFailures( pList->get<int> ("Maximum Number of Stepper Failures"));
820  tsc->setMaxConsecFailures(pList->get<int> ("Maximum Number of Consecutive Stepper Failures"));
821  tsc->setPrintDtChanges( pList->get<bool> ("Print Time Step Changes"));
822  tsc->setNumTimeSteps( pList->get<int> ("Number of Time Steps"));
823 
824 
825  auto tec = rcp(new TimeEventComposite<Scalar>());
826  tec->setName("Time Step Control Events");
827 
828  { // Parse output indices, "Output Index List".
829  std::vector<int> outputIndices;
830  outputIndices.clear();
831  std::string str = pList->get<std::string>("Output Index List");
832  std::string delimiters(",");
833  std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
834  std::string::size_type pos = str.find_first_of(delimiters, lastPos);
835  while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
836  std::string token = str.substr(lastPos,pos-lastPos);
837  outputIndices.push_back(int(std::stoi(token)));
838  if(pos==std::string::npos) break;
839 
840  lastPos = str.find_first_not_of(delimiters, pos);
841  pos = str.find_first_of(delimiters, lastPos);
842  }
843 
844  if (!outputIndices.empty()) {
845  // order output indices
846  std::sort(outputIndices.begin(),outputIndices.end());
847  outputIndices.erase(std::unique(outputIndices.begin(),
848  outputIndices.end() ),
849  outputIndices.end() );
850 
852  outputIndices, "Output Index List")));
853  }
854  }
855 
856  // Parse output index internal
857  tec->add(rcp(new TimeEventRangeIndex<Scalar>(
858  tsc->getInitIndex(), tsc->getFinalIndex(),
859  pList->get<int>("Output Index Interval"),
860  "Output Index Interval")));
861 
862  auto outputExactly = pList->get<bool>("Output Exactly On Output Times",true);
863 
864  { // Parse output times, "Output Time List".
865  std::vector<Scalar> outputTimes;
866  outputTimes.clear();
867  std::string str = pList->get<std::string>("Output Time List");
868  std::string delimiters(",");
869  // Skip delimiters at the beginning
870  std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
871  // Find the first delimiter
872  std::string::size_type pos = str.find_first_of(delimiters, lastPos);
873  while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
874  // Found a token, add it to the vector
875  std::string token = str.substr(lastPos,pos-lastPos);
876  outputTimes.push_back(Scalar(std::stod(token)));
877  if(pos==std::string::npos) break;
878 
879  lastPos = str.find_first_not_of(delimiters, pos); // Skip delimiters
880  pos = str.find_first_of(delimiters, lastPos); // Find next delimiter
881  }
882 
883  if (!outputTimes.empty()) {
884  // order output times
885  std::sort(outputTimes.begin(),outputTimes.end());
886  outputTimes.erase(std::unique(outputTimes.begin(),
887  outputTimes.end() ),
888  outputTimes.end() );
889 
890  tec->add(rcp(new Tempus::TimeEventList<Scalar>(
891  outputTimes, "Output Time List", outputExactly)));
892  }
893  }
894 
895  // Parse output time interval
896  tec->add(rcp(new TimeEventRange<Scalar>(
897  tsc->getInitTime(), tsc->getFinalTime(),
898  pList->get<double>("Output Time Interval"),
899  "Output Time Interval", outputExactly)));
900 
901  if (pList->isSublist("Time Step Control Events")) {
902  RCP<ParameterList> tsctePL =
903  Teuchos::sublist(pList, "Time Step Control Events", true);
904 
905  auto timeEventType = tsctePL->get<std::string>("Type", "Unknown");
906  if (timeEventType == "Range") {
907  tec->add(createTimeEventRange<Scalar>(tsctePL));
908  } else if (timeEventType == "Range Index") {
909  tec->add(createTimeEventRangeIndex<Scalar>(tsctePL));
910  } else if (timeEventType == "List") {
911  tec->add(createTimeEventList<Scalar>(tsctePL));
912  } else if (timeEventType == "List Index") {
913  tec->add(createTimeEventListIndex<Scalar>(tsctePL));
914  } else if (timeEventType == "Composite") {
915  auto tecTmp = createTimeEventComposite<Scalar>(tsctePL);
916  auto timeEvents = tecTmp->getTimeEvents();
917  for(auto& e : timeEvents) tec->add(e);
918  } else {
919  RCP<Teuchos::FancyOStream> out =
920  Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
921  out->setOutputToRootOnly(0);
922  Teuchos::OSTab ostab(out,1, "createTimeStepControl()");
923  *out << "Warning -- Unknown Time Event Type!\n"
924  << "'Type' = '" << timeEventType << "'\n"
925  << "Should call add() with this "
926  << "(app-specific?) Time Event.\n" << std::endl;
927  }
928  }
929 
930  tsc->setTimeEvents(tec);
931 
932  if ( !pList->isParameter("Time Step Control Strategy") ) {
933 
934  tsc->setTimeStepControlStrategy(); // i.e., set default Constant timestep strategy.
935 
936  } else {
937 
938  RCP<ParameterList> tscsPL =
939  Teuchos::sublist(pList, "Time Step Control Strategy", true);
940 
941  auto strategyType = tscsPL->get<std::string>("Strategy Type");
942  if (strategyType == "Constant") {
943  tsc->setTimeStepControlStrategy(
944  createTimeStepControlStrategyConstant<Scalar>(tscsPL));
945  } else if (strategyType == "Basic VS") {
946  tsc->setTimeStepControlStrategy(
947  createTimeStepControlStrategyBasicVS<Scalar>(tscsPL));
948  } else if (strategyType == "Integral Controller") {
949  tsc->setTimeStepControlStrategy(
950  createTimeStepControlStrategyIntegralController<Scalar>(tscsPL));
951  } else if (strategyType == "Composite") {
952  tsc->setTimeStepControlStrategy(
953  createTimeStepControlStrategyComposite<Scalar>(tscsPL));
954  } else {
955  RCP<Teuchos::FancyOStream> out =
956  Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
957  out->setOutputToRootOnly(0);
958  Teuchos::OSTab ostab(out,1, "createTimeStepControl()");
959  *out << "Warning -- Did not find a Tempus strategy to create!\n"
960  << "'Strategy Type' = '" << strategyType << "'\n"
961  << "Should call setTimeStepControlStrategy() with this\n"
962  << "(app-specific?) strategy, and initialize().\n" << std::endl;
963  }
964  }
965 
966  if (runInitialize) tsc->initialize();
967  return tsc;
968 }
969 
970 
971 } // namespace Tempus
972 #endif // Tempus_TimeStepControl_impl_hpp
TimeEventRangeIndex specifies a start, stop and stride index.
Teuchos::RCP< TimeStepControl< Scalar > > createTimeStepControl(Teuchos::RCP< Teuchos::ParameterList > const &pList, bool runInitialize=true)
Nonmember constructor from ParameterList.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
virtual void setNumTimeSteps(int numTimeSteps)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TimeEventListIndex specifies a list of index events.
virtual void setLandOnExactly(bool LOE)
Set if the time event should be landed on exactly.
Ordinal numParams() const
virtual void printDtChanges(int istep, Scalar dt_old, Scalar dt_new, std::string reason) const
virtual void setTimeEvents(Teuchos::RCP< TimeEventComposite< Scalar > > teb=Teuchos::null)
virtual void setOutputTimes(std::vector< Scalar > outputTimes)
bool isParameter(const std::string &name) const
virtual std::string getStepType() const
virtual int getIndexStride() const
Return the stride of the index range.
TimeEventRange specifies a start, stop and stride time.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual void setTimeStepControlStrategy(Teuchos::RCP< TimeStepControlStrategy< Scalar > > tscs=Teuchos::null)
Status
Status for the Integrator, the Stepper and the SolutionState.
bool isSublist(const std::string &name) const
virtual Scalar getTimeStride() const
Return the stride of the time range.
This composite TimeEvent loops over added TimeEvents.
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 mechanisms that effect the time step size a...
TimeStepControl()
Default Constructor.
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
virtual std::vector< int > getIndexList() const
Return a vector of event indices.
virtual void setOutputTimeInterval(Scalar t)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
StepControlStrategy class for TimeStepControl.
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
virtual void setNextTimeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory, Status &integratorStatus)
Determine the time step size.
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
virtual int getOutputIndexInterval() const
virtual Scalar getOutputTimeInterval() const
virtual bool getOutputExactly() const
Return if the output needs to exactly happen on output time.
virtual void setLandOnExactly(bool LOE)
Set if the time event should be landed on exactly.
TimeEventList specifies a list of time events.
TimeStepControlStrategy class for TimeStepControl.
virtual void setOutputIndices(std::vector< int > v)
virtual std::vector< Scalar > getOutputTimes() const
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return ParameterList with current values.
virtual std::vector< int > getOutputIndices() const
virtual std::vector< Scalar > getTimeList() const
Return the list of time events.