Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_StepperExplicitRK_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_StepperExplicitRK_impl_hpp
10 #define Tempus_StepperExplicitRK_impl_hpp
11 
12 #include "Thyra_VectorStdOps.hpp"
13 
14 
16 
17 
18 namespace Tempus {
19 
20 
21 template<class Scalar>
23 {
24  this->setUseEmbedded(false);
25  this->setStageNumber(-1);
26  this->setAppAction(Teuchos::null);
27 }
28 
29 
30 template<class Scalar>
32  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
33  bool useFSAL,
34  std::string ICConsistency,
35  bool ICConsistencyCheck,
36  bool useEmbedded,
37  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
38 {
39  this->setUseFSAL( useFSAL);
40  this->setICConsistency( ICConsistency);
41  this->setICConsistencyCheck( ICConsistencyCheck);
42  this->setUseEmbedded( useEmbedded);
43  this->setStageNumber(-1);
44  this->setErrorNorm();
45 
46  this->setAppAction(stepperRKAppAction);
47 
48  if (appModel != Teuchos::null) {
49  this->setModel(appModel);
50  this->initialize();
51  }
52 }
53 
54 
55 template<class Scalar>
57  const Teuchos::RCP<SolutionHistory<Scalar> >& sh) const
58 {
59 
60  Scalar dt = Scalar(1.0e+99);
61  if (!this->getUseEmbedded()) return dt;
62 
63  Teuchos::RCP<SolutionState<Scalar> > currentState=sh->getCurrentState();
64  const int order = currentState->getOrder();
65  const Scalar time = currentState->getTime();
66  const Scalar errorRel = currentState->getTolRel();
67  const Scalar errorAbs = currentState->getTolAbs();
68 
69  Teuchos::RCP<Thyra::VectorBase<Scalar> > stageX, scratchX;
70  stageX = Thyra::createMember(this->appModel_->get_f_space());
71  scratchX = Thyra::createMember(this->appModel_->get_f_space());
72  Thyra::assign(stageX.ptr(), *(currentState->getX()));
73 
74  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar> > > stageXDot(2);
75  for (int i=0; i<2; ++i) {
76  stageXDot[i] = Thyra::createMember(this->appModel_->get_f_space());
77  assign(stageXDot[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
78  }
79 
80  // A: one function evaluation at F(t_0, X_0)
81  typedef Thyra::ModelEvaluatorBase MEB;
82  MEB::InArgs<Scalar> inArgs = this->appModel_->getNominalValues();
83  MEB::OutArgs<Scalar> outArgs = this->appModel_->createOutArgs();
84  inArgs.set_x(stageX);
85  if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
86  if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
87  outArgs.set_f(stageXDot[0]); // K1
88  this->appModel_->evalModel(inArgs, outArgs);
89 
90  this->stepperErrorNormCalculator_->setRelativeTolerance(errorRel);
91  this->stepperErrorNormCalculator_->setAbsoluteTolerance(errorAbs);
92 
93  Scalar d0 = this->stepperErrorNormCalculator_->errorNorm(stageX);
94  Scalar d1 = this->stepperErrorNormCalculator_->errorNorm(stageXDot[0]);
95 
96  // b) first guess for the step size
97  dt = Teuchos::as<Scalar>(0.01)*(d0/d1);
98 
99  // c) perform one explicit Euler step (X_1)
100  Thyra::Vp_StV(stageX.ptr(), dt, *(stageXDot[0]));
101 
102  // compute F(t_0 + dt, X_1)
103  inArgs.set_x(stageX);
104  if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time + dt);
105  if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
106  outArgs.set_f(stageXDot[1]); // K2
107  this->appModel_->evalModel(inArgs, outArgs);
108 
109  // d) compute estimate of the second derivative of the solution
110  // d2 = || f(t_0 + dt, X_1) - f(t_0, X_0) || / dt
112  errX = Thyra::createMember(this->appModel_->get_f_space());
113  assign(errX.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
114  Thyra::V_VmV(errX.ptr(), *(stageXDot[1]), *(stageXDot[0]));
115  Scalar d2 = this->stepperErrorNormCalculator_->errorNorm(errX) / dt;
116 
117  // e) compute step size h_1 (from m = 0 order Taylor series)
118  Scalar max_d1_d2 = std::max(d1, d2);
119  Scalar h1 = std::pow((0.01/max_d1_d2),(1.0/(order+1)));
120 
121  // f) propose starting step size
122  dt = std::min(100*dt, h1);
123  return dt;
124 }
125 
126 
127 template<class Scalar>
130 {
131  return this->getValidParametersBasicERK();
132 }
133 
134 
135 template<class Scalar>
138 {
139  auto pl = this->getValidParametersBasic();
140  pl->template set<bool>("Use Embedded", this->getUseEmbedded(),
141  "'Whether to use Embedded Stepper (if available) or not\n"
142  " 'true' - Stepper will compute embedded solution and is adaptive.\n"
143  " 'false' - Stepper is not embedded(adaptive).\n");
144  pl->template set<std::string>("Description", this->getDescription());
145 
146  return pl;
147 }
148 
149 
150 template<class Scalar>
152 {
153  TEUCHOS_TEST_FOR_EXCEPTION( this->tableau_ == Teuchos::null, std::logic_error,
154  "Error - Need to set the tableau, before calling "
155  "StepperExplicitRK::initialize()\n");
156 
157  TEUCHOS_TEST_FOR_EXCEPTION( this->appModel_==Teuchos::null, std::logic_error,
158  "Error - Need to set the model, setModel(), before calling "
159  "StepperExplicitRK::initialize()\n");
160 
162 }
163 
164 
165 template<class Scalar>
167  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
168 {
170 
171  // Set the stage vectors
172  int numStages = this->tableau_->numStages();
173  stageXDot_.resize(numStages);
174  for (int i=0; i<numStages; ++i) {
175  stageXDot_[i] = Thyra::createMember(this->appModel_->get_f_space());
176  assign(stageXDot_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
177  }
178 
179  this->setEmbeddedMemory();
180  this->setErrorNorm();
181 
182  this->isInitialized_ = false;
183 }
184 
185 
186 template<class Scalar>
188 {
189  if (this->getModel() == Teuchos::null)
190  return; // Embedded memory will be set when setModel() is called.
191 
192  if ( this->tableau_->isEmbedded() && this->getUseEmbedded() ){
193  this->ee_ = Thyra::createMember(this->appModel_->get_f_space());
194  this->abs_u0 = Thyra::createMember(this->appModel_->get_f_space());
195  this->abs_u = Thyra::createMember(this->appModel_->get_f_space());
196  this->sc = Thyra::createMember(this->appModel_->get_f_space());
197  } else {
198  this->ee_ = Teuchos::null;
199  this->abs_u0 = Teuchos::null;
200  this->abs_u = Teuchos::null;
201  this->sc = Teuchos::null;
202  }
203 }
204 
205 
206 template<class Scalar>
208  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
209 {
210  if (this->getUseFSAL())
211  this->setStepperXDot(stageXDot_.back());
212  else
213  this->setStepperXDot(stageXDot_[0]);
214 
216 
217  auto xDot = solutionHistory->getCurrentState()->getXDot();
218  if (xDot != Teuchos::null && this->getUseFSAL())
219  Thyra::assign(this->getStepperXDot().ptr(), *(xDot));
220 }
221 
222 
223 template<class Scalar>
225  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
226 {
227  this->checkInitialized();
228 
229  using Teuchos::RCP;
230 
231  TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperExplicitRK::takeStep()");
232  {
233  TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
234  std::logic_error,
235  "Error - StepperExplicitRK<Scalar>::takeStep(...)\n"
236  "Need at least two SolutionStates for ExplicitRK.\n"
237  " Number of States = " << solutionHistory->getNumStates() << "\n"
238  "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
239  " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
240 
241  RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
242  RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
243  const Scalar dt = workingState->getTimeStep();
244  const Scalar time = currentState->getTime();
245 
246  const int numStages = this->tableau_->numStages();
247  Teuchos::SerialDenseMatrix<int,Scalar> A = this->tableau_->A();
248  Teuchos::SerialDenseVector<int,Scalar> b = this->tableau_->b();
249  Teuchos::SerialDenseVector<int,Scalar> c = this->tableau_->c();
250 
251  Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
252 
253  RCP<StepperExplicitRK<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
254  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
256 
257  // Compute stage solutions
258  for (int i=0; i < numStages; ++i) {
259  this->setStageNumber(i);
260  Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
261  for (int j=0; j < i; ++j) {
262  if (A(i,j) != Teuchos::ScalarTraits<Scalar>::zero()) {
263  Thyra::Vp_StV(workingState->getX().ptr(), dt*A(i,j), *stageXDot_[j]);
264  }
265  }
266  this->setStepperXDot(stageXDot_[i]);
267 
268  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
270  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
272  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
274  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
276 
277  if ( i == 0 && this->getUseFSAL() &&
278  workingState->getNConsecutiveFailures() == 0 ) {
279  RCP<Thyra::VectorBase<Scalar> > tmp = stageXDot_[0];
280  stageXDot_[0] = stageXDot_.back();
281  stageXDot_.back() = tmp;
282  this->setStepperXDot(stageXDot_[0]);
283 
284  } else {
285  const Scalar ts = time + c(i)*dt;
287 
288  // Evaluate xDot = f(x,t).
289  this->evaluateExplicitODE(stageXDot_[i], workingState->getX(), ts, p);
290  }
291 
292  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
294  }
295  // reset the stage number
296  this->setStageNumber(-1);
297 
298  // Sum for solution: x_n = x_n-1 + Sum{ b(i) * dt*f(i) }
299  Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
300  for (int i=0; i < numStages; ++i) {
301  if (b(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
302  Thyra::Vp_StV((workingState->getX()).ptr(), dt*b(i), *(stageXDot_[i]));
303  }
304  }
305 
306  if (this->getUseFSAL()) {
307  if (numStages == 1) {
308  const Scalar ts = time + dt;
310  // Evaluate xDot = f(x,t).
311  this->evaluateExplicitODE(stageXDot_[0], workingState->getX(), ts, p);
312  }
313  if (workingState->getXDot() != Teuchos::null)
314  Thyra::assign((workingState->getXDot()).ptr(), *(stageXDot_.back()));
315  }
316 
317 
318  // At this point, the stepper has passed.
319  // But when using adaptive time stepping, the embedded method
320  // can change the step status
321  workingState->setSolutionStatus(Status::PASSED);
322 
323  if (this->tableau_->isEmbedded() && this->getUseEmbedded()) {
324 
325  const Scalar tolRel = workingState->getTolRel();
326  const Scalar tolAbs = workingState->getTolAbs();
327 
328  // update the tolerance
329  this->stepperErrorNormCalculator_->setRelativeTolerance(tolRel);
330  this->stepperErrorNormCalculator_->setAbsoluteTolerance(tolAbs);
331 
332  // just compute the error weight vector
333  // (all that is needed is the error, and not the embedded solution)
335  errWght -= this->tableau_->bstar();
336 
337  // Compute local truncation error estimate: | u^{n+1} - \hat{u}^{n+1} |
338  // Sum for solution: ee_n = Sum{ (b(i) - bstar(i)) * dt*f(i) }
339  assign(this->ee_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
340  for (int i=0; i < numStages; ++i) {
341  if (errWght(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
342  Thyra::Vp_StV(this->ee_.ptr(), dt*errWght(i), *(stageXDot_[i]));
343  }
344  }
345 
346 
347  Scalar err = this->stepperErrorNormCalculator_->computeWRMSNorm(currentState->getX(), workingState->getX(), this->ee_);
348  workingState->setErrorRel(err);
349 
350  // Test if step should be rejected
351  if (std::isinf(err) || std::isnan(err) || err > Teuchos::as<Scalar>(1.0))
352  workingState->setSolutionStatus(Status::FAILED);
353  }
354 
355  workingState->setOrder(this->getOrder());
356  workingState->computeNorms(currentState);
357  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
359  }
360  return;
361 }
362 
363 
370 template<class Scalar>
373 {
375  rcp(new StepperState<Scalar>(this->getStepperType()));
376  return stepperState;
377 }
378 
379 
380 template<class Scalar>
383  const Teuchos::EVerbosityLevel verbLevel) const
384 {
385  out.setOutputToRootOnly(0);
386  out << std::endl;
387  Stepper<Scalar>::describe(out, verbLevel);
388  StepperExplicit<Scalar>::describe(out, verbLevel);
389 
390  out << "--- StepperExplicitRK ---\n";
391  if (this->tableau_ != Teuchos::null) this->tableau_->describe(out, verbLevel);
392  out << " tableau_ = " << this->tableau_ << std::endl;
393  out << " stepperRKAppAction_= " << this->stepperRKAppAction_ << std::endl;
394  out << " stageXDot_.size() = " << stageXDot_.size() << std::endl;
395  const int numStages = stageXDot_.size();
396  for (int i=0; i<numStages; ++i)
397  out << " stageXDot_["<<i<<"] = " << stageXDot_[i] << std::endl;
398  out << " useEmbedded_ = "
399  << Teuchos::toString(this->useEmbedded_) << std::endl;
400  out << " ee_ = " << this->ee_ << std::endl;
401  out << " abs_u0 = " << this->abs_u0 << std::endl;
402  out << " abs_u = " << this->abs_u << std::endl;
403  out << " sc = " << this->sc << std::endl;
404  out << "-------------------------" << std::endl;
405 }
406 
407 
408 template<class Scalar>
410 {
411  out.setOutputToRootOnly(0);
412  bool isValidSetup = true;
413 
414  if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
415  if ( !StepperExplicit<Scalar>::isValidSetup(out) ) isValidSetup = false;
416 
417  if (this->tableau_ == Teuchos::null) {
418  isValidSetup = false;
419  out << "The tableau is not set!\n";
420  }
421 
422  if (this->stepperRKAppAction_ == Teuchos::null) {
423  isValidSetup = false;
424  out << "The AppAction is not set!\n";
425  }
426 
427  return isValidSetup;
428 }
429 
430 
431 } // namespace Tempus
432 #endif // Tempus_StepperExplicitRK_impl_hpp
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Set model.
virtual void setupDefault()
Default setup for constructor.
virtual Scalar getInitTimeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void initialize()
Initialize after construction and changing input parameters.
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
Thyra Base interface for time steppers.
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)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Application Action for StepperRKBase.
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Ptr< T > ptr() const
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions, make them consistent, and set needed memory.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Set model.
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasicERK() const
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
virtual void initialize()
Initialize during construction and after changing input parameters.
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
Thyra Base interface for implicit time steppers.
std::string toString(const T &t)
virtual void setup(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Setup for constructor.