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