Tempus  Version of the Day
Time Integration
Tempus_StepperIMEX_RK_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_StepperIMEX_RK_impl_hpp
10 #define Tempus_StepperIMEX_RK_impl_hpp
11 
12 #include "Tempus_RKButcherTableauBuilder.hpp"
13 #include "Tempus_config.hpp"
14 #include "Tempus_StepperFactory.hpp"
15 #include "Tempus_WrapperModelEvaluatorPairIMEX_Basic.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 
27 template<class Scalar>
28 StepperIMEX_RK<Scalar>::StepperIMEX_RK(
29  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
30  std::string stepperType)
31 {
32  this->setTableaus(Teuchos::null, stepperType);
33  this->setParameterList(Teuchos::null);
34  this->setModel(appModel);
35  this->initialize();
36 }
37 
38 
39 template<class Scalar>
40 StepperIMEX_RK<Scalar>::StepperIMEX_RK(
41  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
42  Teuchos::RCP<Teuchos::ParameterList> pList)
43 {
44  this->setTableaus(pList, "IMEX RK SSP2");
45  this->setParameterList(pList);
46  this->setModel(appModel);
47  this->initialize();
48 }
49 
50 
51 template<class Scalar>
52 StepperIMEX_RK<Scalar>::StepperIMEX_RK(
53  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
54  std::string stepperType,
55  Teuchos::RCP<Teuchos::ParameterList> pList)
56 {
57  this->setTableaus(pList, stepperType);
58  this->setParameterList(pList);
59  this->setModel(appModel);
60  this->initialize();
61 }
62 
63 
64 template<class Scalar>
65 void StepperIMEX_RK<Scalar>::setTableaus(
66  Teuchos::RCP<Teuchos::ParameterList> pList,
67  std::string stepperType)
68 {
69  if (stepperType == "") {
70  if (pList == Teuchos::null)
71  stepperType = "IMEX RK SSP2";
72  else
73  stepperType = pList->get<std::string>("Stepper Type", "IMEX RK SSP2");
74  }
75 
76  if (stepperType == "IMEX RK 1st order") {
77  {
78  // Explicit Tableau
79  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
80  pl->setName("IMEX-RK Explicit Stepper");
81  pl->set<std::string>("Stepper Type", "General ERK");
82 
83  // Tableau ParameterList
84  Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
85  tableauPL->set<std::string>("A", "0.0 0.0; 1.0 0.0");
86  tableauPL->set<std::string>("b", "1.0 0.0");
87  tableauPL->set<std::string>("c", "0.0 1.0");
88  tableauPL->set<int>("order", 1);
89  pl->set("Tableau", *tableauPL);
90 
91  this->setExplicitTableau("General ERK", pl);
92  }
93  {
94  // Implicit Tableau
95  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
96  pl->setName("IMEX-RK Implicit Stepper");
97  pl->set<std::string>("Stepper Type", "General DIRK");
98  pl->set("Solver Name", "");
99 
100  // Tableau ParameterList
101  Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
102  tableauPL->set<std::string>("A", "0.0 0.0; 0.0 1.0");
103  tableauPL->set<std::string>("b", "0.0 1.0");
104  tableauPL->set<std::string>("c", "0.0 1.0");
105  tableauPL->set<int>("order", 1);
106  pl->set("Tableau", *tableauPL);
107 
108  this->setImplicitTableau("General DIRK", pl);
109  }
110  description_ = stepperType;
111  order_ = 1;
112 
113  } else if (stepperType == "IMEX RK SSP2") {
114  typedef Teuchos::ScalarTraits<Scalar> ST;
115  // Explicit Tableau
116  this->setExplicitTableau("RK Explicit Trapezoidal", Teuchos::null);
117 
118  // Implicit Tableau
119  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
120  pl->set<std::string>("Stepper Type", "SDIRK 2 Stage 3rd order");
121  pl->set("Solver Name", "");
122  Scalar gamma = 1.0 - 1.0/ST::squareroot(2.0);
123  pl->set<double>("gamma",gamma);
124  this->setImplicitTableau("SDIRK 2 Stage 3rd order", pl);
125 
126  description_ = stepperType;
127  order_ = 2;
128  } else if (stepperType == "IMEX RK ARS 233") {
129  using std::to_string;
130  typedef Teuchos::ScalarTraits<Scalar> ST;
131  const Scalar gammaN = (3.0+ST::squareroot(3.0))/(6.0);
132  std::string gamma = to_string( gammaN);
133  std::string one_gamma = to_string(1.0- gammaN);
134  std::string one_2gamma = to_string(1.0-2.0*gammaN);
135  std::string two_2gamma = to_string(2.0-2.0*gammaN);
136  std::string gamma_one = to_string( gammaN-1.0);
137  {
138  // Explicit Tableau
139  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
140  pl->setName("IMEX-RK Explicit Stepper");
141  pl->set<std::string>("Stepper Type", "General ERK");
142 
143  // Tableau ParameterList
144  Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
145  tableauPL->set<std::string>("A",
146  "0.0 0.0 0.0; "+gamma+" 0.0 0.0; "+gamma_one+" "+two_2gamma+" 0.0");
147  tableauPL->set<std::string>("b", "0.0 0.5 0.5");
148  tableauPL->set<std::string>("c", "0.0 "+gamma+" "+one_gamma);
149  tableauPL->set<int>("order", 2);
150  pl->set("Tableau", *tableauPL);
151 
152  this->setExplicitTableau("General ERK", pl);
153  }
154  {
155  // Implicit Tableau
156  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
157  pl->setName("IMEX-RK Implicit Stepper");
158  pl->set<std::string>("Stepper Type", "General DIRK");
159  pl->set("Solver Name", "");
160 
161  // Tableau ParameterList
162  Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
163  tableauPL->set<std::string>("A",
164  "0.0 0.0 0.0; 0.0 "+gamma+" 0.0; 0.0 "+one_2gamma+" "+gamma);
165  tableauPL->set<std::string>("b", "0.0 0.5 0.5");
166  tableauPL->set<std::string>("c", "0.0 "+gamma+" "+one_gamma);
167  tableauPL->set<int>("order", 3);
168  pl->set("Tableau", *tableauPL);
169 
170  this->setImplicitTableau("General DIRK", pl);
171  }
172  description_ = stepperType;
173  order_ = 3;
174 
175  } else if (stepperType == "General IMEX RK") {
176  Teuchos::RCP<Teuchos::ParameterList> explicitPL = Teuchos::rcp(
177  new Teuchos::ParameterList(pList->sublist("IMEX-RK Explicit Stepper")));
178 
179  Teuchos::RCP<Teuchos::ParameterList> implicitPL = Teuchos::rcp(
180  new Teuchos::ParameterList(pList->sublist("IMEX-RK Implicit Stepper")));
181 
182  // TODO: should probably check the order of the tableau match
183  this->setExplicitTableau("General ERK", explicitPL);
184  this->setImplicitTableau("General DIRK", implicitPL);
185  description_ = stepperType;
186  order_ = pList->get<int>("overall order", 0);
187 
188  } else {
189  TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,
190  "Error - Not a valid StepperIMEX_RK type! Stepper Type = "
191  << stepperType << "\n"
192  << " Current valid types are: " << "\n"
193  << " 'IMEX RK 1st order'" << "\n"
194  << " 'IMEX RK SSP2'" << "\n"
195  << " 'IMEX RK ARS 233'" << "\n"
196  << " 'General IMEX RK'" << "\n");
197  }
198 
199  TEUCHOS_TEST_FOR_EXCEPTION(explicitTableau_==Teuchos::null,
200  std::runtime_error,
201  "Error - StepperIMEX_RK - Explicit tableau is null!");
202  TEUCHOS_TEST_FOR_EXCEPTION(implicitTableau_==Teuchos::null,
203  std::runtime_error,
204  "Error - StepperIMEX_RK - Implicit tableau is null!");
205  TEUCHOS_TEST_FOR_EXCEPTION(
206  explicitTableau_->numStages()!=implicitTableau_->numStages(),
207  std::runtime_error,
208  "Error - StepperIMEX_RK - Number of stages do not match!\n"
209  << " Explicit tableau = " << explicitTableau_->description() << "\n"
210  << " number of stages = " << explicitTableau_->numStages() << "\n"
211  << " Implicit tableau = " << implicitTableau_->description() << "\n"
212  << " number of stages = " << implicitTableau_->numStages() << "\n");
213 }
214 
215 
216 template<class Scalar>
217 void StepperIMEX_RK<Scalar>::setExplicitTableau(
218  std::string stepperType,
219  Teuchos::RCP<Teuchos::ParameterList> pList)
220 {
221  explicitTableau_ = createRKBT<Scalar>(stepperType,pList);
222  TEUCHOS_TEST_FOR_EXCEPTION(explicitTableau_->isImplicit() == true,
223  std::logic_error,
224  "Error - Received an implicit Tableau for setExplicitTableau()!\n"
225  << " Stepper Type = " << stepperType << "\n");
226 }
227 
228 
229 template<class Scalar>
230 void StepperIMEX_RK<Scalar>::setExplicitTableau(
231  Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau)
232 {
233  TEUCHOS_TEST_FOR_EXCEPTION(explicitTableau->isImplicit() == true,
234  std::logic_error,
235  "Error - Received an implicit Tableau for setExplicitTableau()!\n"
236  << " explicitTableau = " << explicitTableau->description() << "\n");
237  explicitTableau_ = explicitTableau;
238 }
239 
240 
241 template<class Scalar>
242 void StepperIMEX_RK<Scalar>::setImplicitTableau(
243  std::string stepperType,
244  Teuchos::RCP<Teuchos::ParameterList> pList)
245 {
246  implicitTableau_ = createRKBT<Scalar>(stepperType,pList);
247  //Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
248  //Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
249  //implicitTableau_->describe(*out,verbLevel);
250  TEUCHOS_TEST_FOR_EXCEPTION( implicitTableau_->isDIRK() != true,
251  std::logic_error,
252  "Error - Did not receive a DIRK Tableau for setImplicitTableau()!\n"
253  << " Stepper Type = " << stepperType << "\n");
254 }
255 
256 
257 template<class Scalar>
258 void StepperIMEX_RK<Scalar>::setImplicitTableau(
259  Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau)
260 {
261  TEUCHOS_TEST_FOR_EXCEPTION( implicitTableau_->isDIRK() != true,
262  std::logic_error,
263  "Error - Did not receive a DIRK Tableau for setImplicitTableau()!\n"
264  << " implicitTableau = " << implicitTableau->description() << "\n");
265  implicitTableau_ = implicitTableau;
266 }
267 
268 template<class Scalar>
269 void StepperIMEX_RK<Scalar>::setModel(
270  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
271 {
272  using Teuchos::RCP;
273  using Teuchos::rcp_const_cast;
274  using Teuchos::rcp_dynamic_cast;
275  RCP<Thyra::ModelEvaluator<Scalar> > ncModel =
276  rcp_const_cast<Thyra::ModelEvaluator<Scalar> > (appModel);
277  RCP<WrapperModelEvaluatorPairIMEX_Basic<Scalar> > modelPairIMEX =
278  rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX_Basic<Scalar> > (ncModel);
279  TEUCHOS_TEST_FOR_EXCEPTION( modelPairIMEX == Teuchos::null, std::logic_error,
280  "Error - StepperIMEX_RK::setModel() was given a ModelEvaluator that\n"
281  " could not be cast to a WrapperModelEvaluatorPairIMEX_Basic!\n"
282  " From: " << appModel << "\n"
283  " To : " << modelPairIMEX << "\n"
284  " Likely have given the wrong ModelEvaluator to this Stepper.\n");
285 
286  setModelPair(modelPairIMEX);
287 }
288 
289 
290 /** \brief Create WrapperModelPairIMEX from user-supplied ModelEvaluator pair.
291  *
292  * The user-supplied ME pair can contain any user-specific IMEX interactions
293  * between explicit and implicit MEs.
294  */
295 template<class Scalar>
296 void StepperIMEX_RK<Scalar>::setModelPair(
297  const Teuchos::RCP<WrapperModelEvaluatorPairIMEX_Basic<Scalar> > &
298  modelPairIMEX)
299 {
300  Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
301  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
302  this->wrapperModel_);
303  this->validExplicitODE (modelPairIMEX->getExplicitModel());
304  this->validImplicitODE_DAE(modelPairIMEX->getImplicitModel());
305  wrapperModelPairIMEX = modelPairIMEX;
306  wrapperModelPairIMEX->initialize();
307  int expXDim = wrapperModelPairIMEX->getExplicitModel()->get_x_space()->dim();
308  int impXDim = wrapperModelPairIMEX->getImplicitModel()->get_x_space()->dim();
309  TEUCHOS_TEST_FOR_EXCEPTION( expXDim != impXDim, std::logic_error,
310  "Error - \n"
311  " Explicit and Implicit x vectors are incompatible!\n"
312  " Explicit vector dim = " << expXDim << "\n"
313  " Implicit vector dim = " << impXDim << "\n");
314 
315  this->wrapperModel_ = wrapperModelPairIMEX;
316 }
317 
318 /** \brief Create WrapperModelPairIMEX from explicit/implicit ModelEvaluators.
319  *
320  * Use the supplied explicit/implicit MEs to create a WrapperModelPairIMEX
321  * with basic IMEX interactions between explicit and implicit MEs.
322  */
323 template<class Scalar>
324 void StepperIMEX_RK<Scalar>::setModelPair(
325  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& explicitModel,
326  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& implicitModel)
327 {
328  this->validExplicitODE (explicitModel);
329  this->validImplicitODE_DAE(implicitModel);
330  this->wrapperModel_ = Teuchos::rcp(
331  new WrapperModelEvaluatorPairIMEX_Basic<Scalar>(
332  explicitModel, implicitModel));
333 }
334 
335 
336 template<class Scalar>
337 void StepperIMEX_RK<Scalar>::setObserver(
338  Teuchos::RCP<StepperObserver<Scalar> > obs)
339 {
340  if (obs == Teuchos::null) {
341  // Create default observer, otherwise keep current observer.
342  if (stepperObserver_ == Teuchos::null) {
343  stepperIMEX_RKObserver_ =
344  Teuchos::rcp(new StepperIMEX_RKObserver<Scalar>());
345  stepperObserver_ =
346  Teuchos::rcp_dynamic_cast<StepperObserver<Scalar> >
347  (stepperIMEX_RKObserver_);
348  }
349  } else {
350  stepperObserver_ = obs;
351  stepperIMEX_RKObserver_ =
352  Teuchos::rcp_dynamic_cast<StepperIMEX_RKObserver<Scalar> >
353  (stepperObserver_);
354  }
355 }
356 
357 
358 template<class Scalar>
359 void StepperIMEX_RK<Scalar>::initialize()
360 {
361  TEUCHOS_TEST_FOR_EXCEPTION(
362  (explicitTableau_ == Teuchos::null) || (implicitTableau_ == Teuchos::null),
363  std::logic_error,
364  "Error - Need to set the Butcher Tableaus, setTableaus(), before calling "
365  "StepperIMEX_RK::initialize()\n");
366 
367  TEUCHOS_TEST_FOR_EXCEPTION( this->wrapperModel_ == Teuchos::null,
368  std::logic_error,
369  "Error - Need to set the model, setModel(), before calling "
370  "StepperIMEX_RK::initialize()\n");
371 
372  this->setTableaus(this->stepperPL_);
373  this->setParameterList(this->stepperPL_);
374  this->setSolver();
375  this->setObserver();
376 
377  // Initialize the stage vectors
378  const int numStages = explicitTableau_->numStages();
379  stageF_.resize(numStages);
380  stageG_.resize(numStages);
381  for(int i=0; i < numStages; i++) {
382  stageF_[i] = Thyra::createMember(this->wrapperModel_->get_f_space());
383  stageG_[i] = Thyra::createMember(this->wrapperModel_->get_f_space());
384  assign(stageF_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
385  assign(stageG_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
386  }
387 
388  xTilde_ = Thyra::createMember(this->wrapperModel_->get_x_space());
389  assign(xTilde_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
390 }
391 
392 
393 template <typename Scalar>
394 void StepperIMEX_RK<Scalar>::evalImplicitModelExplicitly(
395  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & X,
396  Scalar time, Scalar stepSize, Scalar stageNumber,
397  const Teuchos::RCP<Thyra::VectorBase<Scalar> > & G) const
398 {
399  typedef Thyra::ModelEvaluatorBase MEB;
400  Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
401  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
402  this->wrapperModel_);
403  MEB::InArgs<Scalar> inArgs = wrapperModelPairIMEX->getInArgs();
404  inArgs.set_x(X);
405  if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
406  if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(stepSize);
407  if (inArgs.supports(MEB::IN_ARG_stage_number))
408  inArgs.set_stage_number(stageNumber);
409 
410  // For model evaluators whose state function f(x, x_dot, t) describes
411  // an implicit ODE, and which accept an optional x_dot input argument,
412  // make sure the latter is set to null in order to request the evaluation
413  // of a state function corresponding to the explicit ODE formulation
414  // x_dot = f(x, t)
415  if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
416 
417  MEB::OutArgs<Scalar> outArgs = wrapperModelPairIMEX->getOutArgs();
418  outArgs.set_f(G);
419 
420  wrapperModelPairIMEX->getImplicitModel()->evalModel(inArgs,outArgs);
421  Thyra::Vt_S(G.ptr(), -1.0);
422 }
423 
424 
425 template <typename Scalar>
426 void StepperIMEX_RK<Scalar>::evalExplicitModel(
427  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & X,
428  Scalar time, Scalar stepSize, Scalar stageNumber,
429  const Teuchos::RCP<Thyra::VectorBase<Scalar> > & F) const
430 {
431  typedef Thyra::ModelEvaluatorBase MEB;
432 
433  Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
434  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
435  this->wrapperModel_);
436  MEB::InArgs<Scalar> inArgs =
437  wrapperModelPairIMEX->getExplicitModel()->createInArgs();
438  inArgs.set_x(X);
439  if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
440  if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(stepSize);
441  if (inArgs.supports(MEB::IN_ARG_stage_number))
442  inArgs.set_stage_number(stageNumber);
443 
444  // For model evaluators whose state function f(x, x_dot, t) describes
445  // an implicit ODE, and which accept an optional x_dot input argument,
446  // make sure the latter is set to null in order to request the evaluation
447  // of a state function corresponding to the explicit ODE formulation
448  // x_dot = f(x, t)
449  if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
450 
451  MEB::OutArgs<Scalar> outArgs =
452  wrapperModelPairIMEX->getExplicitModel()->createOutArgs();
453  outArgs.set_f(F);
454 
455  wrapperModelPairIMEX->getExplicitModel()->evalModel(inArgs, outArgs);
456  Thyra::Vt_S(F.ptr(), -1.0);
457 }
458 
459 
460 template<class Scalar>
461 void StepperIMEX_RK<Scalar>::takeStep(
462  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
463 {
464  using Teuchos::RCP;
465  using Teuchos::SerialDenseMatrix;
466  using Teuchos::SerialDenseVector;
467 
468  TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperIMEX_RK::takeStep()");
469  {
470  TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
471  std::logic_error,
472  "Error - StepperIMEX_RK<Scalar>::takeStep(...)\n"
473  "Need at least two SolutionStates for IMEX_RK.\n"
474  " Number of States = " << solutionHistory->getNumStates() << "\n"
475  "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
476  " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
477 
478  stepperObserver_->observeBeginTakeStep(solutionHistory, *this);
479  RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
480  RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
481  const Scalar dt = workingState->getTimeStep();
482  const Scalar time = currentState->getTime();
483 
484  const int numStages = explicitTableau_->numStages();
485  const SerialDenseMatrix<int,Scalar> & AHat = explicitTableau_->A();
486  const SerialDenseVector<int,Scalar> & bHat = explicitTableau_->b();
487  const SerialDenseVector<int,Scalar> & cHat = explicitTableau_->c();
488  const SerialDenseMatrix<int,Scalar> & A = implicitTableau_->A();
489  const SerialDenseVector<int,Scalar> & b = implicitTableau_->b();
490  const SerialDenseVector<int,Scalar> & c = implicitTableau_->c();
491 
492  bool pass = true;
493  Thyra::SolveStatus<Scalar> sStatus;
494  stageX_ = workingState->getX();
495  Thyra::assign(stageX_.ptr(), *(currentState->getX()));
496 
497  // Compute stage solutions
498  for (int i = 0; i < numStages; ++i) {
499  if (!Teuchos::is_null(stepperIMEX_RKObserver_))
500  stepperIMEX_RKObserver_->observeBeginStage(solutionHistory, *this);
501  Thyra::assign(xTilde_.ptr(), *(currentState->getX()));
502  for (int j = 0; j < i; ++j) {
503  if (AHat(i,j) != Teuchos::ScalarTraits<Scalar>::zero())
504  Thyra::Vp_StV(xTilde_.ptr(), -dt*AHat(i,j), *(stageF_[j]));
505  if (A (i,j) != Teuchos::ScalarTraits<Scalar>::zero())
506  Thyra::Vp_StV(xTilde_.ptr(), -dt*A (i,j), *(stageG_[j]));
507  }
508 
509  Scalar ts = time + c(i)*dt;
510  Scalar tHats = time + cHat(i)*dt;
511  if (A(i,i) == Teuchos::ScalarTraits<Scalar>::zero()) {
512  // Explicit stage for the ImplicitODE_DAE
513  bool isNeeded = false;
514  for (int k=i+1; k<numStages; ++k) if (A(k,i) != 0.0) isNeeded = true;
515  if (b(i) != 0.0) isNeeded = true;
516  if (isNeeded == false) {
517  // stageG_[i] is not needed.
518  assign(stageG_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
519  } else {
520  Thyra::assign(stageX_.ptr(), *xTilde_);
521  if (!Teuchos::is_null(stepperIMEX_RKObserver_))
522  stepperIMEX_RKObserver_->
523  observeBeforeImplicitExplicitly(solutionHistory, *this);
524  evalImplicitModelExplicitly(stageX_, ts, dt, i, stageG_[i]);
525  }
526  } else {
527  // Implicit stage for the ImplicitODE_DAE
528  Scalar alpha = 1.0/(dt*A(i,i));
529 
530  // Setup TimeDerivative
531  Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
532  Teuchos::rcp(new StepperIMEX_RKTimeDerivative<Scalar>(
533  alpha, xTilde_.getConst()));
534 
535  // Setup InArgs and OutArgs
536  typedef Thyra::ModelEvaluatorBase MEB;
537  MEB::InArgs<Scalar> inArgs = this->wrapperModel_->getInArgs();
538  MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
539  inArgs.set_x(stageX_);
540  if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(stageG_[i]);
541  if (inArgs.supports(MEB::IN_ARG_t )) inArgs.set_t (ts);
542  if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(dt);
543  if (inArgs.supports(MEB::IN_ARG_alpha )) inArgs.set_alpha (alpha);
544  if (inArgs.supports(MEB::IN_ARG_beta )) inArgs.set_beta (1.0);
545  if (inArgs.supports(MEB::IN_ARG_stage_number))
546  inArgs.set_stage_number(i);
547 
548  this->wrapperModel_->setForSolve(timeDer, inArgs, outArgs);
549 
550  if (!Teuchos::is_null(stepperIMEX_RKObserver_))
551  stepperIMEX_RKObserver_->observeBeforeSolve(solutionHistory, *this);
552 
553  sStatus = this->solveImplicitODE(stageX_);
554 
555  if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass = false;
556 
557  if (!Teuchos::is_null(stepperIMEX_RKObserver_))
558  stepperIMEX_RKObserver_->observeAfterSolve(solutionHistory, *this);
559 
560  // Update contributions to stage values
561  Thyra::V_StVpStV(stageG_[i].ptr(), -alpha, *stageX_, alpha, *xTilde_);
562  }
563 
564  if (!Teuchos::is_null(stepperIMEX_RKObserver_))
565  stepperIMEX_RKObserver_->observeBeforeExplicit(solutionHistory, *this);
566  evalExplicitModel(stageX_, tHats, dt, i, stageF_[i]);
567  if (!Teuchos::is_null(stepperIMEX_RKObserver_))
568  stepperIMEX_RKObserver_->observeEndStage(solutionHistory, *this);
569  }
570 
571  // Sum for solution: x_n = x_n-1 - dt*Sum{ bHat(i)*f(i) + b(i)*g(i) }
572  Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
573  for (int i=0; i < numStages; ++i) {
574  if (bHat(i) != Teuchos::ScalarTraits<Scalar>::zero())
575  Thyra::Vp_StV((workingState->getX()).ptr(), -dt*bHat(i), *(stageF_[i]));
576  if (b (i) != Teuchos::ScalarTraits<Scalar>::zero())
577  Thyra::Vp_StV((workingState->getX()).ptr(), -dt*b (i), *(stageG_[i]));
578  }
579 
580  if (pass == true) workingState->setSolutionStatus(Status::PASSED);
581  else workingState->setSolutionStatus(Status::FAILED);
582  workingState->setOrder(this->getOrder());
583  stepperObserver_->observeEndTakeStep(solutionHistory, *this);
584  }
585  return;
586 }
587 
588 /** \brief Provide a StepperState to the SolutionState.
589  * This Stepper does not have any special state data,
590  * so just provide the base class StepperState with the
591  * Stepper description. This can be checked to ensure
592  * that the input StepperState can be used by this Stepper.
593  */
594 template<class Scalar>
595 Teuchos::RCP<Tempus::StepperState<Scalar> >
596 StepperIMEX_RK<Scalar>::
597 getDefaultStepperState()
598 {
599  Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
600  rcp(new StepperState<Scalar>(description()));
601  return stepperState;
602 }
603 
604 
605 template<class Scalar>
606 std::string StepperIMEX_RK<Scalar>::description() const
607 {
608  return(description_);
609 }
610 
611 
612 template<class Scalar>
613 void StepperIMEX_RK<Scalar>::describe(
614  Teuchos::FancyOStream &out,
615  const Teuchos::EVerbosityLevel verbLevel) const
616 {
617  Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
618  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
619  this->wrapperModel_);
620  out << description() << "::describe:" << std::endl
621  << "wrapperModelPairIMEX = " << wrapperModelPairIMEX->description()
622  << std::endl;
623 }
624 
625 
626 template <class Scalar>
627 void StepperIMEX_RK<Scalar>::setParameterList(
628  const Teuchos::RCP<Teuchos::ParameterList> & pList)
629 {
630  if (pList == Teuchos::null) {
631  // Create default parameters if null, otherwise keep current parameters.
632  if (this->stepperPL_ == Teuchos::null) this->stepperPL_ = this->getDefaultParameters();
633  } else {
634  this->stepperPL_ = pList;
635  }
636  // Can not validate because of optional Parameters.
637  //stepperPL_->validateParametersAndSetDefaults(*this->getValidParameters());
638 }
639 
640 
641 template<class Scalar>
642 Teuchos::RCP<const Teuchos::ParameterList>
643 StepperIMEX_RK<Scalar>::getValidParameters() const
644 {
645  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
646  pl->setName("Default Stepper - IMEX RK SSP2");
647  pl->set("Stepper Type", "IMEX RK SSP2");
648  pl->set("Zero Initial Guess", false);
649  pl->set("Solver Name", "",
650  "Name of ParameterList containing the solver specifications.");
651 
652  return pl;
653 }
654 
655 template <class Scalar>
656 Teuchos::RCP<Teuchos::ParameterList>
657 StepperIMEX_RK<Scalar>::getDefaultParameters() const
658 {
659  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
660  pl->setName("Default Stepper - IMEX RK SSP2");
661  pl->set<std::string>("Stepper Type", "IMEX RK SSP2");
662  pl->set<bool> ("Zero Initial Guess", false);
663  pl->set<std::string>("Solver Name", "Default Solver");
664  Teuchos::RCP<Teuchos::ParameterList> solverPL=this->defaultSolverParameters();
665  pl->set("Default Solver", *solverPL);
666 
667  return pl;
668 }
669 
670 template <class Scalar>
671 Teuchos::RCP<Teuchos::ParameterList>
672 StepperIMEX_RK<Scalar>::getNonconstParameterList()
673 {
674  return(this->stepperPL_);
675 }
676 
677 
678 template <class Scalar>
679 Teuchos::RCP<Teuchos::ParameterList>
680 StepperIMEX_RK<Scalar>::unsetParameterList()
681 {
682  Teuchos::RCP<Teuchos::ParameterList> temp_plist = this->stepperPL_;
683  this->stepperPL_ = Teuchos::null;
684  return(temp_plist);
685 }
686 
687 
688 } // namespace Tempus
689 #endif // Tempus_StepperIMEX_RK_impl_hpp
Tempus::StepperState
StepperState is a simple class to hold state information about the stepper.
Definition: Tempus_StepperState.hpp:36
Tempus::RKButcherTableau
Runge-Kutta methods.
Definition: Tempus_RKButcherTableau.hpp:53
Tempus::StepperIMEX_RK::setModel
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Definition: Tempus_StepperIMEX_RK_impl.hpp:269
Tempus::WrapperModelEvaluatorPairIMEX_Basic
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
Definition: Tempus_WrapperModelEvaluatorPairIMEX_Basic_decl.hpp:38
Tempus::StepperIMEX_RK::setObserver
virtual void setObserver(Teuchos::RCP< StepperObserver< Scalar > > obs=Teuchos::null)
Set Observer.
Definition: Tempus_StepperIMEX_RK_impl.hpp:337
Tempus::solutionHistory
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
Definition: Tempus_SolutionHistory_impl.hpp:504
Tempus::StepperIMEX_RK::setParameterList
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl)
Definition: Tempus_StepperIMEX_RK_impl.hpp:627
Tempus::WrapperModelEvaluatorPairIMEX
ModelEvaluator pair for implicit and explicit (IMEX) evaluations.
Definition: Tempus_WrapperModelEvaluatorPairIMEX.hpp:25
Tempus::StepperIMEX_RK::getDefaultParameters
Teuchos::RCP< Teuchos::ParameterList > getDefaultParameters() const
Definition: Tempus_StepperIMEX_RK_impl.hpp:657
Tempus::StepperIMEX_RK::setExplicitTableau
virtual void setExplicitTableau(std::string stepperType, Teuchos::RCP< Teuchos::ParameterList > pList)
Set the explicit tableau from ParameterList.
Definition: Tempus_StepperIMEX_RK_impl.hpp:217
Tempus::WrapperModelEvaluator::getInArgs
virtual Thyra::ModelEvaluatorBase::InArgs< Scalar > getInArgs()=0
Get InArgs the wrapper ModelEvalutor.
Tempus
Definition: Tempus_AdjointAuxSensitivityModelEvaluator_decl.hpp:20
Tempus::StepperIMEX_RK::takeStep
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
Definition: Tempus_StepperIMEX_RK_impl.hpp:461
Tempus::StepperIMEX_RK::initialize
virtual void initialize()
Initialize during construction and after changing input parameters.
Definition: Tempus_StepperIMEX_RK_impl.hpp:359
Tempus::StepperIMEX_RK::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_StepperIMEX_RK_impl.hpp:643
Tempus::StepperIMEX_RK::setImplicitTableau
virtual void setImplicitTableau(std::string stepperType, Teuchos::RCP< Teuchos::ParameterList > pList)
Set the implicit tableau from ParameterList.
Definition: Tempus_StepperIMEX_RK_impl.hpp:242
Tempus::StepperIMEX_RK::getNonconstParameterList
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Definition: Tempus_StepperIMEX_RK_impl.hpp:672
Tempus::StepperObserver
StepperObserver class for Stepper class.
Definition: Tempus_StepperObserver.hpp:38
Tempus::StepperIMEX_RK::getDefaultStepperState
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Provide a StepperState to the SolutionState. This Stepper does not have any special state data,...
Definition: Tempus_StepperIMEX_RK_impl.hpp:597
Tempus::StepperIMEX_RK::setModelPair
virtual void setModelPair(const Teuchos::RCP< WrapperModelEvaluatorPairIMEX_Basic< Scalar > > &mePair)
Create WrapperModelPairIMEX from user-supplied ModelEvaluator pair.
Definition: Tempus_StepperIMEX_RK_impl.hpp:296
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::WrapperModelEvaluatorPairIMEX::getExplicitModel
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getExplicitModel() const =0
Tempus::StepperIMEX_RKTimeDerivative
Time-derivative interface for IMEX RK.
Definition: Tempus_StepperIMEX_RK_decl.hpp:377
Tempus::StepperIMEX_RK::setTableaus
virtual void setTableaus(Teuchos::RCP< Teuchos::ParameterList > pList, std::string stepperType="")
Set both the explicit and implicit tableau from ParameterList.
Definition: Tempus_StepperIMEX_RK_impl.hpp:65
Tempus::FAILED
Definition: Tempus_Types.hpp:18
Tempus::StepperIMEX_RK::evalImplicitModelExplicitly
void evalImplicitModelExplicitly(const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &X, Scalar time, Scalar stepSize, Scalar stageNumber, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &G) const
Definition: Tempus_StepperIMEX_RK_impl.hpp:394
Tempus::StepperIMEX_RK::unsetParameterList
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Definition: Tempus_StepperIMEX_RK_impl.hpp:680
Tempus::StepperIMEX_RK::description
virtual std::string description() const
Definition: Tempus_StepperIMEX_RK_impl.hpp:606
Tempus_StepperFactory.hpp
Tempus::WrapperModelEvaluatorPairIMEX::initialize
virtual void initialize()=0
Initialize after setting member data.
Tempus::StepperIMEX_RK::StepperIMEX_RK
StepperIMEX_RK()
Default Constructor – not allowed.
Tempus::StepperIMEX_RKObserver
StepperIMEX_RKObserver class for StepperIMEX_RK.
Definition: Tempus_StepperIMEX_RKObserver.hpp:35
Tempus::StepperIMEX_RK::evalExplicitModel
void evalExplicitModel(const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &X, Scalar time, Scalar stepSize, Scalar stageNumber, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &F) const
Definition: Tempus_StepperIMEX_RK_impl.hpp:426
Tempus::StepperIMEX_RK::describe
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Definition: Tempus_StepperIMEX_RK_impl.hpp:613