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"
21 template<
class Scalar>
23 : outputScreenIndices_(std::vector<int>()),
24 outputScreenInterval_(1000000),
44 template<
class Scalar>
50 std::vector<int> outputScreenIndices,
51 int outputScreenInterval)
52 : outputScreenIndices_(outputScreenIndices),
53 outputScreenInterval_(outputScreenInterval),
71 template<
class Scalar>
74 this->setIntegratorType (iB->getIntegratorType() );
75 this->setIntegratorName (iB->getIntegratorName() );
76 this->setStepper (iB->getStepper() );
77 this->setSolutionHistory (iB->getNonConstSolutionHistory() );
78 this->setTimeStepControl (iB->getNonConstTimeStepControl() );
79 this->setObserver (iB->getObserver() );
80 this->setScreenOutputIndexList (iB->getScreenOutputIndexList() );
81 this->setScreenOutputIndexInterval(iB->getScreenOutputIndexInterval());
82 this->setStatus (iB->getStatus() );
83 integratorTimer_ = iB->getIntegratorTimer();
84 stepperTimer_ = iB->getStepperTimer();
88 template<
class Scalar>
92 "Error - Integrator Type should be 'Integrator Basic'\n");
94 this->integratorType_ = i;
98 template<
class Scalar>
103 "Error - setModel(), need to set stepper first!\n");
105 stepper_->setModel(model);
109 template<
class Scalar>
113 if (stepper == Teuchos::null)
120 template<
class Scalar>
126 if (solutionHistory_ == Teuchos::null) {
129 solutionHistory_->clear();
133 "Error - initializeSolutionHistory(), need to set stepper first!\n");
135 if (state == Teuchos::null) {
138 "Error - initializeSolutionHistory(), need to set stepper's model first!\n");
141 stepper_->getDefaultStepperState());
143 if (timeStepControl_ != Teuchos::null) {
145 state->setTime (timeStepControl_->getInitTime());
146 state->setIndex (timeStepControl_->getInitIndex());
147 state->setTimeStep(timeStepControl_->getInitTimeStep());
148 state->setTolRel (timeStepControl_->getMaxRelError());
149 state->setTolAbs (timeStepControl_->getMaxAbsError());
151 state->setOrder (stepper_->getOrder());
155 solutionHistory_->addState(state);
157 stepper_->setInitialConditions(solutionHistory_);
161 template<
class Scalar>
171 RCP<Thyra::VectorBase<Scalar> > xdot = x0->clone_v();
172 RCP<Thyra::VectorBase<Scalar> > xdotdot = x0->clone_v();
173 if (xdot0 == Teuchos::null)
176 Thyra::assign(xdot.ptr(), *(xdot0));
177 if (xdotdot0 == Teuchos::null)
180 Thyra::assign(xdotdot.ptr(), *(xdotdot0));
183 "Error - initializeSolutionHistory(), need to set stepper first!\n");
186 state->setStepperState(stepper_->getDefaultStepperState());
189 if (timeStepControl_ != Teuchos::null) {
191 state->setIndex (timeStepControl_->getInitIndex());
192 state->setTimeStep(timeStepControl_->getInitTimeStep());
193 state->setTolRel (timeStepControl_->getMaxRelError());
194 state->setTolAbs (timeStepControl_->getMaxAbsError());
196 state->setOrder (stepper_->getOrder());
199 initializeSolutionHistory(state);
203 template<
class Scalar>
207 if (sh == Teuchos::null) {
209 if (solutionHistory_ == Teuchos::null)
212 solutionHistory_ = sh;
217 template<
class Scalar>
221 if (tsc == Teuchos::null) {
223 if (timeStepControl_ == Teuchos::null) {
228 timeStepControl_ = tsc;
233 template<
class Scalar>
237 if (obs == Teuchos::null)
240 integratorObserver_ = obs;
244 template<
class Scalar>
248 "Error - Need to set the Stepper, setStepper(), before calling "
249 "IntegratorBasic::initialize()\n");
253 "Error - SolutionHistory requires at least one SolutionState.\n"
254 <<
" Supplied SolutionHistory has only "
255 << solutionHistory_->getNumStates() <<
" SolutionStates.\n");
257 stepper_->initialize();
258 solutionHistory_->initialize();
259 timeStepControl_->initialize();
261 isInitialized_ =
true;
265 template<
class Scalar>
268 std::string name =
"Tempus::IntegratorBasic";
273 template<
class Scalar>
278 auto l_out = Teuchos::fancyOStream( out.
getOStream() );
280 l_out->setOutputToRootOnly(0);
282 *l_out <<
"\n--- " << this->description() <<
" ---" << std::endl;
284 if ( solutionHistory_ != Teuchos::null ) {
285 solutionHistory_->describe(*l_out,verbLevel);
287 *l_out <<
"solutionHistory = " << solutionHistory_ << std::endl;
290 if ( timeStepControl_ != Teuchos::null ) {
291 timeStepControl_->describe(out,verbLevel);
293 *l_out <<
"timeStepControl = " << timeStepControl_ << std::endl;
296 if ( stepper_ != Teuchos::null ) {
297 stepper_->describe(out,verbLevel);
299 *l_out <<
"stepper = " << stepper_ << std::endl;
301 *l_out << std::string(this->description().length()+8,
'-') <<std::endl;
305 template <
class Scalar>
308 if (timeStepControl_->timeInRange(timeFinal))
309 timeStepControl_->setFinalTime(timeFinal);
310 bool itgrStatus = advanceTime();
315 template <
class Scalar>
320 if (isInitialized_ ==
false) {
322 *out <<
"Failure - IntegratorBasic is not initialized." << std::endl;
328 auto cs = solutionHistory_->getCurrentState();
329 cs->setTolRel(timeStepControl_->getMaxRelError());
330 cs->setTolAbs(timeStepControl_->getMaxAbsError());
332 integratorTimer_->start();
334 const Scalar initDt =
335 std::min(timeStepControl_->getInitTimeStep(),
336 stepper_->getInitTimeStep(solutionHistory_));
338 timeStepControl_->setInitTimeStep(initDt);
339 timeStepControl_->initialize();
344 template <
class Scalar>
347 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::IntegratorBasic::advanceTime()");
350 integratorObserver_->observeStartIntegrator(*
this);
352 while (integratorStatus_ ==
WORKING &&
353 timeStepControl_->timeInRange (solutionHistory_->getCurrentTime()) &&
354 timeStepControl_->indexInRange(solutionHistory_->getCurrentIndex())){
356 stepperTimer_->reset();
357 stepperTimer_->start();
358 solutionHistory_->initWorkingState();
361 integratorObserver_->observeStartTimeStep(*
this);
363 timeStepControl_->setNextTimeStep(solutionHistory_, integratorStatus_);
364 integratorObserver_->observeNextTimeStep(*
this);
367 solutionHistory_->getWorkingState()->setSolutionStatus(
WORKING);
369 integratorObserver_->observeBeforeTakeStep(*
this);
371 stepper_->takeStep(solutionHistory_);
373 integratorObserver_->observeAfterTakeStep(*
this);
375 stepperTimer_->stop();
377 integratorObserver_->observeAfterCheckTimeStep(*
this);
379 solutionHistory_->promoteWorkingState();
380 integratorObserver_->observeEndTimeStep(*
this);
384 integratorObserver_->observeEndIntegrator(*
this);
391 template <
class Scalar>
394 auto ws = solutionHistory_->getWorkingState();
397 ws->setTolRel(timeStepControl_->getMaxRelError());
398 ws->setTolAbs(timeStepControl_->getMaxAbsError());
401 std::vector<int>::const_iterator it =
402 std::find(outputScreenIndices_.begin(),
403 outputScreenIndices_.end(),
405 if (it == outputScreenIndices_.end())
406 ws->setOutputScreen(
false);
408 ws->setOutputScreen(
true);
410 const int initial = timeStepControl_->getInitIndex();
411 if ( (ws->getIndex() - initial) % outputScreenInterval_ == 0)
412 ws->setOutputScreen(
true);
416 template <
class Scalar>
420 auto ws = solutionHistory_->getWorkingState();
423 if (ws->getNFailures() >= timeStepControl_->getMaxFailures()) {
424 RCP<Teuchos::FancyOStream> out = this->getOStream();
425 out->setOutputToRootOnly(0);
427 *out <<
"Failure - Stepper has failed more than the maximum allowed.\n"
428 <<
" (nFailures = "<<ws->getNFailures()<<
") >= (nFailuresMax = "
429 << timeStepControl_->getMaxFailures()<<
")" << std::endl;
433 if (ws->getNConsecutiveFailures()
434 >= timeStepControl_->getMaxConsecFailures()){
435 RCP<Teuchos::FancyOStream> out = this->getOStream();
436 out->setOutputToRootOnly(0);
438 *out <<
"Failure - Stepper has failed more than the maximum "
439 <<
"consecutive allowed.\n"
440 <<
" (nConsecutiveFailures = "<<ws->getNConsecutiveFailures()
441 <<
") >= (nConsecutiveFailuresMax = "
442 << timeStepControl_->getMaxConsecFailures()
449 if (ws->getTimeStep() <= timeStepControl_->getMinTimeStep() &&
451 RCP<Teuchos::FancyOStream> out = this->getOStream();
452 out->setOutputToRootOnly(0);
454 *out <<
"Failure - Stepper has failed and the time step size is "
455 <<
"at the minimum.\n"
456 <<
" Solution Status = " <<
toString(ws->getSolutionStatus())
458 <<
" (TimeStep = " << ws->getTimeStep()
459 <<
") <= (Minimum TimeStep = "
460 << timeStepControl_->getMinTimeStep()
469 ((timeStepControl_->getStepType() ==
"Constant") &&
470 !
approxEqual(ws->getTimeStep(), timeStepControl_->getInitTimeStep()))
473 RCP<Teuchos::FancyOStream> out = this->getOStream();
474 out->setOutputToRootOnly(0);
476 *out <<std::scientific
477 <<std::setw( 6)<<std::setprecision(3)<<ws->getIndex()
478 <<std::setw(11)<<std::setprecision(3)<<ws->getTime()
479 <<std::setw(11)<<std::setprecision(3)<<ws->getTimeStep()
480 <<
" STEP FAILURE!! - ";
482 *out <<
"Solution Status = " <<
toString(ws->getSolutionStatus())
484 }
else if ((timeStepControl_->getStepType() ==
"Constant") &&
485 (ws->getTimeStep() != timeStepControl_->getInitTimeStep())) {
486 *out <<
"dt != Constant dt (="<<timeStepControl_->getInitTimeStep()<<
")"
490 ws->setNFailures(ws->getNFailures()+1);
491 ws->setNRunningFailures(ws->getNRunningFailures()+1);
492 ws->setNConsecutiveFailures(ws->getNConsecutiveFailures()+1);
503 template <
class Scalar>
506 std::string exitStatus;
507 if (solutionHistory_->getCurrentState()->getSolutionStatus() ==
509 exitStatus =
"Time integration FAILURE!";
512 exitStatus =
"Time integration complete.";
515 integratorTimer_->stop();
519 template <
class Scalar>
523 std::string delimiters(
",");
524 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
525 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
526 while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
527 std::string token = str.substr(lastPos,pos-lastPos);
528 outputScreenIndices_.push_back(
int(std::stoi(token)));
529 if(pos==std::string::npos)
532 lastPos = str.find_first_not_of(delimiters, pos);
533 pos = str.find_first_of(delimiters, lastPos);
537 std::sort(outputScreenIndices_.begin(),outputScreenIndices_.end());
538 outputScreenIndices_.erase(std::unique(outputScreenIndices_.begin(),
539 outputScreenIndices_.end() ),
540 outputScreenIndices_.end() );
545 template <
class Scalar>
548 std::stringstream ss;
549 for(
size_t i = 0; i < outputScreenIndices_.size(); ++i) {
550 if(i != 0) ss <<
", ";
551 ss << outputScreenIndices_[i];
559 template<
class Scalar>
564 Teuchos::parameterList(getIntegratorName());
566 pl->
set(
"Integrator Type", getIntegratorType(),
567 "'Integrator Type' must be 'Integrator Basic'.");
569 pl->
set(
"Screen Output Index List", getScreenOutputIndexListString(),
570 "Screen Output Index List. Required to be in TimeStepControl range "
571 "['Minimum Time Step Index', 'Maximum Time Step Index']");
573 pl->
set(
"Screen Output Index Interval", getScreenOutputIndexInterval(),
574 "Screen Output Index Interval (e.g., every 100 time steps)");
576 pl->
set(
"Stepper Name", stepper_->getStepperName(),
577 "'Stepper Name' selects the Stepper block to construct (Required).");
579 pl->
set(
"Solution History", *solutionHistory_->getValidParameters());
580 pl->
set(
"Time Step Control", *timeStepControl_->getValidParameters());
584 Teuchos::parameterList(
"Tempus");
586 tempusPL->
set(
"Integrator Name", pl->
name());
587 tempusPL->
set(pl->
name(), *pl);
588 tempusPL->
set(stepper_->getStepperName(), *stepper_->getValidParameters());
596 template<
class Scalar>
601 if (tempusPL == Teuchos::null || tempusPL->
numParams() == 0)
604 auto integratorName = tempusPL->
get<std::string>(
"Integrator Name");
605 auto integratorPL = Teuchos::sublist(tempusPL, integratorName,
true);
607 std::string integratorType = integratorPL->get<std::string>(
"Integrator Type");
610 "Error - For IntegratorBasic, 'Integrator Type' should be "
611 <<
"'Integrator Basic'.\n"
612 <<
" Integrator Type = " << integratorType <<
"\n");
614 integrator->setIntegratorName(integratorName);
619 auto vIntegratorName = validPL->template get<std::string>(
"Integrator Name");
620 auto vIntegratorPL = Teuchos::sublist(validPL, vIntegratorName,
true);
621 integratorPL->validateParametersAndSetDefaults(*vIntegratorPL,1);
624 if (integratorPL->isParameter(
"Stepper Name")) {
626 auto stepperName = integratorPL->get<std::string>(
"Stepper Name");
627 auto stepperPL = Teuchos::sublist(tempusPL, stepperName,
true);
628 stepperPL->setName(stepperName);
630 integrator->setStepper(sf->createStepper(stepperPL));
634 integrator->setStepper(stepper);
638 if (integratorPL->isSublist(
"Time Step Control")) {
640 auto tscPL = Teuchos::sublist(integratorPL,
"Time Step Control",
true);
641 integrator->setTimeStepControl(createTimeStepControl<Scalar>(tscPL, runInitialize));
648 if (integratorPL->isSublist(
"Solution History")) {
650 auto shPL = Teuchos::sublist(integratorPL,
"Solution History",
true);
651 auto sh = createSolutionHistoryPL<Scalar>(shPL);
652 integrator->setSolutionHistory(sh);
655 integrator->setSolutionHistory(createSolutionHistory<Scalar>());
659 integrator->setObserver(Teuchos::null);
662 integrator->setScreenOutputIndexInterval(
663 integratorPL->get<
int>(
"Screen Output Index Interval",
664 integrator->getScreenOutputIndexInterval()));
667 auto str = integratorPL->get<std::string>(
"Screen Output Index List",
"");
668 integrator->setScreenOutputIndexList(str);
676 template<
class Scalar>
682 auto integrator = createIntegratorBasic<Scalar>(tempusPL, runInitialize);
683 if ( model == Teuchos::null )
return integrator;
686 integrator->setModel(constModel);
690 integrator->getStepper()->getDefaultStepperState());
691 newState->setTime (integrator->getTimeStepControl()->getInitTime());
692 newState->setIndex (integrator->getTimeStepControl()->getInitIndex());
693 newState->setTimeStep(integrator->getTimeStepControl()->getInitTimeStep());
694 newState->setTolRel (integrator->getTimeStepControl()->getMaxRelError());
695 newState->setTolAbs (integrator->getTimeStepControl()->getMaxAbsError());
696 newState->setOrder (integrator->getStepper()->getOrder());
700 auto sh = integrator->getNonConstSolutionHistory();
701 sh->addState(newState);
702 integrator->getStepper()->setInitialConditions(sh);
704 if (runInitialize) integrator->initialize();
711 template<
class Scalar>
714 std::string stepperType)
719 auto stepper = sf->createStepper(stepperType, model);
720 integrator->setStepper(stepper);
721 integrator->initializeSolutionHistory();
722 integrator->initialize();
729 template<
class Scalar>
737 template<
class Scalar>
743 auto integratorName = tempusPL->
get<std::string>(
"Integrator Name");
744 auto integratorPL = Teuchos::sublist(tempusPL, integratorName,
true);
746 std::string integratorType = integratorPL->get<std::string>(
"Integrator Type");
749 "Error - For IntegratorBasic, 'Integrator Type' should be "
750 <<
"'Integrator Basic'.\n"
751 <<
" Integrator Type = " << integratorType <<
"\n");
754 integrator->setIntegratorName(integratorName);
756 TEUCHOS_TEST_FOR_EXCEPTION( !integratorPL->isParameter(
"Stepper Name"),
758 "Error - Need to set the 'Stepper Name' in 'Integrator Basic'.\n");
760 auto stepperName = integratorPL->get<std::string>(
"Stepper Name");
761 TEUCHOS_TEST_FOR_EXCEPTION( stepperName ==
"Operator Split",
763 "Error - 'Stepper Name' should be 'Operator Split'.\n");
766 auto stepperPL = Teuchos::sublist(tempusPL, stepperName,
true);
767 stepperPL->setName(stepperName);
769 integrator->setStepper(sf->createStepper(stepperPL, models));
772 if (integratorPL->isSublist(
"Time Step Control")) {
774 auto tscPL = Teuchos::sublist(integratorPL,
"Time Step Control",
true);
775 integrator->setTimeStepControl(createTimeStepControl<Scalar>(tscPL, runInitialize));
783 integrator->getStepper()->getDefaultStepperState());
784 newState->setTime (integrator->getTimeStepControl()->getInitTime());
785 newState->setIndex (integrator->getTimeStepControl()->getInitIndex());
786 newState->setTimeStep(integrator->getTimeStepControl()->getInitTimeStep());
787 newState->setTolRel (integrator->getTimeStepControl()->getMaxRelError());
788 newState->setTolAbs (integrator->getTimeStepControl()->getMaxAbsError());
789 newState->setOrder (integrator->getStepper()->getOrder());
793 auto shPL = Teuchos::sublist(integratorPL,
"Solution History",
true);
794 auto sh = createSolutionHistoryPL<Scalar>(shPL);
795 sh->addState(newState);
796 integrator->getStepper()->setInitialConditions(sh);
797 integrator->setSolutionHistory(sh);
800 integrator->setObserver(Teuchos::null);
803 integrator->setScreenOutputIndexInterval(
804 integratorPL->get<
int>(
"Screen Output Index Interval",
805 integrator->getScreenOutputIndexInterval()));
808 auto str = integratorPL->get<std::string>(
"Screen Output Index List",
"");
809 integrator->setScreenOutputIndexList(str);
814 auto vIntegratorName = validPL->template get<std::string>(
"Integrator Name");
815 auto vIntegratorPL = Teuchos::sublist(validPL, vIntegratorName,
true);
816 integratorPL->validateParametersAndSetDefaults(*vIntegratorPL);
819 auto vStepperName = vIntegratorPL->template get<std::string>(
"Stepper Name");
820 auto vStepperPL = Teuchos::sublist(validPL, vStepperName,
true);
821 stepperPL->validateParametersAndSetDefaults(*vStepperPL);
823 integrator->initialize();
830 #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.
IntegratorBasic()
Default constructor (requires calls to setModel and setSolutionHistory for initial conditions before ...
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.