Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_IntegratorBasic_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_IntegratorBasic_impl_hpp
10 #define Tempus_IntegratorBasic_impl_hpp
11 
12 #include "Thyra_VectorStdOps.hpp"
13 
15 #include "Tempus_StepperFactory.hpp"
16 #include "Tempus_StepperForwardEuler.hpp"
17 
18 
19 namespace Tempus {
20 
21 template<class Scalar>
23  : outputScreenIndices_(std::vector<int>()),
24  outputScreenInterval_(1000000),
25  integratorStatus_(WORKING),
26  isInitialized_(false)
27 {
28  setIntegratorName("Integrator Basic");
29  setIntegratorType("Integrator Basic");
30  setStepper(Teuchos::null);
31  setSolutionHistory(Teuchos::null);
32  setTimeStepControl(Teuchos::null);
33  setObserver(Teuchos::null);
34 
35  integratorTimer_ = rcp(new Teuchos::Time("Integrator Timer"));
36  stepperTimer_ = rcp(new Teuchos::Time("Stepper Timer"));
37 
38  // integrator is not initialized. Still requires calls to setModel
39  // and setSolutionHistory for initial conditions before calling
40  // initialize() to be fully constructed.
41 }
42 
43 
44 template<class Scalar>
46  Teuchos::RCP<Stepper<Scalar> > stepper,
47  Teuchos::RCP<SolutionHistory<Scalar> > solutionHistory,
48  Teuchos::RCP<TimeStepControl<Scalar> > timeStepControl,
49  Teuchos::RCP<IntegratorObserver<Scalar> > integratorObserver,
50  std::vector<int> outputScreenIndices,
51  int outputScreenInterval)
52  : outputScreenIndices_(outputScreenIndices),
53  outputScreenInterval_(outputScreenInterval),
54  integratorStatus_(WORKING),
55  isInitialized_(false)
56 {
57  setIntegratorName("Integrator Basic");
58  setIntegratorType("Integrator Basic");
59  setStepper(stepper);
60  setSolutionHistory(solutionHistory);
61  setTimeStepControl(timeStepControl);
62  setObserver(integratorObserver);
63 
64  integratorTimer_ = rcp(new Teuchos::Time("Integrator Timer"));
65  stepperTimer_ = rcp(new Teuchos::Time("Stepper Timer"));
66 
67  initialize();
68 }
69 
70 
71 template<class Scalar>
73 {
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();
85 }
86 
87 
88 template<class Scalar>
90 {
91  TEUCHOS_TEST_FOR_EXCEPTION( i != "Integrator Basic", std::logic_error,
92  "Error - Integrator Type should be 'Integrator Basic'\n");
93 
94  this->integratorType_ = i;
95 }
96 
97 
98 template<class Scalar>
101 {
102  TEUCHOS_TEST_FOR_EXCEPTION( stepper_ == Teuchos::null, std::logic_error,
103  "Error - setModel(), need to set stepper first!\n");
104 
105  stepper_->setModel(model);
106 }
107 
108 
109 template<class Scalar>
111  Teuchos::RCP<Stepper<Scalar> > stepper)
112 {
113  if (stepper == Teuchos::null)
114  stepper_ = Teuchos::rcp(new StepperForwardEuler<Scalar>());
115  else
116  stepper_ = stepper;
117 }
118 
120 template<class Scalar>
123 {
124  using Teuchos::RCP;
125 
126  if (solutionHistory_ == Teuchos::null) {
127  solutionHistory_ = rcp(new SolutionHistory<Scalar>());
128  } else {
129  solutionHistory_->clear();
130  }
131 
132  TEUCHOS_TEST_FOR_EXCEPTION( stepper_ == Teuchos::null, std::logic_error,
133  "Error - initializeSolutionHistory(), need to set stepper first!\n");
134 
135  if (state == Teuchos::null) {
136  TEUCHOS_TEST_FOR_EXCEPTION( stepper_->getModel() == Teuchos::null,
137  std::logic_error,
138  "Error - initializeSolutionHistory(), need to set stepper's model first!\n");
139  // Construct default IC from the application model
140  state = createSolutionStateME( stepper_->getModel(),
141  stepper_->getDefaultStepperState());
142 
143  if (timeStepControl_ != Teuchos::null) {
144  // Set SolutionState from TimeStepControl
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());
150  }
151  state->setOrder (stepper_->getOrder());
152  state->setSolutionStatus(Status::PASSED); // ICs are considered passing.
153  }
154 
155  solutionHistory_->addState(state);
156 
157  stepper_->setInitialConditions(solutionHistory_);
158 }
159 
160 
161 template<class Scalar>
166  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotdot0)
167 {
168  using Teuchos::RCP;
169 
170  // Create and set xdot and xdotdot.
171  RCP<Thyra::VectorBase<Scalar> > xdot = x0->clone_v();
172  RCP<Thyra::VectorBase<Scalar> > xdotdot = x0->clone_v();
173  if (xdot0 == Teuchos::null)
174  Thyra::assign(xdot.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
175  else
176  Thyra::assign(xdot.ptr(), *(xdot0));
177  if (xdotdot0 == Teuchos::null)
178  Thyra::assign(xdotdot.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
179  else
180  Thyra::assign(xdotdot.ptr(), *(xdotdot0));
181 
182  TEUCHOS_TEST_FOR_EXCEPTION( stepper_ == Teuchos::null, std::logic_error,
183  "Error - initializeSolutionHistory(), need to set stepper first!\n");
184 
185  auto state = createSolutionStateX(x0->clone_v(), xdot, xdotdot);
186  state->setStepperState(stepper_->getDefaultStepperState());
187 
188  state->setTime (t0);
189  if (timeStepControl_ != Teuchos::null) {
190  // Set SolutionState from TimeStepControl
191  state->setIndex (timeStepControl_->getInitIndex());
192  state->setTimeStep(timeStepControl_->getInitTimeStep());
193  state->setTolRel (timeStepControl_->getMaxRelError());
194  state->setTolAbs (timeStepControl_->getMaxAbsError());
195  }
196  state->setOrder (stepper_->getOrder());
197  state->setSolutionStatus(Status::PASSED); // ICs are considered passing.
198 
199  initializeSolutionHistory(state);
200 }
201 
202 
203 template<class Scalar>
206 {
207  if (sh == Teuchos::null) {
208  // Create default SolutionHistory, otherwise keep current history.
209  if (solutionHistory_ == Teuchos::null)
210  solutionHistory_ = rcp(new SolutionHistory<Scalar>());
211  } else {
212  solutionHistory_ = sh;
213  }
214 }
215 
216 
217 template<class Scalar>
220 {
221  if (tsc == Teuchos::null) {
222  // Create timeStepControl_ if null, otherwise keep current parameters.
223  if (timeStepControl_ == Teuchos::null) {
224  // Construct default TimeStepControl
225  timeStepControl_ = rcp(new TimeStepControl<Scalar>());
226  }
227  } else {
228  timeStepControl_ = tsc;
229  }
230 }
231 
232 
233 template<class Scalar>
236 {
237  if (obs == Teuchos::null)
238  integratorObserver_ = Teuchos::rcp(new IntegratorObserverBasic<Scalar>());
239  else
240  integratorObserver_ = obs;
241 }
242 
243 
244 template<class Scalar>
246 {
247  TEUCHOS_TEST_FOR_EXCEPTION( stepper_ == Teuchos::null, std::logic_error,
248  "Error - Need to set the Stepper, setStepper(), before calling "
249  "IntegratorBasic::initialize()\n");
250 
251  TEUCHOS_TEST_FOR_EXCEPTION( solutionHistory_->getNumStates() < 1,
252  std::out_of_range,
253  "Error - SolutionHistory requires at least one SolutionState.\n"
254  << " Supplied SolutionHistory has only "
255  << solutionHistory_->getNumStates() << " SolutionStates.\n");
256 
257  stepper_->initialize();
258  solutionHistory_->initialize();
259  timeStepControl_->initialize();
260 
261  isInitialized_ = true;
262 }
263 
264 
265 template<class Scalar>
267 {
268  std::string name = "Tempus::IntegratorBasic";
269  return(name);
270 }
271 
272 
273 template<class Scalar>
276  const Teuchos::EVerbosityLevel verbLevel) const
277 {
278  auto l_out = Teuchos::fancyOStream( out.getOStream() );
279  Teuchos::OSTab ostab(*l_out, 2, this->description());
280  l_out->setOutputToRootOnly(0);
281 
282  *l_out << "\n--- " << this->description() << " ---" << std::endl;
283 
284  if ( solutionHistory_ != Teuchos::null ) {
285  solutionHistory_->describe(*l_out,verbLevel);
286  } else {
287  *l_out << "solutionHistory = " << solutionHistory_ << std::endl;
288  }
289 
290  if ( timeStepControl_ != Teuchos::null ) {
291  timeStepControl_->describe(out,verbLevel);
292  } else {
293  *l_out << "timeStepControl = " << timeStepControl_ << std::endl;
294  }
295 
296  if ( stepper_ != Teuchos::null ) {
297  stepper_->describe(out,verbLevel);
298  } else {
299  *l_out << "stepper = " << stepper_ << std::endl;
300  }
301  *l_out << std::string(this->description().length()+8, '-') <<std::endl;
302 }
303 
304 
305 template <class Scalar>
306 bool IntegratorBasic<Scalar>::advanceTime(const Scalar timeFinal)
307 {
308  if (timeStepControl_->timeInRange(timeFinal))
309  timeStepControl_->setFinalTime(timeFinal);
310  bool itgrStatus = advanceTime();
311  return itgrStatus;
312 }
313 
314 
315 template <class Scalar>
317 {
318  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
319  out->setOutputToRootOnly(0);
320  if (isInitialized_ == false) {
321  Teuchos::OSTab ostab(out,1,"StartIntegrator");
322  *out << "Failure - IntegratorBasic is not initialized." << std::endl;
323  setStatus(Status::FAILED);
324  return;
325  }
326 
327  //set the Abs/Rel tolerance
328  auto cs = solutionHistory_->getCurrentState();
329  cs->setTolRel(timeStepControl_->getMaxRelError());
330  cs->setTolAbs(timeStepControl_->getMaxAbsError());
331 
332  integratorTimer_->start();
333  // get optimal initial time step
334  const Scalar initDt =
335  std::min(timeStepControl_->getInitTimeStep(),
336  stepper_->getInitTimeStep(solutionHistory_));
337  // update initial time step
338  timeStepControl_->setInitTimeStep(initDt);
339  timeStepControl_->initialize();
340  setStatus(WORKING);
341 }
342 
343 
344 template <class Scalar>
346 {
347  TEMPUS_FUNC_TIME_MONITOR("Tempus::IntegratorBasic::advanceTime()");
348  {
349  startIntegrator();
350  integratorObserver_->observeStartIntegrator(*this);
351 
352  while (integratorStatus_ == WORKING &&
353  timeStepControl_->timeInRange (solutionHistory_->getCurrentTime()) &&
354  timeStepControl_->indexInRange(solutionHistory_->getCurrentIndex())){
355 
356  stepperTimer_->reset();
357  stepperTimer_->start();
358  solutionHistory_->initWorkingState();
359 
360  startTimeStep();
361  integratorObserver_->observeStartTimeStep(*this);
362 
363  timeStepControl_->setNextTimeStep(solutionHistory_, integratorStatus_);
364  integratorObserver_->observeNextTimeStep(*this);
365 
366  if (integratorStatus_ == Status::FAILED) break;
367  solutionHistory_->getWorkingState()->setSolutionStatus(WORKING);
368 
369  integratorObserver_->observeBeforeTakeStep(*this);
370 
371  stepper_->takeStep(solutionHistory_);
372 
373  integratorObserver_->observeAfterTakeStep(*this);
374 
375  stepperTimer_->stop();
376  checkTimeStep();
377  integratorObserver_->observeAfterCheckTimeStep(*this);
378 
379  solutionHistory_->promoteWorkingState();
380  integratorObserver_->observeEndTimeStep(*this);
381  }
382 
383  endIntegrator();
384  integratorObserver_->observeEndIntegrator(*this);
385  }
386 
387  return (integratorStatus_ == Status::PASSED);
388 }
389 
390 
391 template <class Scalar>
393 {
394  auto ws = solutionHistory_->getWorkingState();
395 
396  //set the Abs/Rel tolerance
397  ws->setTolRel(timeStepControl_->getMaxRelError());
398  ws->setTolAbs(timeStepControl_->getMaxAbsError());
399 
400  // Check if we need to dump screen output this step
401  std::vector<int>::const_iterator it =
402  std::find(outputScreenIndices_.begin(),
403  outputScreenIndices_.end(),
404  ws->getIndex());
405  if (it == outputScreenIndices_.end())
406  ws->setOutputScreen(false);
407  else
408  ws->setOutputScreen(true);
409 
410  const int initial = timeStepControl_->getInitIndex();
411  if ( (ws->getIndex() - initial) % outputScreenInterval_ == 0)
412  ws->setOutputScreen(true);
413 }
414 
415 
416 template <class Scalar>
418 {
419  using Teuchos::RCP;
420  auto ws = solutionHistory_->getWorkingState();
421 
422  // Too many TimeStep failures, Integrator fails.
423  if (ws->getNFailures() >= timeStepControl_->getMaxFailures()) {
424  RCP<Teuchos::FancyOStream> out = this->getOStream();
425  out->setOutputToRootOnly(0);
426  Teuchos::OSTab ostab(out, 2, "checkTimeStep");
427  *out << "Failure - Stepper has failed more than the maximum allowed.\n"
428  << " (nFailures = "<<ws->getNFailures()<< ") >= (nFailuresMax = "
429  << timeStepControl_->getMaxFailures()<<")" << std::endl;
430  setStatus(Status::FAILED);
431  return;
432  }
433  if (ws->getNConsecutiveFailures()
434  >= timeStepControl_->getMaxConsecFailures()){
435  RCP<Teuchos::FancyOStream> out = this->getOStream();
436  out->setOutputToRootOnly(0);
437  Teuchos::OSTab ostab(out, 1, "checkTimeStep");
438  *out << "Failure - Stepper has failed more than the maximum "
439  << "consecutive allowed.\n"
440  << " (nConsecutiveFailures = "<<ws->getNConsecutiveFailures()
441  << ") >= (nConsecutiveFailuresMax = "
442  << timeStepControl_->getMaxConsecFailures()
443  << ")" << std::endl;
444  setStatus(Status::FAILED);
445  return;
446  }
447 
448  // Timestep size is at the minimum timestep size and the step failed.
449  if (ws->getTimeStep() <= timeStepControl_->getMinTimeStep() &&
450  ws->getSolutionStatus() == Status::FAILED) {
451  RCP<Teuchos::FancyOStream> out = this->getOStream();
452  out->setOutputToRootOnly(0);
453  Teuchos::OSTab ostab(out, 1, "checkTimeStep");
454  *out << "Failure - Stepper has failed and the time step size is "
455  << "at the minimum.\n"
456  << " Solution Status = " << toString(ws->getSolutionStatus())
457  << std::endl
458  << " (TimeStep = " << ws->getTimeStep()
459  << ") <= (Minimum TimeStep = "
460  << timeStepControl_->getMinTimeStep()
461  << ")" << std::endl;
462  setStatus(Status::FAILED);
463  return;
464  }
465 
466  // Check Stepper failure.
467  if (ws->getSolutionStatus() == Status::FAILED ||
468  // Constant time step failure
469  ((timeStepControl_->getStepType() == "Constant") &&
470  !approxEqual(ws->getTimeStep(), timeStepControl_->getInitTimeStep()))
471  )
472  {
473  RCP<Teuchos::FancyOStream> out = this->getOStream();
474  out->setOutputToRootOnly(0);
475  Teuchos::OSTab ostab(out, 0, "checkTimeStep");
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!! - ";
481  if (ws->getSolutionStatus() == Status::FAILED) {
482  *out << "Solution Status = " << toString(ws->getSolutionStatus())
483  << std::endl;
484  } else if ((timeStepControl_->getStepType() == "Constant") &&
485  (ws->getTimeStep() != timeStepControl_->getInitTimeStep())) {
486  *out << "dt != Constant dt (="<<timeStepControl_->getInitTimeStep()<<")"
487  << std::endl;
488  }
489 
490  ws->setNFailures(ws->getNFailures()+1);
491  ws->setNRunningFailures(ws->getNRunningFailures()+1);
492  ws->setNConsecutiveFailures(ws->getNConsecutiveFailures()+1);
493  ws->setSolutionStatus(Status::FAILED);
494  return;
495  }
496 
497  // TIME STEP PASSED basic tests! Ensure it is set as such.
498  ws->setSolutionStatus(Status::PASSED);
499 
500 }
501 
502 
503 template <class Scalar>
505 {
506  std::string exitStatus;
507  if (solutionHistory_->getCurrentState()->getSolutionStatus() ==
508  Status::FAILED || integratorStatus_ == Status::FAILED) {
509  exitStatus = "Time integration FAILURE!";
510  } else {
511  setStatus(Status::PASSED);
512  exitStatus = "Time integration complete.";
513  }
514 
515  integratorTimer_->stop();
516 }
517 
518 
519 template <class Scalar>
521 {
522  // Parse output indices
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)
530  break;
531 
532  lastPos = str.find_first_not_of(delimiters, pos);
533  pos = str.find_first_of(delimiters, lastPos);
534  }
535 
536  // order output indices and remove duplicates.
537  std::sort(outputScreenIndices_.begin(),outputScreenIndices_.end());
538  outputScreenIndices_.erase(std::unique(outputScreenIndices_.begin(),
539  outputScreenIndices_.end() ),
540  outputScreenIndices_.end() );
541  return;
542 }
543 
544 
545 template <class Scalar>
547 {
548  std::stringstream ss;
549  for(size_t i = 0; i < outputScreenIndices_.size(); ++i) {
550  if(i != 0) ss << ", ";
551  ss << outputScreenIndices_[i];
552  }
553  return ss.str();
554 }
555 
556 
559 template<class Scalar>
562 {
564  Teuchos::parameterList(getIntegratorName());
565 
566  pl->set("Integrator Type", getIntegratorType(),
567  "'Integrator Type' must be 'Integrator Basic'.");
568 
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']");
572 
573  pl->set("Screen Output Index Interval", getScreenOutputIndexInterval(),
574  "Screen Output Index Interval (e.g., every 100 time steps)");
575 
576  pl->set("Stepper Name", stepper_->getStepperName(),
577  "'Stepper Name' selects the Stepper block to construct (Required).");
578 
579  pl->set("Solution History", *solutionHistory_->getValidParameters());
580  pl->set("Time Step Control", *timeStepControl_->getValidParameters());
581 
582 
584  Teuchos::parameterList("Tempus");
585 
586  tempusPL->set("Integrator Name", pl->name());
587  tempusPL->set(pl->name(), *pl);
588  tempusPL->set(stepper_->getStepperName(), *stepper_->getValidParameters());
589 
590  return tempusPL;
591 }
592 
593 
594 // Nonmember constructor
595 // ------------------------------------------------------------------------
596 template<class Scalar>
598  Teuchos::RCP<Teuchos::ParameterList> tempusPL, bool runInitialize)
599 {
600  auto integrator = Teuchos::rcp(new IntegratorBasic<Scalar>());
601  if (tempusPL == Teuchos::null || tempusPL->numParams() == 0)
602  return integrator; // integrator is not initialized (missing model and IC).
603 
604  auto integratorName = tempusPL->get<std::string>("Integrator Name");
605  auto integratorPL = Teuchos::sublist(tempusPL, integratorName, true);
606 
607  std::string integratorType = integratorPL->get<std::string>("Integrator Type");
608  TEUCHOS_TEST_FOR_EXCEPTION( integratorType != "Integrator Basic",
609  std::logic_error,
610  "Error - For IntegratorBasic, 'Integrator Type' should be "
611  << "'Integrator Basic'.\n"
612  << " Integrator Type = " << integratorType << "\n");
613 
614  integrator->setIntegratorName(integratorName);
615 
616  // Validate the Integrator ParameterList
617  auto validPL =
618  Teuchos::rcp_const_cast<Teuchos::ParameterList>(integrator->getValidParameters());
619  auto vIntegratorName = validPL->template get<std::string>("Integrator Name");
620  auto vIntegratorPL = Teuchos::sublist(validPL, vIntegratorName, true);
621  integratorPL->validateParametersAndSetDefaults(*vIntegratorPL,1);
622 
623  // Set Stepper
624  if (integratorPL->isParameter("Stepper Name")) {
625  // Construct from Integrator ParameterList
626  auto stepperName = integratorPL->get<std::string>("Stepper Name");
627  auto stepperPL = Teuchos::sublist(tempusPL, stepperName, true);
628  stepperPL->setName(stepperName);
629  auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
630  integrator->setStepper(sf->createStepper(stepperPL));
631  } else {
632  // Construct default Stepper
633  auto stepper = Teuchos::rcp(new StepperForwardEuler<Scalar>());
634  integrator->setStepper(stepper);
635  }
636 
637  // Set TimeStepControl
638  if (integratorPL->isSublist("Time Step Control")) {
639  // Construct from Integrator ParameterList
640  auto tscPL = Teuchos::sublist(integratorPL, "Time Step Control", true);
641  integrator->setTimeStepControl(createTimeStepControl<Scalar>(tscPL, runInitialize));
642  } else {
643  // Construct default TimeStepControl
644  integrator->setTimeStepControl(rcp(new TimeStepControl<Scalar>()));
645  }
646 
647  // Set SolutionHistory
648  if (integratorPL->isSublist("Solution History")) {
649  // Construct from Integrator ParameterList
650  auto shPL = Teuchos::sublist(integratorPL, "Solution History", true);
651  auto sh = createSolutionHistoryPL<Scalar>(shPL);
652  integrator->setSolutionHistory(sh);
653  } else {
654  // Construct default SolutionHistory
655  integrator->setSolutionHistory(createSolutionHistory<Scalar>());
656  }
657 
658  // Set Observer to default.
659  integrator->setObserver(Teuchos::null);
660 
661  // Set screen output interval.
662  integrator->setScreenOutputIndexInterval(
663  integratorPL->get<int>("Screen Output Index Interval",
664  integrator->getScreenOutputIndexInterval()));
665 
666  // Parse screen output indices
667  auto str = integratorPL->get<std::string>("Screen Output Index List", "");
668  integrator->setScreenOutputIndexList(str);
669 
670  return integrator; // integrator is not initialized (missing model and IC).
671 }
672 
673 
674 // Nonmember constructor
675 // ------------------------------------------------------------------------
676 template<class Scalar>
680  bool runInitialize)
681 {
682  auto integrator = createIntegratorBasic<Scalar>(tempusPL, runInitialize);
683  if ( model == Teuchos::null ) return integrator;
684 
686  integrator->setModel(constModel);
687 
688  // Construct default IC state from the application model and TimeStepControl
689  auto newState = createSolutionStateME(integrator->getStepper()->getModel(),
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());
697  newState->setSolutionStatus(Status::PASSED); // ICs are considered passing.
698 
699  // Set SolutionHistory IC
700  auto sh = integrator->getNonConstSolutionHistory();
701  sh->addState(newState);
702  integrator->getStepper()->setInitialConditions(sh);
703 
704  if (runInitialize) integrator->initialize();
705 
706  return integrator;
707 }
708 
709 
711 template<class Scalar>
714  std::string stepperType)
715 {
716  auto integrator = Teuchos::rcp(new IntegratorBasic<Scalar>());
717 
718  auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
719  auto stepper = sf->createStepper(stepperType, model);
720  integrator->setStepper(stepper);
721  integrator->initializeSolutionHistory();
722  integrator->initialize();
723 
724  return integrator;
725 }
726 
727 
729 template<class Scalar>
731 {
733 }
734 
735 
737 template<class Scalar>
740  std::vector<Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > > models,
741  bool runInitialize)
742 {
743  auto integratorName = tempusPL->get<std::string>("Integrator Name");
744  auto integratorPL = Teuchos::sublist(tempusPL, integratorName, true);
745 
746  std::string integratorType = integratorPL->get<std::string>("Integrator Type");
747  TEUCHOS_TEST_FOR_EXCEPTION( integratorType != "Integrator Basic",
748  std::logic_error,
749  "Error - For IntegratorBasic, 'Integrator Type' should be "
750  << "'Integrator Basic'.\n"
751  << " Integrator Type = " << integratorType << "\n");
752 
753  auto integrator = Teuchos::rcp(new IntegratorBasic<Scalar>());
754  integrator->setIntegratorName(integratorName);
755 
756  TEUCHOS_TEST_FOR_EXCEPTION( !integratorPL->isParameter("Stepper Name"),
757  std::logic_error,
758  "Error - Need to set the 'Stepper Name' in 'Integrator Basic'.\n");
759 
760  auto stepperName = integratorPL->get<std::string>("Stepper Name");
761  TEUCHOS_TEST_FOR_EXCEPTION( stepperName == "Operator Split",
762  std::logic_error,
763  "Error - 'Stepper Name' should be 'Operator Split'.\n");
764 
765  // Construct Steppers from Integrator ParameterList
766  auto stepperPL = Teuchos::sublist(tempusPL, stepperName, true);
767  stepperPL->setName(stepperName);
768  auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
769  integrator->setStepper(sf->createStepper(stepperPL, models));
770 
771  // Set TimeStepControl
772  if (integratorPL->isSublist("Time Step Control")) {
773  // Construct from Integrator ParameterList
774  auto tscPL = Teuchos::sublist(integratorPL, "Time Step Control", true);
775  integrator->setTimeStepControl(createTimeStepControl<Scalar>(tscPL, runInitialize));
776  } else {
777  // Construct default TimeStepControl
778  integrator->setTimeStepControl(rcp(new TimeStepControl<Scalar>()));
779  }
780 
781  // Construct default IC state from the application model and TimeStepControl
782  auto newState = createSolutionStateME(integrator->getStepper()->getModel(),
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());
790  newState->setSolutionStatus(Status::PASSED); // ICs are considered passing.
791 
792  // Set SolutionHistory
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);
798 
799  // Set Observer to default.
800  integrator->setObserver(Teuchos::null);
801 
802  // Set screen output interval.
803  integrator->setScreenOutputIndexInterval(
804  integratorPL->get<int>("Screen Output Index Interval",
805  integrator->getScreenOutputIndexInterval()));
806 
807  // Parse screen output indices
808  auto str = integratorPL->get<std::string>("Screen Output Index List", "");
809  integrator->setScreenOutputIndexList(str);
810 
811  auto validPL = Teuchos::rcp_const_cast<Teuchos::ParameterList>(integrator->getValidParameters());
812 
813  // Validate the Integrator ParameterList
814  auto vIntegratorName = validPL->template get<std::string>("Integrator Name");
815  auto vIntegratorPL = Teuchos::sublist(validPL, vIntegratorName, true);
816  integratorPL->validateParametersAndSetDefaults(*vIntegratorPL);
817 
818  // Validate the Stepper ParameterList
819  auto vStepperName = vIntegratorPL->template get<std::string>("Stepper Name");
820  auto vStepperPL = Teuchos::sublist(validPL, vStepperName, true);
821  stepperPL->validateParametersAndSetDefaults(*vStepperPL);
822 
823  integrator->initialize();
824 
825  return integrator;
826 }
827 
828 
829 } // namespace Tempus
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.
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.