Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
19namespace Tempus {
20
21template<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
44template<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
71template<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
88template<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
98template<class Scalar>
100 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > model)
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
109template<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
120template<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
161template<class Scalar>
164 Teuchos::RCP<const Thyra::VectorBase<Scalar> > x0,
165 Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdot0,
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
203template<class Scalar>
205 Teuchos::RCP<SolutionHistory<Scalar> > sh)
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
217template<class Scalar>
219 Teuchos::RCP<TimeStepControl<Scalar> > tsc)
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
233template<class Scalar>
235 Teuchos::RCP<IntegratorObserver<Scalar> > obs)
236{
237 if (obs == Teuchos::null)
238 integratorObserver_ = Teuchos::rcp(new IntegratorObserverBasic<Scalar>());
239 else
240 integratorObserver_ = obs;
241}
242
243
244template<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
265template<class Scalar>
267{
268 std::string name = "Tempus::IntegratorBasic";
269 return(name);
270}
271
272
273template<class Scalar>
275 Teuchos::FancyOStream &out,
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
305template <class Scalar>
306bool 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
315template <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
344template <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
391template <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
416template <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
503template <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
519template <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
545template <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
559template<class Scalar>
560Teuchos::RCP<const Teuchos::ParameterList>
562{
563 Teuchos::RCP<Teuchos::ParameterList> pl =
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
583 Teuchos::RCP<Teuchos::ParameterList> tempusPL =
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// ------------------------------------------------------------------------
596template<class Scalar>
597Teuchos::RCP<IntegratorBasic<Scalar> > createIntegratorBasic(
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// ------------------------------------------------------------------------
676template<class Scalar>
677Teuchos::RCP<IntegratorBasic<Scalar> > createIntegratorBasic(
678 Teuchos::RCP<Teuchos::ParameterList> tempusPL,
679 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
680 bool runInitialize)
681{
682 auto integrator = createIntegratorBasic<Scalar>(tempusPL, runInitialize);
683 if ( model == Teuchos::null ) return integrator;
684
685 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > constModel = model;
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
711template<class Scalar>
712Teuchos::RCP<IntegratorBasic<Scalar> > createIntegratorBasic(
713 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
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
729template<class Scalar>
730Teuchos::RCP<IntegratorBasic<Scalar> > createIntegratorBasic()
731{
732 return Teuchos::rcp(new IntegratorBasic<Scalar>());
733}
734
735
737template<class Scalar>
738Teuchos::RCP<IntegratorBasic<Scalar> > createIntegratorBasic(
739 Teuchos::RCP<Teuchos::ParameterList> tempusPL,
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 void checkTimeStep()
Check if time step has passed or failed.
virtual void initializeSolutionHistory(Teuchos::RCP< SolutionState< Scalar > > state=Teuchos::null)
Set the initial state which has the initial conditions.
virtual void setSolutionHistory(Teuchos::RCP< SolutionHistory< Scalar > > sh=Teuchos::null)
Set the SolutionHistory.
virtual void startIntegrator()
Perform tasks before start of integrator.
std::string description() const override
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Create valid IntegratorBasic ParameterList.
void setIntegratorName(std::string i)
Set the Integrator Name.
virtual void setStepper(Teuchos::RCP< Stepper< Scalar > > stepper)
Set the Stepper.
virtual void setTimeStepControl(Teuchos::RCP< TimeStepControl< Scalar > > tsc=Teuchos::null)
Set the TimeStepControl.
virtual void initialize()
Initializes the Integrator after set* function calls.
void setIntegratorType(std::string i)
Set the Integrator Type.
virtual void endIntegrator()
Perform tasks after end of integrator.
IntegratorBasic()
Default constructor (requires calls to setModel and setSolutionHistory for initial conditions before ...
virtual void startTimeStep()
Start time step.
virtual void copy(Teuchos::RCP< IntegratorBasic< Scalar > > iB)
Copy (a shallow copy)
Teuchos::RCP< Teuchos::Time > integratorTimer_
virtual void setScreenOutputIndexList(std::vector< int > indices)
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void setModel(Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > model)
Set the model on the stepper.
Teuchos::RCP< Teuchos::Time > stepperTimer_
virtual bool advanceTime()
Advance the solution to timeMax, and return true if successful.
virtual void setObserver(Teuchos::RCP< IntegratorObserver< Scalar > > obs=Teuchos::null)
Set the Observer.
virtual std::string getScreenOutputIndexListString() const
IntegratorObserverBasic class for time integrators. This basic class has simple no-op functions,...
IntegratorObserver class for time integrators.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Thyra Base interface for time steppers.
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
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.
bool approxEqual(Scalar a, Scalar b, Scalar relTol=numericalTol< Scalar >())
Test if values are approximately equal within the relative tolerance.
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 toString(const Status status)
Convert Status to string.
Teuchos::RCP< IntegratorBasic< Scalar > > createIntegratorBasic()
Nonmember constructor.