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