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