9 #ifndef Tempus_IntegratorBasic_impl_hpp
10 #define Tempus_IntegratorBasic_impl_hpp
12 #include "Thyra_VectorStdOps.hpp"
15 #include "Tempus_StepperFactory.hpp"
16 #include "Tempus_StepperForwardEuler.hpp"
20 template <
class Scalar>
22 : outputScreenIndices_(std::vector<int>()),
23 outputScreenInterval_(1000000),
42 template <
class Scalar>
48 std::vector<int> outputScreenIndices,
int outputScreenInterval)
49 : outputScreenIndices_(outputScreenIndices),
50 outputScreenInterval_(outputScreenInterval),
67 template <
class Scalar>
70 this->setIntegratorType(iB->getIntegratorType());
71 this->setIntegratorName(iB->getIntegratorName());
72 this->setStepper(iB->getStepper());
73 this->setSolutionHistory(iB->getNonConstSolutionHistory());
74 this->setTimeStepControl(iB->getNonConstTimeStepControl());
75 this->setObserver(iB->getObserver());
76 this->setScreenOutputIndexList(iB->getScreenOutputIndexList());
77 this->setScreenOutputIndexInterval(iB->getScreenOutputIndexInterval());
78 this->setStatus(iB->getStatus());
79 integratorTimer_ = iB->getIntegratorTimer();
80 stepperTimer_ = iB->getStepperTimer();
83 template <
class Scalar>
87 i !=
"Integrator Basic", std::logic_error,
88 "Error - Integrator Type should be 'Integrator Basic'\n");
90 this->integratorType_ = i;
93 template <
class Scalar>
98 stepper_ == Teuchos::null, std::logic_error,
99 "Error - setModel(), need to set stepper first!\n");
101 stepper_->setModel(model);
104 template <
class Scalar>
107 if (stepper == Teuchos::null)
114 template <
class Scalar>
120 if (solutionHistory_ == Teuchos::null) {
124 solutionHistory_->clear();
128 stepper_ == Teuchos::null, std::logic_error,
129 "Error - initializeSolutionHistory(), need to set stepper first!\n");
131 if (state == Teuchos::null) {
134 "Error - initializeSolutionHistory(), need to "
135 "set stepper's model first!\n");
138 stepper_->getDefaultStepperState());
140 if (timeStepControl_ != Teuchos::null) {
142 state->setTime(timeStepControl_->getInitTime());
143 state->setIndex(timeStepControl_->getInitIndex());
144 state->setTimeStep(timeStepControl_->getInitTimeStep());
145 state->setTolRel(timeStepControl_->getMaxRelError());
146 state->setTolAbs(timeStepControl_->getMaxAbsError());
148 state->setOrder(stepper_->getOrder());
152 solutionHistory_->addState(state);
154 stepper_->setInitialConditions(solutionHistory_);
157 template <
class Scalar>
166 RCP<Thyra::VectorBase<Scalar> > xdot = x0->clone_v();
167 RCP<Thyra::VectorBase<Scalar> > xdotdot = x0->clone_v();
168 if (xdot0 == Teuchos::null)
171 Thyra::assign(xdot.ptr(), *(xdot0));
172 if (xdotdot0 == Teuchos::null)
175 Thyra::assign(xdotdot.ptr(), *(xdotdot0));
178 stepper_ == Teuchos::null, std::logic_error,
179 "Error - initializeSolutionHistory(), need to set stepper first!\n");
182 state->setStepperState(stepper_->getDefaultStepperState());
185 if (timeStepControl_ != Teuchos::null) {
187 state->setIndex(timeStepControl_->getInitIndex());
188 state->setTimeStep(timeStepControl_->getInitTimeStep());
189 state->setTolRel(timeStepControl_->getMaxRelError());
190 state->setTolAbs(timeStepControl_->getMaxAbsError());
192 state->setOrder(stepper_->getOrder());
195 initializeSolutionHistory(state);
198 template <
class Scalar>
202 if (sh == Teuchos::null) {
204 if (solutionHistory_ == Teuchos::null)
208 solutionHistory_ = sh;
212 template <
class Scalar>
216 if (tsc == Teuchos::null) {
218 if (timeStepControl_ == Teuchos::null) {
224 timeStepControl_ = tsc;
228 template <
class Scalar>
232 if (obs == Teuchos::null)
235 integratorObserver_ = obs;
238 template <
class Scalar>
242 stepper_ == Teuchos::null, std::logic_error,
243 "Error - Need to set the Stepper, setStepper(), before calling "
244 "IntegratorBasic::initialize()\n");
247 solutionHistory_->getNumStates() < 1, std::out_of_range,
248 "Error - SolutionHistory requires at least one SolutionState.\n"
249 <<
" Supplied SolutionHistory has only "
250 << solutionHistory_->getNumStates() <<
" SolutionStates.\n");
252 stepper_->initialize();
253 solutionHistory_->initialize();
254 timeStepControl_->initialize();
256 isInitialized_ =
true;
259 template <
class Scalar>
262 std::string name =
"Tempus::IntegratorBasic";
266 template <
class Scalar>
270 auto l_out = Teuchos::fancyOStream(out.
getOStream());
272 l_out->setOutputToRootOnly(0);
274 *l_out <<
"\n--- " << this->description() <<
" ---" << std::endl;
276 if (solutionHistory_ != Teuchos::null) {
277 solutionHistory_->describe(*l_out, verbLevel);
280 *l_out <<
"solutionHistory = " << solutionHistory_ << std::endl;
283 if (timeStepControl_ != Teuchos::null) {
284 timeStepControl_->describe(out, verbLevel);
287 *l_out <<
"timeStepControl = " << timeStepControl_ << std::endl;
290 if (stepper_ != Teuchos::null) {
291 stepper_->describe(out, verbLevel);
294 *l_out <<
"stepper = " << stepper_ << std::endl;
296 *l_out << std::string(this->description().length() + 8,
'-') << std::endl;
299 template <
class Scalar>
302 if (timeStepControl_->timeInRange(timeFinal))
303 timeStepControl_->setFinalTime(timeFinal);
304 bool itgrStatus = advanceTime();
308 template <
class Scalar>
313 if (isInitialized_ ==
false) {
315 *out <<
"Failure - IntegratorBasic is not initialized." << std::endl;
321 auto cs = solutionHistory_->getCurrentState();
322 cs->setTolRel(timeStepControl_->getMaxRelError());
323 cs->setTolAbs(timeStepControl_->getMaxAbsError());
325 integratorTimer_->start();
327 const Scalar initDt = std::min(timeStepControl_->getInitTimeStep(),
328 stepper_->getInitTimeStep(solutionHistory_));
330 timeStepControl_->setInitTimeStep(initDt);
331 timeStepControl_->initialize();
335 template <
class Scalar>
338 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::IntegratorBasic::advanceTime()");
341 integratorObserver_->observeStartIntegrator(*
this);
344 integratorStatus_ ==
WORKING &&
345 timeStepControl_->timeInRange(solutionHistory_->getCurrentTime()) &&
346 timeStepControl_->indexInRange(solutionHistory_->getCurrentIndex())) {
347 stepperTimer_->reset();
348 stepperTimer_->start();
349 solutionHistory_->initWorkingState();
352 integratorObserver_->observeStartTimeStep(*
this);
354 timeStepControl_->setNextTimeStep(solutionHistory_, integratorStatus_);
355 integratorObserver_->observeNextTimeStep(*
this);
358 solutionHistory_->getWorkingState()->setSolutionStatus(
WORKING);
360 integratorObserver_->observeBeforeTakeStep(*
this);
362 stepper_->takeStep(solutionHistory_);
364 integratorObserver_->observeAfterTakeStep(*
this);
366 stepperTimer_->stop();
368 integratorObserver_->observeAfterCheckTimeStep(*
this);
370 solutionHistory_->promoteWorkingState();
371 integratorObserver_->observeEndTimeStep(*
this);
375 integratorObserver_->observeEndIntegrator(*
this);
381 template <
class Scalar>
384 auto ws = solutionHistory_->getWorkingState();
387 ws->setTolRel(timeStepControl_->getMaxRelError());
388 ws->setTolAbs(timeStepControl_->getMaxAbsError());
391 std::vector<int>::const_iterator it = std::find(
392 outputScreenIndices_.begin(), outputScreenIndices_.end(), ws->getIndex());
393 if (it == outputScreenIndices_.end())
394 ws->setOutputScreen(
false);
396 ws->setOutputScreen(
true);
398 const int initial = timeStepControl_->getInitIndex();
399 if ((ws->getIndex() - initial) % outputScreenInterval_ == 0)
400 ws->setOutputScreen(
true);
403 template <
class Scalar>
407 auto ws = solutionHistory_->getWorkingState();
410 if (ws->getNFailures() >= timeStepControl_->getMaxFailures()) {
411 RCP<Teuchos::FancyOStream> out = this->getOStream();
412 out->setOutputToRootOnly(0);
414 *out <<
"Failure - Stepper has failed more than the maximum allowed.\n"
415 <<
" (nFailures = " << ws->getNFailures()
416 <<
") >= (nFailuresMax = " << timeStepControl_->getMaxFailures() <<
")"
421 if (ws->getNConsecutiveFailures() >=
422 timeStepControl_->getMaxConsecFailures()) {
423 RCP<Teuchos::FancyOStream> out = this->getOStream();
424 out->setOutputToRootOnly(0);
426 *out <<
"Failure - Stepper has failed more than the maximum "
427 <<
"consecutive allowed.\n"
428 <<
" (nConsecutiveFailures = " << ws->getNConsecutiveFailures()
429 <<
") >= (nConsecutiveFailuresMax = "
430 << timeStepControl_->getMaxConsecFailures() <<
")" << std::endl;
436 if (ws->getTimeStep() <= timeStepControl_->getMinTimeStep() &&
438 RCP<Teuchos::FancyOStream> out = this->getOStream();
439 out->setOutputToRootOnly(0);
441 *out <<
"Failure - Stepper has failed and the time step size is "
442 <<
"at the minimum.\n"
443 <<
" Solution Status = " <<
toString(ws->getSolutionStatus())
445 <<
" (TimeStep = " << ws->getTimeStep()
446 <<
") <= (Minimum TimeStep = " << timeStepControl_->getMinTimeStep()
455 ((timeStepControl_->getStepType() ==
"Constant") &&
456 !
approxEqual(ws->getTimeStep(), timeStepControl_->getInitTimeStep()))) {
457 RCP<Teuchos::FancyOStream> out = this->getOStream();
458 out->setOutputToRootOnly(0);
460 *out << std::scientific << std::setw(6) << std::setprecision(3)
461 << ws->getIndex() << std::setw(11) << std::setprecision(3)
462 << ws->getTime() << std::setw(11) << std::setprecision(3)
463 << ws->getTimeStep() <<
" STEP FAILURE!! - ";
465 *out <<
"Solution Status = " <<
toString(ws->getSolutionStatus())
468 else if ((timeStepControl_->getStepType() ==
"Constant") &&
469 (ws->getTimeStep() != timeStepControl_->getInitTimeStep())) {
470 *out <<
"dt != Constant dt (=" << timeStepControl_->getInitTimeStep()
474 ws->setNFailures(ws->getNFailures() + 1);
475 ws->setNRunningFailures(ws->getNRunningFailures() + 1);
476 ws->setNConsecutiveFailures(ws->getNConsecutiveFailures() + 1);
485 template <
class Scalar>
488 std::string exitStatus;
489 if (solutionHistory_->getCurrentState()->getSolutionStatus() ==
492 exitStatus =
"Time integration FAILURE!";
496 exitStatus =
"Time integration complete.";
499 integratorTimer_->stop();
502 template <
class Scalar>
506 std::string delimiters(
",");
507 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
508 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
509 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
510 std::string token = str.substr(lastPos, pos - lastPos);
511 outputScreenIndices_.push_back(
int(std::stoi(token)));
512 if (pos == std::string::npos)
break;
514 lastPos = str.find_first_not_of(delimiters, pos);
515 pos = str.find_first_of(delimiters, lastPos);
519 std::sort(outputScreenIndices_.begin(), outputScreenIndices_.end());
520 outputScreenIndices_.erase(
521 std::unique(outputScreenIndices_.begin(), outputScreenIndices_.end()),
522 outputScreenIndices_.end());
526 template <
class Scalar>
529 std::stringstream ss;
530 for (
size_t i = 0; i < outputScreenIndices_.size(); ++i) {
531 if (i != 0) ss <<
", ";
532 ss << outputScreenIndices_[i];
539 template <
class Scalar>
544 Teuchos::parameterList(getIntegratorName());
546 pl->
set(
"Integrator Type", getIntegratorType(),
547 "'Integrator Type' must be 'Integrator Basic'.");
549 pl->
set(
"Screen Output Index List", getScreenOutputIndexListString(),
550 "Screen Output Index List. Required to be in TimeStepControl range "
551 "['Minimum Time Step Index', 'Maximum Time Step Index']");
553 pl->
set(
"Screen Output Index Interval", getScreenOutputIndexInterval(),
554 "Screen Output Index Interval (e.g., every 100 time steps)");
556 pl->
set(
"Stepper Name", stepper_->getStepperName(),
557 "'Stepper Name' selects the Stepper block to construct (Required).");
559 pl->
set(
"Solution History", *solutionHistory_->getValidParameters());
560 pl->
set(
"Time Step Control", *timeStepControl_->getValidParameters());
563 Teuchos::parameterList(
"Tempus");
565 tempusPL->
set(
"Integrator Name", pl->
name());
566 tempusPL->
set(pl->
name(), *pl);
567 tempusPL->
set(stepper_->getStepperName(), *stepper_->getValidParameters());
574 template <
class Scalar>
579 if (tempusPL == Teuchos::null || tempusPL->
numParams() == 0)
582 auto integratorName = tempusPL->
get<std::string>(
"Integrator Name");
583 auto integratorPL = Teuchos::sublist(tempusPL, integratorName,
true);
585 std::string integratorType =
586 integratorPL->get<std::string>(
"Integrator Type");
588 integratorType !=
"Integrator Basic", std::logic_error,
589 "Error - For IntegratorBasic, 'Integrator Type' should be "
590 <<
"'Integrator Basic'.\n"
591 <<
" Integrator Type = " << integratorType <<
"\n");
593 integrator->setIntegratorName(integratorName);
597 integrator->getValidParameters());
598 auto vIntegratorName = validPL->template get<std::string>(
"Integrator Name");
599 auto vIntegratorPL = Teuchos::sublist(validPL, vIntegratorName,
true);
600 integratorPL->validateParametersAndSetDefaults(*vIntegratorPL, 1);
603 if (integratorPL->isParameter(
"Stepper Name")) {
605 auto stepperName = integratorPL->get<std::string>(
"Stepper Name");
606 auto stepperPL = Teuchos::sublist(tempusPL, stepperName,
true);
607 stepperPL->setName(stepperName);
609 integrator->setStepper(sf->createStepper(stepperPL));
614 integrator->setStepper(stepper);
618 if (integratorPL->isSublist(
"Time Step Control")) {
620 auto tscPL = Teuchos::sublist(integratorPL,
"Time Step Control",
true);
621 integrator->setTimeStepControl(
622 createTimeStepControl<Scalar>(tscPL, runInitialize));
630 if (integratorPL->isSublist(
"Solution History")) {
632 auto shPL = Teuchos::sublist(integratorPL,
"Solution History",
true);
633 auto sh = createSolutionHistoryPL<Scalar>(shPL);
634 integrator->setSolutionHistory(sh);
638 integrator->setSolutionHistory(createSolutionHistory<Scalar>());
642 integrator->setObserver(Teuchos::null);
645 integrator->setScreenOutputIndexInterval(
646 integratorPL->get<
int>(
"Screen Output Index Interval",
647 integrator->getScreenOutputIndexInterval()));
650 auto str = integratorPL->get<std::string>(
"Screen Output Index List",
"");
651 integrator->setScreenOutputIndexList(str);
658 template <
class Scalar>
664 auto integrator = createIntegratorBasic<Scalar>(tempusPL, runInitialize);
665 if (model == Teuchos::null)
return integrator;
668 integrator->setModel(constModel);
673 integrator->getStepper()->getDefaultStepperState());
674 newState->setTime(integrator->getTimeStepControl()->getInitTime());
675 newState->setIndex(integrator->getTimeStepControl()->getInitIndex());
676 newState->setTimeStep(integrator->getTimeStepControl()->getInitTimeStep());
677 newState->setTolRel(integrator->getTimeStepControl()->getMaxRelError());
678 newState->setTolAbs(integrator->getTimeStepControl()->getMaxAbsError());
679 newState->setOrder(integrator->getStepper()->getOrder());
683 auto sh = integrator->getNonConstSolutionHistory();
684 sh->addState(newState);
685 integrator->getStepper()->setInitialConditions(sh);
687 if (runInitialize) integrator->initialize();
693 template <
class Scalar>
696 std::string stepperType)
701 auto stepper = sf->createStepper(stepperType, model);
702 integrator->setStepper(stepper);
703 integrator->initializeSolutionHistory();
704 integrator->initialize();
710 template <
class Scalar>
717 template <
class Scalar>
723 auto integratorName = tempusPL->
get<std::string>(
"Integrator Name");
724 auto integratorPL = Teuchos::sublist(tempusPL, integratorName,
true);
726 std::string integratorType =
727 integratorPL->get<std::string>(
"Integrator Type");
729 integratorType !=
"Integrator Basic", std::logic_error,
730 "Error - For IntegratorBasic, 'Integrator Type' should be "
731 <<
"'Integrator Basic'.\n"
732 <<
" Integrator Type = " << integratorType <<
"\n");
735 integrator->setIntegratorName(integratorName);
737 TEUCHOS_TEST_FOR_EXCEPTION(
738 !integratorPL->isParameter(
"Stepper Name"), std::logic_error,
739 "Error - Need to set the 'Stepper Name' in 'Integrator Basic'.\n");
741 auto stepperName = integratorPL->get<std::string>(
"Stepper Name");
742 TEUCHOS_TEST_FOR_EXCEPTION(
743 stepperName ==
"Operator Split", std::logic_error,
744 "Error - 'Stepper Name' should be 'Operator Split'.\n");
747 auto stepperPL = Teuchos::sublist(tempusPL, stepperName,
true);
748 stepperPL->setName(stepperName);
750 integrator->setStepper(sf->createStepper(stepperPL, models));
753 if (integratorPL->isSublist(
"Time Step Control")) {
755 auto tscPL = Teuchos::sublist(integratorPL,
"Time Step Control",
true);
756 integrator->setTimeStepControl(
757 createTimeStepControl<Scalar>(tscPL, runInitialize));
767 integrator->getStepper()->getDefaultStepperState());
768 newState->setTime(integrator->getTimeStepControl()->getInitTime());
769 newState->setIndex(integrator->getTimeStepControl()->getInitIndex());
770 newState->setTimeStep(integrator->getTimeStepControl()->getInitTimeStep());
771 newState->setTolRel(integrator->getTimeStepControl()->getMaxRelError());
772 newState->setTolAbs(integrator->getTimeStepControl()->getMaxAbsError());
773 newState->setOrder(integrator->getStepper()->getOrder());
777 auto shPL = Teuchos::sublist(integratorPL,
"Solution History",
true);
778 auto sh = createSolutionHistoryPL<Scalar>(shPL);
779 sh->addState(newState);
780 integrator->getStepper()->setInitialConditions(sh);
781 integrator->setSolutionHistory(sh);
784 integrator->setObserver(Teuchos::null);
787 integrator->setScreenOutputIndexInterval(
788 integratorPL->get<
int>(
"Screen Output Index Interval",
789 integrator->getScreenOutputIndexInterval()));
792 auto str = integratorPL->get<std::string>(
"Screen Output Index List",
"");
793 integrator->setScreenOutputIndexList(str);
796 integrator->getValidParameters());
799 auto vIntegratorName = validPL->template get<std::string>(
"Integrator Name");
800 auto vIntegratorPL = Teuchos::sublist(validPL, vIntegratorName,
true);
801 integratorPL->validateParametersAndSetDefaults(*vIntegratorPL);
804 auto vStepperName = vIntegratorPL->template get<std::string>(
"Stepper Name");
805 auto vStepperPL = Teuchos::sublist(validPL, vStepperName,
true);
806 stepperPL->validateParametersAndSetDefaults(*vStepperPL);
808 integrator->initialize();
814 #endif // Tempus_IntegratorBasic_impl_hpp
virtual bool advanceTime()
Advance the solution to timeMax, and return true if successful.
Teuchos::RCP< SolutionState< Scalar > > createSolutionStateX(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdot=Teuchos::null, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdotdot=Teuchos::null)
Nonmember constructor from non-const solution vectors, x.
const std::string & name() const
virtual void setModel(Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > model)
Set the model on the stepper.
const std::string toString(const Status status)
Convert Status to string.
Teuchos::RCP< Teuchos::Time > integratorTimer_
T & get(const std::string &name, T def_value)
std::string description() const override
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
virtual void setTimeStepControl(Teuchos::RCP< TimeStepControl< Scalar > > tsc=Teuchos::null)
Set the TimeStepControl.
Forward Euler time stepper.
virtual void setObserver(Teuchos::RCP< IntegratorObserver< Scalar > > obs=Teuchos::null)
Set the Observer.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< Teuchos::Time > stepperTimer_
virtual std::string getScreenOutputIndexListString() const
IntegratorObserverBasic class for time integrators. This basic class has simple no-op functions...
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Create valid IntegratorBasic ParameterList.
Ordinal numParams() const
virtual void setStepper(Teuchos::RCP< Stepper< Scalar > > stepper)
Set the Stepper.
Thyra Base interface for time steppers.
virtual void checkTimeStep()
Check if time step has passed or failed.
virtual void setSolutionHistory(Teuchos::RCP< SolutionHistory< Scalar > > sh=Teuchos::null)
Set the SolutionHistory.
bool approxEqual(Scalar a, Scalar b, Scalar relTol=numericalTol< Scalar >())
Test if values are approximately equal within the relative tolerance.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
IntegratorObserver class for time integrators.
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
virtual void initializeSolutionHistory(Teuchos::RCP< SolutionState< Scalar > > state=Teuchos::null)
Set the initial state which has the initial conditions.
virtual void setScreenOutputIndexList(std::vector< int > indices)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
void setIntegratorName(std::string i)
Set the Integrator Name.
Teuchos::RCP< SolutionState< Scalar > > createSolutionStateME(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< StepperState< Scalar > > &stepperState=Teuchos::null, const Teuchos::RCP< PhysicsState< Scalar > > &physicsState=Teuchos::null)
Nonmember constructor from Thyra ModelEvaluator.
virtual void initialize()
Initializes the Integrator after set* function calls.
virtual void endIntegrator()
Perform tasks after end of integrator.
virtual void startIntegrator()
Perform tasks before start of integrator.
void setIntegratorType(std::string i)
Set the Integrator Type.
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
virtual void copy(Teuchos::RCP< IntegratorBasic< Scalar > > iB)
Copy (a shallow copy)
virtual void startTimeStep()
Start time step.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
Teuchos::RCP< IntegratorBasic< Scalar > > createIntegratorBasic(Teuchos::RCP< Teuchos::ParameterList > pList, bool runInitialize=true)
Nonmember constructor.
Solution state for integrators and steppers.