Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_StepperIMEX_RK_impl.hpp
Go to the documentation of this file.
1 //@HEADER
2 // *****************************************************************************
3 // Tempus: Time Integration and Sensitivity Analysis Package
4 //
5 // Copyright 2017 NTESS and the Tempus contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 //@HEADER
9 
10 #ifndef Tempus_StepperIMEX_RK_impl_hpp
11 #define Tempus_StepperIMEX_RK_impl_hpp
12 
13 #include "Thyra_VectorStdOps.hpp"
14 
15 #include "Tempus_StepperFactory.hpp"
17 
18 namespace Tempus {
19 
20 template <class Scalar>
22 {
24  stepperType != "IMEX RK 1st order" && stepperType != "SSP1_111" &&
25  stepperType != "IMEX RK SSP2" && stepperType != "IMEX RK SSP3" &&
26  stepperType != "SSP3_332" && stepperType != "SSP2_222" &&
27  stepperType != "SSP2_222_L" && stepperType != "SSP2_222_A" &&
28  stepperType != "IMEX RK ARS 233" && stepperType != "ARS 233" &&
29  stepperType != "General IMEX RK",
30  std::logic_error,
31  " 'Stepper Type' (='"
32  << stepperType
33  << "')\n"
34  << " does not match one of the types for this Stepper:\n"
35  << " 'IMEX RK 1st order'\n"
36  << " 'SSP1_111'\n"
37  << " 'IMEX RK SSP2'\n"
38  << " 'IMEX RK SSP3'\n"
39  << " 'SSP3_332'\n"
40  << " 'SSP2_222'\n"
41  << " 'SSP2_222_L'\n"
42  << " 'SSP2_222_A'\n"
43  << " 'IMEX RK ARS 233'\n"
44  << " 'ARS 233'\n"
45  << " 'General IMEX RK'\n");
46 
47  this->setStepperName(stepperType);
48  this->setStepperType(stepperType);
49  this->setUseFSAL(false);
50  this->setICConsistency("None");
51  this->setICConsistencyCheck(false);
52  this->setZeroInitialGuess(false);
53 
54  this->setStageNumber(-1);
55 
56  this->setTableaus(stepperType);
57  this->setAppAction(Teuchos::null);
58  this->setDefaultSolver();
59 }
60 
61 template <class Scalar>
63  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
65  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
66  bool zeroInitialGuess,
67  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
68  std::string stepperType,
69  Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau,
70  Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau, Scalar order)
71 {
73  stepperType != "IMEX RK 1st order" && stepperType != "SSP1_111" &&
74  stepperType != "IMEX RK SSP2" && stepperType != "IMEX RK SSP3" &&
75  stepperType != "SSP3_332" && stepperType != "SSP2_222" &&
76  stepperType != "SSP2_222_L" && stepperType != "SSP2_222_A" &&
77  stepperType != "IMEX RK ARS 233" && stepperType != "ARS 233" &&
78  stepperType != "General IMEX RK",
79  std::logic_error,
80  " 'Stepper Type' (='"
81  << stepperType
82  << "')\n"
83  << " does not match one of the types for this Stepper:\n"
84  << " 'IMEX RK 1st order'\n"
85  << " 'SSP1_111'\n"
86  << " 'IMEX RK SSP2'\n"
87  << " 'IMEX RK SSP3'\n"
88  << " 'SSP3_332'\n"
89  << " 'SSP2_222'\n"
90  << " 'SSP2_222_L'\n"
91  << " 'SSP2_222_A'\n"
92  << " 'IMEX RK ARS 233'\n"
93  << " 'ARS 233'\n"
94  << " 'General IMEX RK'\n");
95 
96  this->setStepperName(stepperType);
97  this->setStepperType(stepperType);
98  this->setUseFSAL(useFSAL);
99  this->setICConsistency(ICConsistency);
100  this->setICConsistencyCheck(ICConsistencyCheck);
101  this->setZeroInitialGuess(zeroInitialGuess);
102 
103  this->setStageNumber(-1);
104 
105  if (stepperType == "General IMEX RK") {
106  this->setExplicitTableau(explicitTableau);
107  this->setImplicitTableau(implicitTableau);
108  }
109  else {
110  this->setTableaus(stepperType);
111  }
112  this->setOrder(order);
113  this->setAppAction(stepperRKAppAction);
114  this->setSolver(solver);
115 
116  if (appModel != Teuchos::null) {
117  this->setModel(appModel);
118  this->initialize();
119  }
120 }
121 
122 template <class Scalar>
124  std::string stepperType,
125  Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau,
126  Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau)
127 {
128  if (stepperType == "") stepperType = "IMEX RK SSP2";
129 
130  if (stepperType == "IMEX RK 1st order") {
131  {
132  // Explicit Tableau
134  int NumStages = 2;
135  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
138  const Scalar one = ST::one();
139  const Scalar zero = ST::zero();
140 
141  // Fill A:
142  A(0, 0) = zero;
143  A(0, 1) = zero;
144  A(1, 0) = one;
145  A(1, 1) = zero;
146 
147  // Fill b:
148  b(0) = one;
149  b(1) = zero;
150 
151  // Fill c:
152  c(0) = zero;
153  c(1) = one;
154 
155  int order = 1;
156 
157  auto expTableau = Teuchos::rcp(
158  new RKButcherTableau<Scalar>("Explicit Tableau - IMEX RK 1st order",
159  A, b, c, order, order, order));
160  expTableau->setTVD(true);
161  expTableau->setTVDCoeff(2.0);
162 
163  this->setExplicitTableau(expTableau);
164  }
165  {
166  // Implicit Tableau
168  int NumStages = 2;
169  const Scalar sspcoef = std::numeric_limits<Scalar>::max();
170  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
173  const Scalar one = ST::one();
174  const Scalar zero = ST::zero();
175 
176  // Fill A:
177  A(0, 0) = zero;
178  A(0, 1) = zero;
179  A(1, 0) = zero;
180  A(1, 1) = one;
181 
182  // Fill b:
183  b(0) = zero;
184  b(1) = one;
185 
186  // Fill c:
187  c(0) = zero;
188  c(1) = one;
189 
190  int order = 1;
191 
192  auto impTableau = Teuchos::rcp(
193  new RKButcherTableau<Scalar>("Implicit Tableau - IMEX RK 1st order",
194  A, b, c, order, order, order));
195  impTableau->setTVD(true);
196  impTableau->setTVDCoeff(sspcoef);
197 
198  this->setImplicitTableau(impTableau);
199  }
200  this->setStepperName("IMEX RK 1st order");
201  this->setStepperType("IMEX RK 1st order");
202  this->setOrder(1);
203  }
204  else if (stepperType == "SSP1_111") {
205  {
206  // Explicit Tableau
208  const int NumStages = 1;
209  const int order = 1;
210  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
213  const Scalar one = ST::one();
214  const Scalar zero = ST::zero();
215 
216  // Fill A:
217  A(0, 0) = zero;
218 
219  // Fill b:
220  b(0) = one;
221 
222  // Fill c:
223  c(0) = zero;
224 
225  auto expTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
226  "Explicit Tableau - SSP1_111", A, b, c, order, order, order));
227  expTableau->setTVD(true);
228  expTableau->setTVDCoeff(1.0);
229 
230  this->setExplicitTableau(expTableau);
231  }
232  {
233  // Implicit Tableau
235  const int NumStages = 1;
236  const int order = 1;
237  const Scalar sspcoef = std::numeric_limits<Scalar>::max();
238  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
241  const Scalar one = ST::one();
242 
243  // Fill A:
244  A(0, 0) = one;
245 
246  // Fill b:
247  b(0) = one;
248 
249  // Fill c:
250  c(0) = one;
251 
252  auto impTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
253  "Implicit Tableau - SSP1_111", A, b, c, order, order, order));
254  impTableau->setTVD(true);
255  impTableau->setTVDCoeff(sspcoef);
256 
257  this->setImplicitTableau(impTableau);
258  }
259  this->setStepperName("SSP1_111");
260  this->setStepperType("SSP1_111");
261  this->setOrder(1);
262  }
263  else if (stepperType == "IMEX RK SSP2" || stepperType == "SSP2_222_L") {
264  // Explicit Tableau
265  auto stepperERK = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
266  this->setExplicitTableau(stepperERK->getTableau());
267 
268  // Implicit Tableau
269  auto stepperSDIRK = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
270  stepperSDIRK->setGammaType("2nd Order L-stable");
271  this->setImplicitTableau(stepperSDIRK->getTableau());
272 
273  this->setStepperName("IMEX RK SSP2");
274  this->setStepperType("IMEX RK SSP2");
275  this->setOrder(2);
276  }
277  else if (stepperType == "SSP2_222" || stepperType == "SSP2_222_A") {
278  // Explicit Tableau
279  auto stepperERK = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
280  this->setExplicitTableau(stepperERK->getTableau());
281 
282  // Implicit Tableau
283  auto stepperSDIRK = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
284  stepperSDIRK->setGammaType("gamma");
285  stepperSDIRK->setGamma(0.5);
286  this->setImplicitTableau(stepperSDIRK->getTableau());
287 
288  this->setStepperName("SSP2_222");
289  this->setStepperType("SSP2_222");
290  this->setOrder(2);
291  }
292  else if (stepperType == "IMEX RK SSP3" || stepperType == "SSP3_332") {
293  // Explicit Tableau
294  auto stepperERK = Teuchos::rcp(new StepperERK_3Stage3rdOrderTVD<Scalar>());
295  this->setExplicitTableau(stepperERK->getTableau());
296 
297  // Implicit Tableau
298  auto stepperSDIRK = Teuchos::rcp(new StepperSDIRK_3Stage2ndOrder<Scalar>());
299  this->setImplicitTableau(stepperSDIRK->getTableau());
300 
301  this->setStepperName("IMEX RK SSP3");
302  this->setStepperType("IMEX RK SSP3");
303  this->setOrder(2);
304  }
305  else if (stepperType == "IMEX RK ARS 233" || stepperType == "ARS 233") {
307  int NumStages = 3;
308  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
311  const Scalar one = ST::one();
312  const Scalar zero = ST::zero();
313  const Scalar onehalf = ST::one() / (2 * ST::one());
314  const Scalar gamma = (3 * one + ST::squareroot(3 * one)) / (6 * one);
315  {
316  // Explicit Tableau
317  // Fill A:
318  A(0, 0) = zero;
319  A(0, 1) = zero;
320  A(0, 2) = zero;
321  A(1, 0) = gamma;
322  A(1, 1) = zero;
323  A(1, 2) = zero;
324  A(2, 0) = (gamma - 1.0);
325  A(2, 1) = (2.0 - 2.0 * gamma);
326  A(2, 2) = zero;
327 
328  // Fill b:
329  b(0) = zero;
330  b(1) = onehalf;
331  b(2) = onehalf;
332 
333  // Fill c:
334  c(0) = zero;
335  c(1) = gamma;
336  c(2) = one - gamma;
337 
338  int order = 2;
339 
340  auto expTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
341  "Partition IMEX-RK Explicit Stepper", A, b, c, order, order, order));
342 
343  this->setExplicitTableau(expTableau);
344  }
345  {
346  // Implicit Tableau
347  // Fill A:
348  A(0, 0) = zero;
349  A(0, 1) = zero;
350  A(0, 2) = zero;
351  A(1, 0) = zero;
352  A(1, 1) = gamma;
353  A(1, 2) = zero;
354  A(2, 0) = zero;
355  A(2, 1) = (1.0 - 2.0 * gamma);
356  A(2, 2) = gamma;
357 
358  // Fill b:
359  b(0) = zero;
360  b(1) = onehalf;
361  b(2) = onehalf;
362 
363  // Fill c:
364  c(0) = zero;
365  c(1) = gamma;
366  c(2) = one - gamma;
367 
368  int order = 3;
369 
370  auto impTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
371  "Partition IMEX-RK Implicit Stepper", A, b, c, order, order, order));
372 
373  this->setImplicitTableau(impTableau);
374  }
375  this->setStepperName("IMEX RK ARS 233");
376  this->setStepperType("IMEX RK ARS 233");
377  this->setOrder(3);
378  }
379  else if (stepperType == "General IMEX RK") {
380  if (explicitTableau == Teuchos::null) {
381  // Default Explicit Tableau (i.e., IMEX RK SSP2)
382  auto stepperERK = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
383  this->setExplicitTableau(stepperERK->getTableau());
384  }
385  else {
386  this->setExplicitTableau(explicitTableau);
387  }
388 
389  if (explicitTableau == Teuchos::null) {
390  // Default Implicit Tableau (i.e., IMEX RK SSP2)
391  auto stepperSDIRK =
393  stepperSDIRK->setGammaType("2nd Order L-stable");
394  this->setImplicitTableau(stepperSDIRK->getTableau());
395  }
396  else {
397  this->setImplicitTableau(implicitTableau);
398  }
399 
400  this->setStepperName("General IMEX RK");
401  this->setStepperType("General IMEX RK");
402  this->setOrder(1);
403  }
404  else {
406  true, std::logic_error,
407  "Error - Not a valid StepperIMEX_RK type! Stepper Type = "
408  << stepperType << "\n"
409  << " Current valid types are: \n"
410  << " 'IMEX RK 1st order\n"
411  << " 'SSP1_111\n"
412  << " 'IMEX RK SSP2' ('SSP2_222_L')\n"
413  << " 'SSP2_222' ('SSP2_222_A')\n"
414  << " 'IMEX RK SSP3' ('SSP3_332')\n"
415  << " 'IMEX RK ARS 233' ('ARS 233')\n"
416  << " 'General IMEX RK'\n");
417  }
418 
420  explicitTableau_ == Teuchos::null, std::runtime_error,
421  "Error - StepperIMEX_RK - Explicit tableau is null!");
423  implicitTableau_ == Teuchos::null, std::runtime_error,
424  "Error - StepperIMEX_RK - Implicit tableau is null!");
426  explicitTableau_->numStages() != implicitTableau_->numStages(),
427  std::runtime_error,
428  "Error - StepperIMEX_RK - Number of stages do not match!\n"
429  << " Explicit tableau = " << explicitTableau_->description() << "\n"
430  << " number of stages = " << explicitTableau_->numStages() << "\n"
431  << " Implicit tableau = " << implicitTableau_->description() << "\n"
432  << " number of stages = " << implicitTableau_->numStages() << "\n");
433 
434  this->isInitialized_ = false;
435 }
436 
437 template <class Scalar>
439  Teuchos::RCP<Teuchos::ParameterList> pl, std::string stepperType)
440 {
441  using Teuchos::RCP;
442  if (stepperType == "") {
443  if (pl == Teuchos::null)
444  stepperType = "IMEX RK SSP2";
445  else
446  stepperType = pl->get<std::string>("Stepper Type", "IMEX RK SSP2");
447  }
448 
449  if (stepperType != "General IMEX RK") {
450  this->setTableaus(stepperType);
451  }
452  else {
453  if (pl != Teuchos::null) {
456  if (pl->isSublist("IMEX-RK Explicit Stepper")) {
457  RCP<Teuchos::ParameterList> explicitPL =
459  pl->sublist("IMEX-RK Explicit Stepper")));
460  auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
461  auto stepperTemp = sf->createStepper(explicitPL, Teuchos::null);
462  auto stepperERK = Teuchos::rcp_dynamic_cast<StepperExplicitRK<Scalar> >(
463  stepperTemp, true);
465  stepperERK == Teuchos::null, std::logic_error,
466  "Error - The explicit component of a general IMEX RK stepper was "
467  "not specified as an ExplicitRK stepper");
468  explicitTableau = stepperERK->getTableau();
469  }
470 
471  if (pl->isSublist("IMEX-RK Implicit Stepper")) {
472  RCP<Teuchos::ParameterList> implicitPL =
474  pl->sublist("IMEX-RK Implicit Stepper")));
475  auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
476  auto stepperTemp = sf->createStepper(implicitPL, Teuchos::null);
477  auto stepperDIRK =
478  Teuchos::rcp_dynamic_cast<StepperDIRK<Scalar> >(stepperTemp, true);
480  stepperDIRK == Teuchos::null, std::logic_error,
481  "Error - The implicit component of a general IMEX RK stepper was "
482  "not specified as an DIRK stepper");
483  implicitTableau = stepperDIRK->getTableau();
484  }
485 
487  !(explicitTableau != Teuchos::null &&
488  implicitTableau != Teuchos::null),
489  std::logic_error,
490  "Error - A parameter list was used to setup a general IMEX RK "
491  "stepper, but did not "
492  "specify both an explicit and an implicit tableau!\n");
493 
494  this->setTableaus(stepperType, explicitTableau, implicitTableau);
495  this->setOrder(pl->get<int>("overall order", 1));
496  }
497  }
498 }
499 
500 template <class Scalar>
502  Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau)
503 {
505  explicitTableau->isImplicit() == true, std::logic_error,
506  "Error - Received an implicit Tableau for setExplicitTableau()!\n"
507  << " Tableau = " << explicitTableau->description() << "\n");
508  explicitTableau_ = explicitTableau;
509 
510  this->isInitialized_ = false;
511 }
512 
513 template <class Scalar>
515  Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau)
516 {
518  implicitTableau->isDIRK() != true, std::logic_error,
519  "Error - Did not receive a DIRK Tableau for setImplicitTableau()!\n"
520  << " Tableau = " << implicitTableau->description() << "\n");
521  implicitTableau_ = implicitTableau;
522 
523  this->isInitialized_ = false;
524 }
525 
526 template <class Scalar>
528  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
529 {
530  using Teuchos::RCP;
531  using Teuchos::rcp_const_cast;
532  using Teuchos::rcp_dynamic_cast;
533  RCP<Thyra::ModelEvaluator<Scalar> > ncModel =
534  rcp_const_cast<Thyra::ModelEvaluator<Scalar> >(appModel);
535  RCP<WrapperModelEvaluatorPairIMEX_Basic<Scalar> > modelPairIMEX =
536  rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX_Basic<Scalar> >(ncModel);
538  modelPairIMEX == Teuchos::null, std::logic_error,
539  "Error - StepperIMEX_RK::setModel() was given a ModelEvaluator that\n"
540  << " could not be cast to a WrapperModelEvaluatorPairIMEX_Basic!\n"
541  << " From: " << appModel << "\n"
542  << " To : " << modelPairIMEX << "\n"
543  << " Likely have given the wrong ModelEvaluator to this Stepper.\n");
544 
545  setModelPair(modelPairIMEX);
546 
547  TEUCHOS_TEST_FOR_EXCEPTION(this->solver_ == Teuchos::null, std::logic_error,
548  "Error - Solver is not set!\n");
549  this->solver_->setModel(modelPairIMEX);
550 
551  this->isInitialized_ = false;
552 }
553 
559 template <class Scalar>
562  modelPairIMEX)
563 {
565  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
566  this->wrapperModel_);
567  validExplicitODE(modelPairIMEX->getExplicitModel());
568  validImplicitODE_DAE(modelPairIMEX->getImplicitModel());
569  wrapperModelPairIMEX = modelPairIMEX;
570  wrapperModelPairIMEX->initialize();
571  int expXDim = wrapperModelPairIMEX->getExplicitModel()->get_x_space()->dim();
572  int impXDim = wrapperModelPairIMEX->getImplicitModel()->get_x_space()->dim();
574  expXDim != impXDim, std::logic_error,
575  "Error - \n"
576  << " Explicit and Implicit x vectors are incompatible!\n"
577  << " Explicit vector dim = " << expXDim << "\n"
578  << " Implicit vector dim = " << impXDim << "\n");
579 
580  this->wrapperModel_ = wrapperModelPairIMEX;
581 
582  this->isInitialized_ = false;
583 }
584 
590 template <class Scalar>
592  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& explicitModel,
593  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& implicitModel)
594 {
595  validExplicitODE(explicitModel);
596  validImplicitODE_DAE(implicitModel);
597  this->wrapperModel_ =
599  explicitModel, implicitModel));
600 
601  this->isInitialized_ = false;
602 }
603 
604 template <class Scalar>
606 {
608  this->wrapperModel_ == Teuchos::null, std::logic_error,
609  "Error - Need to set the model, setModel(), before calling "
610  "StepperIMEX_RK::initialize()\n");
611 
612  // Initialize the stage vectors
613  const int numStages = explicitTableau_->numStages();
614  stageF_.resize(numStages);
615  stageG_.resize(numStages);
616  for (int i = 0; i < numStages; i++) {
617  stageF_[i] = Thyra::createMember(this->wrapperModel_->get_f_space());
618  stageG_[i] = Thyra::createMember(this->wrapperModel_->get_f_space());
619  assign(stageF_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
620  assign(stageG_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
621  }
622 
623  xTilde_ = Thyra::createMember(this->wrapperModel_->get_x_space());
624  assign(xTilde_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
625 
627 }
628 
629 template <class Scalar>
631  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
632 {
633  using Teuchos::RCP;
634 
635  int numStates = solutionHistory->getNumStates();
636 
638  numStates < 1, std::logic_error,
639  "Error - setInitialConditions() needs at least one SolutionState\n"
640  " to set the initial condition. Number of States = "
641  << numStates);
642 
643  if (numStates > 1) {
644  RCP<Teuchos::FancyOStream> out = this->getOStream();
645  Teuchos::OSTab ostab(out, 1, "StepperIMEX_RK::setInitialConditions()");
646  *out << "Warning -- SolutionHistory has more than one state!\n"
647  << "Setting the initial conditions on the currentState.\n"
648  << std::endl;
649  }
650 
651  RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
652  RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
653 
654  // Use x from inArgs as ICs, if needed.
655  auto inArgs = this->wrapperModel_->getNominalValues();
656  if (x == Teuchos::null) {
658  (x == Teuchos::null) && (inArgs.get_x() == Teuchos::null),
659  std::logic_error,
660  "Error - setInitialConditions() needs the ICs from the "
661  << "SolutionHistory\n"
662  << " or getNominalValues()!\n");
663 
664  x = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
665  initialState->setX(x);
666  }
667 
668  // Perform IC Consistency
669  std::string icConsistency = this->getICConsistency();
671  icConsistency != "None", std::logic_error,
672  "Error - setInitialConditions() requested a consistency of '"
673  << icConsistency
674  << "'.\n But only 'None' is available for IMEX-RK!\n");
675 
677  this->getUseFSAL(), std::logic_error,
678  "Error - The First-Same-As-Last (FSAL) principle is not "
679  << "available for IMEX-RK. Set useFSAL=false.\n");
680 }
681 
682 template <typename Scalar>
684  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& X, Scalar time,
685  Scalar stepSize, Scalar stageNumber,
686  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& G) const
687 {
688  typedef Thyra::ModelEvaluatorBase MEB;
690  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
691  this->wrapperModel_);
692  MEB::InArgs<Scalar> inArgs = wrapperModelPairIMEX->createInArgs();
693  inArgs.set_x(X);
694  if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
695  if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(stepSize);
696  if (inArgs.supports(MEB::IN_ARG_stage_number))
697  inArgs.set_stage_number(stageNumber);
698 
699  // For model evaluators whose state function f(x, x_dot, t) describes
700  // an implicit ODE, and which accept an optional x_dot input argument,
701  // make sure the latter is set to null in order to request the evaluation
702  // of a state function corresponding to the explicit ODE formulation
703  // x_dot = f(x, t)
704  if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
705 
706  MEB::OutArgs<Scalar> outArgs = wrapperModelPairIMEX->createOutArgs();
707  outArgs.set_f(G);
708 
709  wrapperModelPairIMEX->getImplicitModel()->evalModel(inArgs, outArgs);
710  Thyra::Vt_S(G.ptr(), -1.0);
711 }
712 
713 template <typename Scalar>
715  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& X, Scalar time,
716  Scalar stepSize, Scalar stageNumber,
717  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& F) const
718 {
719  typedef Thyra::ModelEvaluatorBase MEB;
720 
722  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
723  this->wrapperModel_);
724  MEB::InArgs<Scalar> inArgs =
725  wrapperModelPairIMEX->getExplicitModel()->createInArgs();
726  inArgs.set_x(X);
727  if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
728  if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(stepSize);
729  if (inArgs.supports(MEB::IN_ARG_stage_number))
730  inArgs.set_stage_number(stageNumber);
731 
732  // For model evaluators whose state function f(x, x_dot, t) describes
733  // an implicit ODE, and which accept an optional x_dot input argument,
734  // make sure the latter is set to null in order to request the evaluation
735  // of a state function corresponding to the explicit ODE formulation
736  // x_dot = f(x, t)
737  if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
738 
739  MEB::OutArgs<Scalar> outArgs =
740  wrapperModelPairIMEX->getExplicitModel()->createOutArgs();
741  outArgs.set_f(F);
742 
743  wrapperModelPairIMEX->getExplicitModel()->evalModel(inArgs, outArgs);
744  Thyra::Vt_S(F.ptr(), -1.0);
745 }
746 
747 template <class Scalar>
749  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
750 {
751  this->checkInitialized();
752 
753  using Teuchos::RCP;
756 
757  TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperIMEX_RK::takeStep()");
758  {
760  solutionHistory->getNumStates() < 2, std::logic_error,
761  "Error - StepperIMEX_RK<Scalar>::takeStep(...)\n"
762  << "Need at least two SolutionStates for IMEX_RK.\n"
763  << " Number of States = " << solutionHistory->getNumStates()
764  << "\nTry setting in \"Solution History\" "
765  << "\"Storage Type\" = \"Undo\"\n"
766  << " or \"Storage Type\" = \"Static\" and "
767  << "\"Storage Limit\" = \"2\"\n");
768 
769  RCP<SolutionState<Scalar> > currentState =
770  solutionHistory->getCurrentState();
771  RCP<SolutionState<Scalar> > workingState =
772  solutionHistory->getWorkingState();
773  const Scalar dt = workingState->getTimeStep();
774  const Scalar time = currentState->getTime();
775 
776  const int numStages = explicitTableau_->numStages();
777  const SerialDenseMatrix<int, Scalar>& AHat = explicitTableau_->A();
778  const SerialDenseVector<int, Scalar>& bHat = explicitTableau_->b();
779  const SerialDenseVector<int, Scalar>& cHat = explicitTableau_->c();
780  const SerialDenseMatrix<int, Scalar>& A = implicitTableau_->A();
781  const SerialDenseVector<int, Scalar>& b = implicitTableau_->b();
782  const SerialDenseVector<int, Scalar>& c = implicitTableau_->c();
783 
784  bool pass = true;
785  Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
786 
787  RCP<StepperIMEX_RK<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
788  this->stepperRKAppAction_->execute(
789  solutionHistory, thisStepper,
791 
792  // Compute stage solutions
793  for (int i = 0; i < numStages; ++i) {
794  this->setStageNumber(i);
795  Thyra::assign(xTilde_.ptr(), *(currentState->getX()));
796  for (int j = 0; j < i; ++j) {
797  if (AHat(i, j) != Teuchos::ScalarTraits<Scalar>::zero())
798  Thyra::Vp_StV(xTilde_.ptr(), -dt * AHat(i, j), *(stageF_[j]));
799  if (A(i, j) != Teuchos::ScalarTraits<Scalar>::zero())
800  Thyra::Vp_StV(xTilde_.ptr(), -dt * A(i, j), *(stageG_[j]));
801  }
802 
803  this->stepperRKAppAction_->execute(
804  solutionHistory, thisStepper,
806 
807  Scalar ts = time + c(i) * dt;
808  Scalar tHats = time + cHat(i) * dt;
809  if (A(i, i) == Teuchos::ScalarTraits<Scalar>::zero()) {
810  // Explicit stage for the ImplicitODE_DAE
811  bool isNeeded = false;
812  for (int k = i + 1; k < numStages; ++k)
813  if (A(k, i) != 0.0) isNeeded = true;
814  if (b(i) != 0.0) isNeeded = true;
815  if (isNeeded == false) {
816  // stageG_[i] is not needed.
817  assign(stageG_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
818  }
819  else {
820  Thyra::assign(workingState->getX().ptr(), *xTilde_);
821  evalImplicitModelExplicitly(workingState->getX(), ts, dt, i,
822  stageG_[i]);
823  }
824  }
825  else {
826  // Implicit stage for the ImplicitODE_DAE
827  const Scalar alpha = Scalar(1.0) / (dt * A(i, i));
828  const Scalar beta = Scalar(1.0);
829 
830  // Setup TimeDerivative
833  alpha, xTilde_.getConst()));
834 
836  timeDer, dt, alpha, beta, SOLVE_FOR_X, i));
837 
838  this->stepperRKAppAction_->execute(
839  solutionHistory, thisStepper,
841 
842  const Thyra::SolveStatus<Scalar> sStatus =
843  this->solveImplicitODE(workingState->getX(), stageG_[i], ts, p);
844 
845  if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass = false;
846 
847  this->stepperRKAppAction_->execute(
848  solutionHistory, thisStepper,
850 
851  // Update contributions to stage values
852  Thyra::V_StVpStV(stageG_[i].ptr(), -alpha, *workingState->getX(), alpha,
853  *xTilde_);
854  }
855 
856  this->stepperRKAppAction_->execute(
857  solutionHistory, thisStepper,
859  evalExplicitModel(workingState->getX(), tHats, dt, i, stageF_[i]);
860  this->stepperRKAppAction_->execute(
861  solutionHistory, thisStepper,
863  }
864 
865  // reset the stage number
866  this->setStageNumber(-1);
867 
868  // Sum for solution: x_n = x_n-1 - dt*Sum{ bHat(i)*f(i) + b(i)*g(i) }
869  Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
870  for (int i = 0; i < numStages; ++i) {
871  if (bHat(i) != Teuchos::ScalarTraits<Scalar>::zero())
872  Thyra::Vp_StV((workingState->getX()).ptr(), -dt * bHat(i),
873  *(stageF_[i]));
875  Thyra::Vp_StV((workingState->getX()).ptr(), -dt * b(i), *(stageG_[i]));
876  }
877 
878  if (pass == true)
879  workingState->setSolutionStatus(Status::PASSED);
880  else
881  workingState->setSolutionStatus(Status::FAILED);
882  workingState->setOrder(this->getOrder());
883  workingState->computeNorms(currentState);
884  this->stepperRKAppAction_->execute(
885  solutionHistory, thisStepper,
887  }
888  return;
889 }
890 
897 template <class Scalar>
900 {
902  rcp(new StepperState<Scalar>(this->getStepperType()));
903  return stepperState;
904 }
905 
906 template <class Scalar>
908  Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel) const
909 {
910  out.setOutputToRootOnly(0);
911 
912  out << std::endl;
913  Stepper<Scalar>::describe(out, verbLevel);
914  StepperImplicit<Scalar>::describe(out, verbLevel);
915 
916  out << "--- StepperIMEX_RK_Partition ---\n";
917  out << " explicitTableau_ = " << explicitTableau_ << std::endl;
918  if (verbLevel == Teuchos::VERB_HIGH)
919  explicitTableau_->describe(out, verbLevel);
920  out << " implicitTableau_ = " << implicitTableau_ << std::endl;
921  if (verbLevel == Teuchos::VERB_HIGH)
922  implicitTableau_->describe(out, verbLevel);
923  out << " xTilde_ = " << xTilde_ << std::endl;
924  out << " stageF_.size() = " << stageF_.size() << std::endl;
925  int numStages = stageF_.size();
926  for (int i = 0; i < numStages; ++i)
927  out << " stageF_[" << i << "] = " << stageF_[i] << std::endl;
928  out << " stageG_.size() = " << stageG_.size() << std::endl;
929  numStages = stageG_.size();
930  for (int i = 0; i < numStages; ++i)
931  out << " stageG_[" << i << "] = " << stageG_[i] << std::endl;
932  out << " stepperRKAppAction_= " << this->stepperRKAppAction_ << std::endl;
933  out << " order_ = " << order_ << std::endl;
934  out << "--------------------------------" << std::endl;
935 }
936 
937 template <class Scalar>
939 {
940  out.setOutputToRootOnly(0);
941  bool isValidSetup = true;
942 
943  if (!Stepper<Scalar>::isValidSetup(out)) isValidSetup = false;
944 
946  Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
947  this->wrapperModel_);
948 
949  if (wrapperModelPairIMEX->getExplicitModel() == Teuchos::null) {
950  isValidSetup = false;
951  out << "The explicit ModelEvaluator is not set!\n";
952  }
953 
954  if (wrapperModelPairIMEX->getImplicitModel() == Teuchos::null) {
955  isValidSetup = false;
956  out << "The implicit ModelEvaluator is not set!\n";
957  }
958 
959  if (this->wrapperModel_ == Teuchos::null) {
960  isValidSetup = false;
961  out << "The wrapper ModelEvaluator is not set!\n";
962  }
963 
964  if (this->stepperRKAppAction_ == Teuchos::null) {
965  isValidSetup = false;
966  out << "The AppAction is not set!\n";
967  }
968 
969  if (explicitTableau_ == Teuchos::null) {
970  isValidSetup = false;
971  out << "The explicit tableau is not set!\n";
972  }
973 
974  if (implicitTableau_ == Teuchos::null) {
975  isValidSetup = false;
976  out << "The implicit tableau is not set!\n";
977  }
978 
979  return isValidSetup;
980 }
981 
982 template <class Scalar>
985 {
986  auto pl = this->getValidParametersBasicImplicit();
987  pl->template set<int>("overall order", this->getOrder());
988 
989  auto explicitStepper = Teuchos::rcp(new StepperERK_General<Scalar>());
990  explicitStepper->setTableau(
991  explicitTableau_->A(), explicitTableau_->b(), explicitTableau_->c(),
992  explicitTableau_->order(), explicitTableau_->orderMin(),
993  explicitTableau_->orderMax(), explicitTableau_->bstar());
994  pl->set("IMEX-RK Explicit Stepper", *explicitStepper->getValidParameters());
995 
996  auto implicitStepper = Teuchos::rcp(new StepperERK_General<Scalar>());
997  implicitStepper->setTableau(
998  implicitTableau_->A(), implicitTableau_->b(), implicitTableau_->c(),
999  implicitTableau_->order(), implicitTableau_->orderMin(),
1000  implicitTableau_->orderMax(), implicitTableau_->bstar());
1001  pl->set("IMEX-RK Implicit Stepper", *implicitStepper->getValidParameters());
1002 
1003  return pl;
1004 }
1005 
1006 // Nonmember constructor - ModelEvaluator and ParameterList
1007 // ------------------------------------------------------------------------
1008 template <class Scalar>
1010  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1011  std::string stepperType, Teuchos::RCP<Teuchos::ParameterList> pl)
1012 {
1013  auto stepper = Teuchos::rcp(new StepperIMEX_RK<Scalar>(stepperType));
1014  stepper->setStepperImplicitValues(pl);
1015  stepper->setTableaus(pl, stepperType);
1016 
1017  if (model != Teuchos::null) {
1018  stepper->setModel(model);
1019  stepper->initialize();
1020  }
1021 
1022  return stepper;
1023 }
1024 
1025 } // namespace Tempus
1026 #endif // Tempus_StepperIMEX_RK_impl_hpp
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void setImplicitTableau(Teuchos::RCP< const RKButcherTableau< Scalar > > implicitTableau)
Set the implicit tableau from tableau.
Explicit Runge-Kutta time stepper.
T & get(const std::string &name, T def_value)
General Explicit Runge-Kutta Butcher Tableau.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Implicit-Explicit Runge-Kutta (IMEX-RK) time stepper.
virtual void setExplicitTableau(Teuchos::RCP< const RKButcherTableau< Scalar > > explicitTableau)
Set the explicit tableau from tableau.
ModelEvaluator pair for implicit and explicit (IMEX) evaluations.
void validExplicitODE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Validate that the model supports explicit ODE evaluation, f(x,t) [=xdot].
virtual void initialize()
Initialize after construction and changing input parameters.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
StepperIMEX_RK(std::string stepperType="IMEX RK SSP2")
Default constructor.
Thyra Base interface for time steppers.
virtual void setModelPair(const Teuchos::RCP< WrapperModelEvaluatorPairIMEX_Basic< Scalar > > &mePair)
Create WrapperModelPairIMEX from user-supplied ModelEvaluator pair.
Diagonally Implicit Runge-Kutta (DIRK) time stepper.
StepperState is a simple class to hold state information about the stepper.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Application Action for StepperRKBase.
bool isSublist(const std::string &name) const
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
virtual void initialize()
Initialize during construction and after changing input parameters.
void evalExplicitModel(const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &X, Scalar time, Scalar stepSize, Scalar stageNumber, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &F) const
virtual void setTableaus(std::string stepperType="", Teuchos::RCP< const RKButcherTableau< Scalar > > explicitTableau=Teuchos::null, Teuchos::RCP< const RKButcherTableau< Scalar > > implicitTableau=Teuchos::null)
Set both the explicit and implicit tableau from ParameterList.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
ESolveStatus solveStatus
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
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 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...
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Time-derivative interface for IMEX RK.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Set the model.
Teuchos::RCP< StepperIMEX_RK< Scalar > > createStepperIMEX_RK(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, std::string stepperType, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
void validImplicitODE_DAE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Validate ME supports implicit ODE/DAE evaluation, f(xdot,x,t) [= 0].
Solve for x and determine xDot from x.