Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
17 namespace Tempus {
18 
19 
20 template<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 
36 template<class Scalar>
38  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
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 
65 template<class Scalar>
66 void 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 
87 template<class Scalar>
89  Teuchos::RCP<Stepper<Scalar> > predictorStepper)
90 {
91  predictorStepper_ = predictorStepper;
92  if (predictorStepper_ == Teuchos::null) return;
93 
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 
108 template<class Scalar>
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 
123 template<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 
141 template<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 
159 template<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
200  Scalar(1.0)/dt,xOld));
201 
202  const Scalar alpha = Scalar(1.0)/dt;
203  const Scalar beta = Scalar(1.0);
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 
228 template<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 
252 template<class Scalar>
256 {
258  rcp(new StepperState<Scalar>(this->getStepperType()));
259  return stepperState;
260 }
261 
262 
263 template<class Scalar>
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 
288 template<class Scalar>
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 
313 template<class Scalar>
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 
326 template <class Scalar>
327 int
329 {
330  return 2;
331 }
332 
333 
334 template <class Scalar>
335 void
337  Thyra::VectorBase<Scalar>& residual,
339  const Teuchos::Array<Scalar>& t,
340  const Thyra::VectorBase<Scalar>& p,
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 
349 template <class Scalar>
350 void
352  Thyra::LinearOpBase<Scalar>& jacobian,
354  const Teuchos::Array<Scalar>& t,
355  const Thyra::VectorBase<Scalar>& p,
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 
366 template <class Scalar>
367 void
371  const Teuchos::Array<Scalar>& t,
372  const Thyra::VectorBase<Scalar>& p,
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 
383 template <class Scalar>
384 void
386  Thyra::LinearOpWithSolveBase<Scalar>& jacobian_solver,
388  const Teuchos::Array<Scalar>& t,
389  const Thyra::VectorBase<Scalar>& p,
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 
399 template <class Scalar>
400 void
404  const Teuchos::Array<Scalar>& t,
405  const Thyra::VectorBase<Scalar>& p,
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();
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 // ------------------------------------------------------------------------
451 template<class Scalar>
454  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
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
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.
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.
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.
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.
size_type size() const
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.