Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
24namespace Tempus {
25
26template<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{
47 this->initialize();
48}
49
50
51template<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,
71 Teuchos::RCP<TimeEventComposite<Scalar> > timeEvent,
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
124template<class Scalar>
126{
127 TEUCHOS_TEST_FOR_EXCEPTION(
128 (getInitTime() > getFinalTime() ), std::logic_error,
129 "Error - Inconsistent time range.\n"
130 " (timeMin = "<<getInitTime()<<") > (timeMax = "<<getFinalTime()<<")\n");
131
132 TEUCHOS_TEST_FOR_EXCEPTION(
133 (getMinTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
134 std::logic_error,
135 "Error - Negative minimum time step. dtMin = "<<getMinTimeStep()<<")\n");
136 TEUCHOS_TEST_FOR_EXCEPTION(
137 (getMaxTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
138 std::logic_error,
139 "Error - Negative maximum time step. dtMax = "<<getMaxTimeStep()<<")\n");
140 TEUCHOS_TEST_FOR_EXCEPTION(
141 (getMinTimeStep() > getMaxTimeStep() ), std::logic_error,
142 "Error - Inconsistent time step range.\n"
143 " (dtMin = "<<getMinTimeStep()<<") > (dtMax = "<<getMaxTimeStep()<<")\n");
144 TEUCHOS_TEST_FOR_EXCEPTION(
145 (getInitTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
146 std::logic_error,
147 "Error - Negative initial time step. dtInit = "<<getInitTimeStep()<<")\n");
148 TEUCHOS_TEST_FOR_EXCEPTION(
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
157 TEUCHOS_TEST_FOR_EXCEPTION(
158 (getInitIndex() > getFinalIndex() ), std::logic_error,
159 "Error - Inconsistent time index range.\n"
160 " (iStepMin = "<<getInitIndex()<<") > (iStepMax = "
161 <<getFinalIndex()<<")\n");
162
163 TEUCHOS_TEST_FOR_EXCEPTION(
164 (getMaxAbsError() < Teuchos::ScalarTraits<Scalar>::zero() ),
165 std::logic_error,
166 "Error - Negative maximum time step. errorMaxAbs = "
167 <<getMaxAbsError()<<")\n");
168 TEUCHOS_TEST_FOR_EXCEPTION(
169 (getMaxRelError() < Teuchos::ScalarTraits<Scalar>::zero() ),
170 std::logic_error,
171 "Error - Negative maximum time step. errorMaxRel = "
172 <<getMaxRelError()<<")\n");
173
174 TEUCHOS_TEST_FOR_EXCEPTION(
175 (stepControlStrategy_ == Teuchos::null), std::logic_error,
176 "Error - Strategy is unset!\n");
177
178 stepControlStrategy_->initialize();
179
180 TEUCHOS_TEST_FOR_EXCEPTION(
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
192template<class Scalar>
194printDtChanges(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
212template<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
223template<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.
385 TEUCHOS_TEST_FOR_EXCEPTION(
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") {
393 TEUCHOS_TEST_FOR_EXCEPTION(
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
410template<class Scalar>
411bool 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
430template<class Scalar>
431bool TimeStepControl<Scalar>::indexInRange(const int iStep) const
432{
433 bool iir = (getInitIndex() <= iStep && iStep < getFinalIndex());
434 return iir;
435}
436
437
438template<class Scalar>
440{
441 TEUCHOS_TEST_FOR_EXCEPTION(
442 (stepControlStrategy_ == Teuchos::null), std::logic_error,
443 "Error - getStepType() - Strategy is unset!\n");
444 return stepControlStrategy_->getStepType();
445}
446
447
448template<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
461template<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
473template<class Scalar>
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
485template<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
496template<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
507template<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
539template<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
555template<class Scalar>
556void TimeStepControl<Scalar>::setOutputIndices(std::vector<int> outputIndices)
557{
558 timeEvent_->add(Teuchos::rcp(new Tempus::TimeEventListIndex<Scalar>(
559 outputIndices, "Output Index List")));
560 isInitialized_ = false;
561}
562
563
564template<class Scalar>
565void 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
573template<class Scalar>
575{
576 timeEvent_->add(Teuchos::rcp(new TimeEventRangeIndex<Scalar>(
577 getInitIndex(), getFinalIndex(), i, "Output Index Interval")));
578 isInitialized_ = false;
579}
580
581
582template<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
591template<class Scalar>
593 Teuchos::RCP<TimeEventComposite<Scalar> > timeEvent)
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
612template<class Scalar>
614 Teuchos::RCP<TimeStepControlStrategy<Scalar> > tscs)
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
628template<class Scalar>
630{
631 std::string name = "Tempus::TimeStepControl";
632 return(name);
633}
634
635
636template<class Scalar>
638 Teuchos::FancyOStream &out,
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
692template<class Scalar>
693Teuchos::RCP<const Teuchos::ParameterList>
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// ------------------------------------------------------------------------
789template <class Scalar>
790Teuchos::RCP<TimeStepControl<Scalar> > createTimeStepControl(
791 Teuchos::RCP<Teuchos::ParameterList> const& pList, bool runInitialize)
792{
793 using Teuchos::RCP;
794 using Teuchos::rcp;
795 using Teuchos::ParameterList;
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
851 tec->add(rcp(new Tempus::TimeEventListIndex<Scalar>(
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
TimeEventListIndex specifies a list of index events.
TimeEventList specifies a list of time events.
virtual void setNextTimeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory, Status &integratorStatus)
Determine the time step size.
virtual std::vector< Scalar > getOutputTimes() const
virtual void setOutputTimeInterval(Scalar t)
virtual bool timeInRange(const Scalar time) const
Check if time is within minimum and maximum time.
virtual void setNumTimeSteps(int numTimeSteps)
virtual void setTimeEvents(Teuchos::RCP< TimeEventComposite< Scalar > > teb=Teuchos::null)
virtual bool getOutputExactly() const
Return if the output needs to exactly happen on output time.
virtual bool indexInRange(const int iStep) const
Check if time step index is within minimum and maximum index.
virtual void setOutputIndices(std::vector< int > v)
virtual std::string getStepType() const
virtual std::vector< int > getOutputIndices() const
virtual void setTimeStepControlStrategy(Teuchos::RCP< TimeStepControlStrategy< Scalar > > tscs=Teuchos::null)
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void printDtChanges(int istep, Scalar dt_old, Scalar dt_new, std::string reason) const
virtual void setOutputTimes(std::vector< Scalar > outputTimes)
virtual Scalar getOutputTimeInterval() const
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return ParameterList with current values.
Status
Status for the Integrator, the Stepper and the SolutionState.
Teuchos::RCP< TimeStepControl< Scalar > > createTimeStepControl(Teuchos::RCP< Teuchos::ParameterList > const &pList, bool runInitialize=true)
Nonmember constructor from ParameterList.