9 #ifndef Tempus_StepperBackwardEuler_impl_hpp
10 #define Tempus_StepperBackwardEuler_impl_hpp
13 #include "Tempus_WrapperModelEvaluatorBasic.hpp"
14 #include "Tempus_StepperFactory.hpp"
20 template<
class Scalar>
23 this->setStepperName(
"Backward Euler");
24 this->setStepperType(
"Backward Euler");
25 this->setUseFSAL(
false);
26 this->setICConsistency(
"None");
27 this->setICConsistencyCheck(
false);
28 this->setZeroInitialGuess(
false);
30 this->setAppAction(Teuchos::null);
31 this->setDefaultSolver();
36 template<
class Scalar>
42 std::string ICConsistency,
43 bool ICConsistencyCheck,
44 bool zeroInitialGuess,
47 this->setStepperName(
"Backward Euler");
48 this->setStepperType(
"Backward Euler");
49 this->setUseFSAL( useFSAL);
50 this->setICConsistency( ICConsistency);
51 this->setICConsistencyCheck( ICConsistencyCheck);
52 this->setZeroInitialGuess( zeroInitialGuess);
54 this->setAppAction(stepperBEAppAction);
55 this->setSolver(solver);
56 this->setPredictor(predictorStepper);
58 if (appModel != Teuchos::null) {
59 this->setModel(appModel);
65 template<
class Scalar>
68 if (predictorType ==
"None") {
69 predictorStepper_ = Teuchos::null;
74 if (this->wrapperModel_ != Teuchos::null &&
75 this->wrapperModel_->getAppModel() != Teuchos::null) {
76 predictorStepper_ = sf->createStepper(predictorType,
77 this->wrapperModel_->getAppModel());
79 predictorStepper_ = sf->createStepper(predictorType);
82 this->isInitialized_ =
false;
87 template<
class Scalar>
91 predictorStepper_ = predictorStepper;
92 if (predictorStepper_ == Teuchos::null)
return;
95 predictorStepper_->getModel() == Teuchos::null &&
96 this->wrapperModel_->getAppModel() == Teuchos::null, std::logic_error,
97 "Error - Need to set the model, setModel(), before calling "
98 "StepperBackwardEuler::setPredictor()\n");
100 if (predictorStepper_->getModel() == Teuchos::null)
101 predictorStepper_->setModel(this->wrapperModel_->getAppModel());
102 predictorStepper_->initialize();
104 this->isInitialized_ =
false;
108 template<
class Scalar>
112 if (appAction == Teuchos::null) {
114 stepperBEAppAction_ =
117 stepperBEAppAction_ = appAction;
119 this->isInitialized_ =
false;
123 template<
class Scalar>
129 if (predictorStepper_ != Teuchos::null) {
131 if (predictorStepper_->getModel() == Teuchos::null) {
132 predictorStepper_->setModel(appModel);
133 predictorStepper_->initialize();
137 this->isInitialized_ =
false;
141 template<
class Scalar>
147 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
150 if (initialState->getXDot() == Teuchos::null)
151 this->setStepperXDot(initialState->getX()->clone_v());
153 this->setStepperXDot(initialState->getXDot());
159 template<
class Scalar>
163 this->checkInitialized();
167 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperBackwardEuler::takeStep()");
171 "Error - StepperBackwardEuler<Scalar>::takeStep(...)\n"
172 "Need at least two SolutionStates for Backward Euler.\n"
173 " Number of States = " << solutionHistory->getNumStates() <<
"\n"
174 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
175 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
177 RCP<StepperBackwardEuler<Scalar> > thisStepper = Teuchos::rcpFromRef(*
this);
178 stepperBEAppAction_->execute(solutionHistory, thisStepper,
181 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
182 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
184 RCP<const Thyra::VectorBase<Scalar> > xOld = currentState->getX();
185 RCP<Thyra::VectorBase<Scalar> > x = workingState->getX();
186 if (workingState->getXDot() != Teuchos::null)
187 this->setStepperXDot(workingState->getXDot());
188 RCP<Thyra::VectorBase<Scalar> > xDot = this->getStepperXDot();
190 computePredictor(solutionHistory);
194 const Scalar time = workingState->getTime();
195 const Scalar dt = workingState->getTimeStep();
200 Scalar(1.0)/dt,xOld));
202 const Scalar alpha = Scalar(1.0)/dt;
203 const Scalar beta = Scalar(1.0);
205 timeDer, dt, alpha, beta));
207 stepperBEAppAction_->execute(solutionHistory, thisStepper,
211 this->solveImplicitODE(x, xDot, time, p);
213 stepperBEAppAction_->execute(solutionHistory, thisStepper,
216 if (workingState->getXDot() != Teuchos::null)
217 timeDer->compute(x, xDot);
219 workingState->setSolutionStatus(sStatus);
220 workingState->setOrder(this->getOrder());
221 workingState->computeNorms(currentState);
222 stepperBEAppAction_->execute(solutionHistory, thisStepper,
228 template<
class Scalar>
232 if (predictorStepper_ == Teuchos::null)
return;
233 predictorStepper_->takeStep(solutionHistory);
235 if (solutionHistory->getWorkingState()->getSolutionStatus()==
Status::FAILED) {
237 Teuchos::OSTab ostab(out,1,
"StepperBackwardEuler::computePredictor");
238 *out <<
"Warning - predictorStepper has failed." << std::endl;
241 solutionHistory->getWorkingState()->setSolutionStatus(
Status::WORKING);
252 template<
class Scalar>
263 template<
class Scalar>
273 out <<
"--- StepperBackwardEuler ---\n";
274 out <<
" predictorStepper_ = "
275 << predictorStepper_ << std::endl;
276 if (predictorStepper_ != Teuchos::null) {
277 out <<
" predictorStepper_->isInitialized() = "
279 out <<
" predictor stepper type = "
280 << predictorStepper_->description() << std::endl;
282 out <<
" stepperBEAppAction_ = "
283 << stepperBEAppAction_ << std::endl;
284 out <<
"----------------------------" << std::endl;
288 template<
class Scalar>
292 bool isValidSetup =
true;
297 if (predictorStepper_ != Teuchos::null) {
298 if ( !predictorStepper_->isInitialized() ) {
299 isValidSetup =
false;
300 out <<
"The predictor stepper is not initialized!\n";
304 if (stepperBEAppAction_ == Teuchos::null) {
305 isValidSetup =
false;
306 out <<
"The Backward Euler AppAction is not set!\n";
313 template<
class Scalar>
317 auto pl = this->getValidParametersBasicImplicit();
318 if (predictorStepper_ == Teuchos::null)
319 pl->set(
"Predictor Stepper Type",
"None");
321 pl->set(
"Predictor Stepper Type", predictorStepper_->getStepperType());
326 template <
class Scalar>
334 template <
class Scalar>
341 const int param_index)
const
344 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
345 outArgs.set_f(Teuchos::rcpFromRef(residual));
346 computeStepResidDerivImpl(outArgs, x, t, p, param_index);
349 template <
class Scalar>
356 const int param_index,
357 const int deriv_index)
const
360 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
362 outArgs.set_W_op(Teuchos::rcpFromRef(jacobian));
363 computeStepResidDerivImpl(outArgs, x, t, p, param_index, deriv_index);
366 template <
class Scalar>
373 const int param_index)
const
376 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
377 TEUCHOS_ASSERT(outArgs.supports(MEB::OUT_ARG_DfDp, param_index).supports(MEB::DERIV_LINEAR_OP));
378 outArgs.set_DfDp(param_index,
379 MEB::Derivative<Scalar>(Teuchos::rcpFromRef(deriv)));
380 computeStepResidDerivImpl(outArgs, x, t, p, param_index);
383 template <
class Scalar>
390 const int param_index)
const
393 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
395 outArgs.set_W(Teuchos::rcpFromRef(jacobian_solver));
396 computeStepResidDerivImpl(outArgs, x, t, p, param_index, 0);
399 template <
class Scalar>
406 const int param_index,
407 const int deriv_index)
const
414 RCP<const Thyra::VectorBase<Scalar> > xn = x[0];
415 RCP<const Thyra::VectorBase<Scalar> > xo = x[1];
416 const Scalar tn = t[0];
417 const Scalar to = t[1];
418 const Scalar dt = tn-to;
421 RCP<Thyra::VectorBase<Scalar> > x_dot = xn->clone_v();
424 timeDer->compute(xn, x_dot);
427 MEB::InArgs<Scalar> inArgs = this->wrapperModel_->getInArgs();
429 if (inArgs.supports(MEB::IN_ARG_x_dot )) inArgs.set_x_dot (x_dot);
430 if (inArgs.supports(MEB::IN_ARG_t )) inArgs.set_t (tn);
431 if (inArgs.supports(MEB::IN_ARG_step_size )) inArgs.set_step_size(dt);
432 inArgs.set_p(param_index, Teuchos::rcpFromRef(p));
435 if (deriv_index == 0) {
437 inArgs.set_alpha(Scalar(1.0)/dt);
438 inArgs.set_beta(Scalar(1.0));
440 else if (deriv_index == 1) {
442 inArgs.set_alpha(Scalar(-1.0)/dt);
443 inArgs.set_beta(Scalar(0.0));
445 this->wrapperModel_->getAppModel()->evalModel(inArgs, outArgs);
451 template<
class Scalar>
459 stepper->setStepperImplicitValues(pl);
461 if (pl != Teuchos::null) {
462 std::string predictorName =
463 pl->
get<std::string>(
"Predictor Stepper Type",
"None");
464 stepper->setPredictor(predictorName);
467 if (model != Teuchos::null) {
468 stepper->setModel(model);
469 stepper->initialize();
477 #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.