Tempus  Version of the Day
Time Integration
Tempus_StepperIMEX_RK_Partition_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_Partition_impl_hpp
10 #define Tempus_StepperIMEX_RK_Partition_impl_hpp
11 
12 #include "Tempus_RKButcherTableauBuilder.hpp"
13 #include "Tempus_config.hpp"
14 #include "Tempus_StepperFactory.hpp"
15 #include "Tempus_WrapperModelEvaluatorPairPartIMEX_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_Partition<Scalar>::StepperIMEX_RK_Partition(
29  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
30  std::string stepperType)
31 {
32  this->setTableaus(Teuchos::null, stepperType);
33  this->setModel(appModel);
34  this->initialize();
35 }
36 
37 
38 template<class Scalar>
39 StepperIMEX_RK_Partition<Scalar>::StepperIMEX_RK_Partition(
40  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
41  Teuchos::RCP<Teuchos::ParameterList> pList)
42 {
43  this->setTableaus(pList, "Partitioned IMEX RK SSP2");
44  this->setParameterList(pList);
45  this->setModel(appModel);
46  this->initialize();
47 }
48 
49 
50 template<class Scalar>
51 StepperIMEX_RK_Partition<Scalar>::StepperIMEX_RK_Partition(
52  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
53  std::string stepperType,
54  Teuchos::RCP<Teuchos::ParameterList> pList)
55 {
56  this->setTableaus(pList, stepperType);
57  this->setParameterList(pList);
58  this->setModel(appModel);
59  this->initialize();
60 }
61 
62 
63 template<class Scalar>
64 void StepperIMEX_RK_Partition<Scalar>::setTableaus(
65  Teuchos::RCP<Teuchos::ParameterList> pList,
66  std::string stepperType)
67 {
68  if (stepperType == "") {
69  if (pList == Teuchos::null)
70  stepperType = "Partitioned IMEX RK SSP2";
71  else
72  stepperType = pList->get<std::string>("Stepper Type",
73  "Partitioned IMEX RK SSP2");
74  }
75 
76  if (stepperType == "Partitioned 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 == "Partitioned 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 == "Partitioned 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 Partitioned 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_Partition type! Stepper Type = "
191  << stepperType << "\n"
192  << " Current valid types are: " << "\n"
193  << " 'Partitioned IMEX RK 1st order'" << "\n"
194  << " 'Partitioned IMEX RK SSP2'" << "\n"
195  << " 'Partitioned IMEX RK ARS 233'" << "\n"
196  << " 'General Partitioned IMEX RK'" << "\n");
197  }
198 
199  TEUCHOS_TEST_FOR_EXCEPTION(explicitTableau_==Teuchos::null,
200  std::runtime_error,
201  "Error - StepperIMEX_RK_Partition - Explicit tableau is null!");
202  TEUCHOS_TEST_FOR_EXCEPTION(implicitTableau_==Teuchos::null,
203  std::runtime_error,
204  "Error - StepperIMEX_RK_Partition - Implicit tableau is null!");
205  TEUCHOS_TEST_FOR_EXCEPTION(
206  explicitTableau_->numStages()!=implicitTableau_->numStages(),
207  std::runtime_error,
208  "Error - StepperIMEX_RK_Partition - 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_Partition<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_Partition<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_Partition<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_Partition<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_Partition<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<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> > modelPairIMEX =
278  rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_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 WrapperModelEvaluatorPairPartIMEX_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 /** \brief Create WrapperModelPairIMEX from user-supplied ModelEvaluator pair.
290  *
291  * The user-supplied ME pair can contain any user-specific IMEX interactions
292  * between explicit and implicit MEs.
293  */
294 template<class Scalar>
295 void StepperIMEX_RK_Partition<Scalar>::setModelPair(
296  const Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> > &
297  mePairIMEX)
298 {
299  Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
300  wrapperModelPairIMEX =
301  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
302  (this->wrapperModel_);
303  this->validExplicitODE (mePairIMEX->getExplicitModel());
304  this->validImplicitODE_DAE(mePairIMEX->getImplicitModel());
305  wrapperModelPairIMEX = mePairIMEX;
306  wrapperModelPairIMEX->initialize();
307 
308  this->wrapperModel_ = wrapperModelPairIMEX;
309 }
310 
311 /** \brief Create WrapperModelPairIMEX from explicit/implicit ModelEvaluators.
312  *
313  * Use the supplied explicit/implicit MEs to create a WrapperModelPairIMEX
314  * with basic IMEX interactions between explicit and implicit MEs.
315  */
316 template<class Scalar>
317 void StepperIMEX_RK_Partition<Scalar>::setModelPair(
318  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& explicitModel,
319  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& implicitModel)
320 {
321  this->validExplicitODE (explicitModel);
322  this->validImplicitODE_DAE(implicitModel);
323  this->wrapperModel_ = Teuchos::rcp(
324  new WrapperModelEvaluatorPairPartIMEX_Basic<Scalar>(
325  explicitModel, implicitModel));
326 }
327 
328 
329 template<class Scalar>
330 void StepperIMEX_RK_Partition<Scalar>::setObserver(
331  Teuchos::RCP<StepperObserver<Scalar> > obs)
332 {
333  if (obs == Teuchos::null) {
334  // Create default observer, otherwise keep current observer.
335  if (stepperObserver_ == Teuchos::null) {
336  stepperIMEX_RKPartObserver_ =
337  Teuchos::rcp(new StepperIMEX_RKPartObserver<Scalar>());
338  stepperObserver_ =
339  Teuchos::rcp_dynamic_cast<StepperObserver<Scalar> >
340  (stepperIMEX_RKPartObserver_);
341  }
342  } else {
343  stepperObserver_ = obs;
344  stepperIMEX_RKPartObserver_ =
345  Teuchos::rcp_dynamic_cast<StepperIMEX_RKPartObserver<Scalar> >
346  (stepperObserver_);
347  }
348 }
349 
350 
351 template<class Scalar>
352 void StepperIMEX_RK_Partition<Scalar>::initialize()
353 {
354  TEUCHOS_TEST_FOR_EXCEPTION(
355  (explicitTableau_ == Teuchos::null) || (implicitTableau_ == Teuchos::null),
356  std::logic_error,
357  "Error - Need to set the Butcher Tableaus, setTableaus(), before calling "
358  "StepperIMEX_RK_Partition::initialize()\n");
359 
360  Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
361  wrapperModelPairIMEX =
362  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
363  (this->wrapperModel_);
364  TEUCHOS_TEST_FOR_EXCEPTION( wrapperModelPairIMEX == Teuchos::null,
365  std::logic_error,
366  "Error - Need to set the model, setModel(), before calling "
367  "StepperIMEX_RK_Partition::initialize()\n");
368 
369  this->setTableaus(this->stepperPL_);
370  this->setParameterList(this->stepperPL_);
371  this->setSolver();
372  this->setObserver();
373 
374  // Initialize the stage vectors
375  const int numStages = explicitTableau_->numStages();
376  stageF_.resize(numStages);
377  stageGx_.resize(numStages);
378  for(int i=0; i < numStages; i++) {
379  stageF_[i] = Thyra::createMember(wrapperModelPairIMEX->
380  getExplicitModel()->get_f_space());
381  stageGx_[i] = Thyra::createMember(wrapperModelPairIMEX->
382  getImplicitModel()->get_f_space());
383  assign(stageF_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
384  assign(stageGx_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
385  }
386 
387  xTilde_ = Thyra::createMember(wrapperModelPairIMEX->
388  getImplicitModel()->get_x_space());
389  assign(xTilde_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
390 }
391 
392 
393 template <typename Scalar>
394 void StepperIMEX_RK_Partition<Scalar>::evalImplicitModelExplicitly(
395  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & X,
396  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & Y,
397  Scalar time, Scalar stepSize, Scalar stageNumber,
398  const Teuchos::RCP<Thyra::VectorBase<Scalar> > & G) const
399 {
400  typedef Thyra::ModelEvaluatorBase MEB;
401  Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
402  wrapperModelPairIMEX =
403  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
404  (this->wrapperModel_);
405  MEB::InArgs<Scalar> inArgs = wrapperModelPairIMEX->getInArgs();
406  inArgs.set_x(X);
407  inArgs.set_p(wrapperModelPairIMEX->getParameterIndex(), Y);
408  if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
409  if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(stepSize);
410  if (inArgs.supports(MEB::IN_ARG_stage_number))
411  inArgs.set_stage_number(stageNumber);
412 
413  // For model evaluators whose state function f(x, x_dot, t) describes
414  // an implicit ODE, and which accept an optional x_dot input argument,
415  // make sure the latter is set to null in order to request the evaluation
416  // of a state function corresponding to the explicit ODE formulation
417  // x_dot = f(x, t)
418  if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
419 
420  MEB::OutArgs<Scalar> outArgs = wrapperModelPairIMEX->getOutArgs();
421  outArgs.set_f(G);
422 
423  wrapperModelPairIMEX->getImplicitModel()->evalModel(inArgs,outArgs);
424  Thyra::Vt_S(G.ptr(), -1.0);
425 }
426 
427 
428 template <typename Scalar>
429 void StepperIMEX_RK_Partition<Scalar>::evalExplicitModel(
430  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & Z,
431  Scalar time, Scalar stepSize, Scalar stageNumber,
432  const Teuchos::RCP<Thyra::VectorBase<Scalar> > & F) const
433 {
434  typedef Thyra::ModelEvaluatorBase MEB;
435 
436  Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
437  wrapperModelPairIMEX =
438  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
439  (this->wrapperModel_);
440  MEB::InArgs<Scalar> inArgs =
441  wrapperModelPairIMEX->getExplicitModel()->createInArgs();
442  inArgs.set_x(Z);
443  if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
444  if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(stepSize);
445  if (inArgs.supports(MEB::IN_ARG_stage_number))
446  inArgs.set_stage_number(stageNumber);
447 
448  // For model evaluators whose state function f(x, x_dot, t) describes
449  // an implicit ODE, and which accept an optional x_dot input argument,
450  // make sure the latter is set to null in order to request the evaluation
451  // of a state function corresponding to the explicit ODE formulation
452  // x_dot = f(x, t)
453  if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
454 
455  MEB::OutArgs<Scalar> outArgs =
456  wrapperModelPairIMEX->getExplicitModel()->createOutArgs();
457  outArgs.set_f(F);
458 
459  wrapperModelPairIMEX->getExplicitModel()->evalModel(inArgs, outArgs);
460  Thyra::Vt_S(F.ptr(), -1.0);
461 }
462 
463 
464 template<class Scalar>
465 void StepperIMEX_RK_Partition<Scalar>::takeStep(
466  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
467 {
468  using Teuchos::RCP;
469  using Teuchos::SerialDenseMatrix;
470  using Teuchos::SerialDenseVector;
471 
472  TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperIMEX_RK_Partition::takeStep()");
473  {
474  TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
475  std::logic_error,
476  "Error - StepperIMEX_RK_Partition<Scalar>::takeStep(...)\n"
477  "Need at least two SolutionStates for IMEX_RK_Partition.\n"
478  " Number of States = " << solutionHistory->getNumStates() << "\n"
479  "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
480  " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
481 
482  stepperObserver_->observeBeginTakeStep(solutionHistory, *this);
483  RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
484  RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
485  const Scalar dt = workingState->getTimeStep();
486  const Scalar time = currentState->getTime();
487 
488  const int numStages = explicitTableau_->numStages();
489  const SerialDenseMatrix<int,Scalar> & AHat = explicitTableau_->A();
490  const SerialDenseVector<int,Scalar> & bHat = explicitTableau_->b();
491  const SerialDenseVector<int,Scalar> & cHat = explicitTableau_->c();
492  const SerialDenseMatrix<int,Scalar> & A = implicitTableau_->A();
493  const SerialDenseVector<int,Scalar> & b = implicitTableau_->b();
494  const SerialDenseVector<int,Scalar> & c = implicitTableau_->c();
495 
496  Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
497  wrapperModelPairIMEX =
498  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
499  (this->wrapperModel_);
500 
501  bool pass = true;
502  Thyra::SolveStatus<Scalar> sStatus;
503  stageZ_ = workingState->getX();
504  Thyra::assign(stageZ_.ptr(), *(currentState->getX()));
505  RCP<Thyra::VectorBase<Scalar> > stageY =
506  wrapperModelPairIMEX->getExplicitOnlyVector(stageZ_);
507  RCP<Thyra::VectorBase<Scalar> > stageX =
508  wrapperModelPairIMEX->getIMEXVector(stageZ_);
509 
510  // Compute stage solutions
511  for (int i = 0; i < numStages; ++i) {
512  if (!Teuchos::is_null(stepperIMEX_RKPartObserver_))
513  stepperIMEX_RKPartObserver_->observeBeginStage(solutionHistory, *this);
514 
515  Thyra::assign(stageY.ptr(),
516  *(wrapperModelPairIMEX->getExplicitOnlyVector(currentState->getX())));
517  Thyra::assign(xTilde_.ptr(),
518  *(wrapperModelPairIMEX->getIMEXVector(currentState->getX())));
519  for (int j = 0; j < i; ++j) {
520  if (AHat(i,j) != Teuchos::ScalarTraits<Scalar>::zero()) {
521  RCP<Thyra::VectorBase<Scalar> > stageFy =
522  wrapperModelPairIMEX->getExplicitOnlyVector(stageF_[j]);
523  RCP<Thyra::VectorBase<Scalar> > stageFx =
524  wrapperModelPairIMEX->getIMEXVector(stageF_[j]);
525  Thyra::Vp_StV(stageY.ptr(), -dt*AHat(i,j), *stageFy);
526  Thyra::Vp_StV(xTilde_.ptr(), -dt*AHat(i,j), *stageFx);
527  }
528  if (A (i,j) != Teuchos::ScalarTraits<Scalar>::zero())
529  Thyra::Vp_StV(xTilde_.ptr(), -dt*A (i,j), *(stageGx_[j]));
530  }
531 
532  Scalar ts = time + c(i)*dt;
533  Scalar tHats = time + cHat(i)*dt;
534  if (A(i,i) == Teuchos::ScalarTraits<Scalar>::zero()) {
535  // Explicit stage for the ImplicitODE_DAE
536  bool isNeeded = false;
537  for (int k=i+1; k<numStages; ++k) if (A(k,i) != 0.0) isNeeded = true;
538  if (b(i) != 0.0) isNeeded = true;
539  if (isNeeded == false) {
540  // stageGx_[i] is not needed.
541  assign(stageGx_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
542  } else {
543  Thyra::assign(stageX.ptr(), *xTilde_);
544  if (!Teuchos::is_null(stepperIMEX_RKPartObserver_))
545  stepperIMEX_RKPartObserver_->
546  observeBeforeImplicitExplicitly(solutionHistory, *this);
547  evalImplicitModelExplicitly(stageX, stageY, ts, dt, i, stageGx_[i]);
548  }
549  } else {
550  // Implicit stage for the ImplicitODE_DAE
551  Scalar alpha = 1.0/(dt*A(i,i));
552 
553  // Setup TimeDerivative
554  RCP<TimeDerivative<Scalar> > timeDer =
555  Teuchos::rcp(new StepperIMEX_RKPartTimeDerivative<Scalar>(
556  alpha, xTilde_.getConst()));
557 
558  // Setup InArgs and OutArgs
559  typedef Thyra::ModelEvaluatorBase MEB;
560  //MEB::InArgs<Scalar> inArgs = wrapperModelPairIMEX->getInArgs();
561  //MEB::OutArgs<Scalar> outArgs = wrapperModelPairIMEX->getOutArgs();
562  wrapperModelPairIMEX->setUseImplicitModel(true);
563  MEB::InArgs<Scalar> inArgs = wrapperModelPairIMEX->createInArgs();
564  MEB::OutArgs<Scalar> outArgs = wrapperModelPairIMEX->createOutArgs();
565  inArgs.set_x(stageX);
566  if (wrapperModelPairIMEX->getParameterIndex() >= 0)
567  inArgs.set_p(wrapperModelPairIMEX->getParameterIndex(), stageY);
568  if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(stageGx_[i]);
569  if (inArgs.supports(MEB::IN_ARG_t )) inArgs.set_t (ts);
570  if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(dt);
571  if (inArgs.supports(MEB::IN_ARG_alpha )) inArgs.set_alpha (alpha);
572  if (inArgs.supports(MEB::IN_ARG_beta )) inArgs.set_beta (1.0);
573  if (inArgs.supports(MEB::IN_ARG_stage_number))
574  inArgs.set_stage_number(i);
575 
576  wrapperModelPairIMEX->setForSolve(timeDer, inArgs, outArgs);
577 
578  if (!Teuchos::is_null(stepperIMEX_RKPartObserver_))
579  stepperIMEX_RKPartObserver_->observeBeforeSolve(solutionHistory, *this);
580 
581  this->solver_->setModel(wrapperModelPairIMEX);
582  sStatus = this->solveImplicitODE(stageX);
583  if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass = false;
584 
585  wrapperModelPairIMEX->setUseImplicitModel(false);
586 
587  if (!Teuchos::is_null(stepperIMEX_RKPartObserver_))
588  stepperIMEX_RKPartObserver_->observeAfterSolve(solutionHistory, *this);
589 
590  // Update contributions to stage values
591  Thyra::V_StVpStV(stageGx_[i].ptr(), -alpha, *stageX, alpha, *xTilde_);
592  }
593 
594  if (!Teuchos::is_null(stepperIMEX_RKPartObserver_))
595  stepperIMEX_RKPartObserver_->observeBeforeExplicit(solutionHistory,*this);
596  evalExplicitModel(stageZ_, tHats, dt, i, stageF_[i]);
597  if (!Teuchos::is_null(stepperIMEX_RKPartObserver_))
598  stepperIMEX_RKPartObserver_->observeEndStage(solutionHistory, *this);
599  }
600 
601  // Sum for solution: y_n = y_n-1 - dt*Sum{ bHat(i)*fy(i) }
602  // Sum for solution: x_n = x_n-1 - dt*Sum{ bHat(i)*fx(i) + b(i)*gx(i) }
603  Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
604  RCP<Thyra::VectorBase<Scalar> > Z = workingState->getX();
605  RCP<Thyra::VectorBase<Scalar> > X = wrapperModelPairIMEX->getIMEXVector(Z);
606  for (int i=0; i < numStages; ++i) {
607  if (bHat(i) != Teuchos::ScalarTraits<Scalar>::zero())
608  Thyra::Vp_StV(Z.ptr(), -dt*bHat(i), *(stageF_[i]));
609  if (b (i) != Teuchos::ScalarTraits<Scalar>::zero())
610  Thyra::Vp_StV(X.ptr(), -dt*b (i), *(stageGx_[i]));
611  }
612 
613  if (pass == true) workingState->setSolutionStatus(Status::PASSED);
614  else workingState->setSolutionStatus(Status::FAILED);
615  workingState->setOrder(this->getOrder());
616  stepperObserver_->observeEndTakeStep(solutionHistory, *this);
617  }
618  return;
619 }
620 
621 /** \brief Provide a StepperState to the SolutionState.
622  * This Stepper does not have any special state data,
623  * so just provide the base class StepperState with the
624  * Stepper description. This can be checked to ensure
625  * that the input StepperState can be used by this Stepper.
626  */
627 template<class Scalar>
628 Teuchos::RCP<Tempus::StepperState<Scalar> >
629 StepperIMEX_RK_Partition<Scalar>::
630 getDefaultStepperState()
631 {
632  Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
633  rcp(new StepperState<Scalar>(description()));
634  return stepperState;
635 }
636 
637 
638 template<class Scalar>
639 std::string StepperIMEX_RK_Partition<Scalar>::description() const
640 {
641  return(description_);
642 }
643 
644 
645 template<class Scalar>
646 void StepperIMEX_RK_Partition<Scalar>::describe(
647  Teuchos::FancyOStream &out,
648  const Teuchos::EVerbosityLevel verbLevel) const
649 {
650  out << description() << "::describe:" << std::endl
651  << "wrapperModelPairIMEX = " << this->wrapperModel_->description()
652  << std::endl;
653 }
654 
655 
656 template <class Scalar>
657 void StepperIMEX_RK_Partition<Scalar>::setParameterList(
658  const Teuchos::RCP<Teuchos::ParameterList> & pList)
659 {
660  if (pList == Teuchos::null) {
661  // Create default parameters if null, otherwise keep current parameters.
662  if (this->stepperPL_ == Teuchos::null) this->stepperPL_ = this->getDefaultParameters();
663  } else {
664  this->stepperPL_ = pList;
665  }
666  // Can not validate because of optional Parameters.
667  //stepperPL_->validateParametersAndSetDefaults(*this->getValidParameters());
668 }
669 
670 
671 template<class Scalar>
672 Teuchos::RCP<const Teuchos::ParameterList>
673 StepperIMEX_RK_Partition<Scalar>::getValidParameters() const
674 {
675  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
676  pl->setName("Default Stepper - IMEX RK SSP2");
677  pl->set("Stepper Type", "IMEX RK SSP2");
678  pl->set("Zero Initial Guess", false);
679  pl->set("Solver Name", "",
680  "Name of ParameterList containing the solver specifications.");
681 
682  return pl;
683 }
684 
685 template <class Scalar>
686 Teuchos::RCP<Teuchos::ParameterList>
687 StepperIMEX_RK_Partition<Scalar>::getDefaultParameters() const
688 {
689  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
690  pl->setName("Default Stepper - IMEX RK SSP2");
691  pl->set<std::string>("Stepper Type", "IMEX RK SSP2");
692  pl->set<bool> ("Zero Initial Guess", false);
693  pl->set<std::string>("Solver Name", "Default Solver");
694  Teuchos::RCP<Teuchos::ParameterList> solverPL=this->defaultSolverParameters();
695  pl->set("Default Solver", *solverPL);
696 
697  return pl;
698 }
699 
700 template <class Scalar>
701 Teuchos::RCP<Teuchos::ParameterList>
702 StepperIMEX_RK_Partition<Scalar>::getNonconstParameterList()
703 {
704  return(this->stepperPL_);
705 }
706 
707 
708 template <class Scalar>
709 Teuchos::RCP<Teuchos::ParameterList>
710 StepperIMEX_RK_Partition<Scalar>::unsetParameterList()
711 {
712  Teuchos::RCP<Teuchos::ParameterList> temp_plist = this->stepperPL_;
713  this->stepperPL_ = Teuchos::null;
714  return(temp_plist);
715 }
716 
717 
718 } // namespace Tempus
719 #endif // Tempus_StepperIMEX_RK_Partition_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_Partition::description
virtual std::string description() const
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:639
Tempus::StepperIMEX_RK_Partition::setParameterList
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl)
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:657
Tempus::StepperIMEX_RK_Partition::describe
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:646
Tempus::StepperIMEX_RK_Partition::getDefaultParameters
Teuchos::RCP< Teuchos::ParameterList > getDefaultParameters() const
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:687
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_decl.hpp:37
Tempus::solutionHistory
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
Definition: Tempus_SolutionHistory_impl.hpp:504
Tempus::StepperIMEX_RKPartObserver
StepperIMEX_RKPartObserver class for StepperIMEX_RK_Partition.
Definition: Tempus_StepperIMEX_RKPartObserver.hpp:35
Tempus::StepperIMEX_RKPartTimeDerivative
Time-derivative interface for Partitioned IMEX RK.
Definition: Tempus_StepperIMEX_RK_Partition_decl.hpp:378
Tempus::StepperIMEX_RK_Partition::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_Partition_impl.hpp:64
Tempus::StepperIMEX_RK_Partition::setObserver
virtual void setObserver(Teuchos::RCP< StepperObserver< Scalar > > obs=Teuchos::null)
Set Observer.
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:330
Tempus
Definition: Tempus_AdjointAuxSensitivityModelEvaluator_decl.hpp:20
Tempus::StepperIMEX_RK_Partition::setModelPair
virtual void setModelPair(const Teuchos::RCP< WrapperModelEvaluatorPairPartIMEX_Basic< Scalar > > &modelPair)
Create WrapperModelPairIMEX from user-supplied ModelEvaluator pair.
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:295
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::initialize
virtual void initialize()
Initialize after setting member data.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_impl.hpp:61
Tempus::StepperIMEX_RK_Partition::evalImplicitModelExplicitly
void evalImplicitModelExplicitly(const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &X, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &Y, Scalar time, Scalar stepSize, Scalar stageNumber, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &G) const
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:394
Tempus::StepperIMEX_RK_Partition::getNonconstParameterList
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:702
Tempus::StepperIMEX_RK_Partition::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:673
Tempus::StepperIMEX_RK_Partition::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_Partition_impl.hpp:429
Tempus::StepperIMEX_RK_Partition::setModel
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:269
Tempus::StepperObserver
StepperObserver class for Stepper class.
Definition: Tempus_StepperObserver.hpp:38
Tempus::StepperIMEX_RK_Partition::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_Partition_impl.hpp:630
Tempus::StepperIMEX_RK_Partition::unsetParameterList
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:710
Tempus::StepperIMEX_RK_Partition::setImplicitTableau
virtual void setImplicitTableau(std::string stepperType, Teuchos::RCP< Teuchos::ParameterList > pList)
Set the implicit tableau from ParameterList.
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:242
Tempus::PASSED
Definition: Tempus_Types.hpp:17
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::getInArgs
virtual Thyra::ModelEvaluatorBase::InArgs< Scalar > getInArgs()
Get InArgs the wrapper ModelEvalutor.
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_decl.hpp:69
Tempus::StepperIMEX_RK_Partition::takeStep
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:465
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::StepperIMEX_RK_Partition::StepperIMEX_RK_Partition
StepperIMEX_RK_Partition()
Default Constructor – not allowed.
Tempus_StepperFactory.hpp
Tempus::StepperIMEX_RK_Partition::initialize
virtual void initialize()
Initialize during construction and after changing input parameters.
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:352
Tempus::WrapperModelEvaluatorPairPartIMEX_Basic::getExplicitModel
virtual Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getExplicitModel() const
Definition: Tempus_WrapperModelEvaluatorPairPartIMEX_Basic_decl.hpp:117
Tempus::StepperIMEX_RK_Partition::setExplicitTableau
virtual void setExplicitTableau(std::string stepperType, Teuchos::RCP< Teuchos::ParameterList > pList)
Set the explicit tableau from ParameterList.
Definition: Tempus_StepperIMEX_RK_Partition_impl.hpp:217