Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_StepperRKButcherTableau.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_StepperRKButcherTableau_hpp
10 #define Tempus_StepperRKButcherTableau_hpp
11 
12 // disable clang warnings
13 #ifdef __clang__
14 #pragma clang system_header
15 #endif
16 
17 #include "Tempus_config.hpp"
18 #include "Tempus_StepperExplicitRK.hpp"
19 #include "Tempus_StepperDIRK.hpp"
21 
22 namespace Tempus {
23 
24 // ----------------------------------------------------------------------------
40 template <class Scalar>
41 class StepperERK_ForwardEuler : virtual public StepperExplicitRK<Scalar> {
42  public:
49  {
50  this->setStepperName("RK Forward Euler");
51  this->setStepperType("RK Forward Euler");
52  this->setupTableau();
53  this->setupDefault();
54  this->setUseFSAL(true);
55  this->setICConsistency("Consistent");
56  this->setICConsistencyCheck(false);
57  }
58 
60  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
61  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
62  bool useEmbedded,
63  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
64  {
65  this->setStepperName("RK Forward Euler");
66  this->setStepperType("RK Forward Euler");
67  this->setupTableau();
68  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
69  useEmbedded, stepperRKAppAction);
70  }
71 
72  std::string getDescription() const
73  {
74  std::ostringstream Description;
75  Description << this->getStepperType() << "\n"
76  << "c = [ 0 ]'\n"
77  << "A = [ 0 ]\n"
78  << "b = [ 1 ]'";
79  return Description.str();
80  }
81 
82  void setUseFSAL(bool a)
83  {
84  this->useFSAL_ = a;
85  this->isInitialized_ = false;
86  }
87 
88  protected:
89  void setupTableau()
90  {
95  A(0, 0) = ST::zero();
96  b(0) = ST::one();
97  c(0) = ST::zero();
98  int order = 1;
99 
101  this->getStepperType(), A, b, c, order, order, order));
102  this->tableau_->setTVD(true);
103  this->tableau_->setTVDCoeff(1.0);
104  }
105 };
106 
108 // ------------------------------------------------------------------------
109 template <class Scalar>
111  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
113 {
114  auto stepper = Teuchos::rcp(new StepperERK_ForwardEuler<Scalar>());
115 
116  // Test for aliases.
117  if (pl != Teuchos::null) {
118  auto stepperType =
119  pl->get<std::string>("Stepper Type", stepper->getStepperType());
120 
122  stepperType != stepper->getStepperType() && stepperType != "RK1",
123  std::logic_error,
124  " ParameterList 'Stepper Type' (='"
125  << stepperType << "')\n"
126  << " does not match type for this Stepper (='"
127  << stepper->getStepperType()
128  << "')\n or one of its aliases ('RK1').\n");
129 
130  // Reset default StepperType.
131  pl->set<std::string>("Stepper Type", stepper->getStepperType());
132  }
133 
134  stepper->setStepperRKValues(pl);
135 
136  if (model != Teuchos::null) {
137  stepper->setModel(model);
138  stepper->initialize();
139  }
140 
141  return stepper;
142 }
143 
144 // ----------------------------------------------------------------------------
163 template <class Scalar>
164 class StepperERK_4Stage4thOrder : virtual public StepperExplicitRK<Scalar> {
165  public:
172  {
173  this->setStepperName("RK Explicit 4 Stage");
174  this->setStepperType("RK Explicit 4 Stage");
175  this->setupTableau();
176  this->setupDefault();
177  this->setUseFSAL(false);
178  this->setICConsistency("None");
179  this->setICConsistencyCheck(false);
180  }
181 
183  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
184  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
185  bool useEmbedded,
186  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
187  {
188  this->setStepperName("RK Explicit 4 Stage");
189  this->setStepperType("RK Explicit 4 Stage");
190  this->setupTableau();
191  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
192  useEmbedded, stepperRKAppAction);
193  }
194 
195  std::string getDescription() const
196  {
197  std::ostringstream Description;
198  Description << this->getStepperType() << "\n"
199  << "\"The\" Runge-Kutta Method (explicit):\n"
200  << "Solving Ordinary Differential Equations I:\n"
201  << "Nonstiff Problems, 2nd Revised Edition\n"
202  << "E. Hairer, S.P. Norsett, G. Wanner\n"
203  << "Table 1.2, pg 138\n"
204  << "c = [ 0 1/2 1/2 1 ]'\n"
205  << "A = [ 0 ] \n"
206  << " [ 1/2 0 ]\n"
207  << " [ 0 1/2 0 ]\n"
208  << " [ 0 0 1 0 ]\n"
209  << "b = [ 1/6 1/3 1/3 1/6 ]'";
210  return Description.str();
211  }
212 
214  {
216  const Scalar one = ST::one();
217  const Scalar zero = ST::zero();
218  const Scalar onehalf = one / (2 * one);
219  const Scalar onesixth = one / (6 * one);
220  const Scalar onethird = one / (3 * one);
221 
222  int NumStages = 4;
223  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
226 
227  // Fill A:
228  A(0, 0) = zero;
229  A(0, 1) = zero;
230  A(0, 2) = zero;
231  A(0, 3) = zero;
232  A(1, 0) = onehalf;
233  A(1, 1) = zero;
234  A(1, 2) = zero;
235  A(1, 3) = zero;
236  A(2, 0) = zero;
237  A(2, 1) = onehalf;
238  A(2, 2) = zero;
239  A(2, 3) = zero;
240  A(3, 0) = zero;
241  A(3, 1) = zero;
242  A(3, 2) = one;
243  A(3, 3) = zero;
244 
245  // Fill b:
246  b(0) = onesixth;
247  b(1) = onethird;
248  b(2) = onethird;
249  b(3) = onesixth;
250 
251  // fill c:
252  c(0) = zero;
253  c(1) = onehalf;
254  c(2) = onehalf;
255  c(3) = one;
256 
257  int order = 4;
258 
260  this->getStepperType(), A, b, c, order, order, order));
261  }
262 };
263 
265 // ------------------------------------------------------------------------
266 template <class Scalar>
269  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
271 {
272  auto stepper = Teuchos::rcp(new StepperERK_4Stage4thOrder<Scalar>());
273  stepper->setStepperRKValues(pl);
274 
275  if (model != Teuchos::null) {
276  stepper->setModel(model);
277  stepper->initialize();
278  }
279 
280  return stepper;
281 }
282 
283 // ----------------------------------------------------------------------------
307 template <class Scalar>
308 class StepperERK_BogackiShampine32 : virtual public StepperExplicitRK<Scalar> {
309  public:
316  {
317  this->setStepperName("Bogacki-Shampine 3(2) Pair");
318  this->setStepperType("Bogacki-Shampine 3(2) Pair");
319  this->setupTableau();
320  this->setupDefault();
321  this->setUseFSAL(true);
322  this->setICConsistency("Consistent");
323  this->setICConsistencyCheck(false);
324  }
325 
327  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
328  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
329  bool useEmbedded,
330  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
331  {
332  this->setStepperName("Bogacki-Shampine 3(2) Pair");
333  this->setStepperType("Bogacki-Shampine 3(2) Pair");
334  this->setupTableau();
335  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
336  useEmbedded, stepperRKAppAction);
337  }
338 
339  std::string getDescription() const
340  {
341  std::ostringstream Description;
342  Description << this->getStepperType() << "\n"
343  << "P. Bogacki and L.F. Shampine.\n"
344  << "A 3(2) pair of Runge–Kutta formulas.\n"
345  << "Applied Mathematics Letters, 2(4):321 – 325, 1989.\n"
346  << "c = [ 0 1/2 3/4 1 ]'\n"
347  << "A = [ 0 ]\n"
348  << " [ 1/2 0 ]\n"
349  << " [ 0 3/4 0 ]\n"
350  << " [ 2/9 1/3 4/9 0 ]\n"
351  << "b = [ 2/9 1/3 4/9 0 ]'\n"
352  << "bstar = [ 7/24 1/4 1/3 1/8 ]'";
353  return Description.str();
354  }
355 
356  void setUseFSAL(bool a)
357  {
358  this->useFSAL_ = a;
359  this->isInitialized_ = false;
360  }
361 
362  protected:
364  {
366  using Teuchos::as;
367  int NumStages = 4;
368  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
372 
373  const Scalar one = ST::one();
374  const Scalar zero = ST::zero();
375  const Scalar onehalf = one / (2 * one);
376  const Scalar onethird = one / (3 * one);
377  const Scalar threefourths = (3 * one) / (4 * one);
378  const Scalar twoninths = (2 * one) / (9 * one);
379  const Scalar fourninths = (4 * one) / (9 * one);
380 
381  // Fill A:
382  A(0, 0) = zero;
383  A(0, 1) = zero;
384  A(0, 2) = zero;
385  A(0, 3) = zero;
386  A(1, 0) = onehalf;
387  A(1, 1) = zero;
388  A(1, 2) = zero;
389  A(1, 3) = zero;
390  A(2, 0) = zero;
391  A(2, 1) = threefourths;
392  A(2, 2) = zero;
393  A(2, 3) = zero;
394  A(3, 0) = twoninths;
395  A(3, 1) = onethird;
396  A(3, 2) = fourninths;
397  A(3, 3) = zero;
398 
399  // Fill b:
400  b(0) = A(3, 0);
401  b(1) = A(3, 1);
402  b(2) = A(3, 2);
403  b(3) = A(3, 3);
404 
405  // Fill c:
406  c(0) = zero;
407  c(1) = onehalf;
408  c(2) = threefourths;
409  c(3) = one;
410 
411  // Fill bstar
412  bstar(0) = as<Scalar>(7 * one / (24 * one));
413  bstar(1) = as<Scalar>(1 * one / (4 * one));
414  bstar(2) = as<Scalar>(1 * one / (3 * one));
415  bstar(3) = as<Scalar>(1 * one / (8 * one));
416  int order = 3;
417 
419  this->getStepperType(), A, b, c, order, order, order, bstar));
420  }
421 };
422 
424 // ------------------------------------------------------------------------
425 template <class Scalar>
428  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
430 {
432  stepper->setStepperRKValues(pl);
433 
434  if (model != Teuchos::null) {
435  stepper->setModel(model);
436  stepper->initialize();
437  }
438 
439  return stepper;
440 }
441 
442 // ----------------------------------------------------------------------------
468 template <class Scalar>
469 class StepperERK_Merson45 : virtual public StepperExplicitRK<Scalar> {
470  public:
477  {
478  this->setStepperName("Merson 4(5) Pair");
479  this->setStepperType("Merson 4(5) Pair");
480  this->setupTableau();
481  this->setupDefault();
482  this->setUseFSAL(false);
483  this->setICConsistency("None");
484  this->setICConsistencyCheck(false);
485  }
486 
488  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
489  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
490  bool useEmbedded,
491  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
492  {
493  this->setStepperName("Merson 4(5) Pair");
494  this->setStepperType("Merson 4(5) Pair");
495  this->setupTableau();
496  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
497  useEmbedded, stepperRKAppAction);
498  }
499 
500  std::string getDescription() const
501  {
502  std::ostringstream Description;
503  Description << this->getStepperType() << "\n"
504  << "Solving Ordinary Differential Equations I:\n"
505  << "Nonstiff Problems, 2nd Revised Edition\n"
506  << "E. Hairer, S.P. Norsett, G. Wanner\n"
507  << "Table 4.1, pg 167\n"
508  << "c = [ 0 1/3 1/3 1/2 1 ]'\n"
509  << "A = [ 0 ]\n"
510  << " [ 1/3 0 ]\n"
511  << " [ 1/6 1/6 0 ]\n"
512  << " [ 1/8 0 3/8 0 ]\n"
513  << " [ 1/2 0 -3/2 2 0 ]\n"
514  << "b = [ 1/6 0 0 2/3 1/6 ]'\n"
515  << "bstar = [ 1/10 0 3/10 2/5 1/5 ]'";
516  return Description.str();
517  }
518 
519  protected:
521  {
523  using Teuchos::as;
524  int NumStages = 5;
525  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages, true);
526  Teuchos::SerialDenseVector<int, Scalar> b(NumStages, true);
527  Teuchos::SerialDenseVector<int, Scalar> c(NumStages, true);
528  Teuchos::SerialDenseVector<int, Scalar> bstar(NumStages, true);
529 
530  const Scalar one = ST::one();
531  const Scalar zero = ST::zero();
532 
533  // Fill A:
534  A(1, 0) = as<Scalar>(one / (3 * one));
535 
536  A(2, 0) = as<Scalar>(one / (6 * one));
537  A(2, 1) = as<Scalar>(one / (6 * one));
538 
539  A(3, 0) = as<Scalar>(one / (8 * one));
540  A(3, 2) = as<Scalar>(3 * one / (8 * one));
541 
542  A(4, 0) = as<Scalar>(one / (2 * one));
543  A(4, 2) = as<Scalar>(-3 * one / (2 * one));
544  A(4, 3) = 2 * one;
545 
546  // Fill b:
547  b(0) = as<Scalar>(one / (6 * one));
548  b(3) = as<Scalar>(2 * one / (3 * one));
549  b(4) = as<Scalar>(one / (6 * one));
550 
551  // Fill c:
552  c(0) = zero;
553  c(1) = as<Scalar>(1 * one / (3 * one));
554  c(2) = as<Scalar>(1 * one / (3 * one));
555  c(3) = as<Scalar>(1 * one / (2 * one));
556  c(4) = one;
557 
558  // Fill bstar
559  bstar(0) = as<Scalar>(1 * one / (10 * one));
560  bstar(2) = as<Scalar>(3 * one / (10 * one));
561  bstar(3) = as<Scalar>(2 * one / (5 * one));
562  bstar(4) = as<Scalar>(1 * one / (5 * one));
563  int order = 4;
564 
566  this->getStepperType(), A, b, c, order, order, order, bstar));
567  }
568 };
569 
571 // ------------------------------------------------------------------------
572 template <class Scalar>
574  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
576 {
577  auto stepper = Teuchos::rcp(new StepperERK_Merson45<Scalar>());
578  stepper->setStepperRKValues(pl);
579 
580  if (model != Teuchos::null) {
581  stepper->setModel(model);
582  stepper->initialize();
583  }
584 
585  return stepper;
586 }
587 
588 // ----------------------------------------------------------------------------
611 template <class Scalar>
612 class StepperERK_3_8Rule : virtual public StepperExplicitRK<Scalar> {
613  public:
615  {
616  this->setStepperName("RK Explicit 3/8 Rule");
617  this->setStepperType("RK Explicit 3/8 Rule");
618  this->setupTableau();
619  this->setupDefault();
620  this->setUseFSAL(false);
621  this->setICConsistency("None");
622  this->setICConsistencyCheck(false);
623  }
624 
626  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
627  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
628  bool useEmbedded,
629  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
630  {
631  this->setStepperName("RK Explicit 3/8 Rule");
632  this->setStepperType("RK Explicit 3/8 Rule");
633  this->setupTableau();
634  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
635  useEmbedded, stepperRKAppAction);
636  }
637 
638  std::string getDescription() const
639  {
640  std::ostringstream Description;
641  Description << this->getStepperType() << "\n"
642  << "Solving Ordinary Differential Equations I:\n"
643  << "Nonstiff Problems, 2nd Revised Edition\n"
644  << "E. Hairer, S.P. Norsett, G. Wanner\n"
645  << "Table 1.2, pg 138\n"
646  << "c = [ 0 1/3 2/3 1 ]'\n"
647  << "A = [ 0 ]\n"
648  << " [ 1/3 0 ]\n"
649  << " [-1/3 1 0 ]\n"
650  << " [ 1 -1 1 0 ]\n"
651  << "b = [ 1/8 3/8 3/8 1/8 ]'";
652  return Description.str();
653  }
654 
655  protected:
657  {
659  using Teuchos::as;
660  int NumStages = 4;
661  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
664 
665  const Scalar one = ST::one();
666  const Scalar zero = ST::zero();
667  const Scalar onethird = as<Scalar>(one / (3 * one));
668  const Scalar twothirds = as<Scalar>(2 * one / (3 * one));
669  const Scalar oneeighth = as<Scalar>(one / (8 * one));
670  const Scalar threeeighths = as<Scalar>(3 * one / (8 * one));
671 
672  // Fill A:
673  A(0, 0) = zero;
674  A(0, 1) = zero;
675  A(0, 2) = zero;
676  A(0, 3) = zero;
677  A(1, 0) = onethird;
678  A(1, 1) = zero;
679  A(1, 2) = zero;
680  A(1, 3) = zero;
681  A(2, 0) = -onethird;
682  A(2, 1) = one;
683  A(2, 2) = zero;
684  A(2, 3) = zero;
685  A(3, 0) = one;
686  A(3, 1) = -one;
687  A(3, 2) = one;
688  A(3, 3) = zero;
689 
690  // Fill b:
691  b(0) = oneeighth;
692  b(1) = threeeighths;
693  b(2) = threeeighths;
694  b(3) = oneeighth;
695 
696  // Fill c:
697  c(0) = zero;
698  c(1) = onethird;
699  c(2) = twothirds;
700  c(3) = one;
701 
702  int order = 4;
703 
705  this->getStepperType(), A, b, c, order, order, order));
706  }
707 };
708 
710 // ------------------------------------------------------------------------
711 template <class Scalar>
713  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
715 {
716  auto stepper = Teuchos::rcp(new StepperERK_3_8Rule<Scalar>());
717  stepper->setStepperRKValues(pl);
718 
719  if (model != Teuchos::null) {
720  stepper->setModel(model);
721  stepper->initialize();
722  }
723 
724  return stepper;
725 }
726 
727 // ----------------------------------------------------------------------------
750 template <class Scalar>
752  : virtual public StepperExplicitRK<Scalar> {
753  public:
760  {
761  this->setStepperName("RK Explicit 4 Stage 3rd order by Runge");
762  this->setStepperType("RK Explicit 4 Stage 3rd order by Runge");
763  this->setupTableau();
764  this->setupDefault();
765  this->setUseFSAL(false);
766  this->setICConsistency("None");
767  this->setICConsistencyCheck(false);
768  }
769 
771  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
772  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
773  bool useEmbedded,
774  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
775  {
776  this->setStepperName("RK Explicit 4 Stage 3rd order by Runge");
777  this->setStepperType("RK Explicit 4 Stage 3rd order by Runge");
778  this->setupTableau();
779  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
780  useEmbedded, stepperRKAppAction);
781  }
782 
783  std::string getDescription() const
784  {
785  std::ostringstream Description;
786  Description << this->getStepperType() << "\n"
787  << "Solving Ordinary Differential Equations I:\n"
788  << "Nonstiff Problems, 2nd Revised Edition\n"
789  << "E. Hairer, S.P. Norsett, G. Wanner\n"
790  << "Table 1.1, pg 135\n"
791  << "c = [ 0 1/2 1 1 ]'\n"
792  << "A = [ 0 ]\n"
793  << " [ 1/2 0 ]\n"
794  << " [ 0 1 0 ]\n"
795  << " [ 0 0 1 0 ]\n"
796  << "b = [ 1/6 2/3 0 1/6 ]'";
797  return Description.str();
798  }
799 
800  protected:
802  {
804  int NumStages = 4;
805  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
808 
809  const Scalar one = ST::one();
810  const Scalar onehalf = one / (2 * one);
811  const Scalar onesixth = one / (6 * one);
812  const Scalar twothirds = 2 * one / (3 * one);
813  const Scalar zero = ST::zero();
814 
815  // Fill A:
816  A(0, 0) = zero;
817  A(0, 1) = zero;
818  A(0, 2) = zero;
819  A(0, 3) = zero;
820  A(1, 0) = onehalf;
821  A(1, 1) = zero;
822  A(1, 2) = zero;
823  A(1, 3) = zero;
824  A(2, 0) = zero;
825  A(2, 1) = one;
826  A(2, 2) = zero;
827  A(2, 3) = zero;
828  A(3, 0) = zero;
829  A(3, 1) = zero;
830  A(3, 2) = one;
831  A(3, 3) = zero;
832 
833  // Fill b:
834  b(0) = onesixth;
835  b(1) = twothirds;
836  b(2) = zero;
837  b(3) = onesixth;
838 
839  // Fill c:
840  c(0) = zero;
841  c(1) = onehalf;
842  c(2) = one;
843  c(3) = one;
844 
845  int order = 3;
846 
848  this->getStepperType(), A, b, c, order, order, order));
849  }
850 };
851 
853 // ------------------------------------------------------------------------
854 template <class Scalar>
857  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
859 {
861  stepper->setStepperRKValues(pl);
862 
863  if (model != Teuchos::null) {
864  stepper->setModel(model);
865  stepper->initialize();
866  }
867 
868  return stepper;
869 }
870 
871 // ----------------------------------------------------------------------------
893 template <class Scalar>
895  : virtual public StepperExplicitRK<Scalar> {
896  public:
903  {
904  this->setStepperName("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
905  this->setStepperType("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
906  this->setupTableau();
907  this->setupDefault();
908  this->setUseFSAL(false);
909  this->setICConsistency("None");
910  this->setICConsistencyCheck(false);
911  }
912 
914  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
915  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
916  bool useEmbedded,
917  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
918  {
919  this->setStepperName("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
920  this->setStepperType("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
921  this->setupTableau();
922  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
923  useEmbedded, stepperRKAppAction);
924  }
925 
926  std::string getDescription() const
927  {
928  std::ostringstream Description;
929  Description << this->getStepperType() << "\n"
930  << "Kinnmark & Gray 5 stage, 3rd order scheme \n"
931  << "Modified by P. Ullrich. From the prim_advance_mod.F90 \n"
932  << "routine in the HOMME atmosphere model code.\n"
933  << "c = [ 0 1/5 1/5 1/3 2/3 ]'\n"
934  << "A = [ 0 ]\n"
935  << " [ 1/5 0 ]\n"
936  << " [ 0 1/5 0 ]\n"
937  << " [ 0 0 1/3 0 ]\n"
938  << " [ 0 0 0 2/3 0 ]\n"
939  << "b = [ 1/4 0 0 0 3/4 ]'";
940  return Description.str();
941  }
942 
943  protected:
945  {
947  int NumStages = 5;
948  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
951 
952  const Scalar one = ST::one();
953  const Scalar onefifth = one / (5 * one);
954  const Scalar onefourth = one / (4 * one);
955  const Scalar onethird = one / (3 * one);
956  const Scalar twothirds = 2 * one / (3 * one);
957  const Scalar threefourths = 3 * one / (4 * one);
958  const Scalar zero = ST::zero();
959 
960  // Fill A:
961  A(0, 0) = zero;
962  A(0, 1) = zero;
963  A(0, 2) = zero;
964  A(0, 3) = zero;
965  A(0, 4) = zero;
966  A(1, 0) = onefifth;
967  A(1, 1) = zero;
968  A(1, 2) = zero;
969  A(1, 3) = zero;
970  A(1, 4) = zero;
971  A(2, 0) = zero;
972  A(2, 1) = onefifth;
973  A(2, 2) = zero;
974  A(2, 3) = zero;
975  A(2, 4) = zero;
976  A(3, 0) = zero;
977  A(3, 1) = zero;
978  A(3, 2) = onethird;
979  A(3, 3) = zero;
980  A(3, 4) = zero;
981  A(4, 0) = zero;
982  A(4, 1) = zero;
983  A(4, 2) = zero;
984  A(4, 3) = twothirds;
985  A(4, 4) = zero;
986 
987  // Fill b:
988  b(0) = onefourth;
989  b(1) = zero;
990  b(2) = zero;
991  b(3) = zero;
992  b(4) = threefourths;
993 
994  // Fill c:
995  c(0) = zero;
996  c(1) = onefifth;
997  c(2) = onefifth;
998  c(3) = onethird;
999  c(4) = twothirds;
1000 
1001  int order = 3;
1002 
1004  this->getStepperType(), A, b, c, order, order, order));
1005  }
1006 };
1007 
1009 // ------------------------------------------------------------------------
1010 template <class Scalar>
1013  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1015 {
1017  stepper->setStepperRKValues(pl);
1018 
1019  if (model != Teuchos::null) {
1020  stepper->setModel(model);
1021  stepper->initialize();
1022  }
1023 
1024  return stepper;
1025 }
1026 
1027 // ----------------------------------------------------------------------------
1045 template <class Scalar>
1046 class StepperERK_3Stage3rdOrder : virtual public StepperExplicitRK<Scalar> {
1047  public:
1054  {
1055  this->setStepperName("RK Explicit 3 Stage 3rd order");
1056  this->setStepperType("RK Explicit 3 Stage 3rd order");
1057  this->setupTableau();
1058  this->setupDefault();
1059  this->setUseFSAL(false);
1060  this->setICConsistency("None");
1061  this->setICConsistencyCheck(false);
1062  }
1063 
1065  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1066  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
1067  bool useEmbedded,
1068  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1069  {
1070  this->setStepperName("RK Explicit 3 Stage 3rd order");
1071  this->setStepperType("RK Explicit 3 Stage 3rd order");
1072  this->setupTableau();
1073  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
1074  useEmbedded, stepperRKAppAction);
1075  }
1076 
1077  std::string getDescription() const
1078  {
1079  std::ostringstream Description;
1080  Description << this->getStepperType() << "\n"
1081  << "c = [ 0 1/2 1 ]'\n"
1082  << "A = [ 0 ]\n"
1083  << " [ 1/2 0 ]\n"
1084  << " [ -1 2 0 ]\n"
1085  << "b = [ 1/6 4/6 1/6 ]'";
1086  return Description.str();
1087  }
1088 
1089  protected:
1091  {
1092  typedef Teuchos::ScalarTraits<Scalar> ST;
1093  const Scalar one = ST::one();
1094  const Scalar two = Teuchos::as<Scalar>(2 * one);
1095  const Scalar zero = ST::zero();
1096  const Scalar onehalf = one / (2 * one);
1097  const Scalar onesixth = one / (6 * one);
1098  const Scalar foursixth = 4 * one / (6 * one);
1099 
1100  int NumStages = 3;
1101  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
1104 
1105  // Fill A:
1106  A(0, 0) = zero;
1107  A(0, 1) = zero;
1108  A(0, 2) = zero;
1109  A(1, 0) = onehalf;
1110  A(1, 1) = zero;
1111  A(1, 2) = zero;
1112  A(2, 0) = -one;
1113  A(2, 1) = two;
1114  A(2, 2) = zero;
1115 
1116  // Fill b:
1117  b(0) = onesixth;
1118  b(1) = foursixth;
1119  b(2) = onesixth;
1120 
1121  // fill c:
1122  c(0) = zero;
1123  c(1) = onehalf;
1124  c(2) = one;
1125 
1126  int order = 3;
1127 
1129  this->getStepperType(), A, b, c, order, order, order));
1130  }
1131 };
1132 
1134 // ------------------------------------------------------------------------
1135 template <class Scalar>
1138  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1140 {
1141  auto stepper = Teuchos::rcp(new StepperERK_3Stage3rdOrder<Scalar>());
1142  stepper->setStepperRKValues(pl);
1143 
1144  if (model != Teuchos::null) {
1145  stepper->setModel(model);
1146  stepper->initialize();
1147  }
1148 
1149  return stepper;
1150 }
1151 
1152 // ----------------------------------------------------------------------------
1181 template <class Scalar>
1182 class StepperERK_3Stage3rdOrderTVD : virtual public StepperExplicitRK<Scalar> {
1183  public:
1190  {
1191  this->setStepperName("RK Explicit 3 Stage 3rd order TVD");
1192  this->setStepperType("RK Explicit 3 Stage 3rd order TVD");
1193  this->setupTableau();
1194  this->setupDefault();
1195  this->setUseFSAL(false);
1196  this->setICConsistency("None");
1197  this->setICConsistencyCheck(false);
1198  }
1199 
1201  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1202  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
1203  bool useEmbedded,
1204  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1205  {
1206  this->setStepperName("RK Explicit 3 Stage 3rd order TVD");
1207  this->setStepperType("RK Explicit 3 Stage 3rd order TVD");
1208  this->setupTableau();
1209  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
1210  useEmbedded, stepperRKAppAction);
1211  }
1212 
1213  std::string getDescription() const
1214  {
1215  std::ostringstream Description;
1216  Description << this->getStepperType() << "\n"
1217  << "This Stepper is known as 'RK Explicit 3 Stage 3rd order "
1218  "TVD' or 'SSPERK33'.\n"
1219  << "Sigal Gottlieb and Chi-Wang Shu\n"
1220  << "`Total Variation Diminishing Runge-Kutta Schemes'\n"
1221  << "Mathematics of Computation\n"
1222  << "Volume 67, Number 221, January 1998, pp. 73-85\n"
1223  << "c = [ 0 1 1/2 ]'\n"
1224  << "A = [ 0 ]\n"
1225  << " [ 1 0 ]\n"
1226  << " [ 1/4 1/4 0 ]\n"
1227  << "b = [ 1/6 1/6 4/6 ]'\n"
1228  << "This is also written in the following set of updates.\n"
1229  << "u1 = u^n + dt L(u^n)\n"
1230  << "u2 = 3 u^n/4 + u1/4 + dt L(u1)/4\n"
1231  << "u^(n+1) = u^n/3 + 2 u2/2 + 2 dt L(u2)/3";
1232  return Description.str();
1233  }
1234 
1235  protected:
1237  {
1238  typedef Teuchos::ScalarTraits<Scalar> ST;
1239  using Teuchos::as;
1240  const Scalar one = ST::one();
1241  const Scalar zero = ST::zero();
1242  const Scalar onehalf = one / (2 * one);
1243  const Scalar onefourth = one / (4 * one);
1244  const Scalar onesixth = one / (6 * one);
1245  const Scalar foursixth = 4 * one / (6 * one);
1246 
1247  int NumStages = 3;
1248  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
1251  Teuchos::SerialDenseVector<int, Scalar> bstar(NumStages);
1252 
1253  // Fill A:
1254  A(0, 0) = zero;
1255  A(0, 1) = zero;
1256  A(0, 2) = zero;
1257  A(1, 0) = one;
1258  A(1, 1) = zero;
1259  A(1, 2) = zero;
1260  A(2, 0) = onefourth;
1261  A(2, 1) = onefourth;
1262  A(2, 2) = zero;
1263 
1264  // Fill b:
1265  b(0) = onesixth;
1266  b(1) = onesixth;
1267  b(2) = foursixth;
1268 
1269  // fill c:
1270  c(0) = zero;
1271  c(1) = one;
1272  c(2) = onehalf;
1273 
1274  // Fill bstar:
1275  bstar(0) = as<Scalar>(0.291485418878409);
1276  bstar(1) = as<Scalar>(0.291485418878409);
1277  bstar(2) = as<Scalar>(0.417029162243181);
1278 
1279  int order = 3;
1280 
1282  this->getStepperType(), A, b, c, order, order, order, bstar));
1283  this->tableau_->setTVD(true);
1284  this->tableau_->setTVDCoeff(1.0);
1285  }
1286 };
1287 
1289 // ------------------------------------------------------------------------
1290 template <class Scalar>
1293  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1295 {
1296  auto stepper = Teuchos::rcp(new StepperERK_3Stage3rdOrderTVD<Scalar>());
1297 
1298  // Test for aliases.
1299  if (pl != Teuchos::null) {
1300  auto stepperType =
1301  pl->get<std::string>("Stepper Type", stepper->getStepperType());
1302 
1304  stepperType != stepper->getStepperType() && stepperType != "SSPERK33" &&
1305  stepperType != "SSPRK3",
1306  std::logic_error,
1307  " ParameterList 'Stepper Type' (='"
1308  << stepperType
1309  << "')\n does not match type for this Stepper (='"
1310  << stepper->getStepperType()
1311  << "')\n or one of its aliases ('SSPERK33' or 'SSPRK3').\n");
1312 
1313  // Reset default StepperType.
1314  pl->set<std::string>("Stepper Type", stepper->getStepperType());
1315  }
1316 
1317  stepper->setStepperRKValues(pl);
1318 
1319  if (model != Teuchos::null) {
1320  stepper->setModel(model);
1321  stepper->initialize();
1322  }
1323 
1324  return stepper;
1325 }
1326 
1327 // ----------------------------------------------------------------------------
1349 template <class Scalar>
1350 class StepperERK_3Stage3rdOrderHeun : virtual public StepperExplicitRK<Scalar> {
1351  public:
1358  {
1359  this->setStepperName("RK Explicit 3 Stage 3rd order by Heun");
1360  this->setStepperType("RK Explicit 3 Stage 3rd order by Heun");
1361  this->setupTableau();
1362  this->setupDefault();
1363  this->setUseFSAL(false);
1364  this->setICConsistency("None");
1365  this->setICConsistencyCheck(false);
1366  }
1367 
1369  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1370  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
1371  bool useEmbedded,
1372  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1373  {
1374  this->setStepperName("RK Explicit 3 Stage 3rd order by Heun");
1375  this->setStepperType("RK Explicit 3 Stage 3rd order by Heun");
1376  this->setupTableau();
1377  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
1378  useEmbedded, stepperRKAppAction);
1379  }
1380 
1381  std::string getDescription() const
1382  {
1383  std::ostringstream Description;
1384  Description << this->getStepperType() << "\n"
1385  << "Solving Ordinary Differential Equations I:\n"
1386  << "Nonstiff Problems, 2nd Revised Edition\n"
1387  << "E. Hairer, S.P. Norsett, G. Wanner\n"
1388  << "Table 1.1, pg 135\n"
1389  << "c = [ 0 1/3 2/3 ]'\n"
1390  << "A = [ 0 ] \n"
1391  << " [ 1/3 0 ]\n"
1392  << " [ 0 2/3 0 ]\n"
1393  << "b = [ 1/4 0 3/4 ]'";
1394  return Description.str();
1395  }
1396 
1397  protected:
1399  {
1400  typedef Teuchos::ScalarTraits<Scalar> ST;
1401  const Scalar one = ST::one();
1402  const Scalar zero = ST::zero();
1403  const Scalar onethird = one / (3 * one);
1404  const Scalar twothirds = 2 * one / (3 * one);
1405  const Scalar onefourth = one / (4 * one);
1406  const Scalar threefourths = 3 * one / (4 * one);
1407 
1408  int NumStages = 3;
1409  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
1412 
1413  // Fill A:
1414  A(0, 0) = zero;
1415  A(0, 1) = zero;
1416  A(0, 2) = zero;
1417  A(1, 0) = onethird;
1418  A(1, 1) = zero;
1419  A(1, 2) = zero;
1420  A(2, 0) = zero;
1421  A(2, 1) = twothirds;
1422  A(2, 2) = zero;
1423 
1424  // Fill b:
1425  b(0) = onefourth;
1426  b(1) = zero;
1427  b(2) = threefourths;
1428 
1429  // fill c:
1430  c(0) = zero;
1431  c(1) = onethird;
1432  c(2) = twothirds;
1433 
1434  int order = 3;
1435 
1437  this->getStepperType(), A, b, c, order, order, order));
1438  }
1439 };
1440 
1442 // ------------------------------------------------------------------------
1443 template <class Scalar>
1446  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1448 {
1450  stepper->setStepperRKValues(pl);
1451 
1452  if (model != Teuchos::null) {
1453  stepper->setModel(model);
1454  stepper->initialize();
1455  }
1456 
1457  return stepper;
1458 }
1459 
1460 // ----------------------------------------------------------------------------
1481 template <class Scalar>
1482 class StepperERK_Midpoint : virtual public StepperExplicitRK<Scalar> {
1483  public:
1490  {
1491  this->setStepperName("RK Explicit Midpoint");
1492  this->setStepperType("RK Explicit Midpoint");
1493  this->setupTableau();
1494  this->setupDefault();
1495  this->setUseFSAL(false);
1496  this->setICConsistency("None");
1497  this->setICConsistencyCheck(false);
1498  }
1499 
1501  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1502  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
1503  bool useEmbedded,
1504  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1505  {
1506  this->setStepperName("RK Explicit Midpoint");
1507  this->setStepperType("RK Explicit Midpoint");
1508  this->setupTableau();
1509  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
1510  useEmbedded, stepperRKAppAction);
1511  }
1512 
1513  std::string getDescription() const
1514  {
1515  std::ostringstream Description;
1516  Description << this->getStepperType() << "\n"
1517  << "Solving Ordinary Differential Equations I:\n"
1518  << "Nonstiff Problems, 2nd Revised Edition\n"
1519  << "E. Hairer, S.P. Norsett, G. Wanner\n"
1520  << "Table 1.1, pg 135\n"
1521  << "c = [ 0 1/2 ]'\n"
1522  << "A = [ 0 ]\n"
1523  << " [ 1/2 0 ]\n"
1524  << "b = [ 0 1 ]'";
1525  return Description.str();
1526  }
1527 
1528  protected:
1530  {
1531  typedef Teuchos::ScalarTraits<Scalar> ST;
1532  const Scalar one = ST::one();
1533  const Scalar zero = ST::zero();
1534  const Scalar onehalf = one / (2 * one);
1535 
1536  int NumStages = 2;
1537  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
1540 
1541  // Fill A:
1542  A(0, 0) = zero;
1543  A(0, 1) = zero;
1544  A(1, 0) = onehalf;
1545  A(1, 1) = zero;
1546 
1547  // Fill b:
1548  b(0) = zero;
1549  b(1) = one;
1550 
1551  // fill c:
1552  c(0) = zero;
1553  c(1) = onehalf;
1554 
1555  int order = 2;
1556 
1558  this->getStepperType(), A, b, c, order, order, order));
1559  }
1560 };
1561 
1563 // ------------------------------------------------------------------------
1564 template <class Scalar>
1566  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1568 {
1569  auto stepper = Teuchos::rcp(new StepperERK_Midpoint<Scalar>());
1570  stepper->setStepperRKValues(pl);
1571 
1572  if (model != Teuchos::null) {
1573  stepper->setModel(model);
1574  stepper->initialize();
1575  }
1576 
1577  return stepper;
1578 }
1579 
1580 // ----------------------------------------------------------------------------
1597 template <class Scalar>
1598 class StepperERK_Ralston : virtual public StepperExplicitRK<Scalar> {
1599  public:
1606  {
1607  this->setStepperName("RK Explicit Ralston");
1608  this->setStepperType("RK Explicit Ralston");
1609  this->setupTableau();
1610  this->setupDefault();
1611  this->setUseFSAL(false);
1612  this->setICConsistency("None");
1613  this->setICConsistencyCheck(false);
1614  }
1615 
1617  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1618  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
1619  bool useEmbedded,
1620  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1621  {
1622  this->setStepperName("RK Explicit Ralston");
1623  this->setStepperType("RK Explicit Ralston");
1624  this->setupTableau();
1625  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
1626  useEmbedded, stepperRKAppAction);
1627  }
1628 
1629  std::string getDescription() const
1630  {
1631  std::ostringstream Description;
1632  Description << this->getStepperType() << "\n"
1633  << "This Stepper is known as 'RK Explicit Ralston' or 'RK2'.\n"
1634  << "c = [ 0 2/3 ]'\n"
1635  << "A = [ 0 ]\n"
1636  << " [ 2/3 0 ]\n"
1637  << "b = [ 1/4 3/4 ]'";
1638  return Description.str();
1639  }
1640 
1641  protected:
1643  {
1644  typedef Teuchos::ScalarTraits<Scalar> ST;
1645  const Scalar one = ST::one();
1646  const Scalar zero = ST::zero();
1647 
1648  const int NumStages = 2;
1649  const int order = 2;
1650  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
1653 
1654  // Fill A:
1655  A(0, 0) = zero;
1656  A(0, 1) = zero;
1657  A(1, 1) = zero;
1658  A(1, 0) = (2 * one) / (3 * one);
1659 
1660  // Fill b:
1661  b(0) = (1 * one) / (4 * one);
1662  b(1) = (3 * one) / (4 * one);
1663 
1664  // fill c:
1665  c(0) = zero;
1666  c(1) = (2 * one) / (3 * one);
1667 
1669  this->getStepperType(), A, b, c, order, order, order));
1670  this->tableau_->setTVD(true);
1671  this->tableau_->setTVDCoeff(0.5);
1672  }
1673 };
1674 
1676 // ------------------------------------------------------------------------
1677 template <class Scalar>
1679  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1681 {
1682  auto stepper = Teuchos::rcp(new StepperERK_Ralston<Scalar>());
1683 
1684  // Test for aliases.
1685  if (pl != Teuchos::null) {
1686  auto stepperType =
1687  pl->get<std::string>("Stepper Type", stepper->getStepperType());
1688 
1690  stepperType != stepper->getStepperType() && stepperType != "RK2",
1691  std::logic_error,
1692  " ParameterList 'Stepper Type' (='"
1693  << stepperType
1694  << "')\n does not match type for this Stepper (='"
1695  << stepper->getStepperType()
1696  << "')\n or one of its aliases ('RK2').\n");
1697 
1698  // Reset default StepperType.
1699  pl->set<std::string>("Stepper Type", stepper->getStepperType());
1700  }
1701 
1702  stepper->setStepperRKValues(pl);
1703 
1704  if (model != Teuchos::null) {
1705  stepper->setModel(model);
1706  stepper->initialize();
1707  }
1708 
1709  return stepper;
1710 }
1711 
1712 // ----------------------------------------------------------------------------
1730 template <class Scalar>
1731 class StepperERK_Trapezoidal : virtual public StepperExplicitRK<Scalar> {
1732  public:
1739  {
1740  this->setStepperName("RK Explicit Trapezoidal");
1741  this->setStepperType("RK Explicit Trapezoidal");
1742  this->setupTableau();
1743  this->setupDefault();
1744  this->setUseFSAL(false);
1745  this->setICConsistency("None");
1746  this->setICConsistencyCheck(false);
1747  }
1748 
1750  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1751  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
1752  bool useEmbedded,
1753  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1754  {
1755  this->setStepperName("RK Explicit Trapezoidal");
1756  this->setStepperType("RK Explicit Trapezoidal");
1757  this->setupTableau();
1758  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
1759  useEmbedded, stepperRKAppAction);
1760  }
1761 
1762  std::string getDescription() const
1763  {
1764  std::ostringstream Description;
1765  Description << this->getStepperType() << "\n"
1766  << "This Stepper is known as 'RK Explicit Trapezoidal' or "
1767  "'Heuns Method' or 'SSPERK22'.\n"
1768  << "c = [ 0 1 ]'\n"
1769  << "A = [ 0 ]\n"
1770  << " [ 1 0 ]\n"
1771  << "b = [ 1/2 1/2 ]\n"
1772  << "bstart = [ 3/4 1/4 ]'";
1773  return Description.str();
1774  }
1775 
1776  protected:
1778  {
1779  typedef Teuchos::ScalarTraits<Scalar> ST;
1780  using Teuchos::as;
1781  const Scalar one = ST::one();
1782  const Scalar zero = ST::zero();
1783  const Scalar onehalf = one / (2 * one);
1784 
1785  int NumStages = 2;
1786  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
1789  Teuchos::SerialDenseVector<int, Scalar> bstar(NumStages);
1790 
1791  // Fill A:
1792  A(0, 0) = zero;
1793  A(0, 1) = zero;
1794  A(1, 0) = one;
1795  A(1, 1) = zero;
1796 
1797  // Fill b:
1798  b(0) = onehalf;
1799  b(1) = onehalf;
1800 
1801  // fill c:
1802  c(0) = zero;
1803  c(1) = one;
1804 
1805  // Fill bstar
1806  bstar(0) = as<Scalar>(3 * one / (4 * one));
1807  bstar(1) = as<Scalar>(1 * one / (4 * one));
1808 
1809  int order = 2;
1810 
1812  this->getStepperType(), A, b, c, order, order, order, bstar));
1813  this->tableau_->setTVD(true);
1814  this->tableau_->setTVDCoeff(1.0);
1815  }
1816 };
1817 
1819 // ------------------------------------------------------------------------
1820 template <class Scalar>
1822  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1824 {
1825  auto stepper = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
1826 
1827  // Test for aliases.
1828  if (pl != Teuchos::null) {
1829  auto stepperType =
1830  pl->get<std::string>("Stepper Type", stepper->getStepperType());
1831 
1833  stepperType != stepper->getStepperType() &&
1834  stepperType != "Heuns Method" && stepperType != "SSPERK22" &&
1835  stepperType != "SSPRK2",
1836  std::logic_error,
1837  " ParameterList 'Stepper Type' (='"
1838  << stepperType
1839  << "')\n does not match type for this Stepper (='"
1840  << stepper->getStepperType()
1841  << "')\n or one of its aliases ('Heuns Method' or 'SSPERK22' or "
1842  << "'SSPRK2').\n");
1843 
1844  // Reset default StepperType.
1845  pl->set<std::string>("Stepper Type", stepper->getStepperType());
1846  }
1847 
1848  stepper->setStepperRKValues(pl);
1849 
1850  if (model != Teuchos::null) {
1851  stepper->setModel(model);
1852  stepper->initialize();
1853  }
1854 
1855  return stepper;
1856 }
1857 
1858 // ----------------------------------------------------------------------------
1875 template <class Scalar>
1876 class StepperERK_SSPERK54 : virtual public StepperExplicitRK<Scalar> {
1877  public:
1879  {
1880  this->setStepperName("SSPERK54");
1881  this->setStepperType("SSPERK54");
1882  this->setupTableau();
1883  this->setupDefault();
1884  this->setUseFSAL(false);
1885  this->setICConsistency("None");
1886  this->setICConsistencyCheck(false);
1887  }
1888 
1890  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1891  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
1892  bool useEmbedded,
1893  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1894  {
1895  this->setStepperName("SSPERK54");
1896  this->setStepperType("SSPERK54");
1897  this->setupTableau();
1898  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
1899  useEmbedded, stepperRKAppAction);
1900  }
1901 
1902  std::string getDescription() const
1903  {
1904  std::ostringstream Description;
1905  Description
1906  << this->getStepperType() << "\n"
1907  << "Strong Stability Preserving Explicit RK (stage=5, order=4)\n"
1908  << std::endl;
1909  return Description.str();
1910  }
1911 
1912  protected:
1914  {
1915  typedef Teuchos::ScalarTraits<Scalar> ST;
1916  using Teuchos::as;
1917  const int NumStages = 5;
1918  const int order = 4;
1919  const Scalar sspcoef = 1.5082;
1920  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
1923  Teuchos::SerialDenseVector<int, Scalar> bstar(NumStages);
1924  const Scalar zero = ST::zero();
1925 
1926  // Fill A:
1927  A(0, 0) = A(0, 1) = A(0, 2) = A(0, 3) = A(0, 4) = zero;
1928 
1929  A(1, 0) = as<Scalar>(0.391752226571889);
1930  A(1, 1) = A(1, 2) = A(1, 3) = A(0, 4) = zero;
1931 
1932  A(2, 0) = as<Scalar>(0.217669096261169);
1933  A(2, 1) = as<Scalar>(0.368410593050372);
1934  A(2, 2) = A(2, 3) = A(2, 4) = zero;
1935 
1936  A(3, 0) = as<Scalar>(0.082692086657811);
1937  A(3, 1) = as<Scalar>(0.139958502191896);
1938  A(3, 2) = as<Scalar>(0.251891774271693);
1939  A(3, 3) = A(2, 4) = zero;
1940 
1941  A(4, 0) = as<Scalar>(0.067966283637115);
1942  A(4, 1) = as<Scalar>(0.115034698504632);
1943  A(4, 2) = as<Scalar>(0.207034898597385);
1944  A(4, 3) = as<Scalar>(0.544974750228520);
1945  A(4, 4) = zero;
1946 
1947  // Fill b:
1948  b(0) = as<Scalar>(0.146811876084786);
1949  b(1) = as<Scalar>(0.248482909444976);
1950  b(2) = as<Scalar>(0.104258830331980);
1951  b(3) = as<Scalar>(0.274438900901350);
1952  b(4) = as<Scalar>(0.226007483236908);
1953 
1954  // fill c:
1955  c(0) = zero;
1956  c(1) = A(1, 0);
1957  c(2) = A(2, 0) + A(2, 1);
1958  c(3) = A(3, 0) + A(3, 1) + A(3, 2);
1959  c(4) = A(4, 0) + A(4, 1) + A(4, 2) + A(4, 3);
1960 
1961  // Fill bstar:
1962  bstar(0) = as<Scalar>(0.130649104813131);
1963  bstar(1) = as<Scalar>(0.317716031201302);
1964  bstar(2) = as<Scalar>(0.000000869337261);
1965  bstar(3) = as<Scalar>(0.304581512634772);
1966  bstar(4) = as<Scalar>(0.247052482013534);
1967 
1969  this->getStepperType(), A, b, c, order, order, order, bstar));
1970  this->tableau_->setTVD(true);
1971  this->tableau_->setTVDCoeff(sspcoef);
1972  }
1973 };
1974 
1976 // ------------------------------------------------------------------------
1977 template <class Scalar>
1979  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1981 {
1982  auto stepper = Teuchos::rcp(new StepperERK_SSPERK54<Scalar>());
1983  stepper->setStepperRKValues(pl);
1984 
1985  if (model != Teuchos::null) {
1986  stepper->setModel(model);
1987  stepper->initialize();
1988  }
1989 
1990  return stepper;
1991 }
1992 
1993 // ----------------------------------------------------------------------------
2023 template <class Scalar>
2024 class StepperERK_General : virtual public StepperExplicitRK<Scalar> {
2025  public:
2027  {
2028  this->setStepperName("General ERK");
2029  this->setStepperType("General ERK");
2030  this->setupTableau();
2031  this->setupDefault();
2032  this->setUseFSAL(false);
2033  this->setICConsistency("None");
2034  this->setICConsistencyCheck(false);
2035  }
2036 
2038  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2039  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
2040  bool useEmbedded, const Teuchos::SerialDenseMatrix<int, Scalar>& A,
2042  const Teuchos::SerialDenseVector<int, Scalar>& c, const int order,
2043  const int orderMin, const int orderMax,
2045  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction =
2046  Teuchos::null)
2047  {
2048  this->setStepperName("General ERK");
2049  this->setStepperType("General ERK");
2050  this->setTableau(A, b, c, order, orderMin, orderMax, bstar);
2051 
2053  this->tableau_->isImplicit() == true, std::logic_error,
2054  "Error - General ERK received an implicit Butcher Tableau!\n");
2055 
2056  this->setup(appModel, useFSAL, ICConsistency, ICConsistencyCheck,
2057  useEmbedded, stepperRKAppAction);
2058  }
2059 
2060  virtual std::string getDescription() const
2061  {
2062  std::stringstream Description;
2063  Description
2064  << this->getStepperType() << "\n"
2065  << "The format of the Butcher Tableau parameter list is\n"
2066  << " <Parameter name=\"A\" type=\"string\" value=\"# # # ;\n"
2067  << " # # # ;\n"
2068  << " # # #\"/>\n"
2069  << " <Parameter name=\"b\" type=\"string\" value=\"# # #\"/>\n"
2070  << " <Parameter name=\"c\" type=\"string\" value=\"# # #\"/>\n\n"
2071  << "Note the number of stages is implicit in the number of entries.\n"
2072  << "The number of stages must be consistent.\n"
2073  << "\n"
2074  << "Default tableau is RK4 (order=4):\n"
2075  << "c = [ 0 1/2 1/2 1 ]'\n"
2076  << "A = [ 0 ]\n"
2077  << " [ 1/2 0 ]\n"
2078  << " [ 0 1/2 0 ]\n"
2079  << " [ 0 0 1 0 ]\n"
2080  << "b = [ 1/6 1/3 1/3 1/6 ]'";
2081  return Description.str();
2082  }
2083 
2085  {
2086  if (this->tableau_ == Teuchos::null) {
2087  // Set tableau to the default if null, otherwise keep current tableau.
2088  auto stepper = Teuchos::rcp(new StepperERK_4Stage4thOrder<Scalar>());
2089  auto t = stepper->getTableau();
2091  this->getStepperType(), t->A(), t->b(), t->c(), t->order(),
2092  t->orderMin(), t->orderMax(), t->bstar()));
2093  }
2094  }
2095 
2099  const int order, const int orderMin, const int orderMax,
2102  {
2104  this->getStepperType(), A, b, c, order, orderMin, orderMax, bstar));
2105  this->isInitialized_ = false;
2106  }
2107 
2108  virtual std::string getDefaultICConsistency() const { return "Consistent"; }
2109 
2111  {
2112  auto pl = this->getValidParametersBasicERK();
2113 
2114  // Tableau ParameterList
2118  Teuchos::SerialDenseVector<int, Scalar> bstar = this->tableau_->bstar();
2119 
2120  Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
2121 
2122  std::ostringstream Astream;
2123  Astream.precision(15);
2124  for (int i = 0; i < A.numRows(); i++) {
2125  for (int j = 0; j < A.numCols() - 1; j++) {
2126  Astream << A(i, j) << " ";
2127  }
2128  Astream << A(i, A.numCols() - 1);
2129  if (i != A.numRows() - 1) Astream << "; ";
2130  }
2131  tableauPL->set<std::string>("A", Astream.str());
2132 
2133  std::ostringstream bstream;
2134  bstream.precision(15);
2135  for (int i = 0; i < b.length() - 1; i++) {
2136  bstream << b(i) << " ";
2137  }
2138  bstream << b(b.length() - 1);
2139  tableauPL->set<std::string>("b", bstream.str());
2140 
2141  std::ostringstream cstream;
2142  cstream.precision(15);
2143  for (int i = 0; i < c.length() - 1; i++) {
2144  cstream << c(i) << " ";
2145  }
2146  cstream << c(c.length() - 1);
2147  tableauPL->set<std::string>("c", cstream.str());
2148 
2149  tableauPL->set<int>("order", this->tableau_->order());
2150 
2151  if (bstar.length() == 0) {
2152  tableauPL->set("bstar", "");
2153  }
2154  else {
2155  std::ostringstream bstarstream;
2156  bstarstream.precision(15);
2157  for (int i = 0; i < bstar.length() - 1; i++) {
2158  bstarstream << bstar(i) << " ";
2159  }
2160  bstarstream << bstar(bstar.length() - 1);
2161  tableauPL->set<std::string>("bstar", bstarstream.str());
2162  }
2163 
2164  pl->set("Tableau", *tableauPL);
2165 
2166  return pl;
2167  }
2168 };
2169 
2171 // ------------------------------------------------------------------------
2172 template <class Scalar>
2174  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2176 {
2177  auto stepper = Teuchos::rcp(new StepperERK_General<Scalar>());
2178  stepper->setStepperRKValues(pl);
2179 
2180  if (pl != Teuchos::null) {
2181  if (pl->isParameter("Tableau")) {
2182  auto t = stepper->createTableau(pl);
2183  stepper->setTableau(t->A(), t->b(), t->c(), t->order(), t->orderMin(),
2184  t->orderMax(), t->bstar());
2185  }
2186  }
2188  stepper->getTableau()->isImplicit() == true, std::logic_error,
2189  "Error - General ERK received an implicit Butcher Tableau!\n");
2190 
2191  if (model != Teuchos::null) {
2192  stepper->setModel(model);
2193  stepper->initialize();
2194  }
2195 
2196  return stepper;
2197 }
2198 
2199 // ----------------------------------------------------------------------------
2215 template <class Scalar>
2216 class StepperDIRK_BackwardEuler : virtual public StepperDIRK<Scalar> {
2217  public:
2224  {
2225  this->setStepperName("RK Backward Euler");
2226  this->setStepperType("RK Backward Euler");
2227  this->setupTableau();
2228  this->setupDefault();
2229  this->setUseFSAL(false);
2230  this->setICConsistency("None");
2231  this->setICConsistencyCheck(false);
2232  }
2233 
2235  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2237  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
2238  bool useEmbedded, bool zeroInitialGuess,
2239  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
2240  {
2241  this->setStepperName("RK Backward Euler");
2242  this->setStepperType("RK Backward Euler");
2243  this->setupTableau();
2244  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
2245  useEmbedded, zeroInitialGuess, stepperRKAppAction);
2246  }
2247 
2248  std::string getDescription() const
2249  {
2250  std::ostringstream Description;
2251  Description << this->getStepperType() << "\n"
2252  << "c = [ 1 ]'\n"
2253  << "A = [ 1 ]\n"
2254  << "b = [ 1 ]'";
2255  return Description.str();
2256  }
2257 
2258  protected:
2260  {
2261  typedef Teuchos::ScalarTraits<Scalar> ST;
2262  const Scalar sspcoef = std::numeric_limits<Scalar>::max();
2263  int NumStages = 1;
2264  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
2267 
2268  // Fill A:
2269  A(0, 0) = ST::one();
2270 
2271  // Fill b:
2272  b(0) = ST::one();
2273 
2274  // Fill c:
2275  c(0) = ST::one();
2276 
2277  int order = 1;
2278 
2280  this->getStepperType(), A, b, c, order, order, order));
2281  this->tableau_->setTVD(true);
2282  this->tableau_->setTVDCoeff(sspcoef);
2283  }
2284 };
2285 
2287 // ------------------------------------------------------------------------
2288 template <class Scalar>
2291  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2293 {
2294  auto stepper = Teuchos::rcp(new StepperDIRK_BackwardEuler<Scalar>());
2295  stepper->setStepperDIRKValues(pl);
2296 
2297  if (model != Teuchos::null) {
2298  stepper->setModel(model);
2299  stepper->initialize();
2300  }
2301 
2302  return stepper;
2303 }
2304 
2305 // ----------------------------------------------------------------------------
2330 template <class Scalar>
2331 class StepperSDIRK_2Stage2ndOrder : virtual public StepperDIRK<Scalar> {
2332  public:
2339  {
2340  typedef Teuchos::ScalarTraits<Scalar> ST;
2341  const Scalar one = ST::one();
2342  gammaDefault_ =
2343  Teuchos::as<Scalar>((2 * one - ST::squareroot(2 * one)) / (2 * one));
2344 
2345  this->setStepperName("SDIRK 2 Stage 2nd order");
2346  this->setStepperType("SDIRK 2 Stage 2nd order");
2347  this->setGamma(gammaDefault_);
2348  this->setupTableau();
2349  this->setupDefault();
2350  this->setUseFSAL(false);
2351  this->setICConsistency("None");
2352  this->setICConsistencyCheck(false);
2353  }
2354 
2356  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2358  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
2359  bool useEmbedded, bool zeroInitialGuess,
2360  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
2361  Scalar gamma = Scalar(0.2928932188134524))
2362  {
2363  typedef Teuchos::ScalarTraits<Scalar> ST;
2364  const Scalar one = ST::one();
2365  gammaDefault_ =
2366  Teuchos::as<Scalar>((2 * one - ST::squareroot(2 * one)) / (2 * one));
2367 
2368  this->setStepperName("SDIRK 2 Stage 2nd order");
2369  this->setStepperType("SDIRK 2 Stage 2nd order");
2370  this->setGamma(gamma);
2371  this->setupTableau();
2372  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
2373  useEmbedded, zeroInitialGuess, stepperRKAppAction);
2374  }
2375 
2376  void setGamma(Scalar gamma)
2377  {
2378  gamma_ = gamma;
2379  this->isInitialized_ = false;
2380  this->setupTableau();
2381  }
2382 
2383  Scalar getGamma() const { return gamma_; }
2384 
2385  std::string getDescription() const
2386  {
2387  std::ostringstream Description;
2388  Description << this->getStepperType() << "\n"
2389  << "Computer Methods for ODEs and DAEs\n"
2390  << "U. M. Ascher and L. R. Petzold\n"
2391  << "p. 106\n"
2392  << "gamma = (2+-sqrt(2))/2\n"
2393  << "c = [ gamma 1 ]'\n"
2394  << "A = [ gamma 0 ]\n"
2395  << " [ 1-gamma gamma ]\n"
2396  << "b = [ 1-gamma gamma ]'";
2397  return Description.str();
2398  }
2399 
2401  {
2402  auto pl = this->getValidParametersBasicDIRK();
2403  pl->template set<double>(
2404  "gamma", this->getGamma(),
2405  "The default value is gamma = (2-sqrt(2))/2. "
2406  "This will produce an L-stable 2nd order method with the stage "
2407  "times within the timestep. Other values of gamma will still "
2408  "produce an L-stable scheme, but will only be 1st order accurate.");
2409 
2410  return pl;
2411  }
2412 
2413  protected:
2415  {
2416  typedef Teuchos::ScalarTraits<Scalar> ST;
2417  int NumStages = 2;
2418  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
2421 
2422  const Scalar one = ST::one();
2423  const Scalar zero = ST::zero();
2424 
2425  // Fill A:
2426  A(0, 0) = gamma_;
2427  A(0, 1) = zero;
2428  A(1, 0) = Teuchos::as<Scalar>(one - gamma_);
2429  A(1, 1) = gamma_;
2430 
2431  // Fill b:
2432  b(0) = Teuchos::as<Scalar>(one - gamma_);
2433  b(1) = gamma_;
2434 
2435  // Fill c:
2436  c(0) = gamma_;
2437  c(1) = one;
2438 
2439  int order = 1;
2440  if (std::abs((gamma_ - gammaDefault_) / gamma_) < 1.0e-08) order = 2;
2441 
2443  this->getStepperType(), A, b, c, order, order, order));
2444  }
2445 
2446  private:
2448  Scalar gamma_;
2449 };
2450 
2452 // ------------------------------------------------------------------------
2453 template <class Scalar>
2456  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2458 {
2459  auto stepper = Teuchos::rcp(new StepperSDIRK_2Stage2ndOrder<Scalar>());
2460  stepper->setStepperDIRKValues(pl);
2461 
2462  if (pl != Teuchos::null)
2463  stepper->setGamma(pl->get<double>("gamma", 0.2928932188134524));
2464 
2465  if (model != Teuchos::null) {
2466  stepper->setModel(model);
2467  stepper->initialize();
2468  }
2469 
2470  return stepper;
2471 }
2472 
2473 // ----------------------------------------------------------------------------
2501 template <class Scalar>
2502 class StepperSDIRK_3Stage2ndOrder : virtual public StepperDIRK<Scalar> {
2503  public:
2510  {
2511  this->setStepperName("SDIRK 3 Stage 2nd order");
2512  this->setStepperType("SDIRK 3 Stage 2nd order");
2513  this->setupTableau();
2514  this->setupDefault();
2515  this->setUseFSAL(false);
2516  this->setICConsistency("None");
2517  this->setICConsistencyCheck(false);
2518  }
2519 
2521  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2523  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
2524  bool useEmbedded, bool zeroInitialGuess,
2525  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
2526  {
2527  this->setStepperName("SDIRK 3 Stage 2nd order");
2528  this->setStepperType("SDIRK 3 Stage 2nd order");
2529  this->setupTableau();
2530  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
2531  useEmbedded, zeroInitialGuess, stepperRKAppAction);
2532  }
2533 
2534  std::string getDescription() const
2535  {
2536  std::ostringstream Description;
2537  Description << this->getStepperType() << "\n"
2538  << "Implicit-explicit Runge-Kutta schemes and applications to\n"
2539  << "hyperbolic systems with relaxation\n"
2540  << "L Pareschi, G Russo\n"
2541  << "Journal of Scientific computing, 2005 - Springer\n"
2542  << "Table 5\n"
2543  << "gamma = 1/(2+sqrt(2))\n"
2544  << "c = [ gamma (1-gamma) 1/2 ]'\n"
2545  << "A = [ gamma 0 0 ]\n"
2546  << " [ 1-2gamma gamma 0 ]\n"
2547  << " [ 1/2-gamma 0 gamma ]\n"
2548  << "b = [ 1/6 1/6 2/3 ]'";
2549  return Description.str();
2550  }
2551 
2552  protected:
2554  {
2555  typedef Teuchos::ScalarTraits<Scalar> ST;
2556  using Teuchos::as;
2557  const int NumStages = 3;
2558  const int order = 2;
2559  const Scalar sspcoef = 1.0529;
2560  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
2563  const Scalar one = ST::one();
2564  const Scalar zero = ST::zero();
2565  const Scalar gamma = as<Scalar>(one - (one / ST::squareroot(2 * one)));
2566 
2567  // Fill A:
2568  A(0, 0) = A(1, 1) = A(2, 2) = gamma;
2569  A(0, 1) = A(0, 2) = A(1, 2) = A(2, 1) = zero;
2570  A(1, 0) = as<Scalar>(one - 2 * gamma);
2571  A(2, 0) = as<Scalar>((one / (2. * one)) - gamma);
2572 
2573  // Fill b:
2574  b(0) = b(1) = (one / (6 * one));
2575  b(2) = (2 * one) / (3 * one);
2576 
2577  // Fill c:
2578  c(0) = gamma;
2579  c(1) = one - gamma;
2580  c(2) = one / (2 * one);
2581 
2583  this->getStepperType(), A, b, c, order, order, order));
2584  this->tableau_->setTVD(true);
2585  this->tableau_->setTVDCoeff(sspcoef);
2586  }
2587 };
2588 
2590 // ------------------------------------------------------------------------
2591 template <class Scalar>
2594  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2596 {
2597  auto stepper = Teuchos::rcp(new StepperSDIRK_3Stage2ndOrder<Scalar>());
2598  stepper->setStepperDIRKValues(pl);
2599 
2600  if (model != Teuchos::null) {
2601  stepper->setModel(model);
2602  stepper->initialize();
2603  }
2604 
2605  return stepper;
2606 }
2607 
2608 // ----------------------------------------------------------------------------
2637 template <class Scalar>
2638 class StepperSDIRK_2Stage3rdOrder : virtual public StepperDIRK<Scalar> {
2639  public:
2646  {
2647  typedef Teuchos::ScalarTraits<Scalar> ST;
2648  using Teuchos::as;
2649  const Scalar one = ST::one();
2650  gammaDefault_ = as<Scalar>((3 * one + ST::squareroot(3 * one)) / (6 * one));
2651  gammaTypeDefault_ = "3rd Order A-stable";
2652 
2653  this->setStepperName("SDIRK 2 Stage 3rd order");
2654  this->setStepperType("SDIRK 2 Stage 3rd order");
2655  this->setGammaType(gammaTypeDefault_);
2656  this->setGamma(gammaDefault_);
2657  this->setupTableau();
2658  this->setupDefault();
2659  this->setUseFSAL(false);
2660  this->setICConsistency("None");
2661  this->setICConsistencyCheck(false);
2662  }
2663 
2665  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2667  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
2668  bool useEmbedded, bool zeroInitialGuess,
2669  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
2670  std::string gammaType = "3rd Order A-stable",
2671  Scalar gamma = Scalar(0.7886751345948128))
2672  {
2673  typedef Teuchos::ScalarTraits<Scalar> ST;
2674  using Teuchos::as;
2675  const Scalar one = ST::one();
2676  gammaDefault_ = as<Scalar>((3 * one + ST::squareroot(3 * one)) / (6 * one));
2677  gammaTypeDefault_ = "3rd Order A-stable";
2678 
2679  this->setStepperName("SDIRK 2 Stage 3rd order");
2680  this->setStepperType("SDIRK 2 Stage 3rd order");
2681  this->setGammaType(gammaType);
2682  this->setGamma(gamma);
2683  this->setupTableau();
2684  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
2685  useEmbedded, zeroInitialGuess, stepperRKAppAction);
2686  }
2687 
2688  void setGammaType(std::string gammaType)
2689  {
2691  !(gammaType == "3rd Order A-stable" ||
2692  gammaType == "2nd Order L-stable" || gammaType == "gamma"),
2693  std::logic_error,
2694  "gammaType needs to be '3rd Order A-stable', '2nd Order L-stable' "
2695  "or 'gamma'.");
2696 
2697  gammaType_ = gammaType;
2698  this->isInitialized_ = false;
2699  this->setupTableau();
2700  }
2701 
2702  std::string getGammaType() const { return gammaType_; }
2703 
2704  void setGamma(Scalar gamma)
2705  {
2706  if (gammaType_ == "gamma") {
2707  gamma_ = gamma;
2708  this->setupTableau();
2709  }
2710  this->isInitialized_ = false;
2711  }
2712 
2713  Scalar getGamma() const { return gamma_; }
2714 
2715  std::string getDescription() const
2716  {
2717  std::ostringstream Description;
2718  Description << this->getStepperType() << "\n"
2719  << "Solving Ordinary Differential Equations I:\n"
2720  << "Nonstiff Problems, 2nd Revised Edition\n"
2721  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
2722  << "Table 7.2, pg 207\n"
2723  << "gamma = (3+sqrt(3))/6 -> 3rd order and A-stable\n"
2724  << "gamma = (2-sqrt(2))/2 -> 2nd order and L-stable\n"
2725  << "c = [ gamma 1-gamma ]'\n"
2726  << "A = [ gamma 0 ]\n"
2727  << " [ 1-2*gamma gamma ]\n"
2728  << "b = [ 1/2 1/2 ]'";
2729  return Description.str();
2730  }
2731 
2733  {
2734  auto pl = this->getValidParametersBasicDIRK();
2735 
2736  pl->template set<std::string>(
2737  "Gamma Type", this->getGammaType(),
2738  "Valid values are '3rd Order A-stable' ((3+sqrt(3))/6.), "
2739  "'2nd Order L-stable' ((2-sqrt(2))/2) and 'gamma' (user defined). "
2740  "The default value is '3rd Order A-stable'.");
2741  pl->template set<double>(
2742  "gamma", this->getGamma(),
2743  "Equal to (3+sqrt(3))/6 if 'Gamma Type' = '3rd Order A-stable', or "
2744  "(2-sqrt(2))/2 if 'Gamma Type' = '2nd Order L-stable', or "
2745  "user-defined gamma value if 'Gamma Type = 'gamma'. "
2746  "The default value is gamma = (3+sqrt(3))/6, which matches "
2747  "the default 'Gamma Type' = '3rd Order A-stable'.");
2748 
2749  return pl;
2750  }
2751 
2752  protected:
2754  {
2755  typedef Teuchos::ScalarTraits<Scalar> ST;
2756  using Teuchos::as;
2757  int NumStages = 2;
2758  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
2761  const Scalar one = ST::one();
2762  const Scalar zero = ST::zero();
2763 
2764  int order = 0;
2765  if (gammaType_ == "3rd Order A-stable") {
2766  order = 3;
2768  }
2769  else if (gammaType_ == "2nd Order L-stable") {
2770  order = 2;
2771  gamma_ = as<Scalar>((2 * one - ST::squareroot(2 * one)) / (2 * one));
2772  }
2773  else if (gammaType_ == "gamma") {
2774  order = 2;
2775  }
2776 
2777  // Fill A:
2778  A(0, 0) = gamma_;
2779  A(0, 1) = zero;
2780  A(1, 0) = as<Scalar>(one - 2 * gamma_);
2781  A(1, 1) = gamma_;
2782 
2783  // Fill b:
2784  b(0) = as<Scalar>(one / (2 * one));
2785  b(1) = as<Scalar>(one / (2 * one));
2786 
2787  // Fill c:
2788  c(0) = gamma_;
2789  c(1) = as<Scalar>(one - gamma_);
2790 
2792  this->getStepperType(), A, b, c, order, 2, 3));
2793  }
2794 
2795  private:
2796  std::string gammaTypeDefault_;
2797  std::string gammaType_;
2799  Scalar gamma_;
2800 };
2801 
2803 // ------------------------------------------------------------------------
2804 template <class Scalar>
2807  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2809 {
2810  auto stepper = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
2811  stepper->setStepperDIRKValues(pl);
2812 
2813  if (pl != Teuchos::null) {
2814  stepper->setGammaType(
2815  pl->get<std::string>("Gamma Type", "3rd Order A-stable"));
2816  stepper->setGamma(pl->get<double>("gamma", 0.7886751345948128));
2817  }
2818 
2819  if (model != Teuchos::null) {
2820  stepper->setModel(model);
2821  stepper->initialize();
2822  }
2823 
2824  return stepper;
2825 }
2826 
2827 // ----------------------------------------------------------------------------
2848 template <class Scalar>
2849 class StepperEDIRK_2Stage3rdOrder : virtual public StepperDIRK<Scalar> {
2850  public:
2857  {
2858  this->setStepperName("EDIRK 2 Stage 3rd order");
2859  this->setStepperType("EDIRK 2 Stage 3rd order");
2860  this->setupTableau();
2861  this->setupDefault();
2862  this->setUseFSAL(false);
2863  this->setICConsistency("None");
2864  this->setICConsistencyCheck(false);
2865  }
2866 
2868  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2870  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
2871  bool useEmbedded, bool zeroInitialGuess,
2872  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
2873  {
2874  this->setStepperName("EDIRK 2 Stage 3rd order");
2875  this->setStepperType("EDIRK 2 Stage 3rd order");
2876  this->setupTableau();
2877  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
2878  useEmbedded, zeroInitialGuess, stepperRKAppAction);
2879  }
2880 
2881  std::string getDescription() const
2882  {
2883  std::ostringstream Description;
2884  Description << this->getStepperType() << "\n"
2885  << "Hammer & Hollingsworth method\n"
2886  << "Solving Ordinary Differential Equations I:\n"
2887  << "Nonstiff Problems, 2nd Revised Edition\n"
2888  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
2889  << "Table 7.1, pg 205\n"
2890  << "c = [ 0 2/3 ]'\n"
2891  << "A = [ 0 0 ]\n"
2892  << " [ 1/3 1/3 ]\n"
2893  << "b = [ 1/4 3/4 ]'";
2894  return Description.str();
2895  }
2896 
2897  protected:
2899  {
2900  typedef Teuchos::ScalarTraits<Scalar> ST;
2901  using Teuchos::as;
2902  int NumStages = 2;
2903  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
2906  const Scalar one = ST::one();
2907  const Scalar zero = ST::zero();
2908 
2909  // Fill A:
2910  A(0, 0) = zero;
2911  A(0, 1) = zero;
2912  A(1, 0) = as<Scalar>(one / (3 * one));
2913  A(1, 1) = as<Scalar>(one / (3 * one));
2914 
2915  // Fill b:
2916  b(0) = as<Scalar>(one / (4 * one));
2917  b(1) = as<Scalar>(3 * one / (4 * one));
2918 
2919  // Fill c:
2920  c(0) = zero;
2921  c(1) = as<Scalar>(2 * one / (3 * one));
2922  int order = 3;
2923 
2925  this->getStepperType(), A, b, c, order, order, order));
2926  this->tableau_->setTVD(true);
2927  this->tableau_->setTVDCoeff(1.5);
2928  }
2929 };
2930 
2932 template <class Scalar>
2935  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2937 {
2938  auto stepper = Teuchos::rcp(new StepperEDIRK_2Stage3rdOrder<Scalar>());
2939  stepper->setStepperDIRKValues(pl);
2940 
2941  if (model != Teuchos::null) {
2942  stepper->setModel(model);
2943  stepper->initialize();
2944  }
2945 
2946  return stepper;
2947 }
2948 
2949 // ----------------------------------------------------------------------------
2977 template <class Scalar>
2978 class StepperDIRK_1StageTheta : virtual public StepperDIRK<Scalar> {
2979  public:
2986  {
2987  typedef Teuchos::ScalarTraits<Scalar> ST;
2988  thetaDefault_ = ST::one() / (2 * ST::one());
2989 
2990  this->setStepperName("DIRK 1 Stage Theta Method");
2991  this->setStepperType("DIRK 1 Stage Theta Method");
2992  this->setTheta(thetaDefault_);
2993  this->setupTableau();
2994  this->setupDefault();
2995  this->setUseFSAL(false);
2996  this->setICConsistency("None");
2997  this->setICConsistencyCheck(false);
2998  }
2999 
3001  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3003  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
3004  bool useEmbedded, bool zeroInitialGuess,
3005  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
3006  Scalar theta = Scalar(0.5))
3007  {
3008  typedef Teuchos::ScalarTraits<Scalar> ST;
3009  thetaDefault_ = ST::one() / (2 * ST::one());
3010 
3011  this->setStepperName("DIRK 1 Stage Theta Method");
3012  this->setStepperType("DIRK 1 Stage Theta Method");
3013  this->setTheta(theta);
3014  this->setupTableau();
3015  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3016  useEmbedded, zeroInitialGuess, stepperRKAppAction);
3017  }
3018 
3019  void setTheta(Scalar theta)
3020  {
3022  theta == Teuchos::ScalarTraits<Scalar>::zero(), std::logic_error,
3023  "'theta' can not be zero, as it makes this stepper explicit. \n"
3024  "Try using the 'RK Forward Euler' stepper.\n");
3025  theta_ = theta;
3026  this->setupTableau();
3027  this->isInitialized_ = false;
3028  }
3029 
3030  Scalar getTheta() const { return theta_; }
3031 
3032  std::string getDescription() const
3033  {
3034  std::ostringstream Description;
3035  Description << this->getStepperType() << "\n"
3036  << "Non-standard finite-difference methods\n"
3037  << "in dynamical systems, P. Kama,\n"
3038  << "Dissertation, University of Pretoria, pg. 49.\n"
3039  << "Comment: Generalized Implicit Midpoint Method\n"
3040  << "c = [ theta ]'\n"
3041  << "A = [ theta ]\n"
3042  << "b = [ 1 ]'";
3043  return Description.str();
3044  }
3045 
3047  {
3048  auto pl = this->getValidParametersBasicDIRK();
3049 
3050  pl->template set<double>(
3051  "theta", getTheta(),
3052  "Valid values are 0 <= theta <= 1, where theta = 0 "
3053  "implies Forward Euler, theta = 1/2 implies implicit midpoint "
3054  "method (default), and theta = 1 implies Backward Euler. "
3055  "For theta != 1/2, this method is first-order accurate, "
3056  "and with theta = 1/2, it is second-order accurate. "
3057  "This method is A-stable, but becomes L-stable with theta=1.");
3058 
3059  return pl;
3060  }
3061 
3062  protected:
3064  {
3065  typedef Teuchos::ScalarTraits<Scalar> ST;
3066  int NumStages = 1;
3067  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
3070  A(0, 0) = theta_;
3071  b(0) = ST::one();
3072  c(0) = theta_;
3073 
3074  int order = 1;
3075  if (std::abs((theta_ - thetaDefault_) / theta_) < 1.0e-08) order = 2;
3076 
3078  this->getStepperType(), A, b, c, order, 1, 2));
3079  this->tableau_->setTVD(true);
3080  this->tableau_->setTVDCoeff(2.0);
3081  }
3082 
3083  private:
3085  Scalar theta_;
3086 };
3087 
3089 // ------------------------------------------------------------------------
3090 template <class Scalar>
3092  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3094 {
3095  auto stepper = Teuchos::rcp(new StepperDIRK_1StageTheta<Scalar>());
3096  stepper->setStepperDIRKValues(pl);
3097 
3098  if (pl != Teuchos::null) {
3099  stepper->setTheta(pl->get<double>("theta", 0.5));
3100  }
3101 
3102  if (model != Teuchos::null) {
3103  stepper->setModel(model);
3104  stepper->initialize();
3105  }
3106 
3107  return stepper;
3108 }
3109 
3110 // ----------------------------------------------------------------------------
3137 template <class Scalar>
3138 class StepperEDIRK_2StageTheta : virtual public StepperDIRK<Scalar> {
3139  public:
3146  {
3147  typedef Teuchos::ScalarTraits<Scalar> ST;
3148  thetaDefault_ = ST::one() / (2 * ST::one());
3149 
3150  this->setStepperName("EDIRK 2 Stage Theta Method");
3151  this->setStepperType("EDIRK 2 Stage Theta Method");
3152  this->setTheta(thetaDefault_);
3153  this->setupTableau();
3154  this->setupDefault();
3155  this->setUseFSAL(true);
3156  this->setICConsistency("Consistent");
3157  this->setICConsistencyCheck(false);
3158  }
3159 
3161  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3163  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
3164  bool useEmbedded, bool zeroInitialGuess,
3165  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
3166  Scalar theta = Scalar(0.5))
3167  {
3168  typedef Teuchos::ScalarTraits<Scalar> ST;
3169  thetaDefault_ = ST::one() / (2 * ST::one());
3170 
3171  this->setStepperName("EDIRK 2 Stage Theta Method");
3172  this->setStepperType("EDIRK 2 Stage Theta Method");
3173  this->setTheta(theta);
3174  this->setupTableau();
3175  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3176  useEmbedded, zeroInitialGuess, stepperRKAppAction);
3177  }
3178 
3179  void setTheta(Scalar theta)
3180  {
3182  theta == Teuchos::ScalarTraits<Scalar>::zero(), std::logic_error,
3183  "'theta' can not be zero, as it makes this stepper explicit. \n"
3184  "Try using the 'RK Forward Euler' stepper.\n");
3185  theta_ = theta;
3186  this->isInitialized_ = false;
3187  this->setupTableau();
3188  }
3189 
3190  Scalar getTheta() const { return theta_; }
3191 
3192  std::string getDescription() const
3193  {
3194  std::ostringstream Description;
3195  Description << this->getStepperType() << "\n"
3196  << "Computer Methods for ODEs and DAEs\n"
3197  << "U. M. Ascher and L. R. Petzold\n"
3198  << "p. 113\n"
3199  << "c = [ 0 1 ]'\n"
3200  << "A = [ 0 0 ]\n"
3201  << " [ 1-theta theta ]\n"
3202  << "b = [ 1-theta theta ]'";
3203  return Description.str();
3204  }
3205 
3207  {
3208  auto pl = this->getValidParametersBasicDIRK();
3209 
3210  pl->template set<double>(
3211  "theta", getTheta(),
3212  "Valid values are 0 < theta <= 1, where theta = 0 "
3213  "implies Forward Euler, theta = 1/2 implies trapezoidal "
3214  "method (default), and theta = 1 implies Backward Euler. "
3215  "For theta != 1/2, this method is first-order accurate, "
3216  "and with theta = 1/2, it is second-order accurate. "
3217  "This method is A-stable, but becomes L-stable with theta=1.");
3218 
3219  return pl;
3220  }
3221 
3222  void setUseFSAL(bool a)
3223  {
3224  this->useFSAL_ = a;
3225  this->isInitialized_ = false;
3226  }
3227 
3228  protected:
3230  {
3231  typedef Teuchos::ScalarTraits<Scalar> ST;
3232  const Scalar one = ST::one();
3233  const Scalar zero = ST::zero();
3234 
3235  int NumStages = 2;
3236  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
3239 
3240  // Fill A:
3241  A(0, 0) = zero;
3242  A(0, 1) = zero;
3243  A(1, 0) = Teuchos::as<Scalar>(one - theta_);
3244  A(1, 1) = theta_;
3245 
3246  // Fill b:
3247  b(0) = Teuchos::as<Scalar>(one - theta_);
3248  b(1) = theta_;
3249 
3250  // Fill c:
3251  c(0) = zero;
3252  c(1) = one;
3253 
3254  int order = 1;
3255  if (std::abs((theta_ - thetaDefault_) / theta_) < 1.0e-08) order = 2;
3256 
3258  this->getStepperType(), A, b, c, order, 1, 2));
3259  this->tableau_->setTVD(true);
3260  this->tableau_->setTVDCoeff(2.0);
3261  }
3262 
3263  private:
3265  Scalar theta_;
3266 };
3267 
3269 // ------------------------------------------------------------------------
3270 template <class Scalar>
3272  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3274 {
3275  auto stepper = Teuchos::rcp(new StepperEDIRK_2StageTheta<Scalar>());
3276  stepper->setStepperDIRKValues(pl);
3277 
3278  if (pl != Teuchos::null) {
3279  stepper->setTheta(pl->get<double>("theta", 0.5));
3280  }
3281 
3282  if (model != Teuchos::null) {
3283  stepper->setModel(model);
3284  stepper->initialize();
3285  }
3286 
3287  return stepper;
3288 }
3289 
3290 // ----------------------------------------------------------------------------
3311 template <class Scalar>
3312 class StepperEDIRK_TrapezoidalRule : virtual public StepperDIRK<Scalar> {
3313  public:
3320  {
3321  this->setStepperName("RK Trapezoidal Rule");
3322  this->setStepperType("RK Trapezoidal Rule");
3323  this->setupTableau();
3324  this->setupDefault();
3325  this->setUseFSAL(true);
3326  this->setICConsistency("Consistent");
3327  this->setICConsistencyCheck(false);
3328  }
3329 
3331  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3333  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
3334  bool useEmbedded, bool zeroInitialGuess,
3335  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3336  {
3337  this->setStepperName("RK Trapezoidal Rule");
3338  this->setStepperType("RK Trapezoidal Rule");
3339  this->setupTableau();
3340  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3341  useEmbedded, zeroInitialGuess, stepperRKAppAction);
3342  }
3343 
3344  std::string getDescription() const
3345  {
3346  std::ostringstream Description;
3347  Description << this->getStepperType() << "\n"
3348  << "Also known as Crank-Nicolson Method.\n"
3349  << "c = [ 0 1 ]'\n"
3350  << "A = [ 0 0 ]\n"
3351  << " [ 1/2 1/2 ]\n"
3352  << "b = [ 1/2 1/2 ]'";
3353  return Description.str();
3354  }
3355 
3356  void setUseFSAL(bool a)
3357  {
3358  this->useFSAL_ = a;
3359  this->isInitialized_ = false;
3360  }
3361 
3362  protected:
3364  {
3365  typedef Teuchos::ScalarTraits<Scalar> ST;
3366  const Scalar one = ST::one();
3367  const Scalar zero = ST::zero();
3368  const Scalar onehalf = ST::one() / (2 * ST::one());
3369 
3370  int NumStages = 2;
3371  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
3374 
3375  // Fill A:
3376  A(0, 0) = zero;
3377  A(0, 1) = zero;
3378  A(1, 0) = onehalf;
3379  A(1, 1) = onehalf;
3380 
3381  // Fill b:
3382  b(0) = onehalf;
3383  b(1) = onehalf;
3384 
3385  // Fill c:
3386  c(0) = zero;
3387  c(1) = one;
3388 
3389  int order = 2;
3390 
3392  this->getStepperType(), A, b, c, order, order, order));
3393  this->tableau_->setTVD(true);
3394  this->tableau_->setTVDCoeff(2.0);
3395  }
3396 };
3397 
3399 // ------------------------------------------------------------------------
3400 template <class Scalar>
3403  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3405 {
3406  auto stepper = Teuchos::rcp(new StepperEDIRK_TrapezoidalRule<Scalar>());
3407 
3408  // Test for aliases.
3409  if (pl != Teuchos::null) {
3410  auto stepperType =
3411  pl->get<std::string>("Stepper Type", stepper->getStepperType());
3412 
3414  stepperType != stepper->getStepperType() &&
3415  stepperType != "RK Crank-Nicolson",
3416  std::logic_error,
3417  " ParameterList 'Stepper Type' (='" + stepperType
3418  << "')\n does not match type for this Stepper (='"
3419  << stepper->getStepperType()
3420  << "')\n or one of its aliases ('RK Crank-Nicolson').\n");
3421 
3422  // Reset default StepperType.
3423  pl->set<std::string>("Stepper Type", stepper->getStepperType());
3424  }
3425 
3426  stepper->setStepperDIRKValues(pl);
3427 
3428  if (model != Teuchos::null) {
3429  stepper->setModel(model);
3430  stepper->initialize();
3431  }
3432 
3433  return stepper;
3434 }
3435 
3436 // ----------------------------------------------------------------------------
3463 template <class Scalar>
3464 class StepperSDIRK_ImplicitMidpoint : virtual public StepperDIRK<Scalar> {
3465  public:
3472  {
3473  this->setStepperName("RK Implicit Midpoint");
3474  this->setStepperType("RK Implicit Midpoint");
3475  this->setupTableau();
3476  this->setupDefault();
3477  this->setUseFSAL(false);
3478  this->setICConsistency("None");
3479  this->setICConsistencyCheck(false);
3480  }
3481 
3483  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3485  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
3486  bool useEmbedded, bool zeroInitialGuess,
3487  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3488  {
3489  this->setStepperName("RK Implicit Midpoint");
3490  this->setStepperType("RK Implicit Midpoint");
3491  this->setupTableau();
3492  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3493  useEmbedded, zeroInitialGuess, stepperRKAppAction);
3494  }
3495 
3496  std::string getDescription() const
3497  {
3498  std::ostringstream Description;
3499  Description << this->getStepperType() << "\n"
3500  << "A-stable\n"
3501  << "Solving Ordinary Differential Equations II:\n"
3502  << "Stiff and Differential-Algebraic Problems,\n"
3503  << "2nd Revised Edition\n"
3504  << "E. Hairer and G. Wanner\n"
3505  << "Table 5.2, pg 72\n"
3506  << "Solving Ordinary Differential Equations I:\n"
3507  << "Nonstiff Problems, 2nd Revised Edition\n"
3508  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
3509  << "Table 7.1, pg 205\n"
3510  << "c = [ 1/2 ]'\n"
3511  << "A = [ 1/2 ]\n"
3512  << "b = [ 1 ]'";
3513  return Description.str();
3514  }
3515 
3516  protected:
3518  {
3519  typedef Teuchos::ScalarTraits<Scalar> ST;
3520  int NumStages = 1;
3521  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
3524  const Scalar onehalf = ST::one() / (2 * ST::one());
3525  const Scalar one = ST::one();
3526 
3527  // Fill A:
3528  A(0, 0) = onehalf;
3529 
3530  // Fill b:
3531  b(0) = one;
3532 
3533  // Fill c:
3534  c(0) = onehalf;
3535 
3536  int order = 2;
3537 
3539  this->getStepperType(), A, b, c, order, order, order));
3540  this->tableau_->setTVD(true);
3541  this->tableau_->setTVDCoeff(2.0);
3542  }
3543 };
3544 
3546 // ------------------------------------------------------------------------
3547 template <class Scalar>
3550  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3552 {
3554  stepper->setStepperDIRKValues(pl);
3555 
3556  if (model != Teuchos::null) {
3557  stepper->setModel(model);
3558  stepper->initialize();
3559  }
3560 
3561  return stepper;
3562 }
3563 
3564 // ----------------------------------------------------------------------------
3585 template <class Scalar>
3586 class StepperSDIRK_SSPDIRK22 : virtual public StepperDIRK<Scalar> {
3587  public:
3589  {
3590  this->setStepperName("SSPDIRK22");
3591  this->setStepperType("SSPDIRK22");
3592  this->setupTableau();
3593  this->setupDefault();
3594  this->setUseFSAL(false);
3595  this->setICConsistency("None");
3596  this->setICConsistencyCheck(false);
3597  }
3598 
3600  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3602  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
3603  bool useEmbedded, bool zeroInitialGuess,
3604  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3605  {
3606  this->setStepperName("SSPDIRK22");
3607  this->setStepperType("SSPDIRK22");
3608  this->setupTableau();
3609  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3610  useEmbedded, zeroInitialGuess, stepperRKAppAction);
3611  }
3612 
3613  std::string getDescription() const
3614  {
3615  std::ostringstream Description;
3616  Description << this->getStepperType() << "\n"
3617  << "Strong Stability Preserving Diagonally-Implicit RK "
3618  "(stage=2, order=2)\n"
3619  << "SSP-Coef = 4\n"
3620  << "c = [ 1/4 3/4 ]'\n"
3621  << "A = [ 1/4 ]\n"
3622  << " [ 1/2 1/4 ]\n"
3623  << "b = [ 1/2 1/2 ]\n"
3624  << std::endl;
3625  return Description.str();
3626  }
3627 
3628  protected:
3630  {
3631  typedef Teuchos::ScalarTraits<Scalar> ST;
3632  using Teuchos::as;
3633  const int NumStages = 2;
3634  const int order = 2;
3635  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
3638 
3639  const Scalar one = ST::one();
3640  const Scalar zero = ST::zero();
3641  const Scalar onehalf = one / (2 * one);
3642  const Scalar onefourth = one / (4 * one);
3643 
3644  // Fill A:
3645  A(0, 0) = A(1, 1) = onefourth;
3646  A(0, 1) = zero;
3647  A(1, 0) = onehalf;
3648 
3649  // Fill b:
3650  b(0) = b(1) = onehalf;
3651 
3652  // Fill c:
3653  c(0) = A(0, 0);
3654  c(1) = A(1, 0) + A(1, 1);
3655 
3657  this->getStepperType(), A, b, c, order, order, order));
3658  this->tableau_->setTVD(true);
3659  this->tableau_->setTVDCoeff(4.0);
3660  }
3661 };
3662 
3664 // ------------------------------------------------------------------------
3665 template <class Scalar>
3667  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3669 {
3670  auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK22<Scalar>());
3671  stepper->setStepperDIRKValues(pl);
3672 
3673  if (model != Teuchos::null) {
3674  stepper->setModel(model);
3675  stepper->initialize();
3676  }
3677 
3678  return stepper;
3679 }
3680 
3681 // ----------------------------------------------------------------------------
3703 template <class Scalar>
3704 class StepperSDIRK_SSPDIRK32 : virtual public StepperDIRK<Scalar> {
3705  public:
3707  {
3708  this->setStepperName("SSPDIRK32");
3709  this->setStepperType("SSPDIRK32");
3710  this->setupTableau();
3711  this->setupDefault();
3712  this->setUseFSAL(false);
3713  this->setICConsistency("None");
3714  this->setICConsistencyCheck(false);
3715  }
3716 
3718  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3720  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
3721  bool useEmbedded, bool zeroInitialGuess,
3722  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3723  {
3724  this->setStepperName("SSPDIRK32");
3725  this->setStepperType("SSPDIRK32");
3726  this->setupTableau();
3727  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3728  useEmbedded, zeroInitialGuess, stepperRKAppAction);
3729  }
3730 
3731  std::string getDescription() const
3732  {
3733  std::ostringstream Description;
3734  Description << this->getStepperType() << "\n"
3735  << "Strong Stability Preserving Diagonally-Implicit RK "
3736  "(stage=3, order=2)\n"
3737  << "SSP-Coef = 6\n"
3738  << "c = [ 1/6 1/2 5/6 ]'\n"
3739  << "A = [ 1/6 ]\n"
3740  << " [ 1/3 1/6 ]\n"
3741  << " [ 1/3 1/3 1/6 ]\n"
3742  << "b = [ 1/3 1/3 1/3 ]\n"
3743  << std::endl;
3744  return Description.str();
3745  }
3746 
3747  protected:
3749  {
3750  typedef Teuchos::ScalarTraits<Scalar> ST;
3751  using Teuchos::as;
3752  const int NumStages = 3;
3753  const int order = 2;
3754  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
3757 
3758  const Scalar one = ST::one();
3759  const Scalar zero = ST::zero();
3760  const Scalar onethird = one / (3 * one);
3761  const Scalar onesixth = one / (6 * one);
3762 
3763  // Fill A:
3764  A(0, 0) = A(1, 1) = A(2, 2) = onesixth;
3765  A(1, 0) = A(2, 0) = A(2, 1) = onethird;
3766  A(0, 1) = A(0, 2) = A(1, 2) = zero;
3767 
3768  // Fill b:
3769  b(0) = b(1) = b(2) = onethird;
3770 
3771  // Fill c:
3772  c(0) = A(0, 0);
3773  c(1) = A(1, 0) + A(1, 1);
3774  c(2) = A(2, 0) + A(2, 1) + A(2, 2);
3775 
3777  this->getStepperType(), A, b, c, order, order, order));
3778  this->tableau_->setTVD(true);
3779  this->tableau_->setTVDCoeff(6.0);
3780  }
3781 };
3782 
3784 // ------------------------------------------------------------------------
3785 template <class Scalar>
3787  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3789 {
3790  auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK32<Scalar>());
3791  stepper->setStepperDIRKValues(pl);
3792 
3793  if (model != Teuchos::null) {
3794  stepper->setModel(model);
3795  stepper->initialize();
3796  }
3797 
3798  return stepper;
3799 }
3800 
3801 // ----------------------------------------------------------------------------
3820 template <class Scalar>
3821 class StepperSDIRK_SSPDIRK23 : virtual public StepperDIRK<Scalar> {
3822  public:
3824  {
3825  this->setStepperName("SSPDIRK23");
3826  this->setStepperType("SSPDIRK23");
3827  this->setupTableau();
3828  this->setupDefault();
3829  this->setUseFSAL(false);
3830  this->setICConsistency("None");
3831  this->setICConsistencyCheck(false);
3832  }
3833 
3835  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3837  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
3838  bool useEmbedded, bool zeroInitialGuess,
3839  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3840  {
3841  this->setStepperName("SSPDIRK23");
3842  this->setStepperType("SSPDIRK23");
3843  this->setupTableau();
3844  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3845  useEmbedded, zeroInitialGuess, stepperRKAppAction);
3846  }
3847 
3848  std::string getDescription() const
3849  {
3850  std::ostringstream Description;
3851  Description << this->getStepperType() << "\n"
3852  << "Strong Stability Preserving Diagonally-Implicit RK "
3853  "(stage=2, order=3)\n"
3854  << "SSP-Coef = 1 + sqrt( 3 )\n"
3855  << "c = [ 1/(3 + sqrt( 3 )) (1/6)(3 + sqrt( 3 )) ] '\n"
3856  << "A = [ 1/(3 + sqrt( 3 )) ] \n"
3857  << " [ 1/sqrt( 3 ) 1/(3 + sqrt( 3 )) ] \n"
3858  << "b = [ 1/2 1/2 ] \n"
3859  << std::endl;
3860  return Description.str();
3861  }
3862 
3863  protected:
3865  {
3866  typedef Teuchos::ScalarTraits<Scalar> ST;
3867  using Teuchos::as;
3868  const int NumStages = 2;
3869  const int order = 3;
3870  const Scalar sspcoef = 2.7321;
3871  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
3874 
3875  const Scalar one = ST::one();
3876  const Scalar zero = ST::zero();
3877  const Scalar onehalf = one / (2 * one);
3878  const Scalar rootthree = ST::squareroot(3 * one);
3879 
3880  // Fill A:
3881  A(0, 0) = A(1, 1) = one / (3 * one + rootthree);
3882  A(1, 0) = one / rootthree;
3883  A(0, 1) = zero;
3884 
3885  // Fill b:
3886  b(0) = b(1) = onehalf;
3887 
3888  // Fill c:
3889  c(0) = A(0, 0);
3890  c(1) = A(1, 0) + A(1, 1);
3891 
3893  this->getStepperType(), A, b, c, order, order, order));
3894  this->tableau_->setTVD(true);
3895  this->tableau_->setTVDCoeff(sspcoef);
3896  }
3897 };
3898 
3900 // ------------------------------------------------------------------------
3901 template <class Scalar>
3903  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3905 {
3906  auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK23<Scalar>());
3907  stepper->setStepperDIRKValues(pl);
3908 
3909  if (model != Teuchos::null) {
3910  stepper->setModel(model);
3911  stepper->initialize();
3912  }
3913 
3914  return stepper;
3915 }
3916 
3917 // ----------------------------------------------------------------------------
3939 template <class Scalar>
3940 class StepperSDIRK_SSPDIRK33 : virtual public StepperDIRK<Scalar> {
3941  public:
3943  {
3944  this->setStepperName("SSPDIRK33");
3945  this->setStepperType("SSPDIRK33");
3946  this->setupTableau();
3947  this->setupDefault();
3948  this->setUseFSAL(false);
3949  this->setICConsistency("None");
3950  this->setICConsistencyCheck(false);
3951  }
3952 
3954  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3956  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
3957  bool useEmbedded, bool zeroInitialGuess,
3958  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3959  {
3960  this->setStepperName("SSPDIRK33");
3961  this->setStepperType("SSPDIRK33");
3962  this->setupTableau();
3963  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3964  useEmbedded, zeroInitialGuess, stepperRKAppAction);
3965  }
3966 
3967  std::string getDescription() const
3968  {
3969  std::ostringstream Description;
3970  Description << this->getStepperType() << "\n"
3971  << "Strong Stability Preserving Diagonally-Implicit RK "
3972  "(stage=3, order=3)\n"
3973  << "SSP-Coef = 2 + 2 sqrt(2)\n"
3974  << "c = [ 1/( 4 + 2 sqrt(2) 1/2 (1/4)(2 + "
3975  "sqrt(2) ] '\n"
3976  << "A = [ 1/( 4 + 2 sqrt(2) "
3977  " ] \n"
3978  << " [ 1/(2 sqrt(2) 1/( 4 + 2 sqrt(2) "
3979  " ] \n"
3980  << " [ 1/(2 sqrt(2) 1/(2 sqrt(2) 1/( 4 + 2 "
3981  "sqrt(2) ] \n"
3982  << "b = [ 1/3 1/3 1/3 "
3983  " ] \n"
3984  << std::endl;
3985  return Description.str();
3986  }
3987 
3988  protected:
3990  {
3991  typedef Teuchos::ScalarTraits<Scalar> ST;
3992  using Teuchos::as;
3993  const int NumStages = 3;
3994  const int order = 3;
3995  const Scalar sspcoef = 4.8284;
3996  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
3999 
4000  const Scalar one = ST::one();
4001  const Scalar zero = ST::zero();
4002  const Scalar onethird = one / (3 * one);
4003  const Scalar rootwo = ST::squareroot(2 * one);
4004 
4005  // Fill A:
4006  A(0, 0) = A(1, 1) = A(2, 2) = one / (4 * one + 2 * rootwo);
4007  A(1, 0) = A(2, 0) = A(2, 1) = one / (2 * rootwo);
4008  A(0, 1) = A(0, 2) = A(1, 2) = zero;
4009 
4010  // Fill b:
4011  b(0) = b(1) = b(2) = onethird;
4012 
4013  // Fill c:
4014  c(0) = A(0, 0);
4015  c(1) = A(1, 0) + A(1, 1);
4016  c(2) = A(2, 0) + A(2, 1) + A(2, 2);
4017 
4019  this->getStepperType(), A, b, c, order, order, order));
4020  this->tableau_->setTVD(true);
4021  this->tableau_->setTVDCoeff(sspcoef);
4022  }
4023 };
4024 
4026 // ------------------------------------------------------------------------
4027 template <class Scalar>
4029  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4031 {
4032  auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK33<Scalar>());
4033  stepper->setStepperDIRKValues(pl);
4034 
4035  if (model != Teuchos::null) {
4036  stepper->setModel(model);
4037  stepper->initialize();
4038  }
4039 
4040  return stepper;
4041 }
4042 
4043 // ----------------------------------------------------------------------------
4064 template <class Scalar>
4065 class StepperDIRK_1Stage1stOrderRadauIA : virtual public StepperDIRK<Scalar> {
4066  public:
4073  {
4074  this->setStepperName("RK Implicit 1 Stage 1st order Radau IA");
4075  this->setStepperType("RK Implicit 1 Stage 1st order Radau IA");
4076  this->setupTableau();
4077  this->setupDefault();
4078  this->setUseFSAL(false);
4079  this->setICConsistency("None");
4080  this->setICConsistencyCheck(false);
4081  }
4082 
4084  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4086  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
4087  bool useEmbedded, bool zeroInitialGuess,
4088  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4089  {
4090  this->setStepperName("RK Implicit 1 Stage 1st order Radau IA");
4091  this->setStepperType("RK Implicit 1 Stage 1st order Radau IA");
4092  this->setupTableau();
4093  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4094  useEmbedded, zeroInitialGuess, stepperRKAppAction);
4095  }
4096 
4097  std::string getDescription() const
4098  {
4099  std::ostringstream Description;
4100  Description << this->getStepperType() << "\n"
4101  << "A-stable\n"
4102  << "Solving Ordinary Differential Equations II:\n"
4103  << "Stiff and Differential-Algebraic Problems,\n"
4104  << "2nd Revised Edition\n"
4105  << "E. Hairer and G. Wanner\n"
4106  << "Table 5.3, pg 73\n"
4107  << "c = [ 0 ]'\n"
4108  << "A = [ 1 ]\n"
4109  << "b = [ 1 ]'";
4110  return Description.str();
4111  }
4112 
4113  protected:
4115  {
4116  typedef Teuchos::ScalarTraits<Scalar> ST;
4117  int NumStages = 1;
4118  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
4121  const Scalar one = ST::one();
4122  const Scalar zero = ST::zero();
4123  A(0, 0) = one;
4124  b(0) = one;
4125  c(0) = zero;
4126  int order = 1;
4127 
4128  auto emptyBStar = Teuchos::SerialDenseVector<int, Scalar>();
4129  this->tableau_ = Teuchos::rcp(
4130  new RKButcherTableau<Scalar>(this->getStepperType(), A, b, c, order,
4131  order, order, emptyBStar, false));
4132  }
4133 };
4134 
4136 // ------------------------------------------------------------------------
4137 template <class Scalar>
4140  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4142 {
4144  stepper->setStepperDIRKValues(pl);
4145 
4146  if (model != Teuchos::null) {
4147  stepper->setModel(model);
4148  stepper->initialize();
4149  }
4150 
4151  return stepper;
4152 }
4153 
4154 // ----------------------------------------------------------------------------
4177 template <class Scalar>
4179  : virtual public StepperDIRK<Scalar> {
4180  public:
4187  {
4188  this->setStepperName("RK Implicit 2 Stage 2nd order Lobatto IIIB");
4189  this->setStepperType("RK Implicit 2 Stage 2nd order Lobatto IIIB");
4190  this->setupTableau();
4191  this->setupDefault();
4192  this->setUseFSAL(false);
4193  this->setICConsistency("None");
4194  this->setICConsistencyCheck(false);
4195  }
4196 
4198  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4200  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
4201  bool useEmbedded, bool zeroInitialGuess,
4202  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4203  {
4204  this->setStepperName("RK Implicit 2 Stage 2nd order Lobatto IIIB");
4205  this->setStepperType("RK Implicit 2 Stage 2nd order Lobatto IIIB");
4206  this->setupTableau();
4207  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4208  useEmbedded, zeroInitialGuess, stepperRKAppAction);
4209  }
4210 
4211  std::string getDescription() const
4212  {
4213  std::ostringstream Description;
4214  Description << this->getStepperType() << "\n"
4215  << "A-stable\n"
4216  << "Solving Ordinary Differential Equations II:\n"
4217  << "Stiff and Differential-Algebraic Problems,\n"
4218  << "2nd Revised Edition\n"
4219  << "E. Hairer and G. Wanner\n"
4220  << "Table 5.9, pg 76\n"
4221  << "c = [ 0 1 ]'\n"
4222  << "A = [ 1/2 0 ]\n"
4223  << " [ 1/2 0 ]\n"
4224  << "b = [ 1/2 1/2 ]'";
4225  return Description.str();
4226  }
4227 
4228  protected:
4230  {
4231  typedef Teuchos::ScalarTraits<Scalar> ST;
4232  using Teuchos::as;
4233  int NumStages = 2;
4234  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
4237  const Scalar zero = ST::zero();
4238  const Scalar one = ST::one();
4239 
4240  // Fill A:
4241  A(0, 0) = as<Scalar>(one / (2 * one));
4242  A(0, 1) = zero;
4243  A(1, 0) = as<Scalar>(one / (2 * one));
4244  A(1, 1) = zero;
4245 
4246  // Fill b:
4247  b(0) = as<Scalar>(one / (2 * one));
4248  b(1) = as<Scalar>(one / (2 * one));
4249 
4250  // Fill c:
4251  c(0) = zero;
4252  c(1) = one;
4253  int order = 2;
4254 
4255  auto emptyBStar = Teuchos::SerialDenseVector<int, Scalar>();
4256  this->tableau_ = Teuchos::rcp(
4257  new RKButcherTableau<Scalar>(this->getStepperType(), A, b, c, order,
4258  order, order, emptyBStar, false));
4259  this->tableau_->setTVD(true);
4260  this->tableau_->setTVDCoeff(2.0);
4261  }
4262 };
4263 
4265 // ------------------------------------------------------------------------
4266 template <class Scalar>
4269  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4271 {
4272  auto stepper =
4274  stepper->setStepperDIRKValues(pl);
4275 
4276  if (model != Teuchos::null) {
4277  stepper->setModel(model);
4278  stepper->initialize();
4279  }
4280 
4281  return stepper;
4282 }
4283 
4284 // ----------------------------------------------------------------------------
4310 template <class Scalar>
4311 class StepperSDIRK_5Stage4thOrder : virtual public StepperDIRK<Scalar> {
4312  public:
4319  {
4320  this->setStepperName("SDIRK 5 Stage 4th order");
4321  this->setStepperType("SDIRK 5 Stage 4th order");
4322  this->setupTableau();
4323  this->setupDefault();
4324  this->setUseFSAL(false);
4325  this->setICConsistency("None");
4326  this->setICConsistencyCheck(false);
4327  }
4328 
4330  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4332  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
4333  bool useEmbedded, bool zeroInitialGuess,
4334  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4335  {
4336  this->setStepperName("SDIRK 5 Stage 4th order");
4337  this->setStepperType("SDIRK 5 Stage 4th order");
4338  this->setupTableau();
4339  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4340  useEmbedded, zeroInitialGuess, stepperRKAppAction);
4341  }
4342 
4343  std::string getDescription() const
4344  {
4345  std::ostringstream Description;
4346  Description << this->getStepperType() << "\n"
4347  << "L-stable\n"
4348  << "Solving Ordinary Differential Equations II:\n"
4349  << "Stiff and Differential-Algebraic Problems,\n"
4350  << "2nd Revised Edition\n"
4351  << "E. Hairer and G. Wanner\n"
4352  << "pg100 \n"
4353  << "c = [ 1/4 3/4 11/20 1/2 1 ]'\n"
4354  << "A = [ 1/4 ]\n"
4355  << " [ 1/2 1/4 ]\n"
4356  << " [ 17/50 -1/25 1/4 ]\n"
4357  << " [ 371/1360 -137/2720 15/544 1/4 ]\n"
4358  << " [ 25/24 -49/48 125/16 -85/12 1/4 ]\n"
4359  << "b = [ 25/24 -49/48 125/16 -85/12 1/4 ]'";
4360  // << "b = [ 59/48 -17/96 225/32 -85/12 0 ]'";
4361  return Description.str();
4362  }
4363 
4364  protected:
4366  {
4367  typedef Teuchos::ScalarTraits<Scalar> ST;
4368  using Teuchos::as;
4369  int NumStages = 5;
4370  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
4373  const Scalar zero = ST::zero();
4374  const Scalar one = ST::one();
4375  const Scalar onequarter = as<Scalar>(one / (4 * one));
4376 
4377  // Fill A:
4378  A(0, 0) = onequarter;
4379  A(0, 1) = zero;
4380  A(0, 2) = zero;
4381  A(0, 3) = zero;
4382  A(0, 4) = zero;
4383 
4384  A(1, 0) = as<Scalar>(one / (2 * one));
4385  A(1, 1) = onequarter;
4386  A(1, 2) = zero;
4387  A(1, 3) = zero;
4388  A(1, 4) = zero;
4389 
4390  A(2, 0) = as<Scalar>(17 * one / (50 * one));
4391  A(2, 1) = as<Scalar>(-one / (25 * one));
4392  A(2, 2) = onequarter;
4393  A(2, 3) = zero;
4394  A(2, 4) = zero;
4395 
4396  A(3, 0) = as<Scalar>(371 * one / (1360 * one));
4397  A(3, 1) = as<Scalar>(-137 * one / (2720 * one));
4398  A(3, 2) = as<Scalar>(15 * one / (544 * one));
4399  A(3, 3) = onequarter;
4400  A(3, 4) = zero;
4401 
4402  A(4, 0) = as<Scalar>(25 * one / (24 * one));
4403  A(4, 1) = as<Scalar>(-49 * one / (48 * one));
4404  A(4, 2) = as<Scalar>(125 * one / (16 * one));
4405  A(4, 3) = as<Scalar>(-85 * one / (12 * one));
4406  A(4, 4) = onequarter;
4407 
4408  // Fill b:
4409  b(0) = as<Scalar>(25 * one / (24 * one));
4410  b(1) = as<Scalar>(-49 * one / (48 * one));
4411  b(2) = as<Scalar>(125 * one / (16 * one));
4412  b(3) = as<Scalar>(-85 * one / (12 * one));
4413  b(4) = onequarter;
4414 
4415  /*
4416  // Alternate version
4417  b(0) = as<Scalar>( 59*one/(48*one) );
4418  b(1) = as<Scalar>( -17*one/(96*one) );
4419  b(2) = as<Scalar>( 225*one/(32*one) );
4420  b(3) = as<Scalar>( -85*one/(12*one) );
4421  b(4) = zero;
4422  */
4423 
4424  // Fill c:
4425  c(0) = onequarter;
4426  c(1) = as<Scalar>(3 * one / (4 * one));
4427  c(2) = as<Scalar>(11 * one / (20 * one));
4428  c(3) = as<Scalar>(one / (2 * one));
4429  c(4) = one;
4430 
4431  int order = 4;
4432 
4434  this->getStepperType(), A, b, c, order, order, order));
4435  }
4436 };
4437 
4439 // ------------------------------------------------------------------------
4440 template <class Scalar>
4443  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4445 {
4446  auto stepper = Teuchos::rcp(new StepperSDIRK_5Stage4thOrder<Scalar>());
4447  stepper->setStepperDIRKValues(pl);
4448 
4449  if (model != Teuchos::null) {
4450  stepper->setModel(model);
4451  stepper->initialize();
4452  }
4453 
4454  return stepper;
4455 }
4456 
4457 // ----------------------------------------------------------------------------
4482 template <class Scalar>
4483 class StepperSDIRK_3Stage4thOrder : virtual public StepperDIRK<Scalar> {
4484  public:
4491  {
4492  this->setStepperName("SDIRK 3 Stage 4th order");
4493  this->setStepperType("SDIRK 3 Stage 4th order");
4494  this->setupTableau();
4495  this->setupDefault();
4496  this->setUseFSAL(false);
4497  this->setICConsistency("None");
4498  this->setICConsistencyCheck(false);
4499  }
4500 
4502  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4504  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
4505  bool useEmbedded, bool zeroInitialGuess,
4506  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4507  {
4508  this->setStepperName("SDIRK 3 Stage 4th order");
4509  this->setStepperType("SDIRK 3 Stage 4th order");
4510  this->setupTableau();
4511  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4512  useEmbedded, zeroInitialGuess, stepperRKAppAction);
4513  }
4514 
4515  std::string getDescription() const
4516  {
4517  std::ostringstream Description;
4518  Description << this->getStepperType() << "\n"
4519  << "A-stable\n"
4520  << "Solving Ordinary Differential Equations II:\n"
4521  << "Stiff and Differential-Algebraic Problems,\n"
4522  << "2nd Revised Edition\n"
4523  << "E. Hairer and G. Wanner\n"
4524  << "p. 100 \n"
4525  << "gamma = (1/sqrt(3))*cos(pi/18)+1/2\n"
4526  << "delta = 1/(6*(2*gamma-1)^2)\n"
4527  << "c = [ gamma 1/2 1-gamma ]'\n"
4528  << "A = [ gamma ]\n"
4529  << " [ 1/2-gamma gamma ]\n"
4530  << " [ 2*gamma 1-4*gamma gamma ]\n"
4531  << "b = [ delta 1-2*delta delta ]'";
4532  return Description.str();
4533  }
4534 
4535  protected:
4537  {
4538  typedef Teuchos::ScalarTraits<Scalar> ST;
4539  using Teuchos::as;
4540  int NumStages = 3;
4541  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
4544  const Scalar zero = ST::zero();
4545  const Scalar one = ST::one();
4546  const Scalar pi = as<Scalar>(4 * one) * std::atan(one);
4547  const Scalar gamma =
4548  as<Scalar>(one / ST::squareroot(3 * one) * std::cos(pi / (18 * one)) +
4549  one / (2 * one));
4550  const Scalar delta =
4551  as<Scalar>(one / (6 * one * std::pow(2 * gamma - one, 2 * one)));
4552 
4553  // Fill A:
4554  A(0, 0) = gamma;
4555  A(0, 1) = zero;
4556  A(0, 2) = zero;
4557 
4558  A(1, 0) = as<Scalar>(one / (2 * one) - gamma);
4559  A(1, 1) = gamma;
4560  A(1, 2) = zero;
4561 
4562  A(2, 0) = as<Scalar>(2 * gamma);
4563  A(2, 1) = as<Scalar>(one - 4 * gamma);
4564  A(2, 2) = gamma;
4565 
4566  // Fill b:
4567  b(0) = delta;
4568  b(1) = as<Scalar>(one - 2 * delta);
4569  b(2) = delta;
4570 
4571  // Fill c:
4572  c(0) = gamma;
4573  c(1) = as<Scalar>(one / (2 * one));
4574  c(2) = as<Scalar>(one - gamma);
4575 
4576  int order = 4;
4577 
4579  this->getStepperType(), A, b, c, order, order, order));
4580  }
4581 };
4582 
4584 // ------------------------------------------------------------------------
4585 template <class Scalar>
4588  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4590 {
4591  auto stepper = Teuchos::rcp(new StepperSDIRK_3Stage4thOrder<Scalar>());
4592  stepper->setStepperDIRKValues(pl);
4593 
4594  if (model != Teuchos::null) {
4595  stepper->setModel(model);
4596  stepper->initialize();
4597  }
4598 
4599  return stepper;
4600 }
4601 
4602 // ----------------------------------------------------------------------------
4631 template <class Scalar>
4632 class StepperSDIRK_5Stage5thOrder : virtual public StepperDIRK<Scalar> {
4633  public:
4640  {
4641  this->setStepperName("SDIRK 5 Stage 5th order");
4642  this->setStepperType("SDIRK 5 Stage 5th order");
4643  this->setupTableau();
4644  this->setupDefault();
4645  this->setUseFSAL(false);
4646  this->setICConsistency("None");
4647  this->setICConsistencyCheck(false);
4648  }
4649 
4651  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4653  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
4654  bool useEmbedded, bool zeroInitialGuess,
4655  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4656  {
4657  this->setStepperName("SDIRK 5 Stage 5th order");
4658  this->setStepperType("SDIRK 5 Stage 5th order");
4659  this->setupTableau();
4660  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4661  useEmbedded, zeroInitialGuess, stepperRKAppAction);
4662  }
4663 
4664  std::string getDescription() const
4665  {
4666  std::ostringstream Description;
4667  Description << this->getStepperType() << "\n"
4668  << "Solving Ordinary Differential Equations II:\n"
4669  << "Stiff and Differential-Algebraic Problems,\n"
4670  << "2nd Revised Edition\n"
4671  << "E. Hairer and G. Wanner\n"
4672  << "pg101 \n"
4673  << "c = [ (6-sqrt(6))/10 ]\n"
4674  << " [ (6+9*sqrt(6))/35 ]\n"
4675  << " [ 1 ]\n"
4676  << " [ (4-sqrt(6))/10 ]\n"
4677  << " [ (4+sqrt(6))/10 ]\n"
4678  << "A = [ A1 A2 A3 A4 A5 ]\n"
4679  << " A1 = [ (6-sqrt(6))/10 ]\n"
4680  << " [ (-6+5*sqrt(6))/14 ]\n"
4681  << " [ (888+607*sqrt(6))/2850 ]\n"
4682  << " [ (3153-3082*sqrt(6))/14250 ]\n"
4683  << " [ (-32583+14638*sqrt(6))/71250 ]\n"
4684  << " A2 = [ 0 ]\n"
4685  << " [ (6-sqrt(6))/10 ]\n"
4686  << " [ (126-161*sqrt(6))/1425 ]\n"
4687  << " [ (3213+1148*sqrt(6))/28500 ]\n"
4688  << " [ (-17199+364*sqrt(6))/142500 ]\n"
4689  << " A3 = [ 0 ]\n"
4690  << " [ 0 ]\n"
4691  << " [ (6-sqrt(6))/10 ]\n"
4692  << " [ (-267+88*sqrt(6))/500 ]\n"
4693  << " [ (1329-544*sqrt(6))/2500 ]\n"
4694  << " A4 = [ 0 ]\n"
4695  << " [ 0 ]\n"
4696  << " [ 0 ]\n"
4697  << " [ (6-sqrt(6))/10 ]\n"
4698  << " [ (-96+131*sqrt(6))/625 ]\n"
4699  << " A5 = [ 0 ]\n"
4700  << " [ 0 ]\n"
4701  << " [ 0 ]\n"
4702  << " [ 0 ]\n"
4703  << " [ (6-sqrt(6))/10 ]\n"
4704  << "b = [ 0 ]\n"
4705  << " [ 0 ]\n"
4706  << " [ 1/9 ]\n"
4707  << " [ (16-sqrt(6))/36 ]\n"
4708  << " [ (16+sqrt(6))/36 ]'";
4709  return Description.str();
4710  }
4711 
4712  protected:
4714  {
4715  typedef Teuchos::ScalarTraits<Scalar> ST;
4716  using Teuchos::as;
4717  int NumStages = 5;
4718  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
4721  const Scalar zero = ST::zero();
4722  const Scalar one = ST::one();
4723  const Scalar sqrt6 = ST::squareroot(as<Scalar>(6 * one));
4724  const Scalar gamma =
4725  as<Scalar>((6 * one - sqrt6) / (10 * one)); // diagonal
4726 
4727  // Fill A:
4728  A(0, 0) = gamma;
4729  A(0, 1) = zero;
4730  A(0, 2) = zero;
4731  A(0, 3) = zero;
4732  A(0, 4) = zero;
4733 
4734  A(1, 0) = as<Scalar>((-6 * one + 5 * one * sqrt6) / (14 * one));
4735  A(1, 1) = gamma;
4736  A(1, 2) = zero;
4737  A(1, 3) = zero;
4738  A(1, 4) = zero;
4739 
4740  A(2, 0) = as<Scalar>((888 * one + 607 * one * sqrt6) / (2850 * one));
4741  A(2, 1) = as<Scalar>((126 * one - 161 * one * sqrt6) / (1425 * one));
4742  A(2, 2) = gamma;
4743  A(2, 3) = zero;
4744  A(2, 4) = zero;
4745 
4746  A(3, 0) = as<Scalar>((3153 * one - 3082 * one * sqrt6) / (14250 * one));
4747  A(3, 1) = as<Scalar>((3213 * one + 1148 * one * sqrt6) / (28500 * one));
4748  A(3, 2) = as<Scalar>((-267 * one + 88 * one * sqrt6) / (500 * one));
4749  A(3, 3) = gamma;
4750  A(3, 4) = zero;
4751 
4752  A(4, 0) = as<Scalar>((-32583 * one + 14638 * one * sqrt6) / (71250 * one));
4753  A(4, 1) = as<Scalar>((-17199 * one + 364 * one * sqrt6) / (142500 * one));
4754  A(4, 2) = as<Scalar>((1329 * one - 544 * one * sqrt6) / (2500 * one));
4755  A(4, 3) = as<Scalar>((-96 * one + 131 * sqrt6) / (625 * one));
4756  A(4, 4) = gamma;
4757 
4758  // Fill b:
4759  b(0) = zero;
4760  b(1) = zero;
4761  b(2) = as<Scalar>(one / (9 * one));
4762  b(3) = as<Scalar>((16 * one - sqrt6) / (36 * one));
4763  b(4) = as<Scalar>((16 * one + sqrt6) / (36 * one));
4764 
4765  // Fill c:
4766  c(0) = gamma;
4767  c(1) = as<Scalar>((6 * one + 9 * one * sqrt6) / (35 * one));
4768  c(2) = one;
4769  c(3) = as<Scalar>((4 * one - sqrt6) / (10 * one));
4770  c(4) = as<Scalar>((4 * one + sqrt6) / (10 * one));
4771 
4772  int order = 5;
4773 
4775  this->getStepperType(), A, b, c, order, order, order));
4776  }
4777 };
4778 
4780 // ------------------------------------------------------------------------
4781 template <class Scalar>
4784  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4786 {
4787  auto stepper = Teuchos::rcp(new StepperSDIRK_5Stage5thOrder<Scalar>());
4788  stepper->setStepperDIRKValues(pl);
4789 
4790  if (model != Teuchos::null) {
4791  stepper->setModel(model);
4792  stepper->initialize();
4793  }
4794 
4795  return stepper;
4796 }
4797 
4798 // ----------------------------------------------------------------------------
4817 template <class Scalar>
4818 class StepperSDIRK_21Pair : virtual public StepperDIRK<Scalar> {
4819  public:
4826  {
4827  this->setStepperName("SDIRK 2(1) Pair");
4828  this->setStepperType("SDIRK 2(1) Pair");
4829  this->setupTableau();
4830  this->setupDefault();
4831  this->setUseFSAL(false);
4832  this->setICConsistency("None");
4833  this->setICConsistencyCheck(false);
4834  }
4835 
4837  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4839  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
4840  bool useEmbedded, bool zeroInitialGuess,
4841  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4842  {
4843  this->setStepperName("SDIRK 2(1) Pair");
4844  this->setStepperType("SDIRK 2(1) Pair");
4845  this->setupTableau();
4846  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4847  useEmbedded, zeroInitialGuess, stepperRKAppAction);
4848  }
4849 
4850  std::string getDescription() const
4851  {
4852  std::ostringstream Description;
4853  Description << this->getStepperType() << "\n"
4854  << "c = [ 1 0 ]'\n"
4855  << "A = [ 1 ]\n"
4856  << " [ -1 1 ]\n"
4857  << "b = [ 1/2 1/2 ]'\n"
4858  << "bstar = [ 1 0 ]'";
4859  return Description.str();
4860  }
4861 
4862  protected:
4864  {
4865  typedef Teuchos::ScalarTraits<Scalar> ST;
4866  using Teuchos::as;
4867  int NumStages = 2;
4868  Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
4871  Teuchos::SerialDenseVector<int, Scalar> bstar(NumStages);
4872 
4873  const Scalar one = ST::one();
4874  const Scalar zero = ST::zero();
4875 
4876  // Fill A:
4877  A(0, 0) = one;
4878  A(0, 1) = zero;
4879  A(1, 0) = -one;
4880  A(1, 1) = one;
4881 
4882  // Fill b:
4883  b(0) = as<Scalar>(one / (2 * one));
4884  b(1) = as<Scalar>(one / (2 * one));
4885 
4886  // Fill c:
4887  c(0) = one;
4888  c(1) = zero;
4889 
4890  // Fill bstar
4891  bstar(0) = one;
4892  bstar(1) = zero;
4893  int order = 2;
4894 
4896  this->getStepperType(), A, b, c, order, order, order, bstar));
4897  }
4898 };
4899 
4901 // ------------------------------------------------------------------------
4902 template <class Scalar>
4904  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4906 {
4907  auto stepper = Teuchos::rcp(new StepperSDIRK_21Pair<Scalar>());
4908  stepper->setStepperDIRKValues(pl);
4909 
4910  if (model != Teuchos::null) {
4911  stepper->setModel(model);
4912  stepper->initialize();
4913  }
4914 
4915  return stepper;
4916 }
4917 
4918 // ----------------------------------------------------------------------------
4951 template <class Scalar>
4952 class StepperDIRK_General : virtual public StepperDIRK<Scalar> {
4953  public:
4960  {
4961  this->setStepperName("General DIRK");
4962  this->setStepperType("General DIRK");
4963  this->setupTableau();
4964  this->setupDefault();
4965  this->setUseFSAL(false);
4966  this->setICConsistency("None");
4967  this->setICConsistencyCheck(false);
4968  }
4969 
4971  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4973  bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
4974  bool useEmbedded, bool zeroInitialGuess,
4975  const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
4978  const Teuchos::SerialDenseVector<int, Scalar>& c, const int order,
4979  const int orderMin, const int orderMax,
4981  {
4982  this->setStepperName("General DIRK");
4983  this->setStepperType("General DIRK");
4984  this->setTableau(A, b, c, order, orderMin, orderMax, bstar);
4985 
4987  this->tableau_->isImplicit() != true, std::logic_error,
4988  "Error - General DIRK did not receive a DIRK Butcher Tableau!\n");
4989 
4990  this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4991  useEmbedded, zeroInitialGuess, stepperRKAppAction);
4992  }
4993 
4994  std::string getDescription() const
4995  {
4996  std::stringstream Description;
4997  Description
4998  << this->getStepperType() << "\n"
4999  << "The format of the Butcher Tableau parameter list is\n"
5000  << " <Parameter name=\"A\" type=\"string\" value=\"# # # ;\n"
5001  << " # # # ;\n"
5002  << " # # #\"/>\n"
5003  << " <Parameter name=\"b\" type=\"string\" value=\"# # #\"/>\n"
5004  << " <Parameter name=\"c\" type=\"string\" value=\"# # #\"/>\n\n"
5005  << "Note the number of stages is implicit in the number of entries.\n"
5006  << "The number of stages must be consistent.\n"
5007  << "\n"
5008  << "Default tableau is 'SDIRK 2 Stage 2nd order':\n"
5009  << " Computer Methods for ODEs and DAEs\n"
5010  << " U. M. Ascher and L. R. Petzold\n"
5011  << " p. 106\n"
5012  << " gamma = (2-sqrt(2))/2\n"
5013  << " c = [ gamma 1 ]'\n"
5014  << " A = [ gamma 0 ]\n"
5015  << " [ 1-gamma gamma ]\n"
5016  << " b = [ 1-gamma gamma ]'";
5017  return Description.str();
5018  }
5019 
5021  {
5022  if (this->tableau_ == Teuchos::null) {
5023  // Set tableau to the default if null, otherwise keep current tableau.
5024  auto stepper = Teuchos::rcp(new StepperSDIRK_2Stage2ndOrder<Scalar>());
5025  auto t = stepper->getTableau();
5027  this->getStepperType(), t->A(), t->b(), t->c(), t->order(),
5028  t->orderMin(), t->orderMax(), t->bstar()));
5029  this->isInitialized_ = false;
5030  }
5031  }
5032 
5036  const int order, const int orderMin, const int orderMax,
5039  {
5041  this->getStepperType(), A, b, c, order, orderMin, orderMax, bstar));
5042  this->isInitialized_ = false;
5043  }
5044 
5046  {
5047  auto pl = this->getValidParametersBasicDIRK();
5048 
5049  // Tableau ParameterList
5053  Teuchos::SerialDenseVector<int, Scalar> bstar = this->tableau_->bstar();
5054 
5055  Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
5056 
5057  std::ostringstream Astream;
5058  Astream.precision(15);
5059  for (int i = 0; i < A.numRows(); i++) {
5060  for (int j = 0; j < A.numCols() - 1; j++) {
5061  Astream << A(i, j) << " ";
5062  }
5063  Astream << A(i, A.numCols() - 1);
5064  if (i != A.numRows() - 1) Astream << "; ";
5065  }
5066  tableauPL->set<std::string>("A", Astream.str());
5067 
5068  std::ostringstream bstream;
5069  bstream.precision(15);
5070  for (int i = 0; i < b.length() - 1; i++) {
5071  bstream << b(i) << " ";
5072  }
5073  bstream << b(b.length() - 1);
5074  tableauPL->set<std::string>("b", bstream.str());
5075 
5076  std::ostringstream cstream;
5077  cstream.precision(15);
5078  for (int i = 0; i < c.length() - 1; i++) {
5079  cstream << c(i) << " ";
5080  }
5081  cstream << c(c.length() - 1);
5082  tableauPL->set<std::string>("c", cstream.str());
5083 
5084  tableauPL->set<int>("order", this->tableau_->order());
5085 
5086  if (bstar.length() == 0) {
5087  tableauPL->set("bstar", "");
5088  }
5089  else {
5090  std::ostringstream bstarstream;
5091  bstarstream.precision(15);
5092  for (int i = 0; i < bstar.length() - 1; i++) {
5093  bstarstream << bstar(i) << " ";
5094  }
5095  bstarstream << bstar(bstar.length() - 1);
5096  tableauPL->set<std::string>("bstar", bstarstream.str());
5097  }
5098 
5099  pl->set("Tableau", *tableauPL);
5100 
5101  return pl;
5102  }
5103 };
5104 
5106 // ------------------------------------------------------------------------
5107 template <class Scalar>
5109  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
5111 {
5112  auto stepper = Teuchos::rcp(new StepperDIRK_General<Scalar>());
5113  stepper->setStepperDIRKValues(pl);
5114 
5115  if (pl != Teuchos::null) {
5116  if (pl->isParameter("Tableau")) {
5117  auto t = stepper->createTableau(pl);
5118  stepper->setTableau(t->A(), t->b(), t->c(), t->order(), t->orderMin(),
5119  t->orderMax(), t->bstar());
5120  }
5121  }
5123  stepper->getTableau()->isDIRK() != true, std::logic_error,
5124  "Error - General DIRK did not receive a DIRK Butcher Tableau!\n");
5125 
5126  if (model != Teuchos::null) {
5127  stepper->setModel(model);
5128  stepper->initialize();
5129  }
5130 
5131  return stepper;
5132 }
5133 
5134 } // namespace Tempus
5135 
5136 #endif // Tempus_StepperRKButcherTableau_hpp
StepperERK_3Stage3rdOrderTVD(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
StepperERK_Merson45(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Teuchos::RCP< StepperEDIRK_TrapezoidalRule< Scalar > > createStepperEDIRK_TrapezoidalRule(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
General Implicit Runge-Kutta Butcher Tableau.
StepperEDIRK_TrapezoidalRule(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
StepperERK_3Stage3rdOrderHeun(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
RK Implicit 2 Stage 2nd order Lobatto IIIB.
virtual std::string getDefaultICConsistency() const
std::string getStepperType() const
Get the stepper type. The stepper type is used as an identifier for the stepper, and can only be set ...
Teuchos::RCP< StepperSDIRK_SSPDIRK33< Scalar > > createStepperSDIRK_SSPDIRK33(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Explicit Runge-Kutta time stepper.
virtual void setupDefault()
Default setup for constructor.
Strong Stability Preserving Diagonally-Implicit RK Butcher Tableau.
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasicDIRK() const
StepperSDIRK_SSPDIRK33(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
StepperDIRK_BackwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Teuchos::RCP< StepperEDIRK_2Stage3rdOrder< Scalar > > createStepperEDIRK_2Stage3rdOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
StepperERK_SSPERK54(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Teuchos::RCP< StepperSDIRK_5Stage4thOrder< Scalar > > createStepperSDIRK_5Stage4thOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
StepperEDIRK_2StageTheta(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction, Scalar theta=Scalar(0.5))
Backward Euler Runge-Kutta Butcher Tableau.
StepperSDIRK_3Stage4thOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Teuchos::RCP< RKButcherTableau< Scalar > > tableau_
T & get(const std::string &name, T def_value)
Teuchos::RCP< StepperERK_SSPERK54< Scalar > > createStepperERK_SSPERK54(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
StepperSDIRK_ImplicitMidpoint(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
General Explicit Runge-Kutta Butcher Tableau.
StepperDIRK_General(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction, const Teuchos::SerialDenseMatrix< int, Scalar > &A, const Teuchos::SerialDenseVector< int, Scalar > &b, const Teuchos::SerialDenseVector< int, Scalar > &c, const int order, const int orderMin, const int orderMax, const Teuchos::SerialDenseVector< int, Scalar > &bstar)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Teuchos::RCP< StepperERK_3Stage3rdOrderHeun< Scalar > > createStepperERK_3Stage3rdOrderHeun(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
RK Explicit 5 Stage 3rd order by Kinnmark and Gray.
Teuchos::RCP< StepperERK_Merson45< Scalar > > createStepperERK_Merson45(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< StepperEDIRK_2StageTheta< Scalar > > createStepperEDIRK_2StageTheta(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
StepperSDIRK_SSPDIRK23(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
StepperSDIRK_3Stage2ndOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
void setStepperName(std::string s)
Set the stepper name.
Strong Stability Preserving Explicit RK Butcher Tableau.
RK Explicit 4 Stage 3rd order by Runge.
StepperERK_5Stage3rdOrderKandG(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
StepperEDIRK_2Stage3rdOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Teuchos::RCP< StepperERK_3_8Rule< Scalar > > createStepperERK_3_8Rule(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< StepperERK_BogackiShampine32< Scalar > > createStepperERK_BogackiShampine32(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< StepperSDIRK_5Stage5thOrder< Scalar > > createStepperSDIRK_5Stage5thOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< StepperERK_3Stage3rdOrder< Scalar > > createStepperERK_3Stage3rdOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
StepperDIRK_2Stage2ndOrderLobattoIIIB(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
RK Trapezoidal Rule (A.K.A. RK Crank-Nicolson)
Teuchos::RCP< StepperDIRK_2Stage2ndOrderLobattoIIIB< Scalar > > createStepperDIRK_2Stage2ndOrderLobattoIIIB(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< StepperERK_ForwardEuler< Scalar > > createStepperERK_ForwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< StepperSDIRK_21Pair< Scalar > > createStepperSDIRK_21Pair(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
bool isInitialized_
True if stepper&#39;s member data is initialized.
Teuchos::RCP< StepperSDIRK_SSPDIRK23< Scalar > > createStepperSDIRK_SSPDIRK23(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
StepperERK_4Stage3rdOrderRunge(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Teuchos::RCP< StepperERK_3Stage3rdOrderTVD< Scalar > > createStepperERK_3Stage3rdOrderTVD(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< StepperSDIRK_3Stage4thOrder< Scalar > > createStepperSDIRK_3Stage4thOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Explicit RK 3/8th Rule Butcher Tableau.
StepperSDIRK_2Stage3rdOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction, std::string gammaType="3rd Order A-stable", Scalar gamma=Scalar(0.7886751345948128))
bool isParameter(const std::string &name) const
StepperERK_Midpoint(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Diagonally Implicit Runge-Kutta (DIRK) time stepper.
StepperERK_4Stage4thOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Strong Stability Preserving Diagonally-Implicit RK Butcher Tableau.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Application Action for StepperRKBase.
Teuchos::RCP< StepperDIRK_1Stage1stOrderRadauIA< Scalar > > createStepperDIRK_1Stage1stOrderRadauIA(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
bool useFSAL_
Use First-Same-As-Last (FSAL) principle.
virtual void setupDefault()
Default setup for constructor.
Explicit RK Bogacki-Shampine Butcher Tableau.
StepperSDIRK_5Stage5thOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
void setICConsistencyCheck(bool c)
Teuchos::RCP< StepperSDIRK_3Stage2ndOrder< Scalar > > createStepperSDIRK_3Stage2ndOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< StepperERK_5Stage3rdOrderKandG< Scalar > > createStepperERK_5Stage3rdOrderKandG(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
OrdinalType length() const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
StepperSDIRK_21Pair(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Teuchos::RCP< StepperDIRK_BackwardEuler< Scalar > > createStepperDIRK_BackwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
StepperERK_Trapezoidal(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
StepperDIRK_1Stage1stOrderRadauIA(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Strong Stability Preserving Diagonally-Implicit RK Butcher Tableau.
Teuchos::RCP< StepperSDIRK_2Stage3rdOrder< Scalar > > createStepperSDIRK_2Stage3rdOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasicERK() const
StepperSDIRK_SSPDIRK22(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
void setTableau(const Teuchos::SerialDenseMatrix< int, Scalar > &A, const Teuchos::SerialDenseVector< int, Scalar > &b, const Teuchos::SerialDenseVector< int, Scalar > &c, const int order, const int orderMin, const int orderMax, const Teuchos::SerialDenseVector< int, Scalar > &bstar=Teuchos::SerialDenseVector< int, Scalar >())
Teuchos::RCP< StepperERK_Midpoint< Scalar > > createStepperERK_Midpoint(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
OrdinalType numCols() const
Teuchos::RCP< StepperSDIRK_SSPDIRK22< Scalar > > createStepperSDIRK_SSPDIRK22(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
StepperERK_ForwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Teuchos::RCP< StepperERK_Trapezoidal< Scalar > > createStepperERK_Trapezoidal(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Forward Euler Runge-Kutta Butcher Tableau.
StepperDIRK_1StageTheta(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction, Scalar theta=Scalar(0.5))
StepperERK_BogackiShampine32(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Teuchos::RCP< StepperSDIRK_ImplicitMidpoint< Scalar > > createStepperSDIRK_ImplicitMidpoint(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
StepperERK_3Stage3rdOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Teuchos::RCP< StepperERK_General< Scalar > > createStepperERK_General(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
virtual std::string getDescription() const
TypeTo as(const TypeFrom &t)
Teuchos::RCP< StepperSDIRK_2Stage2ndOrder< Scalar > > createStepperSDIRK_2Stage2ndOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Explicit RK Merson Butcher Tableau.
Runge-Kutta 4th order Butcher Tableau.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void setTableau(const Teuchos::SerialDenseMatrix< int, Scalar > &A, const Teuchos::SerialDenseVector< int, Scalar > &b, const Teuchos::SerialDenseVector< int, Scalar > &c, const int order, const int orderMin, const int orderMax, const Teuchos::SerialDenseVector< int, Scalar > &bstar=Teuchos::SerialDenseVector< int, Scalar >())
StepperERK_Ralston(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
virtual void setup(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &wrapperModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Setup for constructor.
StepperSDIRK_5Stage4thOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Teuchos::RCP< StepperERK_Ralston< Scalar > > createStepperERK_Ralston(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
virtual void setUseFSAL(bool a)
StepperSDIRK_2Stage2ndOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction, Scalar gamma=Scalar(0.2928932188134524))
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
StepperSDIRK_SSPDIRK32(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, const Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > &solver, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, bool zeroInitialGuess, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
StepperERK_3_8Rule(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
void setStepperType(std::string s)
Set the stepper type.
RK Explicit 3 Stage 3rd order by Heun.
Definition: