9 #ifndef Tempus_StepperBackwardEuler_impl_hpp
10 #define Tempus_StepperBackwardEuler_impl_hpp
13 #include "Tempus_WrapperModelEvaluatorBasic.hpp"
14 #include "Tempus_StepperFactory.hpp"
18 template <
class Scalar>
21 this->setStepperName(
"Backward Euler");
22 this->setStepperType(
"Backward Euler");
23 this->setUseFSAL(
false);
24 this->setICConsistency(
"None");
25 this->setICConsistencyCheck(
false);
26 this->setZeroInitialGuess(
false);
28 this->setAppAction(Teuchos::null);
29 this->setDefaultSolver();
33 template <
class Scalar>
38 std::string ICConsistency,
bool ICConsistencyCheck,
bool zeroInitialGuess,
42 this->setStepperName(
"Backward Euler");
43 this->setStepperType(
"Backward Euler");
44 this->setUseFSAL(useFSAL);
45 this->setICConsistency(ICConsistency);
46 this->setICConsistencyCheck(ICConsistencyCheck);
47 this->setZeroInitialGuess(zeroInitialGuess);
49 this->setAppAction(stepperBEAppAction);
50 this->setSolver(solver);
51 this->setPredictor(predictorStepper);
53 if (appModel != Teuchos::null) {
54 this->setModel(appModel);
59 template <
class Scalar>
62 if (predictorType ==
"None") {
63 predictorStepper_ = Teuchos::null;
68 if (this->wrapperModel_ != Teuchos::null &&
69 this->wrapperModel_->getAppModel() != Teuchos::null) {
71 sf->createStepper(predictorType, this->wrapperModel_->getAppModel());
74 predictorStepper_ = sf->createStepper(predictorType);
77 this->isInitialized_ =
false;
81 template <
class Scalar>
85 predictorStepper_ = predictorStepper;
86 if (predictorStepper_ == Teuchos::null)
return;
89 predictorStepper_->getModel() == Teuchos::null &&
90 this->wrapperModel_->getAppModel() == Teuchos::null,
92 "Error - Need to set the model, setModel(), before calling "
93 "StepperBackwardEuler::setPredictor()\n");
95 if (predictorStepper_->getModel() == Teuchos::null)
96 predictorStepper_->setModel(this->wrapperModel_->getAppModel());
97 predictorStepper_->initialize();
99 this->isInitialized_ =
false;
102 template <
class Scalar>
106 if (appAction == Teuchos::null) {
108 stepperBEAppAction_ =
112 stepperBEAppAction_ = appAction;
114 this->isInitialized_ =
false;
117 template <
class Scalar>
123 if (predictorStepper_ != Teuchos::null) {
125 if (predictorStepper_->getModel() == Teuchos::null) {
126 predictorStepper_->setModel(appModel);
127 predictorStepper_->initialize();
131 this->isInitialized_ =
false;
134 template <
class Scalar>
140 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
143 if (initialState->getXDot() == Teuchos::null)
144 this->setStepperXDot(initialState->getX()->clone_v());
146 this->setStepperXDot(initialState->getXDot());
151 template <
class Scalar>
155 this->checkInitialized();
159 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperBackwardEuler::takeStep()");
162 solutionHistory->getNumStates() < 2, std::logic_error,
163 "Error - StepperBackwardEuler<Scalar>::takeStep(...)\n"
164 <<
"Need at least two SolutionStates for Backward Euler.\n"
165 <<
" Number of States = " << solutionHistory->getNumStates()
166 <<
"\n Try setting in \"Solution History\" \"Storage Type\" = "
167 <<
"\"Undo\"\n or \"Storage Type\" = \"Static\" and \"Storage Limit\" = "
170 RCP<StepperBackwardEuler<Scalar> > thisStepper = Teuchos::rcpFromRef(*
this);
171 stepperBEAppAction_->execute(
172 solutionHistory, thisStepper,
175 RCP<SolutionState<Scalar> > workingState =
176 solutionHistory->getWorkingState();
177 RCP<SolutionState<Scalar> > currentState =
178 solutionHistory->getCurrentState();
180 RCP<const Thyra::VectorBase<Scalar> > xOld = currentState->getX();
181 RCP<Thyra::VectorBase<Scalar> > x = workingState->getX();
182 if (workingState->getXDot() != Teuchos::null)
183 this->setStepperXDot(workingState->getXDot());
184 RCP<Thyra::VectorBase<Scalar> > xDot = this->getStepperXDot();
186 computePredictor(solutionHistory);
189 const Scalar time = workingState->getTime();
190 const Scalar dt = workingState->getTimeStep();
196 const Scalar alpha = Scalar(1.0) / dt;
197 const Scalar beta = Scalar(1.0);
201 stepperBEAppAction_->execute(
202 solutionHistory, thisStepper,
206 this->solveImplicitODE(x, xDot, time, p);
208 stepperBEAppAction_->execute(
209 solutionHistory, thisStepper,
212 if (workingState->getXDot() != Teuchos::null) timeDer->compute(x, xDot);
214 workingState->setSolutionStatus(sStatus);
215 workingState->setOrder(this->getOrder());
216 workingState->computeNorms(currentState);
217 stepperBEAppAction_->execute(
218 solutionHistory, thisStepper,
224 template <
class Scalar>
228 if (predictorStepper_ == Teuchos::null)
return;
229 predictorStepper_->takeStep(solutionHistory);
231 if (solutionHistory->getWorkingState()->getSolutionStatus() ==
234 Teuchos::OSTab ostab(out, 1,
"StepperBackwardEuler::computePredictor");
235 *out <<
"Warning - predictorStepper has failed." << std::endl;
239 solutionHistory->getWorkingState()->setSolutionStatus(
Status::WORKING);
249 template <
class Scalar>
258 template <
class Scalar>
267 out <<
"--- StepperBackwardEuler ---\n";
268 out <<
" predictorStepper_ = " << predictorStepper_
270 if (predictorStepper_ != Teuchos::null) {
271 out <<
" predictorStepper_->isInitialized() = "
273 out <<
" predictor stepper type = "
274 << predictorStepper_->description() << std::endl;
276 out <<
" stepperBEAppAction_ = " << stepperBEAppAction_
278 out <<
"----------------------------" << std::endl;
281 template <
class Scalar>
286 bool isValidSetup =
true;
291 if (predictorStepper_ != Teuchos::null) {
292 if (!predictorStepper_->isInitialized()) {
293 isValidSetup =
false;
294 out <<
"The predictor stepper is not initialized!\n";
298 if (stepperBEAppAction_ == Teuchos::null) {
299 isValidSetup =
false;
300 out <<
"The Backward Euler AppAction is not set!\n";
306 template <
class Scalar>
310 auto pl = this->getValidParametersBasicImplicit();
311 if (predictorStepper_ == Teuchos::null)
312 pl->set(
"Predictor Stepper Type",
"None");
314 pl->set(
"Predictor Stepper Type", predictorStepper_->getStepperType());
318 template <
class Scalar>
324 template <
class Scalar>
329 const int param_index)
const
332 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->createOutArgs();
333 outArgs.set_f(Teuchos::rcpFromRef(residual));
334 computeStepResidDerivImpl(outArgs, x, t, p, param_index);
337 template <
class Scalar>
342 const int param_index,
const int deriv_index)
const
345 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->createOutArgs();
347 outArgs.set_W_op(Teuchos::rcpFromRef(jacobian));
348 computeStepResidDerivImpl(outArgs, x, t, p, param_index, deriv_index);
351 template <
class Scalar>
356 const int param_index)
const
359 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->createOutArgs();
361 .supports(MEB::DERIV_LINEAR_OP));
362 outArgs.set_DfDp(param_index,
363 MEB::Derivative<Scalar>(Teuchos::rcpFromRef(deriv)));
364 computeStepResidDerivImpl(outArgs, x, t, p, param_index);
367 template <
class Scalar>
372 const int param_index)
const
375 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->createOutArgs();
377 outArgs.set_W(Teuchos::rcpFromRef(jacobian_solver));
378 computeStepResidDerivImpl(outArgs, x, t, p, param_index, 0);
381 template <
class Scalar>
386 const int param_index,
const int deriv_index)
const
393 RCP<const Thyra::VectorBase<Scalar> > xn = x[0];
394 RCP<const Thyra::VectorBase<Scalar> > xo = x[1];
395 const Scalar tn = t[0];
396 const Scalar to = t[1];
397 const Scalar dt = tn - to;
400 RCP<Thyra::VectorBase<Scalar> > x_dot = xn->clone_v();
403 timeDer->compute(xn, x_dot);
406 MEB::InArgs<Scalar> inArgs = this->wrapperModel_->createInArgs();
408 if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(x_dot);
409 if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(tn);
410 if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(dt);
411 inArgs.set_p(param_index, Teuchos::rcpFromRef(p));
414 if (deriv_index == 0) {
416 inArgs.set_alpha(Scalar(1.0) / dt);
417 inArgs.set_beta(Scalar(1.0));
419 else if (deriv_index == 1) {
421 inArgs.set_alpha(Scalar(-1.0) / dt);
422 inArgs.set_beta(Scalar(0.0));
424 this->wrapperModel_->getAppModel()->evalModel(inArgs, outArgs);
429 template <
class Scalar>
436 stepper->setStepperImplicitValues(pl);
438 if (pl != Teuchos::null) {
439 std::string predictorName =
440 pl->
get<std::string>(
"Predictor Stepper Type",
"None");
441 stepper->setPredictor(predictorName);
444 if (model != Teuchos::null) {
445 stepper->setModel(model);
446 stepper->initialize();
453 #endif // Tempus_StepperBackwardEuler_impl_hpp
void computeStepResidDerivImpl(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs, const Teuchos::Array< Teuchos::RCP< const Thyra::VectorBase< Scalar > > > &x, const Teuchos::Array< Scalar > &t, const Thyra::VectorBase< Scalar > &p, const int param_index, const int deriv_index=0) const
Implementation of computeStep*() methods.
StepperBackwardEuler()
Default constructor.
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void computePredictor(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Compute predictor given the supplied stepper.
virtual void computeStepParamDeriv(Thyra::LinearOpBase< Scalar > &deriv, const Teuchos::Array< Teuchos::RCP< const Thyra::VectorBase< Scalar > > > &x, const Teuchos::Array< Scalar > &t, const Thyra::VectorBase< Scalar > &p, const int param_index) const override
Compute time step derivative w.r.t. model parameters.
Default modifier for StepperBackwardEuler.
virtual void computeStepResidual(Thyra::VectorBase< Scalar > &residual, const Teuchos::Array< Teuchos::RCP< const Thyra::VectorBase< Scalar > > > &x, const Teuchos::Array< Scalar > &t, const Thyra::VectorBase< Scalar > &p, const int param_index) const override
Compute time step residual.
Application Action for StepperBackwardEuler.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void computeStepSolver(Thyra::LinearOpWithSolveBase< Scalar > &jacobian_solver, const Teuchos::Array< Teuchos::RCP< const Thyra::VectorBase< Scalar > > > &x, const Teuchos::Array< Scalar > &t, const Thyra::VectorBase< Scalar > &p, const int param_index) const override
Compute time step Jacobian solver.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
Thyra Base interface for time steppers.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Set the initial conditions and make them consistent.
StepperState is a simple class to hold state information about the stepper.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Thyra Base interface for implicit time steppers.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Take the specified timestep, dt, and return true if successful.
virtual bool isValidSetup(Teuchos::FancyOStream &out) const override
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Return a valid ParameterList with current settings.
void setPredictor(std::string predictorType="None")
Set the predictor.
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState() override
Get a default (initial) StepperState.
Backward Euler time stepper.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
virtual void setAppAction(Teuchos::RCP< StepperBackwardEulerAppAction< Scalar > > appAction)
virtual void computeStepJacobian(Thyra::LinearOpBase< Scalar > &jacobian, const Teuchos::Array< Teuchos::RCP< const Thyra::VectorBase< Scalar > > > &x, const Teuchos::Array< Scalar > &t, const Thyra::VectorBase< Scalar > &p, const int param_index, const int deriv_index) const override
Compute time step Jacobian.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Set the initial conditions and make them consistent.
Teuchos::RCP< StepperBackwardEuler< Scalar > > createStepperBackwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
Time-derivative interface for Backward Euler.
#define TEUCHOS_ASSERT(assertion_test)
std::string toString(const T &t)
virtual int stencilLength() const override
Return the number of solution vectors in the time step stencil.