Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_StepperBackwardEuler_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_StepperBackwardEuler_impl_hpp
10#define Tempus_StepperBackwardEuler_impl_hpp
11
13#include "Tempus_WrapperModelEvaluatorBasic.hpp"
14#include "Tempus_StepperFactory.hpp"
15
16
17namespace Tempus {
18
19
20template<class Scalar>
22{
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);
29
30 this->setAppAction(Teuchos::null);
31 this->setDefaultSolver();
32 this->setPredictor();
33}
34
35
36template<class Scalar>
38 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
39 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
40 const Teuchos::RCP<Stepper<Scalar> >& predictorStepper,
41 bool useFSAL,
42 std::string ICConsistency,
43 bool ICConsistencyCheck,
44 bool zeroInitialGuess,
45 const Teuchos::RCP<StepperBackwardEulerAppAction<Scalar> >& stepperBEAppAction)
46{
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);
53
54 this->setAppAction(stepperBEAppAction);
55 this->setSolver(solver);
56 this->setPredictor(predictorStepper);
57
58 if (appModel != Teuchos::null) {
59 this->setModel(appModel);
60 this->initialize();
61 }
62}
63
64
65template<class Scalar>
66void StepperBackwardEuler<Scalar>::setPredictor(std::string predictorType)
67{
68 if (predictorType == "None") {
69 predictorStepper_ = Teuchos::null;
70 return;
71 }
72
73 auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
74 if (this->wrapperModel_ != Teuchos::null &&
75 this->wrapperModel_->getAppModel() != Teuchos::null) {
76 predictorStepper_ = sf->createStepper(predictorType,
77 this->wrapperModel_->getAppModel());
78 } else {
79 predictorStepper_ = sf->createStepper(predictorType);
80 }
81
82 this->isInitialized_ = false;
83}
84
85
87template<class Scalar>
89 Teuchos::RCP<Stepper<Scalar> > predictorStepper)
90{
91 predictorStepper_ = predictorStepper;
92 if (predictorStepper_ == Teuchos::null) return;
93
94 TEUCHOS_TEST_FOR_EXCEPTION(
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");
99
100 if (predictorStepper_->getModel() == Teuchos::null)
101 predictorStepper_->setModel(this->wrapperModel_->getAppModel());
102 predictorStepper_->initialize();
103
104 this->isInitialized_ = false;
105}
106
107
108template<class Scalar>
110 Teuchos::RCP<StepperBackwardEulerAppAction<Scalar> > appAction)
111{
112 if (appAction == Teuchos::null) {
113 // Create default appAction
114 stepperBEAppAction_ =
116 } else {
117 stepperBEAppAction_ = appAction;
118 }
119 this->isInitialized_ = false;
120}
121
122
123template<class Scalar>
125 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
126{
128
129 if (predictorStepper_ != Teuchos::null) {
130 // If predictor's model is not set, set it to the stepper model.
131 if (predictorStepper_->getModel() == Teuchos::null) {
132 predictorStepper_->setModel(appModel);
133 predictorStepper_->initialize();
134 }
135 }
136
137 this->isInitialized_ = false;
138}
139
140
141template<class Scalar>
143 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
144{
145 using Teuchos::RCP;
146
147 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
148
149 // Check if we need Stepper storage for xDot
150 if (initialState->getXDot() == Teuchos::null)
151 this->setStepperXDot(initialState->getX()->clone_v());
152 else
153 this->setStepperXDot(initialState->getXDot());
154
156}
157
158
159template<class Scalar>
161 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
162{
163 this->checkInitialized();
164
165 using Teuchos::RCP;
166
167 TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperBackwardEuler::takeStep()");
168 {
169 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
170 std::logic_error,
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");
176
177 RCP<StepperBackwardEuler<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
178 stepperBEAppAction_->execute(solutionHistory, thisStepper,
180
181 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
182 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
183
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();
189
190 computePredictor(solutionHistory);
191 if (workingState->getSolutionStatus() == Status::FAILED)
192 return;
193
194 const Scalar time = workingState->getTime();
195 const Scalar dt = workingState->getTimeStep();
196
197 // Setup TimeDerivative
198 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
200 Scalar(1.0)/dt,xOld));
201
202 const Scalar alpha = Scalar(1.0)/dt;
203 const Scalar beta = Scalar(1.0);
204 auto p = Teuchos::rcp(new ImplicitODEParameters<Scalar>(
205 timeDer, dt, alpha, beta));
206
207 stepperBEAppAction_->execute(solutionHistory, thisStepper,
209
210 const Thyra::SolveStatus<Scalar> sStatus =
211 this->solveImplicitODE(x, xDot, time, p);
212
213 stepperBEAppAction_->execute(solutionHistory, thisStepper,
215
216 if (workingState->getXDot() != Teuchos::null)
217 timeDer->compute(x, xDot);
218
219 workingState->setSolutionStatus(sStatus); // Converged --> pass.
220 workingState->setOrder(this->getOrder());
221 workingState->computeNorms(currentState);
222 stepperBEAppAction_->execute(solutionHistory, thisStepper,
224 }
225 return;
226}
227
228template<class Scalar>
230 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
231{
232 if (predictorStepper_ == Teuchos::null) return;
233 predictorStepper_->takeStep(solutionHistory);
234
235 if (solutionHistory->getWorkingState()->getSolutionStatus()==Status::FAILED) {
236 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
237 Teuchos::OSTab ostab(out,1,"StepperBackwardEuler::computePredictor");
238 *out << "Warning - predictorStepper has failed." << std::endl;
239 } else {
240 // Reset status to WORKING since this is the predictor
241 solutionHistory->getWorkingState()->setSolutionStatus(Status::WORKING);
242 }
243}
244
245
252template<class Scalar>
253Teuchos::RCP<Tempus::StepperState<Scalar> >
256{
257 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
258 rcp(new StepperState<Scalar>(this->getStepperType()));
259 return stepperState;
260}
261
262
263template<class Scalar>
265 Teuchos::FancyOStream &out,
266 const Teuchos::EVerbosityLevel verbLevel) const
267{
268 out.setOutputToRootOnly(0);
269 out << std::endl;
270 Stepper<Scalar>::describe(out, verbLevel);
271 StepperImplicit<Scalar>::describe(out, verbLevel);
272
273 out << "--- StepperBackwardEuler ---\n";
274 out << " predictorStepper_ = "
275 << predictorStepper_ << std::endl;
276 if (predictorStepper_ != Teuchos::null) {
277 out << " predictorStepper_->isInitialized() = "
278 << Teuchos::toString(predictorStepper_->isInitialized()) << std::endl;
279 out << " predictor stepper type = "
280 << predictorStepper_->description() << std::endl;
281 }
282 out << " stepperBEAppAction_ = "
283 << stepperBEAppAction_ << std::endl;
284 out << "----------------------------" << std::endl;
285}
286
287
288template<class Scalar>
289bool StepperBackwardEuler<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
290{
291 out.setOutputToRootOnly(0);
292 bool isValidSetup = true;
293
294 if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
295 if ( !StepperImplicit<Scalar>::isValidSetup(out) ) isValidSetup = false;
296
297 if (predictorStepper_ != Teuchos::null) {
298 if ( !predictorStepper_->isInitialized() ) {
299 isValidSetup = false;
300 out << "The predictor stepper is not initialized!\n";
301 }
302 }
303
304 if (stepperBEAppAction_ == Teuchos::null) {
305 isValidSetup = false;
306 out << "The Backward Euler AppAction is not set!\n";
307 }
308
309 return isValidSetup;
310}
311
312
313template<class Scalar>
314Teuchos::RCP<const Teuchos::ParameterList>
316{
317 auto pl = this->getValidParametersBasicImplicit();
318 if (predictorStepper_ == Teuchos::null)
319 pl->set("Predictor Stepper Type", "None");
320 else
321 pl->set("Predictor Stepper Type", predictorStepper_->getStepperType());
322 return pl;
323}
324
325
326template <class Scalar>
327int
329{
330 return 2;
331}
332
333
334template <class Scalar>
335void
338 const Teuchos::Array< Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x,
339 const Teuchos::Array<Scalar>& t,
341 const int param_index) const
342{
343 typedef Thyra::ModelEvaluatorBase MEB;
344 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
345 outArgs.set_f(Teuchos::rcpFromRef(residual));
346 computeStepResidDerivImpl(outArgs, x, t, p, param_index);
347}
348
349template <class Scalar>
350void
353 const Teuchos::Array< Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x,
354 const Teuchos::Array<Scalar>& t,
356 const int param_index,
357 const int deriv_index) const
358{
359 typedef Thyra::ModelEvaluatorBase MEB;
360 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
361 TEUCHOS_ASSERT(outArgs.supports(MEB::OUT_ARG_W_op));
362 outArgs.set_W_op(Teuchos::rcpFromRef(jacobian));
363 computeStepResidDerivImpl(outArgs, x, t, p, param_index, deriv_index);
364}
365
366template <class Scalar>
367void
370 const Teuchos::Array< Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x,
371 const Teuchos::Array<Scalar>& t,
373 const int param_index) const
374{
375 typedef Thyra::ModelEvaluatorBase MEB;
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);
381}
382
383template <class Scalar>
384void
387 const Teuchos::Array< Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x,
388 const Teuchos::Array<Scalar>& t,
390 const int param_index) const
391{
392 typedef Thyra::ModelEvaluatorBase MEB;
393 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
394 TEUCHOS_ASSERT(outArgs.supports(MEB::OUT_ARG_W));
395 outArgs.set_W(Teuchos::rcpFromRef(jacobian_solver));
396 computeStepResidDerivImpl(outArgs, x, t, p, param_index, 0);
397}
398
399template <class Scalar>
400void
402 const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs,
403 const Teuchos::Array< Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x,
404 const Teuchos::Array<Scalar>& t,
406 const int param_index,
407 const int deriv_index) const
408{
409 using Teuchos::RCP;
410 typedef Thyra::ModelEvaluatorBase MEB;
411
412 TEUCHOS_ASSERT(x.size() == 2);
413 TEUCHOS_ASSERT(t.size() == 2);
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;
419
420 // compute x_dot
421 RCP<Thyra::VectorBase<Scalar> > x_dot = xn->clone_v();
422 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
423 Teuchos::rcp(new StepperBackwardEulerTimeDerivative<Scalar>(Scalar(1.0)/dt,xo));
424 timeDer->compute(xn, x_dot);
425
426 // evaluate model
427 MEB::InArgs<Scalar> inArgs = this->wrapperModel_->getInArgs();
428 inArgs.set_x(xn);
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));
433 TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_alpha));
434 TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_beta));
435 if (deriv_index == 0) {
436 // df/dx_n = df/dx_dot * dx_dot/dx_n + df/dx_n = 1/dt*df/dx_dot + df/dx_n
437 inArgs.set_alpha(Scalar(1.0)/dt);
438 inArgs.set_beta(Scalar(1.0));
439 }
440 else if (deriv_index == 1) {
441 // df/dx_{n-1} = df/dx_dot * dx_dot/dx_{n-1} = -1/dt*df/dx_dot
442 inArgs.set_alpha(Scalar(-1.0)/dt);
443 inArgs.set_beta(Scalar(0.0));
444 }
445 this->wrapperModel_->getAppModel()->evalModel(inArgs, outArgs);
446}
447
448
449// Nonmember constructor - ModelEvaluator and ParameterList
450// ------------------------------------------------------------------------
451template<class Scalar>
452Teuchos::RCP<StepperBackwardEuler<Scalar> >
454 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
455 Teuchos::RCP<Teuchos::ParameterList> pl)
456{
457 auto stepper = Teuchos::rcp(new StepperBackwardEuler<Scalar>());
458
459 stepper->setStepperImplicitValues(pl);
460
461 if (pl != Teuchos::null) {
462 std::string predictorName =
463 pl->get<std::string>("Predictor Stepper Type", "None");
464 stepper->setPredictor(predictorName);
465 }
466
467 if (model != Teuchos::null) {
468 stepper->setModel(model);
469 stepper->initialize();
470 }
471
472 return stepper;
473}
474
475
476} // namespace Tempus
477#endif // Tempus_StepperBackwardEuler_impl_hpp
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Application Action for StepperBackwardEuler.
Time-derivative interface for Backward Euler.
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState() override
Get a default (initial) StepperState.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void setAppAction(Teuchos::RCP< StepperBackwardEulerAppAction< Scalar > > appAction)
virtual int stencilLength() const override
Return the number of solution vectors in the time step stencil.
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.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
virtual void computePredictor(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Compute predictor given the supplied stepper.
void setPredictor(std::string predictorType="None")
Set the predictor.
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 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 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.
virtual bool isValidSetup(Teuchos::FancyOStream &out) const override
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Set the initial conditions and make them consistent.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Return a valid ParameterList with current settings.
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.
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Take the specified timestep, dt, and return true if successful.
Thyra Base interface for implicit time steppers.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Set the initial conditions and make them consistent.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
StepperState is a simple class to hold state information about the stepper.
Thyra Base interface for time steppers.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Teuchos::RCP< StepperBackwardEuler< Scalar > > createStepperBackwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.