39 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
40 const Teuchos::RCP<Stepper<Scalar> >& startUpStepper,
42 std::string ICConsistency,
43 bool ICConsistencyCheck,
44 bool zeroInitialGuess,
45 const Teuchos::RCP<StepperBDF2AppAction<Scalar> >& stepperBDF2AppAction)
47 this->setStepperName(
"BDF2");
48 this->setStepperType(
"BDF2");
49 this->setUseFSAL( useFSAL);
50 this->setICConsistency( ICConsistency);
51 this->setICConsistencyCheck( ICConsistencyCheck);
52 this->setZeroInitialGuess( zeroInitialGuess);
54 this->setAppAction(stepperBDF2AppAction);
55 this->setSolver(solver);
56 this->setStartUpStepper(startUpStepper);
58 if (appModel != Teuchos::null) {
59 this->setModel(appModel);
82 auto sf = Teuchos::rcp(
new StepperFactory<Scalar>());
83 if (this->wrapperModel_ != Teuchos::null &&
84 this->wrapperModel_->getAppModel() != Teuchos::null) {
85 startUpStepper_ = sf->createStepper(startupStepperType,
86 this->wrapperModel_->getAppModel());
88 startUpStepper_ = sf->createStepper(startupStepperType);
91 this->isInitialized_ =
false;
97 Teuchos::RCP<Stepper<Scalar> > startUpStepper)
99 startUpStepper_ = startUpStepper;
101 if (this->wrapperModel_ != Teuchos::null) {
102 TEUCHOS_TEST_FOR_EXCEPTION(
103 this->wrapperModel_->getAppModel() == Teuchos::null, std::logic_error,
104 "Error - Can not set the startUpStepper to Teuchos::null.\n");
106 if (startUpStepper->getModel() == Teuchos::null &&
107 this->wrapperModel_->getAppModel() != Teuchos::null) {
108 startUpStepper_->setModel(this->wrapperModel_->getAppModel());
109 startUpStepper_->initialize();
113 this->isInitialized_ =
false;
159 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
161 this->checkInitialized();
165 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperBDF2::takeStep()");
167 int numStates = solutionHistory->getNumStates();
169 RCP<Thyra::VectorBase<Scalar> > xOld;
170 RCP<Thyra::VectorBase<Scalar> > xOldOld;
175 computeStartUp(solutionHistory);
178 TEUCHOS_TEST_FOR_EXCEPTION( (numStates < 3), std::logic_error,
179 "Error in Tempus::StepperBDF2::takeStep(): numStates after \n"
180 <<
"startup stepper must be at least 3, whereas numStates = "
181 << numStates <<
"!\n" <<
"If running with Storage Type = Static, "
182 <<
"make sure Storage Limit > 2.\n");
187 RCP<StepperBDF2<Scalar> > thisStepper = Teuchos::rcpFromRef(*
this);
188 stepperBDF2AppAction_->execute(solutionHistory, thisStepper,
192 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
193 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
195 RCP<Thyra::VectorBase<Scalar> > x = workingState->getX();
196 if (workingState->getXDot() != Teuchos::null)
197 this->setStepperXDot(workingState->getXDot());
198 RCP<Thyra::VectorBase<Scalar> > xDot = this->getStepperXDot();
201 const Scalar time = workingState->getTime();
202 const Scalar dt = workingState->getTimeStep();
203 const Scalar dtOld = currentState->getTimeStep();
205 xOld = solutionHistory->getStateTimeIndexNM1()->getX();
206 xOldOld = solutionHistory->getStateTimeIndexNM2()->getX();
207 order_ = Scalar(2.0);
210 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
211 Teuchos::rcp(
new StepperBDF2TimeDerivative<Scalar>(dt, dtOld, xOld, xOldOld));
213 const Scalar alpha = getAlpha(dt, dtOld);
214 const Scalar beta = getBeta (dt);
216 auto p = Teuchos::rcp(
new ImplicitODEParameters<Scalar>(
217 timeDer, dt, alpha, beta));
218 stepperBDF2AppAction_->execute(solutionHistory, thisStepper,
221 const Thyra::SolveStatus<Scalar> sStatus =
222 this->solveImplicitODE(x, xDot, time, p);
224 stepperBDF2AppAction_->execute(solutionHistory, thisStepper,
228 if (workingState->getXDot() != Teuchos::null)
229 timeDer->compute(x, xDot);
231 workingState->setSolutionStatus(sStatus);
232 workingState->setOrder(getOrder());
233 workingState->computeNorms(currentState);
234 stepperBDF2AppAction_->execute(solutionHistory, thisStepper,
242 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
244 auto out = Teuchos::fancyOStream( this->getOStream() );
245 out->setOutputToRootOnly(0);
246 Teuchos::OSTab ostab(out,1,
"StepperBDF2::computeStartUp()");
247 *out <<
"Warning -- Taking a startup step for BDF2 using '"
248 << startUpStepper_->getStepperType()<<
"'!" << std::endl;
251 startUpStepper_->takeStep(solutionHistory);
253 order_ = startUpStepper_->getOrder();
275 Teuchos::FancyOStream &out,
276 const Teuchos::EVerbosityLevel verbLevel )
const
278 auto l_out = Teuchos::fancyOStream( out.getOStream() );
279 Teuchos::OSTab ostab(*l_out, 2, this->description());
280 l_out->setOutputToRootOnly(0);
285 *l_out <<
"--- StepperBDF2 ---\n";
286 if (startUpStepper_ != Teuchos::null) {
287 *l_out <<
" startup stepper type = "
288 << startUpStepper_->description() << std::endl;
290 *l_out <<
" startUpStepper_ = "
291 << startUpStepper_ << std::endl;
292 *l_out <<
" startUpStepper_->isInitialized() = "
293 << Teuchos::toString(startUpStepper_->isInitialized()) << std::endl;
294 *l_out <<
" stepperBDF2AppAction_ = "
295 << stepperBDF2AppAction_ << std::endl;
296 *l_out <<
"----------------------------" << std::endl;
297 *l_out <<
" order_ = " << order_ << std::endl;
298 *l_out <<
"-------------------" << std::endl;
340 Teuchos::RCP<Teuchos::ParameterList> pl)
342 auto stepper = Teuchos::rcp(
new StepperBDF2<Scalar>());
343 stepper->setStepperImplicitValues(pl);
345 std::string startUpStepperName =
"DIRK 1 Stage Theta Method";
346 if (pl != Teuchos::null) startUpStepperName =
347 pl->get<std::string>(
"Start Up Stepper Type", startUpStepperName);
348 stepper->setStartUpStepper(startUpStepperName);
350 if (model != Teuchos::null) {
351 stepper->setModel(model);
352 stepper->initialize();