Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_StepperDIRK_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_StepperDIRK_impl_hpp
10 #define Tempus_StepperDIRK_impl_hpp
11 
12 #include "Thyra_VectorStdOps.hpp"
13 
14 #include "Tempus_WrapperModelEvaluatorBasic.hpp"
15 
16 
17 namespace Tempus {
18 
19 
20 template<class Scalar>
22 {
23  this->setUseEmbedded( false);
24  this->setZeroInitialGuess( false);
25  this->setStageNumber(-1);
26 
27  this->setAppAction(Teuchos::null);
28  this->setDefaultSolver();
29 }
30 
31 
32 template<class Scalar>
34  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
36  bool useFSAL,
37  std::string ICConsistency,
38  bool ICConsistencyCheck,
39  bool useEmbedded,
40  bool zeroInitialGuess,
41  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
42 {
43  this->setUseFSAL( useFSAL);
44  this->setICConsistency( ICConsistency);
45  this->setICConsistencyCheck( ICConsistencyCheck);
46  this->setUseEmbedded( useEmbedded);
47  this->setZeroInitialGuess( zeroInitialGuess);
48 
49  this->setStageNumber(-1);
50  this->setErrorNorm();
51  this->setAppAction(stepperRKAppAction);
52  this->setSolver(solver);
53 
54  if (appModel != Teuchos::null) {
55  this->setModel(appModel);
56  this->initialize();
57  }
58 }
59 
60 
61 template<class Scalar>
64 {
65  auto pl = this->getValidParametersBasicImplicit();
66  pl->template set<bool>("Use Embedded", this->getUseEmbedded(),
67  "'Whether to use Embedded Stepper (if available) or not\n"
68  " 'true' - Stepper will compute embedded solution and is adaptive.\n"
69  " 'false' - Stepper is not embedded(adaptive).\n");
70  pl->template set<std::string>("Description", this->getDescription());
71  pl->template set<bool>("Reset Initial Guess", this->getResetInitialGuess());
72 
73  return pl;
74 }
75 
76 
77 template<class Scalar>
79 {
81  this->tableau_ == Teuchos::null, std::logic_error,
82  "Error - Need to set the tableau, before calling "
83  "StepperDIRK::initialize()\n");
84 
86  this->wrapperModel_==Teuchos::null, std::logic_error,
87  "Error - Need to set the model, setModel(), before calling "
88  "StepperDIRK::initialize()\n");
89 
91 }
92 
93 
94 template<class Scalar>
96  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
97 {
99 
100  // Set the stage vectors
101  const int numStages = this->tableau_->numStages();
102  stageXDot_.resize(numStages);
103  for (int i=0; i<numStages; ++i) {
104  stageXDot_[i] = Thyra::createMember(this->wrapperModel_->get_f_space());
105  assign(stageXDot_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
106  }
107  xTilde_ = Thyra::createMember(this->wrapperModel_->get_x_space());
108  assign(xTilde_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
109 
110  this->setEmbeddedMemory();
111  this->setErrorNorm();
112 
113  this->isInitialized_ = false;
114 }
115 
116 
117 template<class Scalar>
119 {
120  if (this->getModel() == Teuchos::null)
121  return; // Embedded memory will be set when setModel() is called.
122 
123  if (this->tableau_->isEmbedded() && this->getUseEmbedded()) {
124  this->ee_ = Thyra::createMember(this->wrapperModel_->get_f_space());
125  this->abs_u0 = Thyra::createMember(this->wrapperModel_->get_f_space());
126  this->abs_u = Thyra::createMember(this->wrapperModel_->get_f_space());
127  this->sc = Thyra::createMember(this->wrapperModel_->get_f_space());
128  } else {
129  this->ee_ = Teuchos::null;
130  this->abs_u0 = Teuchos::null;
131  this->abs_u = Teuchos::null;
132  this->sc = Teuchos::null;
133  }
134 }
135 
136 
137 template<class Scalar>
139  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
140 {
141  this->setStepperXDot(stageXDot_.back());
143 }
144 
145 
146 template<class Scalar>
148  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
149 {
150  this->checkInitialized();
151 
152  using Teuchos::RCP;
153 
154  TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperDIRK::takeStep()");
155  {
156  TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
157  std::logic_error,
158  "Error - StepperDIRK<Scalar>::takeStep(...)\n"
159  "Need at least two SolutionStates for DIRK.\n"
160  " Number of States = " << solutionHistory->getNumStates() << "\n"
161  "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
162  " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
163 
164  RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
165  RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
166  const Scalar dt = workingState->getTimeStep();
167  const Scalar time = currentState->getTime();
168 
169  const int numStages = this->tableau_->numStages();
170  Teuchos::SerialDenseMatrix<int,Scalar> A = this->tableau_->A();
171  Teuchos::SerialDenseVector<int,Scalar> b = this->tableau_->b();
172  Teuchos::SerialDenseVector<int,Scalar> c = this->tableau_->c();
173 
174  // Reset non-zero initial guess.
175  if ( this->getResetInitialGuess() && (!this->getZeroInitialGuess()) )
176  Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
177 
178  RCP<StepperDIRK<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
179  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
181 
182  // Compute stage solutions
183  bool pass = true;
185  for (int i=0; i < numStages; ++i) {
186  this->setStageNumber(i);
187 
188  Thyra::assign(xTilde_.ptr(), *(currentState->getX()));
189  for (int j=0; j < i; ++j) {
190  if (A(i,j) != Teuchos::ScalarTraits<Scalar>::zero()) {
191  Thyra::Vp_StV(xTilde_.ptr(), dt*A(i,j), *(stageXDot_[j]));
192  }
193  }
194  this->setStepperXDot(stageXDot_[i]);
195 
196  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
198 
199  Scalar ts = time + c(i)*dt;
200 
201  // Check if stageXDot_[i] is needed.
202  bool isNeeded = false;
203  for (int k=i+1; k<numStages; ++k) if (A(k,i) != 0.0) isNeeded = true;
204  if (b(i) != 0.0) isNeeded = true;
205  if (this->tableau_->isEmbedded() && this->getUseEmbedded() &&
206  this->tableau_->bstar()(i) != 0.0)
207  isNeeded = true;
208  if (isNeeded == false) {
209  assign(stageXDot_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
210  } else if (A(i,i) == Teuchos::ScalarTraits<Scalar>::zero()) {
211  // Explicit stage for the ImplicitODE_DAE
212  if (i == 0 && this->getUseFSAL() &&
213  workingState->getNConsecutiveFailures() == 0) {
214  // Reuse last evaluation for first step
215  RCP<Thyra::VectorBase<Scalar> > tmp = stageXDot_[0];
216  stageXDot_[0] = stageXDot_.back();
217  stageXDot_.back() = tmp;
218  this->setStepperXDot(stageXDot_[0]);
219  } else {
220  // Calculate explicit stage
221  typedef Thyra::ModelEvaluatorBase MEB;
222  MEB::InArgs<Scalar> inArgs = this->wrapperModel_->getInArgs();
223  MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
224  inArgs.set_x(xTilde_);
225  if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(ts);
226  if (inArgs.supports(MEB::IN_ARG_x_dot))
227  inArgs.set_x_dot(Teuchos::null);
228  outArgs.set_f(stageXDot_[i]);
229 
230  this->wrapperModel_->getAppModel()->evalModel(inArgs,outArgs);
231  }
232  } else {
233  // Implicit stage for the ImplicitODE_DAE
234  const Scalar alpha = 1.0/(dt*A(i,i));
235  const Scalar beta = 1.0;
236 
237  // Setup TimeDerivative
240  alpha,xTilde_.getConst()));
241 
243  timeDer, dt, alpha, beta, SOLVE_FOR_X, i));
244 
245  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
247 
248  sStatus = this->solveImplicitODE(workingState->getX(), stageXDot_[i], ts, p);
249 
250  if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass=false;
251 
252  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
254 
255  timeDer->compute(workingState->getX(), stageXDot_[i]);
256  }
257  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
259  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
261  }
262  // reset the stage number
263  this->setStageNumber(-1);
264 
265  // Sum for solution: x_n = x_n-1 + Sum{ dt*b(i) * f(i) }
266  Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
267  for (int i=0; i < numStages; ++i) {
268  if (b(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
269  Thyra::Vp_StV((workingState->getX()).ptr(), dt*b(i), *(stageXDot_[i]));
270  }
271  }
272 
273  if (this->tableau_->isEmbedded() && this->getUseEmbedded()) {
274  const Scalar tolRel = workingState->getTolRel();
275  const Scalar tolAbs = workingState->getTolAbs();
276 
277  // update the tolerance
278  this->stepperErrorNormCalculator_->setRelativeTolerance(tolRel);
279  this->stepperErrorNormCalculator_->setAbsoluteTolerance(tolAbs);
280 
281  // just compute the error weight vector
282  // (all that is needed is the error, and not the embedded solution)
284  errWght -= this->tableau_->bstar();
285 
286  // compute local truncation error estimate: | u^{n+1} - \hat{u}^{n+1} |
287  // Sum for solution: ee_n = Sum{ (b(i) - bstar(i)) * dt*f(i) }
288  assign(this->ee_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
289  for (int i=0; i < numStages; ++i) {
290  if (errWght(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
291  Thyra::Vp_StV(this->ee_.ptr(), dt*errWght(i), *(stageXDot_[i]));
292  }
293  }
294 
295  Scalar err = this->stepperErrorNormCalculator_->computeWRMSNorm(currentState->getX(), workingState->getX(), this->ee_);
296  workingState->setErrorRel(err);
297 
298  // test if step should be rejected
299  if (std::isinf(err) || std::isnan(err) || err > Teuchos::as<Scalar>(1.0))
300  pass = false;
301  }
302 
303  if (pass) workingState->setSolutionStatus(Status::PASSED);
304  else workingState->setSolutionStatus(Status::FAILED);
305 
306  workingState->setOrder(this->getOrder());
307  workingState->computeNorms(currentState);
308  this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
310  }
311  return;
312 }
313 
320 template<class Scalar>
324 {
326  rcp(new StepperState<Scalar>(this->getStepperType()));
327  return stepperState;
328 }
329 
330 
331 template<class Scalar>
334  const Teuchos::EVerbosityLevel verbLevel) const
335 {
336  out.setOutputToRootOnly(0);
337  out << std::endl;
338  Stepper<Scalar>::describe(out, verbLevel);
339  StepperImplicit<Scalar>::describe(out, verbLevel);
340 
341  out << "--- StepperDIRK ---\n";
342  out << " tableau_ = " << this->tableau_ << std::endl;
343  if (this->tableau_ != Teuchos::null) this->tableau_->describe(out, verbLevel);
344  out << " stepperRKAppAction_ = " << this->stepperRKAppAction_ << std::endl;
345  out << " xTilde_ = " << xTilde_ << std::endl;
346  out << " stageXDot_.size() = " << stageXDot_.size() << std::endl;
347  const int numStages = stageXDot_.size();
348  for (int i=0; i<numStages; ++i)
349  out << " stageXDot_["<<i<<"] = " << stageXDot_[i] << std::endl;
350  out << " useEmbedded_ = "
351  << Teuchos::toString(this->useEmbedded_) << std::endl;
352  out << " ee_ = " << this->ee_ << std::endl;
353  out << " abs_u0 = " << this->abs_u0 << std::endl;
354  out << " abs_u = " << this->abs_u << std::endl;
355  out << " sc = " << this->sc << std::endl;
356  out << "-------------------" << std::endl;
357 }
358 
359 
360 template<class Scalar>
362 {
363  out.setOutputToRootOnly(0);
364  bool isValidSetup = true;
365 
366  if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
367  if ( !StepperImplicit<Scalar>::isValidSetup(out) ) isValidSetup = false;
368 
369  if (this->tableau_ == Teuchos::null) {
370  isValidSetup = false;
371  out << "The tableau is not set!\n";
372  }
373 
374  if (this->stepperRKAppAction_ == Teuchos::null) {
375  isValidSetup = false;
376  out << "The AppAction is not set!\n";
377  }
378 
379  return isValidSetup;
380 }
381 
382 
383 template<class Scalar>
386 {
387  return this->getValidParametersBasicDIRK();
388 }
389 
390 
391 } // namespace Tempus
392 #endif // Tempus_StepperDIRK_impl_hpp
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasicDIRK() const
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void initialize()
Initialize after construction and changing input parameters.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
Thyra Base interface for time steppers.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
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
Application Action for StepperRKBase.
Time-derivative interface for DIRK.
virtual void setupDefault()
Default setup for constructor.
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
virtual bool isValidSetup(Teuchos::FancyOStream &out) const override
ESolveStatus solveStatus
virtual void initialize() override
Initialize after construction and changing input parameters.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Set the initial conditions and make them consistent.
virtual void setEmbeddedMemory() override
virtual void setup(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &wrapperModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Setup for constructor.
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState() override
Get a default (initial) StepperState.
Solve for x and determine xDot from x.
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Take the specified timestep, dt, and return true if successful.
std::string toString(const T &t)