9 #ifndef Tempus_StepperBackwardEuler_impl_hpp
10 #define Tempus_StepperBackwardEuler_impl_hpp
12 #include "Tempus_config.hpp"
14 #include "Tempus_WrapperModelEvaluatorBasic.hpp"
15 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
16 #include "NOX_Thyra.H"
25 template<
class Scalar>
28 this->setParameterList(Teuchos::null);
33 template<
class Scalar>
35 const Teuchos::RCP<
const Thyra::ModelEvaluator<Scalar> >& appModel,
36 Teuchos::RCP<Teuchos::ParameterList> pList)
38 this->setParameterList(pList);
40 if (appModel == Teuchos::null) {
44 this->setModel(appModel);
55 template<
class Scalar>
59 using Teuchos::ParameterList;
61 RCP<ParameterList> predPL =
62 Teuchos::sublist(this->stepperPL_, predictorName,
true);
63 this->stepperPL_->set(
"Predictor Name", predictorName);
64 if (predictorStepper_ != Teuchos::null) predictorStepper_ = Teuchos::null;
74 template<
class Scalar>
76 Teuchos::RCP<Teuchos::ParameterList> predPL)
79 using Teuchos::ParameterList;
81 Teuchos::RCP<Teuchos::ParameterList> stepperPL = this->stepperPL_;
82 std::string predictorName =
83 stepperPL->get<std::string>(
"Predictor Name",
"None");
84 if (is_null(predPL)) {
85 if (predictorName !=
"None") {
86 predPL = Teuchos::sublist(stepperPL, predictorName,
true);
87 RCP<StepperFactory<Scalar> > sf =
90 sf->createStepper(predPL, this->wrapperModel_->getAppModel());
93 TEUCHOS_TEST_FOR_EXCEPTION( predictorName == predPL->name(),
95 "Error - Trying to add a predictor that is already in ParameterList!\n"
96 <<
" Stepper Type = " << stepperPL->get<std::string>(
"Stepper Type")
97 <<
"\n" <<
" Predictor Name = "<<predictorName<<
"\n");
98 predictorName = predPL->name();
99 stepperPL->set(
"Predictor Name", predictorName);
100 stepperPL->set(predictorName, *predPL);
101 if (predictorStepper_ != Teuchos::null) predictorStepper_ = Teuchos::null;
102 RCP<StepperFactory<Scalar> > sf =
105 sf->createStepper(predPL, this->wrapperModel_->getAppModel());
110 template<
class Scalar>
114 if (obs == Teuchos::null) {
116 if (this->stepperObserver_ == Teuchos::null) {
119 this->stepperObserver_ =
123 this->stepperObserver_ = obs;
126 (this->stepperObserver_);
131 template<
class Scalar>
134 TEUCHOS_TEST_FOR_EXCEPTION(
135 this->wrapperModel_ == Teuchos::null, std::logic_error,
136 "Error - Need to set the model, setModel(), before calling "
137 "StepperBackwardEuler::initialize()\n");
139 this->setParameterList(this->stepperPL_);
141 this->setPredictor();
146 template<
class Scalar>
152 RCP<SolutionState<Scalar> > initialState =
solutionHistory->getCurrentState();
155 if (initialState->getXDot() == Teuchos::null)
156 this->setStepperXDot(initialState->getX()->clone_v());
160 if (this->getUseFSAL()) {
161 RCP<Teuchos::FancyOStream> out = this->getOStream();
162 Teuchos::OSTab ostab(out,1,
"StepperBackwardEuler::setInitialConditions()");
163 *out <<
"\nWarning -- The First-Step-As-Last (FSAL) principle is not "
164 <<
"needed with Backward Euler. The default is to set useFSAL=false, "
165 <<
"however useFSAL=true will also work but have no affect "
166 <<
"(i.e., no-op).\n" << std::endl;
171 template<
class Scalar>
177 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperBackwardEuler::takeStep()");
181 "Error - StepperBackwardEuler<Scalar>::takeStep(...)\n"
182 "Need at least two SolutionStates for Backward Euler.\n"
184 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
185 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
188 RCP<SolutionState<Scalar> > workingState=
solutionHistory->getWorkingState();
189 RCP<SolutionState<Scalar> > currentState=
solutionHistory->getCurrentState();
191 RCP<const Thyra::VectorBase<Scalar> > xOld = currentState->getX();
192 RCP<Thyra::VectorBase<Scalar> > x = workingState->getX();
194 RCP<Thyra::VectorBase<Scalar> > xDot = this->getStepperXDot(workingState);
200 const Scalar time = workingState->getTime();
201 const Scalar dt = workingState->getTimeStep();
204 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
206 Scalar(1.0)/dt,xOld));
208 const Scalar alpha = Scalar(1.0)/dt;
209 const Scalar beta = Scalar(1.0);
210 Teuchos::RCP<ImplicitODEParameters<Scalar> > p =
214 if (!Teuchos::is_null(stepperBEObserver_))
217 const Thyra::SolveStatus<Scalar> sStatus =
218 this->solveImplicitODE(x, xDot, time, p);
220 if (!Teuchos::is_null(stepperBEObserver_))
223 if (workingState->getXDot() != Teuchos::null)
224 timeDer->compute(x, xDot);
226 workingState->setSolutionStatus(sStatus);
227 workingState->setOrder(this->getOrder());
233 template<
class Scalar>
237 if (predictorStepper_ == Teuchos::null)
return;
241 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
242 Teuchos::OSTab ostab(out,1,
"StepperBackwardEuler::computePredictor");
243 *out <<
"Warning - predictorStepper has failed." << std::endl;
257 template<
class Scalar>
258 Teuchos::RCP<Tempus::StepperState<Scalar> >
262 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
268 template<
class Scalar>
271 std::string name =
"Backward Euler";
276 template<
class Scalar>
278 Teuchos::FancyOStream &out,
279 const Teuchos::EVerbosityLevel )
const
281 out << description() <<
"::describe:" << std::endl
282 <<
"wrapperModel_ = " << this->wrapperModel_->description() << std::endl;
286 template <
class Scalar>
288 Teuchos::RCP<Teuchos::ParameterList>
const& pList)
290 Teuchos::RCP<Teuchos::ParameterList> stepperPL = this->stepperPL_;
291 if (pList == Teuchos::null) {
293 if (stepperPL == Teuchos::null) stepperPL = this->getDefaultParameters();
297 if (!(stepperPL->isParameter(
"Solver Name"))) {
298 stepperPL->set<std::string>(
"Solver Name",
"Default Solver");
299 Teuchos::RCP<Teuchos::ParameterList> solverPL =
300 this->defaultSolverParameters();
301 stepperPL->set(
"Default Solver", *solverPL);
306 std::string stepperType = stepperPL->get<std::string>(
"Stepper Type");
307 TEUCHOS_TEST_FOR_EXCEPTION( stepperType !=
"Backward Euler",
309 "Error - Stepper Type is not 'Backward Euler'!\n"
310 <<
" Stepper Type = "<<stepperPL->get<std::string>(
"Stepper Type")<<
"\n");
312 this->stepperPL_ = stepperPL;
316 template<
class Scalar>
317 Teuchos::RCP<const Teuchos::ParameterList>
320 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
321 pl->setName(
"Default Stepper - " + this->description());
322 pl->set<std::string>(
"Stepper Type", this->description());
323 this->getValidParametersBasic(pl);
324 pl->set<
bool> (
"Zero Initial Guess",
false);
325 pl->set<std::string>(
"Solver Name",
"",
326 "Name of ParameterList containing the solver specifications.");
332 template<
class Scalar>
333 Teuchos::RCP<Teuchos::ParameterList>
337 using Teuchos::ParameterList;
338 using Teuchos::rcp_const_cast;
340 RCP<ParameterList> pl =
341 rcp_const_cast<ParameterList>(this->getValidParameters());
343 pl->set<std::string>(
"Solver Name",
"Default Solver");
344 RCP<ParameterList> solverPL = this->defaultSolverParameters();
345 pl->set(
"Default Solver", *solverPL);
351 template <
class Scalar>
352 Teuchos::RCP<Teuchos::ParameterList>
355 return(this->stepperPL_);
359 template <
class Scalar>
360 Teuchos::RCP<Teuchos::ParameterList>
363 Teuchos::RCP<Teuchos::ParameterList> temp_plist = this->stepperPL_;
364 this->stepperPL_ = Teuchos::null;
369 template <
class Scalar>
377 template <
class Scalar>
380 Thyra::VectorBase<Scalar>& residual,
381 const Teuchos::Array< Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >& x,
382 const Teuchos::Array<Scalar>& t,
383 const Thyra::VectorBase<Scalar>& p,
384 const int param_index)
const
386 typedef Thyra::ModelEvaluatorBase MEB;
387 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
388 outArgs.set_f(Teuchos::rcpFromRef(residual));
389 computeStepResidDerivImpl(outArgs, x, t, p, param_index);
392 template <
class Scalar>
395 Thyra::LinearOpBase<Scalar>& jacobian,
396 const Teuchos::Array< Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >& x,
397 const Teuchos::Array<Scalar>& t,
398 const Thyra::VectorBase<Scalar>& p,
399 const int param_index,
400 const int deriv_index)
const
402 typedef Thyra::ModelEvaluatorBase MEB;
403 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
404 TEUCHOS_ASSERT(outArgs.supports(MEB::OUT_ARG_W_op));
405 outArgs.set_W_op(Teuchos::rcpFromRef(jacobian));
406 computeStepResidDerivImpl(outArgs, x, t, p, param_index, deriv_index);
409 template <
class Scalar>
412 Thyra::LinearOpBase<Scalar>& deriv,
413 const Teuchos::Array< Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >& x,
414 const Teuchos::Array<Scalar>& t,
415 const Thyra::VectorBase<Scalar>& p,
416 const int param_index)
const
418 typedef Thyra::ModelEvaluatorBase MEB;
419 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
420 TEUCHOS_ASSERT(outArgs.supports(MEB::OUT_ARG_DfDp, param_index).supports(MEB::DERIV_LINEAR_OP));
421 outArgs.set_DfDp(param_index,
422 MEB::Derivative<Scalar>(Teuchos::rcpFromRef(deriv)));
423 computeStepResidDerivImpl(outArgs, x, t, p, param_index);
426 template <
class Scalar>
429 Thyra::LinearOpWithSolveBase<Scalar>& jacobian_solver,
430 const Teuchos::Array< Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >& x,
431 const Teuchos::Array<Scalar>& t,
432 const Thyra::VectorBase<Scalar>& p,
433 const int param_index)
const
435 typedef Thyra::ModelEvaluatorBase MEB;
436 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
437 TEUCHOS_ASSERT(outArgs.supports(MEB::OUT_ARG_W));
438 outArgs.set_W(Teuchos::rcpFromRef(jacobian_solver));
439 computeStepResidDerivImpl(outArgs, x, t, p, param_index, 0);
442 template <
class Scalar>
445 const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs,
446 const Teuchos::Array< Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >& x,
447 const Teuchos::Array<Scalar>& t,
448 const Thyra::VectorBase<Scalar>& p,
449 const int param_index,
450 const int deriv_index)
const
453 typedef Thyra::ModelEvaluatorBase MEB;
455 TEUCHOS_ASSERT(x.size() == 2);
456 TEUCHOS_ASSERT(t.size() == 2);
457 RCP<const Thyra::VectorBase<Scalar> > xn = x[0];
458 RCP<const Thyra::VectorBase<Scalar> > xo = x[1];
459 const Scalar tn = t[0];
460 const Scalar to = t[1];
461 const Scalar dt = tn-to;
464 RCP<Thyra::VectorBase<Scalar> > x_dot = xn->clone_v();
465 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
467 timeDer->compute(xn, x_dot);
470 MEB::InArgs<Scalar> inArgs = this->wrapperModel_->getInArgs();
472 if (inArgs.supports(MEB::IN_ARG_x_dot )) inArgs.set_x_dot (x_dot);
473 if (inArgs.supports(MEB::IN_ARG_t )) inArgs.set_t (tn);
474 if (inArgs.supports(MEB::IN_ARG_step_size )) inArgs.set_step_size(dt);
475 inArgs.set_p(param_index, Teuchos::rcpFromRef(p));
476 TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_alpha));
477 TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_beta));
478 if (deriv_index == 0) {
480 inArgs.set_alpha(Scalar(1.0)/dt);
481 inArgs.set_beta(Scalar(1.0));
483 else if (deriv_index == 1) {
485 inArgs.set_alpha(Scalar(-1.0)/dt);
486 inArgs.set_beta(Scalar(0.0));
488 this->wrapperModel_->getAppModel()->evalModel(inArgs, outArgs);
492 #endif // Tempus_StepperBackwardEuler_impl_hpp
StepperBackwardEulerObserver class for StepperBackwardEuler.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
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.
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
Compute time step residual.
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
virtual void computePredictor(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Compute predictor given the supplied stepper.
virtual int stencilLength() const
Return the number of solution vectors in the time step stencil.
Teuchos::RCP< Teuchos::ParameterList > getDefaultParameters() const
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
void setPredictor(std::string predictorName)
Set the predictor.
virtual void setObserver(Teuchos::RCP< StepperObserver< Scalar > > obs=Teuchos::null)
Set Observer.
StepperState is a simple class to hold state information about the stepper.
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
Compute time step Jacobian solver.
virtual std::string description() const
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
Compute time step Jacobian.
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
Compute time step derivative w.r.t. model parameters.
StepperObserver class for Stepper class.
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl)
Time-derivative interface for Backward Euler.
virtual void initialize()
Initialize during construction and after changing input parameters.
Solve for x and determine xDot from x.