Tempus  Version of the Day
Time Integration
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 "Tempus_RKButcherTableauBuilder.hpp"
13 #include "Tempus_config.hpp"
14 #include "Tempus_StepperFactory.hpp"
15 #include "Tempus_WrapperModelEvaluatorBasic.hpp"
16 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
17 #include "Thyra_VectorStdOps.hpp"
18 #include "NOX_Thyra.H"
19 
20 
21 namespace Tempus {
22 
23 // Forward Declaration for recursive includes (this Stepper <--> StepperFactory)
24 template<class Scalar> class StepperFactory;
25 
26 template<class Scalar>
27 StepperDIRK<Scalar>::StepperDIRK(
28  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
29  std::string stepperType)
30 {
31  this->setTableau(stepperType);
32  this->setModel(appModel);
33  this->initialize();
34 }
35 
36 template<class Scalar>
37 StepperDIRK<Scalar>::StepperDIRK(
38  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
39  Teuchos::RCP<Teuchos::ParameterList> pList)
40 {
41  this->setTableau(pList);
42  this->setParameterList(pList);
43  this->setModel(appModel);
44  this->initialize();
45 }
46 
47 template<class Scalar>
48 StepperDIRK<Scalar>::StepperDIRK(
49  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
50  std::string stepperType,
51  Teuchos::RCP<Teuchos::ParameterList> pList)
52 {
53  this->setTableau(stepperType);
54  this->setParameterList(pList);
55  this->setModel(appModel);
56  this->initialize();
57 }
58 
59 
60 template<class Scalar>
61 void StepperDIRK<Scalar>::setTableau(std::string stepperType)
62 {
63  if (stepperType == "") {
64  this->setTableau();
65  } else {
66  DIRK_ButcherTableau_ = createRKBT<Scalar>(stepperType, this->stepperPL_);
67  }
68 
69  TEUCHOS_TEST_FOR_EXCEPTION( DIRK_ButcherTableau_->isDIRK() != true,
70  std::logic_error,
71  "Error - StepperDIRK did not receive a DIRK Butcher Tableau!\n"
72  << " Stepper Type = " << stepperType << "\n");
73  description_ = DIRK_ButcherTableau_->description();
74 }
75 
76 
77 template<class Scalar>
78 void StepperDIRK<Scalar>::setTableau(Teuchos::RCP<Teuchos::ParameterList> pList)
79 {
80  if (pList == Teuchos::null) {
81  // Create default parameters if null, otherwise keep current parameters.
82  if (this->stepperPL_ == Teuchos::null)
83  this->stepperPL_ = this->getDefaultParameters();
84  } else {
85  this->stepperPL_ = pList;
86  }
87 
88  std::string stepperType =
89  this->stepperPL_->template get<std::string>("Stepper Type",
90  "SDIRK 2 Stage 2nd order");
91  DIRK_ButcherTableau_ = createRKBT<Scalar>(stepperType, this->stepperPL_);
92 
93  TEUCHOS_TEST_FOR_EXCEPTION( DIRK_ButcherTableau_->isDIRK() != true,
94  std::logic_error,
95  "Error - StepperDIRK did not receive a DIRK Butcher Tableau!\n"
96  << " Stepper Type = " << stepperType << "\n");
97  description_ = DIRK_ButcherTableau_->description();
98 }
99 
100 
101 template<class Scalar>
102 void StepperDIRK<Scalar>::setObserver(
103  Teuchos::RCP<StepperObserver<Scalar> > obs)
104 {
105  if (obs == Teuchos::null) {
106  // Create default observer, otherwise keep current observer.
107  if (stepperObserver_ == Teuchos::null) {
108  stepperDIRKObserver_ =
109  Teuchos::rcp(new StepperDIRKObserver<Scalar>());
110  stepperObserver_ =
111  Teuchos::rcp_dynamic_cast<StepperObserver<Scalar> >
112  (stepperDIRKObserver_);
113  }
114  } else {
115  stepperObserver_ = obs;
116  stepperDIRKObserver_ =
117  Teuchos::rcp_dynamic_cast<StepperDIRKObserver<Scalar> >(stepperObserver_);
118  }
119 }
120 
121 
122 template<class Scalar>
123 void StepperDIRK<Scalar>::initialize()
124 {
125  TEUCHOS_TEST_FOR_EXCEPTION( this->wrapperModel_ == Teuchos::null,
126  std::logic_error,
127  "Error - Need to set the model, setModel(), before calling "
128  "StepperDIRK::initialize()\n");
129 
130  this->setTableau(this->stepperPL_);
131  this->setParameterList(this->stepperPL_);
132  this->setSolver();
133  this->setObserver();
134 
135  // Initialize the stage vectors
136  const int numStages = DIRK_ButcherTableau_->numStages();
137  stageX_ = this->wrapperModel_->getNominalValues().get_x()->clone_v();
138  stageXDot_.resize(numStages);
139  for (int i=0; i<numStages; ++i) {
140  stageXDot_[i] = Thyra::createMember(this->wrapperModel_->get_f_space());
141  assign(stageXDot_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
142  }
143  xTilde_ = Thyra::createMember(this->wrapperModel_->get_x_space());
144  assign(xTilde_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
145 
146  if (DIRK_ButcherTableau_->isEmbedded() and this->getEmbedded()) {
147  ee_ = Thyra::createMember(this->wrapperModel_->get_f_space());
148  abs_u0 = Thyra::createMember(this->wrapperModel_->get_f_space());
149  abs_u = Thyra::createMember(this->wrapperModel_->get_f_space());
150  sc = Thyra::createMember(this->wrapperModel_->get_f_space());
151  }
152 }
153 
154 template<class Scalar>
155 void StepperDIRK<Scalar>::takeStep(
156  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
157 {
158  using Teuchos::RCP;
159 
160  TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperDIRK::takeStep()");
161  {
162  TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
163  std::logic_error,
164  "Error - StepperDIRK<Scalar>::takeStep(...)\n"
165  "Need at least two SolutionStates for DIRK.\n"
166  " Number of States = " << solutionHistory->getNumStates() << "\n"
167  "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
168  " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
169 
170  stepperObserver_->observeBeginTakeStep(solutionHistory, *this);
171  RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
172  RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
173  const Scalar dt = workingState->getTimeStep();
174  const Scalar time = currentState->getTime();
175 
176  const int numStages = DIRK_ButcherTableau_->numStages();
177  Teuchos::SerialDenseMatrix<int,Scalar> A = DIRK_ButcherTableau_->A();
178  Teuchos::SerialDenseVector<int,Scalar> b = DIRK_ButcherTableau_->b();
179  Teuchos::SerialDenseVector<int,Scalar> c = DIRK_ButcherTableau_->c();
180 
181  // Compute stage solutions
182  bool pass = true;
183  Thyra::SolveStatus<Scalar> sStatus;
184  for (int i=0; i < numStages; ++i) {
185  if (!Teuchos::is_null(stepperDIRKObserver_))
186  stepperDIRKObserver_->observeBeginStage(solutionHistory, *this);
187  Thyra::assign(xTilde_.ptr(), *(currentState->getX()));
188  for (int j=0; j < i; ++j) {
189  if (A(i,j) != Teuchos::ScalarTraits<Scalar>::zero()) {
190  Thyra::Vp_StV(xTilde_.ptr(), dt*A(i,j), *(stageXDot_[j]));
191  }
192  }
193 
194  Scalar ts = time + c(i)*dt;
195  if (A(i,i) == Teuchos::ScalarTraits<Scalar>::zero()) {
196  // Explicit stage for the ImplicitODE_DAE
197  bool isNeeded = false;
198  for (int k=i+1; k<numStages; ++k) if (A(k,i) != 0.0) isNeeded = true;
199  if (b(i) != 0.0) isNeeded = true;
200  if (isNeeded == false) {
201  // stageXDot_[i] is not needed.
202  assign(stageXDot_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
203  } else {
204  typedef Thyra::ModelEvaluatorBase MEB;
205  MEB::InArgs<Scalar> inArgs = this->wrapperModel_->getInArgs();
206  MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
207  inArgs.set_x(xTilde_);
208  if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(ts);
209  if (inArgs.supports(MEB::IN_ARG_x_dot))
210  inArgs.set_x_dot(Teuchos::null);
211  outArgs.set_f(stageXDot_[i]);
212 
213  if (!Teuchos::is_null(stepperDIRKObserver_))
214  stepperDIRKObserver_->observeBeforeExplicit(solutionHistory,*this);
215  this->wrapperModel_->getAppModel()->evalModel(inArgs,outArgs);
216  }
217  } else {
218  // Implicit stage for the ImplicitODE_DAE
219  Scalar alpha = 1.0/(dt*A(i,i));
220 
221  // Setup TimeDerivative
222  Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
223  Teuchos::rcp(new StepperDIRKTimeDerivative<Scalar>(
224  alpha,xTilde_.getConst()));
225 
226  // Setup InArgs and OutArgs
227  typedef Thyra::ModelEvaluatorBase MEB;
228  MEB::InArgs<Scalar> inArgs = this->wrapperModel_->getInArgs();
229  MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
230  inArgs.set_x(stageX_);
231  if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(stageXDot_[i]);
232  if (inArgs.supports(MEB::IN_ARG_t )) inArgs.set_t (ts);
233  if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(dt);
234  if (inArgs.supports(MEB::IN_ARG_alpha )) inArgs.set_alpha (alpha);
235  if (inArgs.supports(MEB::IN_ARG_beta )) inArgs.set_beta (1.0);
236 
237  this->wrapperModel_->setForSolve(timeDer, inArgs, outArgs);
238 
239  if (!Teuchos::is_null(stepperDIRKObserver_))
240  stepperDIRKObserver_->observeBeforeSolve(solutionHistory, *this);
241 
242  sStatus = this->solveImplicitODE(stageX_);
243 
244  if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED ) pass=false;
245 
246  if (!Teuchos::is_null(stepperDIRKObserver_))
247  stepperDIRKObserver_->observeAfterSolve(solutionHistory, *this);
248 
249  timeDer->compute(stageX_, stageXDot_[i]);
250  }
251  if (!Teuchos::is_null(stepperDIRKObserver_))
252  stepperDIRKObserver_->observeEndStage(solutionHistory, *this);
253  }
254 
255  // Sum for solution: x_n = x_n-1 + Sum{ dt*b(i) * f(i) }
256  Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
257  for (int i=0; i < numStages; ++i) {
258  if (b(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
259  Thyra::Vp_StV((workingState->getX()).ptr(), dt*b(i), *(stageXDot_[i]));
260  }
261  }
262 
263  if (DIRK_ButcherTableau_->isEmbedded() and this->getEmbedded()) {
264  RCP<SolutionStateMetaData<Scalar> > metaData = workingState->getMetaData();
265  const Scalar tolAbs = metaData->getTolRel();
266  const Scalar tolRel = metaData->getTolAbs();
267 
268  // just compute the error weight vector
269  // (all that is needed is the error, and not the embedded solution)
270  Teuchos::SerialDenseVector<int,Scalar> errWght = b ;
271  errWght -= DIRK_ButcherTableau_->bstar();
272 
273  //compute local truncation error estimate: | u^{n+1} - \hat{u}^{n+1} |
274  // Sum for solution: ee_n = Sum{ (b(i) - bstar(i)) * dt*f(i) }
275  assign(ee_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
276  for (int i=0; i < numStages; ++i) {
277  if (errWght(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
278  Thyra::Vp_StV(ee_.ptr(), dt*errWght(i), *(stageXDot_[i]));
279  }
280  }
281 
282  // compute: Atol + max(|u^n|, |u^{n+1}| ) * Rtol
283  Thyra::abs( *(currentState->getX()), abs_u0.ptr());
284  Thyra::abs( *(workingState->getX()), abs_u.ptr());
285  Thyra::pair_wise_max_update(tolRel, *abs_u0, abs_u.ptr());
286  Thyra::add_scalar(tolAbs, abs_u.ptr());
287 
288  //compute: || ee / sc ||
289  assign(sc.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
290  Thyra::ele_wise_divide(Teuchos::as<Scalar>(1.0), *ee_, *abs_u, sc.ptr());
291  Scalar err = Thyra::norm_inf(*sc);
292  metaData->setErrorRel(err);
293 
294  // test if step should be rejected
295  if (err > 1.0) workingState->setSolutionStatus(Status::FAILED);
296  }
297 
298  if (pass == true) workingState->setSolutionStatus(Status::PASSED);
299  else workingState->setSolutionStatus(Status::FAILED);
300  workingState->setOrder(this->getOrder());
301  stepperObserver_->observeEndTakeStep(solutionHistory, *this);
302  }
303  return;
304 }
305 
306 /** \brief Provide a StepperState to the SolutionState.
307  * This Stepper does not have any special state data,
308  * so just provide the base class StepperState with the
309  * Stepper description. This can be checked to ensure
310  * that the input StepperState can be used by this Stepper.
311  */
312 template<class Scalar>
313 Teuchos::RCP<Tempus::StepperState<Scalar> >
314 StepperDIRK<Scalar>::
315 getDefaultStepperState()
316 {
317  Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
318  rcp(new StepperState<Scalar>(description()));
319  return stepperState;
320 }
321 
322 
323 template<class Scalar>
324 std::string StepperDIRK<Scalar>::description() const
325 {
326  return(description_);
327 }
328 
329 
330 template<class Scalar>
331 void StepperDIRK<Scalar>::describe(
332  Teuchos::FancyOStream &out,
333  const Teuchos::EVerbosityLevel verbLevel) const
334 {
335  out << description() << "::describe:" << std::endl
336  << "wrapperModel_ = " << this->wrapperModel_->description() << std::endl;
337 }
338 
339 
340 template <class Scalar>
341 void StepperDIRK<Scalar>::setParameterList(
342  const Teuchos::RCP<Teuchos::ParameterList> & pList)
343 {
344  if (pList == Teuchos::null) {
345  // Create default parameters if null, otherwise keep current parameters.
346  if (this->stepperPL_ == Teuchos::null) this->stepperPL_ = this->getDefaultParameters();
347  } else {
348  this->stepperPL_ = pList;
349  }
350  // Can not validate because of optional Parameters.
351  //stepperPL_->validateParametersAndSetDefaults(*this->getValidParameters());
352 }
353 
354 
355 template<class Scalar>
356 Teuchos::RCP<const Teuchos::ParameterList>
357 StepperDIRK<Scalar>::getValidParameters() const
358 {
359  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
360  *pl = *(DIRK_ButcherTableau_->getValidParameters());
361  pl->set<bool> ("Zero Initial Guess", false);
362  return pl;
363 }
364 
365 template <class Scalar>
366 Teuchos::RCP<Teuchos::ParameterList>
367 StepperDIRK<Scalar>::getDefaultParameters() const
368 {
369  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
370  if (DIRK_ButcherTableau_ == Teuchos::null) {
371  auto DIRK_ButcherTableau =
372  createRKBT<Scalar>("SDIRK 2 Stage 2nd order", Teuchos::null);
373  pl->setParameters(*(DIRK_ButcherTableau->getValidParameters()));
374  } else {
375  pl->setParameters(*(DIRK_ButcherTableau_->getValidParameters()));
376  }
377  pl->set<std::string>("Solver Name", "Default Solver");
378  pl->set<bool> ("Zero Initial Guess", false);
379  Teuchos::RCP<Teuchos::ParameterList> solverPL=this->defaultSolverParameters();
380  pl->set("Default Solver", *solverPL);
381 
382  return pl;
383 }
384 
385 template <class Scalar>
386 Teuchos::RCP<Teuchos::ParameterList>
387 StepperDIRK<Scalar>::getNonconstParameterList()
388 {
389  return(this->stepperPL_);
390 }
391 
392 
393 template <class Scalar>
394 Teuchos::RCP<Teuchos::ParameterList>
395 StepperDIRK<Scalar>::unsetParameterList()
396 {
397  Teuchos::RCP<Teuchos::ParameterList> temp_plist = this->stepperPL_;
398  this->stepperPL_ = Teuchos::null;
399  return(temp_plist);
400 }
401 
402 
403 } // namespace Tempus
404 #endif // Tempus_StepperDIRK_impl_hpp
Tempus::StepperDIRK::StepperDIRK
StepperDIRK()
Default Constructor – not allowed.
Tempus::StepperState
StepperState is a simple class to hold state information about the stepper.
Definition: Tempus_StepperState.hpp:36
Tempus::StepperDIRK::describe
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Definition: Tempus_StepperDIRK_impl.hpp:331
Tempus::StepperDIRK::unsetParameterList
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Definition: Tempus_StepperDIRK_impl.hpp:395
Tempus::solutionHistory
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
Definition: Tempus_SolutionHistory_impl.hpp:504
Tempus::StepperDIRK::takeStep
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
Definition: Tempus_StepperDIRK_impl.hpp:155
Tempus::StepperDIRKTimeDerivative
Time-derivative interface for DIRK.
Definition: Tempus_StepperDIRK_decl.hpp:204
Tempus
Definition: Tempus_AdjointAuxSensitivityModelEvaluator_decl.hpp:20
Tempus::StepperDIRK::description
virtual std::string description() const
Definition: Tempus_StepperDIRK_impl.hpp:324
Tempus::StepperDIRK::getNonconstParameterList
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Definition: Tempus_StepperDIRK_impl.hpp:387
Tempus::StepperDIRKObserver
StepperDIRKObserver class for StepperDIRK.
Definition: Tempus_StepperDIRKObserver.hpp:35
Tempus::StepperObserver
StepperObserver class for Stepper class.
Definition: Tempus_StepperObserver.hpp:38
Tempus::StepperDIRK::getDefaultParameters
Teuchos::RCP< Teuchos::ParameterList > getDefaultParameters() const
Definition: Tempus_StepperDIRK_impl.hpp:367
Tempus::StepperDIRK::getDefaultStepperState
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
Definition: Tempus_StepperDIRK_impl.hpp:315
Tempus::StepperDIRK::initialize
virtual void initialize()
Initialize during construction and after changing input parameters.
Definition: Tempus_StepperDIRK_impl.hpp:123
Tempus::PASSED
Definition: Tempus_Types.hpp:17
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:25
Tempus::FAILED
Definition: Tempus_Types.hpp:18
Tempus::StepperDIRK::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_StepperDIRK_impl.hpp:357
Tempus::StepperDIRK::setParameterList
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl)
Definition: Tempus_StepperDIRK_impl.hpp:341
Tempus_StepperFactory.hpp
Tempus::StepperDIRK::setTableau
virtual void setTableau(std::string stepperType)
Definition: Tempus_StepperDIRK_impl.hpp:61
Tempus::StepperDIRK::setObserver
virtual void setObserver(Teuchos::RCP< StepperObserver< Scalar > > obs=Teuchos::null)
Set Observer.
Definition: Tempus_StepperDIRK_impl.hpp:102