Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
23namespace Tempus {
24
25
26// ----------------------------------------------------------------------------
42template<class Scalar>
44 virtual public StepperExplicitRK<Scalar>
45{
46public:
53 {
54 this->setStepperName("RK Forward Euler");
55 this->setStepperType("RK Forward Euler");
56 this->setupTableau();
57 this->setupDefault();
58 this->setUseFSAL( true);
59 this->setICConsistency( "Consistent");
60 this->setICConsistencyCheck( false);
61 }
62
64 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
65 bool useFSAL,
66 std::string ICConsistency,
67 bool ICConsistencyCheck,
68 bool useEmbedded,
69 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
70 {
71 this->setStepperName("RK Forward Euler");
72 this->setStepperType("RK Forward Euler");
73 this->setupTableau();
74 this->setup(appModel, useFSAL, ICConsistency,
75 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
76 }
77
78 std::string getDescription() const
79 {
80 std::ostringstream Description;
81 Description << this->getStepperType() << "\n"
82 << "c = [ 0 ]'\n"
83 << "A = [ 0 ]\n"
84 << "b = [ 1 ]'";
85 return Description.str();
86 }
87
88 void setUseFSAL(bool a) { this->useFSAL_ = a; this->isInitialized_ = false; }
89
90protected:
91
93 {
94 typedef Teuchos::ScalarTraits<Scalar> ST;
95 Teuchos::SerialDenseMatrix<int,Scalar> A(1,1);
96 Teuchos::SerialDenseVector<int,Scalar> b(1);
97 Teuchos::SerialDenseVector<int,Scalar> c(1);
98 A(0,0) = ST::zero();
99 b(0) = ST::one();
100 c(0) = ST::zero();
101 int order = 1;
102
103 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
104 this->getStepperType(),A,b,c,order,order,order));
105 this->tableau_->setTVD(true);
106 this->tableau_->setTVDCoeff(1.0);
107 }
108};
109
110
112// ------------------------------------------------------------------------
113template<class Scalar>
114Teuchos::RCP<StepperERK_ForwardEuler<Scalar> >
116 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
117 Teuchos::RCP<Teuchos::ParameterList> pl)
118{
119 auto stepper = Teuchos::rcp(new StepperERK_ForwardEuler<Scalar>());
120
121 // Test for aliases.
122 if (pl != Teuchos::null) {
123 auto stepperType =
124 pl->get<std::string>("Stepper Type", stepper->getStepperType());
125
126 TEUCHOS_TEST_FOR_EXCEPTION(
127 stepperType != stepper->getStepperType() &&
128 stepperType != "RK1", std::logic_error,
129 " ParameterList 'Stepper Type' (='" + stepperType +"')\n"
130 " does not match type for this Stepper (='" + stepper->getStepperType() +
131 "')\n or one of its aliases ('RK1').\n");
132
133 // Reset default StepperType.
134 pl->set<std::string>("Stepper Type", stepper->getStepperType());
135 }
136
137 stepper->setStepperRKValues(pl);
138
139 if (model != Teuchos::null) {
140 stepper->setModel(model);
141 stepper->initialize();
142 }
143
144 return stepper;
145}
146
147
148// ----------------------------------------------------------------------------
167template<class Scalar>
169 virtual public StepperExplicitRK<Scalar>
170{
171public:
178 {
179 this->setStepperName("RK Explicit 4 Stage");
180 this->setStepperType("RK Explicit 4 Stage");
181 this->setupTableau();
182 this->setupDefault();
183 this->setUseFSAL( false);
184 this->setICConsistency( "None");
185 this->setICConsistencyCheck( false);
186 }
187
189 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
190 bool useFSAL,
191 std::string ICConsistency,
192 bool ICConsistencyCheck,
193 bool useEmbedded,
194 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
195 {
196 this->setStepperName("RK Explicit 4 Stage");
197 this->setStepperType("RK Explicit 4 Stage");
198 this->setupTableau();
199 this->setup(appModel, useFSAL, ICConsistency,
200 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
201 }
202
203 std::string getDescription() const
204 {
205 std::ostringstream Description;
206 Description << this->getStepperType() << "\n"
207 << "\"The\" Runge-Kutta Method (explicit):\n"
208 << "Solving Ordinary Differential Equations I:\n"
209 << "Nonstiff Problems, 2nd Revised Edition\n"
210 << "E. Hairer, S.P. Norsett, G. Wanner\n"
211 << "Table 1.2, pg 138\n"
212 << "c = [ 0 1/2 1/2 1 ]'\n"
213 << "A = [ 0 ] \n"
214 << " [ 1/2 0 ]\n"
215 << " [ 0 1/2 0 ]\n"
216 << " [ 0 0 1 0 ]\n"
217 << "b = [ 1/6 1/3 1/3 1/6 ]'";
218 return Description.str();
219 }
220
222 {
223 typedef Teuchos::ScalarTraits<Scalar> ST;
224 const Scalar one = ST::one();
225 const Scalar zero = ST::zero();
226 const Scalar onehalf = one/(2*one);
227 const Scalar onesixth = one/(6*one);
228 const Scalar onethird = one/(3*one);
229
230 int NumStages = 4;
231 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
232 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
233 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
234
235 // Fill A:
236 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero; A(0,3) = zero;
237 A(1,0) = onehalf; A(1,1) = zero; A(1,2) = zero; A(1,3) = zero;
238 A(2,0) = zero; A(2,1) = onehalf; A(2,2) = zero; A(2,3) = zero;
239 A(3,0) = zero; A(3,1) = zero; A(3,2) = one; A(3,3) = zero;
240
241 // Fill b:
242 b(0) = onesixth; b(1) = onethird; b(2) = onethird; b(3) = onesixth;
243
244 // fill c:
245 c(0) = zero; c(1) = onehalf; c(2) = onehalf; c(3) = one;
246
247 int order = 4;
248
249 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
250 this->getStepperType(),A,b,c,order,order,order));
251 }
252};
253
254
256// ------------------------------------------------------------------------
257template<class Scalar>
258Teuchos::RCP<StepperERK_4Stage4thOrder<Scalar> >
260 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
261 Teuchos::RCP<Teuchos::ParameterList> pl)
262{
263 auto stepper = Teuchos::rcp(new StepperERK_4Stage4thOrder<Scalar>());
264 stepper->setStepperRKValues(pl);
265
266 if (model != Teuchos::null) {
267 stepper->setModel(model);
268 stepper->initialize();
269 }
270
271 return stepper;
272}
273
274
275// ----------------------------------------------------------------------------
299template<class Scalar>
301 virtual public StepperExplicitRK<Scalar>
302{
303public:
310 {
311 this->setStepperName("Bogacki-Shampine 3(2) Pair");
312 this->setStepperType("Bogacki-Shampine 3(2) Pair");
313 this->setupTableau();
314 this->setupDefault();
315 this->setUseFSAL( true);
316 this->setICConsistency( "Consistent");
317 this->setICConsistencyCheck( false);
318 }
319
321 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
322 bool useFSAL,
323 std::string ICConsistency,
324 bool ICConsistencyCheck,
325 bool useEmbedded,
326 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
327 {
328 this->setStepperName("Bogacki-Shampine 3(2) Pair");
329 this->setStepperType("Bogacki-Shampine 3(2) Pair");
330 this->setupTableau();
331 this->setup(appModel, useFSAL, ICConsistency,
332 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
333 }
334
335 std::string getDescription() const
336 {
337 std::ostringstream Description;
338 Description << this->getStepperType() << "\n"
339 << "P. Bogacki and L.F. Shampine.\n"
340 << "A 3(2) pair of Runge–Kutta formulas.\n"
341 << "Applied Mathematics Letters, 2(4):321 – 325, 1989.\n"
342 << "c = [ 0 1/2 3/4 1 ]'\n"
343 << "A = [ 0 ]\n"
344 << " [ 1/2 0 ]\n"
345 << " [ 0 3/4 0 ]\n"
346 << " [ 2/9 1/3 4/9 0 ]\n"
347 << "b = [ 2/9 1/3 4/9 0 ]'\n"
348 << "bstar = [ 7/24 1/4 1/3 1/8 ]'";
349 return Description.str();
350 }
351
352 void setUseFSAL(bool a) { this->useFSAL_ = a; this->isInitialized_ = false; }
353
354protected:
355
357 {
358 typedef Teuchos::ScalarTraits<Scalar> ST;
359 using Teuchos::as;
360 int NumStages = 4;
361 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
362 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
363 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
364 Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);
365
366 const Scalar one = ST::one();
367 const Scalar zero = ST::zero();
368 const Scalar onehalf = one/(2*one);
369 const Scalar onethird = one/(3*one);
370 const Scalar threefourths = (3*one)/(4*one);
371 const Scalar twoninths = (2*one)/(9*one);
372 const Scalar fourninths = (4*one)/(9*one);
373
374 // Fill A:
375 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero; A(0,3) = zero;
376 A(1,0) = onehalf; A(1,1) = zero; A(1,2) = zero; A(1,3) = zero;
377 A(2,0) = zero; A(2,1) =threefourths; A(2,2) = zero; A(2,3) = zero;
378 A(3,0) =twoninths; A(3,1) = onethird; A(3,2) =fourninths; A(3,3) = zero;
379
380 // Fill b:
381 b(0) = A(3,0); b(1) = A(3,1); b(2) = A(3,2); b(3) = A(3,3);
382
383 // Fill c:
384 c(0) = zero; c(1) = onehalf; c(2) = threefourths; c(3) = one;
385
386 // Fill bstar
387 bstar(0) = as<Scalar>(7*one/(24*one));
388 bstar(1) = as<Scalar>(1*one/(4*one));
389 bstar(2) = as<Scalar>(1*one/(3*one));
390 bstar(3) = as<Scalar>(1*one/(8*one));
391 int order = 3;
392
393 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
394 this->getStepperType(),A,b,c,order,order,order,bstar));
395 }
396};
397
398
400// ------------------------------------------------------------------------
401template<class Scalar>
402Teuchos::RCP<StepperERK_BogackiShampine32<Scalar> >
404 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
405 Teuchos::RCP<Teuchos::ParameterList> pl)
406{
407 auto stepper = Teuchos::rcp(new StepperERK_BogackiShampine32<Scalar>());
408 stepper->setStepperRKValues(pl);
409
410 if (model != Teuchos::null) {
411 stepper->setModel(model);
412 stepper->initialize();
413 }
414
415 return stepper;
416}
417
418
419// ----------------------------------------------------------------------------
445template<class Scalar>
447 virtual public StepperExplicitRK<Scalar>
448{
449public:
456 {
457 this->setStepperName("Merson 4(5) Pair");
458 this->setStepperType("Merson 4(5) Pair");
459 this->setupTableau();
460 this->setupDefault();
461 this->setUseFSAL( false);
462 this->setICConsistency( "None");
463 this->setICConsistencyCheck( false);
464 }
465
467 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
468 bool useFSAL,
469 std::string ICConsistency,
470 bool ICConsistencyCheck,
471 bool useEmbedded,
472 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
473 {
474 this->setStepperName("Merson 4(5) Pair");
475 this->setStepperType("Merson 4(5) Pair");
476 this->setupTableau();
477 this->setup(appModel, useFSAL, ICConsistency,
478 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
479 }
480
481 std::string getDescription() const
482 {
483 std::ostringstream Description;
484 Description << this->getStepperType() << "\n"
485 << "Solving Ordinary Differential Equations I:\n"
486 << "Nonstiff Problems, 2nd Revised Edition\n"
487 << "E. Hairer, S.P. Norsett, G. Wanner\n"
488 << "Table 4.1, pg 167\n"
489 << "c = [ 0 1/3 1/3 1/2 1 ]'\n"
490 << "A = [ 0 ]\n"
491 << " [ 1/3 0 ]\n"
492 << " [ 1/6 1/6 0 ]\n"
493 << " [ 1/8 0 3/8 0 ]\n"
494 << " [ 1/2 0 -3/2 2 0 ]\n"
495 << "b = [ 1/6 0 0 2/3 1/6 ]'\n"
496 << "bstar = [ 1/10 0 3/10 2/5 1/5 ]'";
497 return Description.str();
498 }
499
500
501protected:
502
504 {
505 typedef Teuchos::ScalarTraits<Scalar> ST;
506 using Teuchos::as;
507 int NumStages = 5;
508 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages, true);
509 Teuchos::SerialDenseVector<int,Scalar> b(NumStages, true);
510 Teuchos::SerialDenseVector<int,Scalar> c(NumStages, true);
511 Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages, true);
512
513 const Scalar one = ST::one();
514 const Scalar zero = ST::zero();
515
516 // Fill A:
517 A(1,0) = as<Scalar>(one/(3*one));;
518
519 A(2,0) = as<Scalar>(one/(6*one));;
520 A(2,1) = as<Scalar>(one/(6*one));;
521
522 A(3,0) = as<Scalar>(one/(8*one));;
523 A(3,2) = as<Scalar>(3*one/(8*one));;
524
525 A(4,0) = as<Scalar>(one/(2*one));;
526 A(4,2) = as<Scalar>(-3*one/(2*one));;
527 A(4,3) = 2*one;
528
529 // Fill b:
530 b(0) = as<Scalar>(one/(6*one));
531 b(3) = as<Scalar>(2*one/(3*one));
532 b(4) = as<Scalar>(one/(6*one));
533
534 // Fill c:
535 c(0) = zero;
536 c(1) = as<Scalar>(1*one/(3*one));
537 c(2) = as<Scalar>(1*one/(3*one));
538 c(3) = as<Scalar>(1*one/(2*one));
539 c(4) = one;
540
541 // Fill bstar
542 bstar(0) = as<Scalar>(1*one/(10*one));
543 bstar(2) = as<Scalar>(3*one/(10*one));
544 bstar(3) = as<Scalar>(2*one/(5*one));
545 bstar(4) = as<Scalar>(1*one/(5*one));
546 int order = 4;
547
548 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
549 this->getStepperType(),A,b,c,order,order,order,bstar));
550 }
551};
552
553
555// ------------------------------------------------------------------------
556template<class Scalar>
557Teuchos::RCP<StepperERK_Merson45<Scalar> >
559 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
560 Teuchos::RCP<Teuchos::ParameterList> pl)
561{
562 auto stepper = Teuchos::rcp(new StepperERK_Merson45<Scalar>());
563 stepper->setStepperRKValues(pl);
564
565 if (model != Teuchos::null) {
566 stepper->setModel(model);
567 stepper->initialize();
568 }
569
570 return stepper;
571}
572
573
574// ----------------------------------------------------------------------------
597template<class Scalar>
599 virtual public StepperExplicitRK<Scalar>
600{
601public:
602
604 {
605 this->setStepperName("RK Explicit 3/8 Rule");
606 this->setStepperType("RK Explicit 3/8 Rule");
607 this->setupTableau();
608 this->setupDefault();
609 this->setUseFSAL( false);
610 this->setICConsistency( "None");
611 this->setICConsistencyCheck( false);
612 }
613
615 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
616 bool useFSAL,
617 std::string ICConsistency,
618 bool ICConsistencyCheck,
619 bool useEmbedded,
620 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
621 {
622 this->setStepperName("RK Explicit 3/8 Rule");
623 this->setStepperType("RK Explicit 3/8 Rule");
624 this->setupTableau();
625 this->setup(appModel, useFSAL, ICConsistency,
626 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
627 }
628
629 std::string getDescription() const
630 {
631 std::ostringstream Description;
632 Description << this->getStepperType() << "\n"
633 << "Solving Ordinary Differential Equations I:\n"
634 << "Nonstiff Problems, 2nd Revised Edition\n"
635 << "E. Hairer, S.P. Norsett, G. Wanner\n"
636 << "Table 1.2, pg 138\n"
637 << "c = [ 0 1/3 2/3 1 ]'\n"
638 << "A = [ 0 ]\n"
639 << " [ 1/3 0 ]\n"
640 << " [-1/3 1 0 ]\n"
641 << " [ 1 -1 1 0 ]\n"
642 << "b = [ 1/8 3/8 3/8 1/8 ]'";
643 return Description.str();
644 }
645
646
647protected:
648
650 {
651 typedef Teuchos::ScalarTraits<Scalar> ST;
652 using Teuchos::as;
653 int NumStages = 4;
654 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
655 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
656 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
657
658 const Scalar one = ST::one();
659 const Scalar zero = ST::zero();
660 const Scalar onethird = as<Scalar>(one/(3*one));
661 const Scalar twothirds = as<Scalar>(2*one/(3*one));
662 const Scalar oneeighth = as<Scalar>(one/(8*one));
663 const Scalar threeeighths = as<Scalar>(3*one/(8*one));
664
665 // Fill A:
666 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero; A(0,3) = zero;
667 A(1,0) = onethird; A(1,1) = zero; A(1,2) = zero; A(1,3) = zero;
668 A(2,0) = -onethird; A(2,1) = one; A(2,2) = zero; A(2,3) = zero;
669 A(3,0) = one; A(3,1) = -one; A(3,2) = one; A(3,3) = zero;
670
671 // Fill b:
672 b(0) =oneeighth; b(1) =threeeighths; b(2) =threeeighths; b(3) =oneeighth;
673
674 // Fill c:
675 c(0) = zero; c(1) = onethird; c(2) = twothirds; c(3) = one;
676
677 int order = 4;
678
679 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
680 this->getStepperType(),A,b,c,order,order,order));
681 }
682};
683
684
686// ------------------------------------------------------------------------
687template<class Scalar>
688Teuchos::RCP<StepperERK_3_8Rule<Scalar> >
690 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
691 Teuchos::RCP<Teuchos::ParameterList> pl)
692{
693 auto stepper = Teuchos::rcp(new StepperERK_3_8Rule<Scalar>());
694 stepper->setStepperRKValues(pl);
695
696 if (model != Teuchos::null) {
697 stepper->setModel(model);
698 stepper->initialize();
699 }
700
701 return stepper;
702}
703
704
705// ----------------------------------------------------------------------------
728template<class Scalar>
730 virtual public StepperExplicitRK<Scalar>
731{
732public:
739 {
740 this->setStepperName("RK Explicit 4 Stage 3rd order by Runge");
741 this->setStepperType("RK Explicit 4 Stage 3rd order by Runge");
742 this->setupTableau();
743 this->setupDefault();
744 this->setUseFSAL( false);
745 this->setICConsistency( "None");
746 this->setICConsistencyCheck( false);
747 }
748
750 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
751 bool useFSAL,
752 std::string ICConsistency,
753 bool ICConsistencyCheck,
754 bool useEmbedded,
755 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
756 {
757 this->setStepperName("RK Explicit 4 Stage 3rd order by Runge");
758 this->setStepperType("RK Explicit 4 Stage 3rd order by Runge");
759 this->setupTableau();
760 this->setup(appModel, useFSAL, ICConsistency,
761 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
762 }
763
764 std::string getDescription() const
765 {
766 std::ostringstream Description;
767 Description << this->getStepperType() << "\n"
768 << "Solving Ordinary Differential Equations I:\n"
769 << "Nonstiff Problems, 2nd Revised Edition\n"
770 << "E. Hairer, S.P. Norsett, G. Wanner\n"
771 << "Table 1.1, pg 135\n"
772 << "c = [ 0 1/2 1 1 ]'\n"
773 << "A = [ 0 ]\n"
774 << " [ 1/2 0 ]\n"
775 << " [ 0 1 0 ]\n"
776 << " [ 0 0 1 0 ]\n"
777 << "b = [ 1/6 2/3 0 1/6 ]'";
778 return Description.str();
779 }
780protected:
781
783 {
784 typedef Teuchos::ScalarTraits<Scalar> ST;
785 int NumStages = 4;
786 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
787 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
788 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
789
790 const Scalar one = ST::one();
791 const Scalar onehalf = one/(2*one);
792 const Scalar onesixth = one/(6*one);
793 const Scalar twothirds = 2*one/(3*one);
794 const Scalar zero = ST::zero();
795
796 // Fill A:
797 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero; A(0,3) = zero;
798 A(1,0) = onehalf; A(1,1) = zero; A(1,2) = zero; A(1,3) = zero;
799 A(2,0) = zero; A(2,1) = one; A(2,2) = zero; A(2,3) = zero;
800 A(3,0) = zero; A(3,1) = zero; A(3,2) = one; A(3,3) = zero;
801
802 // Fill b:
803 b(0) = onesixth; b(1) = twothirds; b(2) = zero; b(3) = onesixth;
804
805 // Fill c:
806 c(0) = zero; c(1) = onehalf; c(2) = one; c(3) = one;
807
808 int order = 3;
809
810 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
811 this->getStepperType(),A,b,c,order,order,order));
812 }
813};
814
815
817// ------------------------------------------------------------------------
818template<class Scalar>
819Teuchos::RCP<StepperERK_4Stage3rdOrderRunge<Scalar> >
821 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
822 Teuchos::RCP<Teuchos::ParameterList> pl)
823{
824 auto stepper = Teuchos::rcp(new StepperERK_4Stage3rdOrderRunge<Scalar>());
825 stepper->setStepperRKValues(pl);
826
827 if (model != Teuchos::null) {
828 stepper->setModel(model);
829 stepper->initialize();
830 }
831
832 return stepper;
833}
834
835
836// ----------------------------------------------------------------------------
858template<class Scalar>
860 virtual public StepperExplicitRK<Scalar>
861{
862public:
869 {
870 this->setStepperName("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
871 this->setStepperType("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
872 this->setupTableau();
873 this->setupDefault();
874 this->setUseFSAL( false);
875 this->setICConsistency( "None");
876 this->setICConsistencyCheck( false);
877 }
878
880 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
881 bool useFSAL,
882 std::string ICConsistency,
883 bool ICConsistencyCheck,
884 bool useEmbedded,
885 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
886 {
887 this->setStepperName("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
888 this->setStepperType("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
889 this->setupTableau();
890 this->setup(appModel, useFSAL, ICConsistency,
891 ICConsistencyCheck, useEmbedded,stepperRKAppAction);
892 }
893
894 std::string getDescription() const
895 {
896 std::ostringstream Description;
897 Description << this->getStepperType() << "\n"
898 << "Kinnmark & Gray 5 stage, 3rd order scheme \n"
899 << "Modified by P. Ullrich. From the prim_advance_mod.F90 \n"
900 << "routine in the HOMME atmosphere model code.\n"
901 << "c = [ 0 1/5 1/5 1/3 2/3 ]'\n"
902 << "A = [ 0 ]\n"
903 << " [ 1/5 0 ]\n"
904 << " [ 0 1/5 0 ]\n"
905 << " [ 0 0 1/3 0 ]\n"
906 << " [ 0 0 0 2/3 0 ]\n"
907 << "b = [ 1/4 0 0 0 3/4 ]'";
908 return Description.str();
909 }
910
911protected:
912
914 {
915 typedef Teuchos::ScalarTraits<Scalar> ST;
916 int NumStages = 5;
917 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
918 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
919 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
920
921 const Scalar one = ST::one();
922 const Scalar onefifth = one/(5*one);
923 const Scalar onefourth = one/(4*one);
924 const Scalar onethird = one/(3*one);
925 const Scalar twothirds = 2*one/(3*one);
926 const Scalar threefourths = 3*one/(4*one);
927 const Scalar zero = ST::zero();
928
929 // Fill A:
930 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero; A(0,3) = zero; A(0,4) = zero;
931 A(1,0) = onefifth; A(1,1) = zero; A(1,2) = zero; A(1,3) = zero; A(1,4) = zero;
932 A(2,0) = zero; A(2,1) = onefifth; A(2,2) = zero; A(2,3) = zero; A(2,4) = zero;
933 A(3,0) = zero; A(3,1) = zero; A(3,2) = onethird; A(3,3) = zero; A(3,4) = zero;
934 A(4,0) = zero; A(4,1) = zero; A(4,2) = zero; A(4,3) = twothirds; A(4,4) = zero;
935
936 // Fill b:
937 b(0) =onefourth; b(1) =zero; b(2) =zero; b(3) =zero; b(4) =threefourths;
938
939 // Fill c:
940 c(0) =zero; c(1) =onefifth; c(2) =onefifth; c(3) =onethird; c(4) =twothirds;
941
942 int order = 3;
943
944 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
945 this->getStepperType(),A,b,c,order,order,order));
946 }
947};
948
949
951// ------------------------------------------------------------------------
952template<class Scalar>
953Teuchos::RCP<StepperERK_5Stage3rdOrderKandG<Scalar> >
955 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
956 Teuchos::RCP<Teuchos::ParameterList> pl)
957{
958 auto stepper = Teuchos::rcp(new StepperERK_5Stage3rdOrderKandG<Scalar>());
959 stepper->setStepperRKValues(pl);
960
961 if (model != Teuchos::null) {
962 stepper->setModel(model);
963 stepper->initialize();
964 }
965
966 return stepper;
967}
968
969
970// ----------------------------------------------------------------------------
988template<class Scalar>
990 virtual public StepperExplicitRK<Scalar>
991{
992public:
999 {
1000 this->setStepperName("RK Explicit 3 Stage 3rd order");
1001 this->setStepperType("RK Explicit 3 Stage 3rd order");
1002 this->setupTableau();
1003 this->setupDefault();
1004 this->setUseFSAL( false);
1005 this->setICConsistency( "None");
1006 this->setICConsistencyCheck( false);
1007 }
1008
1010 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1011 bool useFSAL,
1012 std::string ICConsistency,
1013 bool ICConsistencyCheck,
1014 bool useEmbedded,
1015 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1016 {
1017 this->setStepperName("RK Explicit 3 Stage 3rd order");
1018 this->setStepperType("RK Explicit 3 Stage 3rd order");
1019 this->setupTableau();
1020 this->setup(appModel, useFSAL, ICConsistency,
1021 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
1022 }
1023
1024 std::string getDescription() const
1025 {
1026 std::ostringstream Description;
1027 Description << this->getStepperType() << "\n"
1028 << "c = [ 0 1/2 1 ]'\n"
1029 << "A = [ 0 ]\n"
1030 << " [ 1/2 0 ]\n"
1031 << " [ -1 2 0 ]\n"
1032 << "b = [ 1/6 4/6 1/6 ]'";
1033 return Description.str();
1034 }
1035
1036protected:
1037
1039 {
1040 typedef Teuchos::ScalarTraits<Scalar> ST;
1041 const Scalar one = ST::one();
1042 const Scalar two = Teuchos::as<Scalar>(2*one);
1043 const Scalar zero = ST::zero();
1044 const Scalar onehalf = one/(2*one);
1045 const Scalar onesixth = one/(6*one);
1046 const Scalar foursixth = 4*one/(6*one);
1047
1048 int NumStages = 3;
1049 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1050 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1051 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1052
1053 // Fill A:
1054 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero;
1055 A(1,0) = onehalf; A(1,1) = zero; A(1,2) = zero;
1056 A(2,0) = -one; A(2,1) = two; A(2,2) = zero;
1057
1058 // Fill b:
1059 b(0) = onesixth; b(1) = foursixth; b(2) = onesixth;
1060
1061 // fill c:
1062 c(0) = zero; c(1) = onehalf; c(2) = one;
1063
1064 int order = 3;
1065
1066 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
1067 this->getStepperType(),A,b,c,order,order,order));
1068 }
1069};
1070
1071
1073// ------------------------------------------------------------------------
1074template<class Scalar>
1075Teuchos::RCP<StepperERK_3Stage3rdOrder<Scalar> >
1077 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1078 Teuchos::RCP<Teuchos::ParameterList> pl)
1079{
1080 auto stepper = Teuchos::rcp(new StepperERK_3Stage3rdOrder<Scalar>());
1081 stepper->setStepperRKValues(pl);
1082
1083 if (model != Teuchos::null) {
1084 stepper->setModel(model);
1085 stepper->initialize();
1086 }
1087
1088 return stepper;
1089}
1090
1091
1092// ----------------------------------------------------------------------------
1121template<class Scalar>
1123 virtual public StepperExplicitRK<Scalar>
1124{
1125public:
1132 {
1133 this->setStepperName("RK Explicit 3 Stage 3rd order TVD");
1134 this->setStepperType("RK Explicit 3 Stage 3rd order TVD");
1135 this->setupTableau();
1136 this->setupDefault();
1137 this->setUseFSAL( false);
1138 this->setICConsistency( "None");
1139 this->setICConsistencyCheck( false);
1140 }
1141
1143 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1144 bool useFSAL,
1145 std::string ICConsistency,
1146 bool ICConsistencyCheck,
1147 bool useEmbedded,
1148 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1149 {
1150 this->setStepperName("RK Explicit 3 Stage 3rd order TVD");
1151 this->setStepperType("RK Explicit 3 Stage 3rd order TVD");
1152 this->setupTableau();
1153 this->setup(appModel, useFSAL, ICConsistency,
1154 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
1155 }
1156
1157 std::string getDescription() const
1158 {
1159 std::ostringstream Description;
1160 Description << this->getStepperType() << "\n"
1161 << "This Stepper is known as 'RK Explicit 3 Stage 3rd order TVD' or 'SSPERK33'.\n"
1162 << "Sigal Gottlieb and Chi-Wang Shu\n"
1163 << "`Total Variation Diminishing Runge-Kutta Schemes'\n"
1164 << "Mathematics of Computation\n"
1165 << "Volume 67, Number 221, January 1998, pp. 73-85\n"
1166 << "c = [ 0 1 1/2 ]'\n"
1167 << "A = [ 0 ]\n"
1168 << " [ 1 0 ]\n"
1169 << " [ 1/4 1/4 0 ]\n"
1170 << "b = [ 1/6 1/6 4/6 ]'\n"
1171 << "This is also written in the following set of updates.\n"
1172 << "u1 = u^n + dt L(u^n)\n"
1173 << "u2 = 3 u^n/4 + u1/4 + dt L(u1)/4\n"
1174 << "u^(n+1) = u^n/3 + 2 u2/2 + 2 dt L(u2)/3";
1175 return Description.str();
1176 }
1177
1178protected:
1179
1181 {
1182 typedef Teuchos::ScalarTraits<Scalar> ST;
1183 using Teuchos::as;
1184 const Scalar one = ST::one();
1185 const Scalar zero = ST::zero();
1186 const Scalar onehalf = one/(2*one);
1187 const Scalar onefourth = one/(4*one);
1188 const Scalar onesixth = one/(6*one);
1189 const Scalar foursixth = 4*one/(6*one);
1190
1191 int NumStages = 3;
1192 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1193 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1194 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1195 Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);
1196
1197 // Fill A:
1198 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero;
1199 A(1,0) = one; A(1,1) = zero; A(1,2) = zero;
1200 A(2,0) = onefourth; A(2,1) = onefourth; A(2,2) = zero;
1201
1202 // Fill b:
1203 b(0) = onesixth; b(1) = onesixth; b(2) = foursixth;
1204
1205 // fill c:
1206 c(0) = zero; c(1) = one; c(2) = onehalf;
1207
1208 // Fill bstar:
1209 bstar(0) = as<Scalar>(0.291485418878409);
1210 bstar(1) = as<Scalar>(0.291485418878409);
1211 bstar(2) = as<Scalar>(0.417029162243181);
1212
1213 int order = 3;
1214
1215 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
1216 this->getStepperType(),A,b,c,order,order,order,bstar));
1217 this->tableau_->setTVD(true);
1218 this->tableau_->setTVDCoeff(1.0);
1219 }
1220};
1221
1222
1224// ------------------------------------------------------------------------
1225template<class Scalar>
1226Teuchos::RCP<StepperERK_3Stage3rdOrderTVD<Scalar> >
1228 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1229 Teuchos::RCP<Teuchos::ParameterList> pl)
1230{
1231 auto stepper = Teuchos::rcp(new StepperERK_3Stage3rdOrderTVD<Scalar>());
1232
1233 // Test for aliases.
1234 if (pl != Teuchos::null) {
1235 auto stepperType =
1236 pl->get<std::string>("Stepper Type", stepper->getStepperType());
1237
1238 TEUCHOS_TEST_FOR_EXCEPTION(
1239 stepperType != stepper->getStepperType() &&
1240 stepperType != "SSPERK33" && stepperType != "SSPRK3", std::logic_error,
1241 " ParameterList 'Stepper Type' (='" + stepperType +"')\n"
1242 " does not match type for this Stepper (='" + stepper->getStepperType() +
1243 "')\n or one of its aliases ('SSPERK33' or 'SSPRK3').\n");
1244
1245 // Reset default StepperType.
1246 pl->set<std::string>("Stepper Type", stepper->getStepperType());
1247 }
1248
1249 stepper->setStepperRKValues(pl);
1250
1251 if (model != Teuchos::null) {
1252 stepper->setModel(model);
1253 stepper->initialize();
1254 }
1255
1256 return stepper;
1257}
1258
1259
1260// ----------------------------------------------------------------------------
1282template<class Scalar>
1284 virtual public StepperExplicitRK<Scalar>
1285{
1286public:
1293 {
1294 this->setStepperName("RK Explicit 3 Stage 3rd order by Heun");
1295 this->setStepperType("RK Explicit 3 Stage 3rd order by Heun");
1296 this->setupTableau();
1297 this->setupDefault();
1298 this->setUseFSAL( false);
1299 this->setICConsistency( "None");
1300 this->setICConsistencyCheck( false);
1301 }
1302
1304 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1305 bool useFSAL,
1306 std::string ICConsistency,
1307 bool ICConsistencyCheck,
1308 bool useEmbedded,
1309 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1310 {
1311 this->setStepperName("RK Explicit 3 Stage 3rd order by Heun");
1312 this->setStepperType("RK Explicit 3 Stage 3rd order by Heun");
1313 this->setupTableau();
1314 this->setup(appModel, useFSAL, ICConsistency,
1315 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
1316 }
1317
1318 std::string getDescription() const
1319 {
1320 std::ostringstream Description;
1321 Description << this->getStepperType() << "\n"
1322 << "Solving Ordinary Differential Equations I:\n"
1323 << "Nonstiff Problems, 2nd Revised Edition\n"
1324 << "E. Hairer, S.P. Norsett, G. Wanner\n"
1325 << "Table 1.1, pg 135\n"
1326 << "c = [ 0 1/3 2/3 ]'\n"
1327 << "A = [ 0 ] \n"
1328 << " [ 1/3 0 ]\n"
1329 << " [ 0 2/3 0 ]\n"
1330 << "b = [ 1/4 0 3/4 ]'";
1331 return Description.str();
1332 }
1333
1334protected:
1335
1337 {
1338 typedef Teuchos::ScalarTraits<Scalar> ST;
1339 const Scalar one = ST::one();
1340 const Scalar zero = ST::zero();
1341 const Scalar onethird = one/(3*one);
1342 const Scalar twothirds = 2*one/(3*one);
1343 const Scalar onefourth = one/(4*one);
1344 const Scalar threefourths = 3*one/(4*one);
1345
1346 int NumStages = 3;
1347 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1348 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1349 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1350
1351 // Fill A:
1352 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero;
1353 A(1,0) = onethird; A(1,1) = zero; A(1,2) = zero;
1354 A(2,0) = zero; A(2,1) = twothirds; A(2,2) = zero;
1355
1356 // Fill b:
1357 b(0) = onefourth; b(1) = zero; b(2) = threefourths;
1358
1359 // fill c:
1360 c(0) = zero; c(1) = onethird; c(2) = twothirds;
1361
1362 int order = 3;
1363
1364 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
1365 this->getStepperType(),A,b,c,order,order,order));
1366 }
1367};
1368
1369
1371// ------------------------------------------------------------------------
1372template<class Scalar>
1373Teuchos::RCP<StepperERK_3Stage3rdOrderHeun<Scalar> >
1375 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1376 Teuchos::RCP<Teuchos::ParameterList> pl)
1377{
1378 auto stepper = Teuchos::rcp(new StepperERK_3Stage3rdOrderHeun<Scalar>());
1379 stepper->setStepperRKValues(pl);
1380
1381 if (model != Teuchos::null) {
1382 stepper->setModel(model);
1383 stepper->initialize();
1384 }
1385
1386 return stepper;
1387}
1388
1389
1390// ----------------------------------------------------------------------------
1411template<class Scalar>
1413 virtual public StepperExplicitRK<Scalar>
1414{
1415public:
1422 {
1423 this->setStepperName("RK Explicit Midpoint");
1424 this->setStepperType("RK Explicit Midpoint");
1425 this->setupTableau();
1426 this->setupDefault();
1427 this->setUseFSAL( false);
1428 this->setICConsistency( "None");
1429 this->setICConsistencyCheck( false);
1430 }
1431
1433 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1434 bool useFSAL,
1435 std::string ICConsistency,
1436 bool ICConsistencyCheck,
1437 bool useEmbedded,
1438 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1439 {
1440 this->setStepperName("RK Explicit Midpoint");
1441 this->setStepperType("RK Explicit Midpoint");
1442 this->setupTableau();
1443 this->setup(appModel, useFSAL, ICConsistency,
1444 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
1445 }
1446
1447 std::string getDescription() const
1448 {
1449 std::ostringstream Description;
1450 Description << this->getStepperType() << "\n"
1451 << "Solving Ordinary Differential Equations I:\n"
1452 << "Nonstiff Problems, 2nd Revised Edition\n"
1453 << "E. Hairer, S.P. Norsett, G. Wanner\n"
1454 << "Table 1.1, pg 135\n"
1455 << "c = [ 0 1/2 ]'\n"
1456 << "A = [ 0 ]\n"
1457 << " [ 1/2 0 ]\n"
1458 << "b = [ 0 1 ]'";
1459 return Description.str();
1460 }
1461
1462protected:
1463
1465 {
1466 typedef Teuchos::ScalarTraits<Scalar> ST;
1467 const Scalar one = ST::one();
1468 const Scalar zero = ST::zero();
1469 const Scalar onehalf = one/(2*one);
1470
1471 int NumStages = 2;
1472 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1473 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1474 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1475
1476 // Fill A:
1477 A(0,0) = zero; A(0,1) = zero;
1478 A(1,0) = onehalf; A(1,1) = zero;
1479
1480 // Fill b:
1481 b(0) = zero; b(1) = one;
1482
1483 // fill c:
1484 c(0) = zero; c(1) = onehalf;
1485
1486 int order = 2;
1487
1488 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
1489 this->getStepperType(),A,b,c,order,order,order));
1490 }
1491};
1492
1493
1495// ------------------------------------------------------------------------
1496template<class Scalar>
1497Teuchos::RCP<StepperERK_Midpoint<Scalar> >
1499 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1500 Teuchos::RCP<Teuchos::ParameterList> pl)
1501{
1502 auto stepper = Teuchos::rcp(new StepperERK_Midpoint<Scalar>());
1503 stepper->setStepperRKValues(pl);
1504
1505 if (model != Teuchos::null) {
1506 stepper->setModel(model);
1507 stepper->initialize();
1508 }
1509
1510 return stepper;
1511}
1512
1513
1514// ----------------------------------------------------------------------------
1531template<class Scalar>
1533 virtual public StepperExplicitRK<Scalar>
1534{
1535public:
1542 {
1543 this->setStepperName("RK Explicit Ralston");
1544 this->setStepperType("RK Explicit Ralston");
1545 this->setupTableau();
1546 this->setupDefault();
1547 this->setUseFSAL( false);
1548 this->setICConsistency( "None");
1549 this->setICConsistencyCheck( false);
1550 }
1551
1553 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1554 bool useFSAL,
1555 std::string ICConsistency,
1556 bool ICConsistencyCheck,
1557 bool useEmbedded,
1558 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1559 {
1560 this->setStepperName("RK Explicit Ralston");
1561 this->setStepperType("RK Explicit Ralston");
1562 this->setupTableau();
1563 this->setup(appModel, useFSAL, ICConsistency,
1564 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
1565 }
1566
1567 std::string getDescription() const
1568 {
1569 std::ostringstream Description;
1570 Description << this->getStepperType() << "\n"
1571 << "This Stepper is known as 'RK Explicit Ralston' or 'RK2'.\n"
1572 << "c = [ 0 2/3 ]'\n"
1573 << "A = [ 0 ]\n"
1574 << " [ 2/3 0 ]\n"
1575 << "b = [ 1/4 3/4 ]'";
1576 return Description.str();
1577 }
1578
1579protected:
1580
1582 {
1583 typedef Teuchos::ScalarTraits<Scalar> ST;
1584 const Scalar one = ST::one();
1585 const Scalar zero = ST::zero();
1586
1587 const int NumStages = 2;
1588 const int order = 2;
1589 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1590 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1591 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1592
1593 // Fill A:
1594 A(0,0) = zero; A(0,1) = zero; A(1,1) = zero;
1595 A(1,0) = (2*one)/(3*one);
1596
1597 // Fill b:
1598 b(0) = (1*one)/(4*one);
1599 b(1) = (3*one)/(4*one);
1600
1601 // fill c:
1602 c(0) = zero;
1603 c(1) = (2*one)/(3*one);
1604
1605
1606 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
1607 this->getStepperType(),A,b,c,order,order,order));
1608 this->tableau_->setTVD(true);
1609 this->tableau_->setTVDCoeff(0.5);
1610 }
1611};
1612
1613
1615// ------------------------------------------------------------------------
1616template<class Scalar>
1617Teuchos::RCP<StepperERK_Ralston<Scalar> >
1619 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1620 Teuchos::RCP<Teuchos::ParameterList> pl)
1621{
1622 auto stepper = Teuchos::rcp(new StepperERK_Ralston<Scalar>());
1623
1624 // Test for aliases.
1625 if (pl != Teuchos::null) {
1626 auto stepperType =
1627 pl->get<std::string>("Stepper Type", stepper->getStepperType());
1628
1629 TEUCHOS_TEST_FOR_EXCEPTION(
1630 stepperType != stepper->getStepperType() &&
1631 stepperType != "RK2", std::logic_error,
1632 " ParameterList 'Stepper Type' (='" + stepperType +"')\n"
1633 " does not match type for this Stepper (='" + stepper->getStepperType() +
1634 "')\n or one of its aliases ('RK2').\n");
1635
1636 // Reset default StepperType.
1637 pl->set<std::string>("Stepper Type", stepper->getStepperType());
1638 }
1639
1640 stepper->setStepperRKValues(pl);
1641
1642 if (model != Teuchos::null) {
1643 stepper->setModel(model);
1644 stepper->initialize();
1645 }
1646
1647 return stepper;
1648}
1649
1650
1651// ----------------------------------------------------------------------------
1669template<class Scalar>
1671 virtual public StepperExplicitRK<Scalar>
1672{
1673public:
1680 {
1681 this->setStepperName("RK Explicit Trapezoidal");
1682 this->setStepperType("RK Explicit Trapezoidal");
1683 this->setupTableau();
1684 this->setupDefault();
1685 this->setUseFSAL( false);
1686 this->setICConsistency( "None");
1687 this->setICConsistencyCheck( false);
1688 }
1689
1691 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1692 bool useFSAL,
1693 std::string ICConsistency,
1694 bool ICConsistencyCheck,
1695 bool useEmbedded,
1696 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1697 {
1698 this->setStepperName("RK Explicit Trapezoidal");
1699 this->setStepperType("RK Explicit Trapezoidal");
1700 this->setupTableau();
1701 this->setup(appModel, useFSAL, ICConsistency,
1702 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
1703 }
1704
1705 std::string getDescription() const
1706 {
1707 std::ostringstream Description;
1708 Description << this->getStepperType() << "\n"
1709 << "This Stepper is known as 'RK Explicit Trapezoidal' or 'Heuns Method' or 'SSPERK22'.\n"
1710 << "c = [ 0 1 ]'\n"
1711 << "A = [ 0 ]\n"
1712 << " [ 1 0 ]\n"
1713 << "b = [ 1/2 1/2 ]\n"
1714 << "bstart = [ 3/4 1/4 ]'";
1715 return Description.str();
1716 }
1717
1718protected:
1719
1721 {
1722 typedef Teuchos::ScalarTraits<Scalar> ST;
1723 using Teuchos::as;
1724 const Scalar one = ST::one();
1725 const Scalar zero = ST::zero();
1726 const Scalar onehalf = one/(2*one);
1727
1728 int NumStages = 2;
1729 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1730 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1731 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1732 Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);
1733
1734 // Fill A:
1735 A(0,0) = zero; A(0,1) = zero;
1736 A(1,0) = one; A(1,1) = zero;
1737
1738 // Fill b:
1739 b(0) = onehalf; b(1) = onehalf;
1740
1741 // fill c:
1742 c(0) = zero; c(1) = one;
1743
1744 // Fill bstar
1745 bstar(0) = as<Scalar>(3*one/(4*one));
1746 bstar(1) = as<Scalar>(1*one/(4*one));
1747
1748 int order = 2;
1749
1750 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
1751 this->getStepperType(),A,b,c,order,order,order,bstar));
1752 this->tableau_->setTVD(true);
1753 this->tableau_->setTVDCoeff(1.0);
1754 }
1755};
1756
1757
1759// ------------------------------------------------------------------------
1760template<class Scalar>
1761Teuchos::RCP<StepperERK_Trapezoidal<Scalar> >
1763 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1764 Teuchos::RCP<Teuchos::ParameterList> pl)
1765{
1766 auto stepper = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
1767
1768 // Test for aliases.
1769 if (pl != Teuchos::null) {
1770 auto stepperType =
1771 pl->get<std::string>("Stepper Type", stepper->getStepperType());
1772
1773 TEUCHOS_TEST_FOR_EXCEPTION(
1774 stepperType != stepper->getStepperType() &&
1775 stepperType != "Heuns Method" && stepperType != "SSPERK22" &&
1776 stepperType != "SSPRK2", std::logic_error,
1777 " ParameterList 'Stepper Type' (='" + stepperType +"')\n"
1778 " does not match type for this Stepper (='" + stepper->getStepperType() +
1779 "')\n or one of its aliases ('Heuns Method' or 'SSPERK22' or 'SSPRK2').\n");
1780
1781 // Reset default StepperType.
1782 pl->set<std::string>("Stepper Type", stepper->getStepperType());
1783 }
1784
1785 stepper->setStepperRKValues(pl);
1786
1787 if (model != Teuchos::null) {
1788 stepper->setModel(model);
1789 stepper->initialize();
1790 }
1791
1792 return stepper;
1793}
1794
1795
1796// ----------------------------------------------------------------------------
1813template<class Scalar>
1815 virtual public StepperExplicitRK<Scalar>
1816{
1817 public:
1819 {
1820 this->setStepperName("SSPERK54");
1821 this->setStepperType("SSPERK54");
1822 this->setupTableau();
1823 this->setupDefault();
1824 this->setUseFSAL( false);
1825 this->setICConsistency( "None");
1826 this->setICConsistencyCheck( false);
1827 }
1828
1830 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1831 bool useFSAL,
1832 std::string ICConsistency,
1833 bool ICConsistencyCheck,
1834 bool useEmbedded,
1835 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
1836 {
1837 this->setStepperName("SSPERK54");
1838 this->setStepperType("SSPERK54");
1839 this->setupTableau();
1840 this->setup(appModel, useFSAL, ICConsistency,
1841 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
1842 }
1843
1844 std::string getDescription() const
1845 {
1846 std::ostringstream Description;
1847 Description << this->getStepperType() << "\n"
1848 << "Strong Stability Preserving Explicit RK (stage=5, order=4)\n"
1849 << std::endl;
1850 return Description.str();
1851 }
1852
1853protected:
1854
1856 {
1857
1858 typedef Teuchos::ScalarTraits<Scalar> ST;
1859 using Teuchos::as;
1860 const int NumStages = 5;
1861 const int order = 4;
1862 const Scalar sspcoef = 1.5082;
1863 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1864 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1865 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1866 Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);
1867 const Scalar zero = ST::zero();
1868
1869 // Fill A:
1870 A(0,0) = A(0,1) = A(0,2) = A(0,3) = A(0,4) = zero;
1871
1872 A(1,0) = as<Scalar>(0.391752226571889);
1873 A(1,1) = A(1,2) = A(1,3) = A(0,4) = zero;
1874
1875 A(2,0) = as<Scalar>(0.217669096261169);
1876 A(2,1) = as<Scalar>(0.368410593050372);
1877 A(2,2) = A(2,3) = A(2,4) = zero;
1878
1879 A(3,0) = as<Scalar>(0.082692086657811);
1880 A(3,1) = as<Scalar>(0.139958502191896);
1881 A(3,2) = as<Scalar>(0.251891774271693);
1882 A(3,3) = A(2,4) = zero;
1883
1884 A(4,0) = as<Scalar>(0.067966283637115);
1885 A(4,1) = as<Scalar>(0.115034698504632);
1886 A(4,2) = as<Scalar>(0.207034898597385);
1887 A(4,3) = as<Scalar>(0.544974750228520);
1888 A(4,4) = zero;
1889
1890 // Fill b:
1891 b(0) = as<Scalar>(0.146811876084786);
1892 b(1) = as<Scalar>(0.248482909444976);
1893 b(2) = as<Scalar>(0.104258830331980);
1894 b(3) = as<Scalar>(0.274438900901350);
1895 b(4) = as<Scalar>(0.226007483236908);
1896
1897 // fill c:
1898 c(0) = zero;
1899 c(1) = A(1,0);
1900 c(2) = A(2,0) + A(2,1);
1901 c(3) = A(3,0) + A(3,1) + A(3,2);
1902 c(4) = A(4,0) + A(4,1) + A(4,2) + A(4,3);
1903
1904 // Fill bstar:
1905 bstar(0) = as<Scalar>(0.130649104813131);
1906 bstar(1) = as<Scalar>(0.317716031201302);
1907 bstar(2) = as<Scalar>(0.000000869337261);
1908 bstar(3) = as<Scalar>(0.304581512634772);
1909 bstar(4) = as<Scalar>(0.247052482013534);
1910
1911 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
1912 this->getStepperType(),A,b,c,order,order,order,bstar));
1913 this->tableau_->setTVD(true);
1914 this->tableau_->setTVDCoeff(sspcoef);
1915 }
1916};
1917
1918
1920// ------------------------------------------------------------------------
1921template<class Scalar>
1922Teuchos::RCP<StepperERK_SSPERK54<Scalar> >
1924 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1925 Teuchos::RCP<Teuchos::ParameterList> pl)
1926{
1927 auto stepper = Teuchos::rcp(new StepperERK_SSPERK54<Scalar>());
1928 stepper->setStepperRKValues(pl);
1929
1930 if (model != Teuchos::null) {
1931 stepper->setModel(model);
1932 stepper->initialize();
1933 }
1934
1935 return stepper;
1936}
1937
1938
1939// ----------------------------------------------------------------------------
1969template<class Scalar>
1971 virtual public StepperExplicitRK<Scalar>
1972{
1973public:
1975 {
1976 this->setStepperName("General ERK");
1977 this->setStepperType("General ERK");
1978 this->setupTableau();
1979 this->setupDefault();
1980 this->setUseFSAL( false);
1981 this->setICConsistency( "None");
1982 this->setICConsistencyCheck( false);
1983 }
1984
1986 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
1987 bool useFSAL,
1988 std::string ICConsistency,
1989 bool ICConsistencyCheck,
1990 bool useEmbedded,
1991 const Teuchos::SerialDenseMatrix<int,Scalar>& A,
1992 const Teuchos::SerialDenseVector<int,Scalar>& b,
1993 const Teuchos::SerialDenseVector<int,Scalar>& c,
1994 const int order,
1995 const int orderMin,
1996 const int orderMax,
1997 const Teuchos::SerialDenseVector<int,Scalar>& bstar,
1998 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction=Teuchos::null)
1999 {
2000 this->setStepperName("General ERK");
2001 this->setStepperType("General ERK");
2002 this->setTableau(A,b,c,order,orderMin,orderMax,bstar);
2003
2004 TEUCHOS_TEST_FOR_EXCEPTION(
2005 this->tableau_->isImplicit() == true, std::logic_error,
2006 "Error - General ERK received an implicit Butcher Tableau!\n");
2007
2008 this->setup(appModel, useFSAL, ICConsistency,
2009 ICConsistencyCheck, useEmbedded, stepperRKAppAction);
2010 }
2011
2012 virtual std::string getDescription() const
2013 {
2014 std::stringstream Description;
2015 Description << this->getStepperType() << "\n"
2016 << "The format of the Butcher Tableau parameter list is\n"
2017 << " <Parameter name=\"A\" type=\"string\" value=\"# # # ;\n"
2018 << " # # # ;\n"
2019 << " # # #\"/>\n"
2020 << " <Parameter name=\"b\" type=\"string\" value=\"# # #\"/>\n"
2021 << " <Parameter name=\"c\" type=\"string\" value=\"# # #\"/>\n\n"
2022 << "Note the number of stages is implicit in the number of entries.\n"
2023 << "The number of stages must be consistent.\n"
2024 << "\n"
2025 << "Default tableau is RK4 (order=4):\n"
2026 << "c = [ 0 1/2 1/2 1 ]'\n"
2027 << "A = [ 0 ]\n"
2028 << " [ 1/2 0 ]\n"
2029 << " [ 0 1/2 0 ]\n"
2030 << " [ 0 0 1 0 ]\n"
2031 << "b = [ 1/6 1/3 1/3 1/6 ]'";
2032 return Description.str();
2033 }
2034
2036 {
2037 if (this->tableau_ == Teuchos::null) {
2038 // Set tableau to the default if null, otherwise keep current tableau.
2039 auto stepper = Teuchos::rcp(new StepperERK_4Stage4thOrder<Scalar>());
2040 auto t = stepper->getTableau();
2041 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
2042 this->getStepperType(),
2043 t->A(),t->b(),t->c(),
2044 t->order(),t->orderMin(),t->orderMax(),
2045 t->bstar()));
2046 }
2047 }
2048
2049 void setTableau(const Teuchos::SerialDenseMatrix<int,Scalar>& A,
2050 const Teuchos::SerialDenseVector<int,Scalar>& b,
2051 const Teuchos::SerialDenseVector<int,Scalar>& c,
2052 const int order,
2053 const int orderMin,
2054 const int orderMax,
2055 const Teuchos::SerialDenseVector<int,Scalar>&
2056 bstar = Teuchos::SerialDenseVector<int,Scalar>())
2057 {
2058 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
2059 this->getStepperType(),A,b,c,order,orderMin,orderMax,bstar));
2060 this->isInitialized_ = false;
2061 }
2062
2063 virtual std::string getDefaultICConsistency() const { return "Consistent"; }
2064
2065 Teuchos::RCP<const Teuchos::ParameterList>
2067 {
2068 auto pl = this->getValidParametersBasicERK();
2069
2070 // Tableau ParameterList
2071 Teuchos::SerialDenseMatrix<int,Scalar> A = this->tableau_->A();
2072 Teuchos::SerialDenseVector<int,Scalar> b = this->tableau_->b();
2073 Teuchos::SerialDenseVector<int,Scalar> c = this->tableau_->c();
2074 Teuchos::SerialDenseVector<int,Scalar> bstar = this->tableau_->bstar();
2075
2076 Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
2077
2078 std::ostringstream Astream;
2079 Astream.precision(15);
2080 for(int i = 0; i < A.numRows(); i++) {
2081 for(int j = 0; j < A.numCols()-1; j++) {
2082 Astream << A(i,j) << " ";
2083 }
2084 Astream << A(i,A.numCols()-1);
2085 if ( i != A.numRows()-1 ) Astream << "; ";
2086 }
2087 tableauPL->set<std::string>("A", Astream.str());
2088
2089 std::ostringstream bstream;
2090 bstream.precision(15);
2091 for(int i = 0; i < b.length()-1; i++) {
2092 bstream << b(i) << " ";
2093 }
2094 bstream << b(b.length()-1);
2095 tableauPL->set<std::string>("b", bstream.str());
2096
2097 std::ostringstream cstream;
2098 cstream.precision(15);
2099 for(int i = 0; i < c.length()-1; i++) {
2100 cstream << c(i) << " ";
2101 }
2102 cstream << c(c.length()-1);
2103 tableauPL->set<std::string>("c", cstream.str());
2104
2105 tableauPL->set<int>("order", this->tableau_->order());
2106
2107 if ( bstar.length() == 0 ) {
2108 tableauPL->set("bstar", "");
2109 } else {
2110 std::ostringstream bstarstream;
2111 bstarstream.precision(15);
2112 for(int i = 0; i < bstar.length()-1; i++) {
2113 bstarstream << bstar(i) << " ";
2114 }
2115 bstarstream << bstar(bstar.length()-1);
2116 tableauPL->set<std::string>("bstar", bstarstream.str());
2117 }
2118
2119 pl->set("Tableau", *tableauPL);
2120
2121 return pl;
2122 }
2123};
2124
2125
2127// ------------------------------------------------------------------------
2128template<class Scalar>
2129Teuchos::RCP<StepperERK_General<Scalar> >
2131 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2132 Teuchos::RCP<Teuchos::ParameterList> pl)
2133{
2134 auto stepper = Teuchos::rcp(new StepperERK_General<Scalar>());
2135 stepper->setStepperRKValues(pl);
2136
2137 if (pl != Teuchos::null) {
2138 if (pl->isParameter("Tableau")) {
2139 auto t = stepper->createTableau(pl);
2140 stepper->setTableau( t->A(),t->b(),t->c(),
2141 t->order(),t->orderMin(),t->orderMax(),
2142 t->bstar() );
2143 }
2144 }
2145 TEUCHOS_TEST_FOR_EXCEPTION(
2146 stepper->getTableau()->isImplicit() == true, std::logic_error,
2147 "Error - General ERK received an implicit Butcher Tableau!\n");
2148
2149 if (model != Teuchos::null) {
2150 stepper->setModel(model);
2151 stepper->initialize();
2152 }
2153
2154 return stepper;
2155}
2156
2157
2158// ----------------------------------------------------------------------------
2174template<class Scalar>
2176 virtual public StepperDIRK<Scalar>
2177{
2178public:
2185 {
2186 this->setStepperName("RK Backward Euler");
2187 this->setStepperType("RK Backward Euler");
2188 this->setupTableau();
2189 this->setupDefault();
2190 this->setUseFSAL( false);
2191 this->setICConsistency( "None");
2192 this->setICConsistencyCheck( false);
2193 }
2194
2196 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2197 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
2198 bool useFSAL,
2199 std::string ICConsistency,
2200 bool ICConsistencyCheck,
2201 bool useEmbedded,
2202 bool zeroInitialGuess,
2203 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
2204 {
2205 this->setStepperName("RK Backward Euler");
2206 this->setStepperType("RK Backward Euler");
2207 this->setupTableau();
2208 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
2209 useEmbedded, zeroInitialGuess, stepperRKAppAction);
2210 }
2211
2212 std::string getDescription() const
2213 {
2214 std::ostringstream Description;
2215 Description << this->getStepperType() << "\n"
2216 << "c = [ 1 ]'\n"
2217 << "A = [ 1 ]\n"
2218 << "b = [ 1 ]'";
2219 return Description.str();
2220 }
2221
2222protected:
2223
2225 {
2226 typedef Teuchos::ScalarTraits<Scalar> ST;
2227 const Scalar sspcoef = std::numeric_limits<Scalar>::max();
2228 int NumStages = 1;
2229 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2230 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2231 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2232
2233 // Fill A:
2234 A(0,0) = ST::one();
2235
2236 // Fill b:
2237 b(0) = ST::one();
2238
2239 // Fill c:
2240 c(0) = ST::one();
2241
2242 int order = 1;
2243
2244 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
2245 this->getStepperType(),A,b,c,order,order,order));
2246 this->tableau_->setTVD(true);
2247 this->tableau_->setTVDCoeff(sspcoef);
2248 }
2249};
2250
2251
2253// ------------------------------------------------------------------------
2254template<class Scalar>
2255Teuchos::RCP<StepperDIRK_BackwardEuler<Scalar> >
2257 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2258 Teuchos::RCP<Teuchos::ParameterList> pl)
2259{
2260 auto stepper = Teuchos::rcp(new StepperDIRK_BackwardEuler<Scalar>());
2261 stepper->setStepperDIRKValues(pl);
2262
2263 if (model != Teuchos::null) {
2264 stepper->setModel(model);
2265 stepper->initialize();
2266 }
2267
2268 return stepper;
2269}
2270
2271
2272// ----------------------------------------------------------------------------
2297template<class Scalar>
2299 virtual public StepperDIRK<Scalar>
2300{
2301public:
2308 {
2309 typedef Teuchos::ScalarTraits<Scalar> ST;
2310 const Scalar one = ST::one();
2311 gammaDefault_ = Teuchos::as<Scalar>((2*one-ST::squareroot(2*one))/(2*one));
2312
2313 this->setStepperName("SDIRK 2 Stage 2nd order");
2314 this->setStepperType("SDIRK 2 Stage 2nd order");
2315 this->setGamma(gammaDefault_);
2316 this->setupTableau();
2317 this->setupDefault();
2318 this->setUseFSAL( false);
2319 this->setICConsistency( "None");
2320 this->setICConsistencyCheck( false);
2321 }
2322
2324 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2325 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
2326 bool useFSAL,
2327 std::string ICConsistency,
2328 bool ICConsistencyCheck,
2329 bool useEmbedded,
2330 bool zeroInitialGuess,
2331 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
2332 Scalar gamma = Scalar(0.2928932188134524))
2333 {
2334 typedef Teuchos::ScalarTraits<Scalar> ST;
2335 const Scalar one = ST::one();
2336 gammaDefault_ = Teuchos::as<Scalar>((2*one-ST::squareroot(2*one))/(2*one));
2337
2338 this->setStepperName("SDIRK 2 Stage 2nd order");
2339 this->setStepperType("SDIRK 2 Stage 2nd order");
2340 this->setGamma(gamma);
2341 this->setupTableau();
2342 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
2343 useEmbedded, zeroInitialGuess, stepperRKAppAction);
2344 }
2345
2346 void setGamma(Scalar gamma)
2347 {
2348 gamma_ = gamma;
2349 this->isInitialized_ = false;
2350 this->setupTableau();
2351 }
2352
2353 Scalar getGamma() const { return gamma_; }
2354
2355 std::string getDescription() const
2356 {
2357 std::ostringstream Description;
2358 Description << this->getStepperType() << "\n"
2359 << "Computer Methods for ODEs and DAEs\n"
2360 << "U. M. Ascher and L. R. Petzold\n"
2361 << "p. 106\n"
2362 << "gamma = (2+-sqrt(2))/2\n"
2363 << "c = [ gamma 1 ]'\n"
2364 << "A = [ gamma 0 ]\n"
2365 << " [ 1-gamma gamma ]\n"
2366 << "b = [ 1-gamma gamma ]'";
2367 return Description.str();
2368 }
2369
2370 Teuchos::RCP<const Teuchos::ParameterList>
2372 {
2373 auto pl = this->getValidParametersBasicDIRK();
2374 pl->template set<double>("gamma", this->getGamma(),
2375 "The default value is gamma = (2-sqrt(2))/2. "
2376 "This will produce an L-stable 2nd order method with the stage "
2377 "times within the timestep. Other values of gamma will still "
2378 "produce an L-stable scheme, but will only be 1st order accurate.");
2379
2380 return pl;
2381 }
2382
2383protected:
2384
2386 {
2387 typedef Teuchos::ScalarTraits<Scalar> ST;
2388 int NumStages = 2;
2389 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2390 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2391 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2392
2393 const Scalar one = ST::one();
2394 const Scalar zero = ST::zero();
2395
2396 // Fill A:
2397 A(0,0) = gamma_; A(0,1) = zero;
2398 A(1,0) = Teuchos::as<Scalar>( one - gamma_ ); A(1,1) = gamma_;
2399
2400 // Fill b:
2401 b(0) = Teuchos::as<Scalar>( one - gamma_ ); b(1) = gamma_;
2402
2403 // Fill c:
2404 c(0) = gamma_; c(1) = one;
2405
2406 int order = 1;
2407 if ( std::abs((gamma_-gammaDefault_)/gamma_) < 1.0e-08 ) order = 2;
2408
2409 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
2410 this->getStepperType(),A,b,c,order,order,order));
2411 }
2412
2413 private:
2415 Scalar gamma_;
2416};
2417
2418
2420// ------------------------------------------------------------------------
2421template<class Scalar>
2422Teuchos::RCP<StepperSDIRK_2Stage2ndOrder<Scalar> >
2424 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2425 Teuchos::RCP<Teuchos::ParameterList> pl)
2426{
2427 auto stepper = Teuchos::rcp(new StepperSDIRK_2Stage2ndOrder<Scalar>());
2428 stepper->setStepperDIRKValues(pl);
2429
2430 if (pl != Teuchos::null)
2431 stepper->setGamma(pl->get<double>("gamma", 0.2928932188134524));
2432
2433 if (model != Teuchos::null) {
2434 stepper->setModel(model);
2435 stepper->initialize();
2436 }
2437
2438 return stepper;
2439}
2440
2441
2442// ----------------------------------------------------------------------------
2470template<class Scalar>
2472 virtual public StepperDIRK<Scalar>
2473{
2474public:
2481 {
2482 this->setStepperName("SDIRK 3 Stage 2nd order");
2483 this->setStepperType("SDIRK 3 Stage 2nd order");
2484 this->setupTableau();
2485 this->setupDefault();
2486 this->setUseFSAL( false);
2487 this->setICConsistency( "None");
2488 this->setICConsistencyCheck( false);
2489 }
2490
2492 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2493 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
2494 bool useFSAL,
2495 std::string ICConsistency,
2496 bool ICConsistencyCheck,
2497 bool useEmbedded,
2498 bool zeroInitialGuess,
2499 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
2500 {
2501
2502 this->setStepperName("SDIRK 3 Stage 2nd order");
2503 this->setStepperType("SDIRK 3 Stage 2nd order");
2504 this->setupTableau();
2505 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
2506 useEmbedded, zeroInitialGuess, stepperRKAppAction);
2507 }
2508
2509 std::string getDescription() const
2510 {
2511 std::ostringstream Description;
2512 Description << this->getStepperType() << "\n"
2513 << "Implicit-explicit Runge-Kutta schemes and applications to\n"
2514 << "hyperbolic systems with relaxation\n"
2515 << "L Pareschi, G Russo\n"
2516 << "Journal of Scientific computing, 2005 - Springer\n"
2517 << "Table 5\n"
2518 << "gamma = 1/(2+sqrt(2))\n"
2519 << "c = [ gamma (1-gamma) 1/2 ]'\n"
2520 << "A = [ gamma 0 0 ]\n"
2521 << " [ 1-2gamma gamma 0 ]\n"
2522 << " [ 1/2-gamma 0 gamma ]\n"
2523 << "b = [ 1/6 1/6 2/3 ]'";
2524 return Description.str();
2525 }
2526
2527protected:
2528
2530 {
2531 typedef Teuchos::ScalarTraits<Scalar> ST;
2532 using Teuchos::as;
2533 const int NumStages = 3;
2534 const int order = 2;
2535 const Scalar sspcoef = 1.0529;
2536 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2537 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2538 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2539 const Scalar one = ST::one();
2540 const Scalar zero = ST::zero();
2541 const Scalar gamma = as<Scalar>(one - ( one / ST::squareroot(2*one) ) );
2542
2543 // Fill A:
2544 A(0,0) = A(1,1) = A(2,2) = gamma;
2545 A(0,1) = A(0,2) = A(1,2) = A(2,1) = zero;
2546 A(1,0) = as<Scalar>(one - 2*gamma);
2547 A(2,0) = as<Scalar>( ( one/ (2.*one)) - gamma );
2548
2549 // Fill b:
2550 b(0) = b(1) = ( one / (6*one) );
2551 b(2) = (2*one)/(3*one);
2552
2553 // Fill c:
2554 c(0) = gamma;
2555 c(1) = one - gamma;
2556 c(2) = one / (2*one);
2557
2558 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
2559 this->getStepperType(),A,b,c,order,order,order));
2560 this->tableau_->setTVD(true);
2561 this->tableau_->setTVDCoeff(sspcoef);
2562 }
2563
2564};
2565
2566
2568// ------------------------------------------------------------------------
2569template<class Scalar>
2570Teuchos::RCP<StepperSDIRK_3Stage2ndOrder<Scalar> >
2572 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2573 Teuchos::RCP<Teuchos::ParameterList> pl)
2574{
2575 auto stepper = Teuchos::rcp(new StepperSDIRK_3Stage2ndOrder<Scalar>());
2576 stepper->setStepperDIRKValues(pl);
2577
2578 if (model != Teuchos::null) {
2579 stepper->setModel(model);
2580 stepper->initialize();
2581 }
2582
2583 return stepper;
2584}
2585
2586
2587// ----------------------------------------------------------------------------
2616template<class Scalar>
2618 virtual public StepperDIRK<Scalar>
2619{
2620public:
2627 {
2628 typedef Teuchos::ScalarTraits<Scalar> ST;
2629 using Teuchos::as;
2630 const Scalar one = ST::one();
2631 gammaDefault_ = as<Scalar>((3*one+ST::squareroot(3*one))/(6*one));
2632 gammaTypeDefault_ = "3rd Order A-stable";
2633
2634 this->setStepperName("SDIRK 2 Stage 3rd order");
2635 this->setStepperType("SDIRK 2 Stage 3rd order");
2637 this->setGamma(gammaDefault_);
2638 this->setupTableau();
2639 this->setupDefault();
2640 this->setUseFSAL( false);
2641 this->setICConsistency( "None");
2642 this->setICConsistencyCheck( false);
2643 }
2644
2646 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2647 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
2648 bool useFSAL,
2649 std::string ICConsistency,
2650 bool ICConsistencyCheck,
2651 bool useEmbedded,
2652 bool zeroInitialGuess,
2653 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
2654 std::string gammaType = "3rd Order A-stable",
2655 Scalar gamma = Scalar(0.7886751345948128))
2656 {
2657 typedef Teuchos::ScalarTraits<Scalar> ST;
2658 using Teuchos::as;
2659 const Scalar one = ST::one();
2660 gammaDefault_ = as<Scalar>((3*one+ST::squareroot(3*one))/(6*one));
2661 gammaTypeDefault_ = "3rd Order A-stable";
2662
2663 this->setStepperName("SDIRK 2 Stage 3rd order");
2664 this->setStepperType("SDIRK 2 Stage 3rd order");
2665 this->setGammaType(gammaType);
2666 this->setGamma(gamma);
2667 this->setupTableau();
2668 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
2669 useEmbedded, zeroInitialGuess, stepperRKAppAction);
2670 }
2671
2672 void setGammaType(std::string gammaType)
2673 {
2674 TEUCHOS_TEST_FOR_EXCEPTION(
2675 !(gammaType == "3rd Order A-stable" ||
2676 gammaType == "2nd Order L-stable" ||
2677 gammaType == "gamma"), std::logic_error,
2678 "gammaType needs to be '3rd Order A-stable', '2nd Order L-stable' "
2679 "or 'gamma'.");
2680
2681 gammaType_ = gammaType;
2682 this->isInitialized_ = false;
2683 this->setupTableau();
2684 }
2685
2686 std::string getGammaType() const { return gammaType_; }
2687
2688 void setGamma(Scalar gamma)
2689 {
2690 if ( gammaType_ == "gamma" ) {
2691 gamma_ = gamma;
2692 this->setupTableau();
2693 }
2694 this->isInitialized_ = false;
2695 }
2696
2697 Scalar getGamma() const { return gamma_; }
2698
2699 std::string getDescription() const
2700 {
2701 std::ostringstream Description;
2702 Description << this->getStepperType() << "\n"
2703 << "Solving Ordinary Differential Equations I:\n"
2704 << "Nonstiff Problems, 2nd Revised Edition\n"
2705 << "E. Hairer, S. P. Norsett, and G. Wanner\n"
2706 << "Table 7.2, pg 207\n"
2707 << "gamma = (3+sqrt(3))/6 -> 3rd order and A-stable\n"
2708 << "gamma = (2-sqrt(2))/2 -> 2nd order and L-stable\n"
2709 << "c = [ gamma 1-gamma ]'\n"
2710 << "A = [ gamma 0 ]\n"
2711 << " [ 1-2*gamma gamma ]\n"
2712 << "b = [ 1/2 1/2 ]'";
2713 return Description.str();
2714 }
2715
2716 Teuchos::RCP<const Teuchos::ParameterList>
2718 {
2719 auto pl = this->getValidParametersBasicDIRK();
2720
2721 pl->template set<std::string>("Gamma Type", this->getGammaType(),
2722 "Valid values are '3rd Order A-stable' ((3+sqrt(3))/6.), "
2723 "'2nd Order L-stable' ((2-sqrt(2))/2) and 'gamma' (user defined). "
2724 "The default value is '3rd Order A-stable'.");
2725 pl->template set<double>("gamma", this->getGamma(),
2726 "Equal to (3+sqrt(3))/6 if 'Gamma Type' = '3rd Order A-stable', or "
2727 "(2-sqrt(2))/2 if 'Gamma Type' = '2nd Order L-stable', or "
2728 "user-defined gamma value if 'Gamma Type = 'gamma'. "
2729 "The default value is gamma = (3+sqrt(3))/6, which matches "
2730 "the default 'Gamma Type' = '3rd Order A-stable'.");
2731
2732 return pl;
2733 }
2734
2735protected:
2736
2738 {
2739 typedef Teuchos::ScalarTraits<Scalar> ST;
2740 using Teuchos::as;
2741 int NumStages = 2;
2742 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2743 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2744 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2745 const Scalar one = ST::one();
2746 const Scalar zero = ST::zero();
2747
2748 int order = 0;
2749 if (gammaType_ == "3rd Order A-stable") {
2750 order = 3;
2752 } else if (gammaType_ == "2nd Order L-stable") {
2753 order = 2;
2754 gamma_ = as<Scalar>( (2*one - ST::squareroot(2*one))/(2*one) );
2755 } else if (gammaType_ == "gamma") {
2756 order = 2;
2757 }
2758
2759 // Fill A:
2760 A(0,0) = gamma_; A(0,1) = zero;
2761 A(1,0) = as<Scalar>(one - 2*gamma_); A(1,1) = gamma_;
2762
2763 // Fill b:
2764 b(0) = as<Scalar>( one/(2*one) ); b(1) = as<Scalar>( one/(2*one) );
2765
2766 // Fill c:
2767 c(0) = gamma_; c(1) = as<Scalar>( one - gamma_ );
2768
2769 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
2770 this->getStepperType(),A,b,c,order,2,3));
2771 }
2772
2773 private:
2775 std::string gammaType_;
2777 Scalar gamma_;
2778};
2779
2780
2782// ------------------------------------------------------------------------
2783template<class Scalar>
2784Teuchos::RCP<StepperSDIRK_2Stage3rdOrder<Scalar> >
2786 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2787 Teuchos::RCP<Teuchos::ParameterList> pl)
2788{
2789 auto stepper = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
2790 stepper->setStepperDIRKValues(pl);
2791
2792 if (pl != Teuchos::null) {
2793 stepper->setGammaType(
2794 pl->get<std::string>("Gamma Type", "3rd Order A-stable"));
2795 stepper->setGamma(pl->get<double>("gamma", 0.7886751345948128));
2796 }
2797
2798 if (model != Teuchos::null) {
2799 stepper->setModel(model);
2800 stepper->initialize();
2801 }
2802
2803 return stepper;
2804}
2805
2806
2807// ----------------------------------------------------------------------------
2828template<class Scalar>
2830 virtual public StepperDIRK<Scalar>
2831{
2832public:
2839 {
2840 this->setStepperName("EDIRK 2 Stage 3rd order");
2841 this->setStepperType("EDIRK 2 Stage 3rd order");
2842 this->setupTableau();
2843 this->setupDefault();
2844 this->setUseFSAL( false);
2845 this->setICConsistency( "None");
2846 this->setICConsistencyCheck( false);
2847 }
2848
2850 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2851 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
2852 bool useFSAL,
2853 std::string ICConsistency,
2854 bool ICConsistencyCheck,
2855 bool useEmbedded,
2856 bool zeroInitialGuess,
2857 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
2858 {
2859 this->setStepperName("EDIRK 2 Stage 3rd order");
2860 this->setStepperType("EDIRK 2 Stage 3rd order");
2861 this->setupTableau();
2862 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
2863 useEmbedded, zeroInitialGuess, stepperRKAppAction);
2864 }
2865
2866 std::string getDescription() const
2867 {
2868 std::ostringstream Description;
2869 Description << this->getStepperType() << "\n"
2870 << "Hammer & Hollingsworth method\n"
2871 << "Solving Ordinary Differential Equations I:\n"
2872 << "Nonstiff Problems, 2nd Revised Edition\n"
2873 << "E. Hairer, S. P. Norsett, and G. Wanner\n"
2874 << "Table 7.1, pg 205\n"
2875 << "c = [ 0 2/3 ]'\n"
2876 << "A = [ 0 0 ]\n"
2877 << " [ 1/3 1/3 ]\n"
2878 << "b = [ 1/4 3/4 ]'";
2879 return Description.str();
2880 }
2881
2882protected:
2883
2885 {
2886 typedef Teuchos::ScalarTraits<Scalar> ST;
2887 using Teuchos::as;
2888 int NumStages = 2;
2889 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2890 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2891 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2892 const Scalar one = ST::one();
2893 const Scalar zero = ST::zero();
2894
2895 // Fill A:
2896 A(0,0) = zero; A(0,1) = zero;
2897 A(1,0) = as<Scalar>( one/(3*one) ); A(1,1) = as<Scalar>( one/(3*one) );
2898
2899 // Fill b:
2900 b(0) = as<Scalar>( one/(4*one) ); b(1) = as<Scalar>( 3*one/(4*one) );
2901
2902 // Fill c:
2903 c(0) = zero; c(1) = as<Scalar>( 2*one/(3*one) );
2904 int order = 3;
2905
2906 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
2907 this->getStepperType(),A,b,c,order,order,order));
2908 this->tableau_->setTVD(true);
2909 this->tableau_->setTVDCoeff(1.5);
2910 }
2911};
2912
2913
2915template<class Scalar>
2916Teuchos::RCP<StepperEDIRK_2Stage3rdOrder<Scalar> >
2918 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
2919 Teuchos::RCP<Teuchos::ParameterList> pl)
2920{
2921 auto stepper = Teuchos::rcp(new StepperEDIRK_2Stage3rdOrder<Scalar>());
2922 stepper->setStepperDIRKValues(pl);
2923
2924 if (model != Teuchos::null) {
2925 stepper->setModel(model);
2926 stepper->initialize();
2927 }
2928
2929 return stepper;
2930}
2931
2932
2933// ----------------------------------------------------------------------------
2961template<class Scalar>
2963 virtual public StepperDIRK<Scalar>
2964{
2965public:
2972 {
2973 typedef Teuchos::ScalarTraits<Scalar> ST;
2974 thetaDefault_ = ST::one()/(2*ST::one());
2975
2976 this->setStepperName("DIRK 1 Stage Theta Method");
2977 this->setStepperType("DIRK 1 Stage Theta Method");
2978 this->setTheta(thetaDefault_);
2979 this->setupTableau();
2980 this->setupDefault();
2981 this->setUseFSAL( false);
2982 this->setICConsistency( "None");
2983 this->setICConsistencyCheck( false);
2984 }
2985
2987 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
2988 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
2989 bool useFSAL,
2990 std::string ICConsistency,
2991 bool ICConsistencyCheck,
2992 bool useEmbedded,
2993 bool zeroInitialGuess,
2994 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
2995 Scalar theta = Scalar(0.5))
2996 {
2997 typedef Teuchos::ScalarTraits<Scalar> ST;
2998 thetaDefault_ = ST::one()/(2*ST::one());
2999
3000 this->setStepperName("DIRK 1 Stage Theta Method");
3001 this->setStepperType("DIRK 1 Stage Theta Method");
3002 this->setTheta(theta);
3003 this->setupTableau();
3004 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3005 useEmbedded, zeroInitialGuess, stepperRKAppAction);
3006 }
3007
3008 void setTheta(Scalar theta)
3009 {
3010 TEUCHOS_TEST_FOR_EXCEPTION(
3011 theta == Teuchos::ScalarTraits<Scalar>::zero(), std::logic_error,
3012 "'theta' can not be zero, as it makes this stepper explicit. \n"
3013 "Try using the 'RK Forward Euler' stepper.\n");
3014 theta_ = theta;
3015 this->setupTableau();
3016 this->isInitialized_ = false;
3017 }
3018
3019 Scalar getTheta() const { return theta_; }
3020
3021 std::string getDescription() const
3022 {
3023 std::ostringstream Description;
3024 Description << this->getStepperType() << "\n"
3025 << "Non-standard finite-difference methods\n"
3026 << "in dynamical systems, P. Kama,\n"
3027 << "Dissertation, University of Pretoria, pg. 49.\n"
3028 << "Comment: Generalized Implicit Midpoint Method\n"
3029 << "c = [ theta ]'\n"
3030 << "A = [ theta ]\n"
3031 << "b = [ 1 ]'";
3032 return Description.str();
3033 }
3034
3035 Teuchos::RCP<const Teuchos::ParameterList>
3037 {
3038 auto pl = this->getValidParametersBasicDIRK();
3039
3040 pl->template set<double>("theta", getTheta(),
3041 "Valid values are 0 <= theta <= 1, where theta = 0 "
3042 "implies Forward Euler, theta = 1/2 implies implicit midpoint "
3043 "method (default), and theta = 1 implies Backward Euler. "
3044 "For theta != 1/2, this method is first-order accurate, "
3045 "and with theta = 1/2, it is second-order accurate. "
3046 "This method is A-stable, but becomes L-stable with theta=1.");
3047
3048 return pl;
3049 }
3050
3051protected:
3052
3054 {
3055 typedef Teuchos::ScalarTraits<Scalar> ST;
3056 int NumStages = 1;
3057 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3058 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3059 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3060 A(0,0) = theta_;
3061 b(0) = ST::one();
3062 c(0) = theta_;
3063
3064 int order = 1;
3065 if ( std::abs((theta_-thetaDefault_)/theta_) < 1.0e-08 ) order = 2;
3066
3067 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
3068 this->getStepperType(),A,b,c,order,1,2));
3069 this->tableau_->setTVD(true);
3070 this->tableau_->setTVDCoeff(2.0);
3071 }
3072
3073 private:
3075 Scalar theta_;
3076};
3077
3078
3080// ------------------------------------------------------------------------
3081template<class Scalar>
3082Teuchos::RCP<StepperDIRK_1StageTheta<Scalar> >
3084 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3085 Teuchos::RCP<Teuchos::ParameterList> pl)
3086{
3087 auto stepper = Teuchos::rcp(new StepperDIRK_1StageTheta<Scalar>());
3088 stepper->setStepperDIRKValues(pl);
3089
3090 if (pl != Teuchos::null) {
3091 stepper->setTheta(pl->get<double>("theta", 0.5));
3092 }
3093
3094 if (model != Teuchos::null) {
3095 stepper->setModel(model);
3096 stepper->initialize();
3097 }
3098
3099 return stepper;
3100}
3101
3102
3103// ----------------------------------------------------------------------------
3130template<class Scalar>
3132 virtual public StepperDIRK<Scalar>
3133{
3134public:
3141 {
3142 typedef Teuchos::ScalarTraits<Scalar> ST;
3143 thetaDefault_ = ST::one()/(2*ST::one());
3144
3145 this->setStepperName("EDIRK 2 Stage Theta Method");
3146 this->setStepperType("EDIRK 2 Stage Theta Method");
3147 this->setTheta(thetaDefault_);
3148 this->setupTableau();
3149 this->setupDefault();
3150 this->setUseFSAL( true);
3151 this->setICConsistency( "Consistent");
3152 this->setICConsistencyCheck( false);
3153 }
3154
3156 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3157 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
3158 bool useFSAL,
3159 std::string ICConsistency,
3160 bool ICConsistencyCheck,
3161 bool useEmbedded,
3162 bool zeroInitialGuess,
3163 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
3164 Scalar theta = Scalar(0.5))
3165 {
3166 typedef Teuchos::ScalarTraits<Scalar> ST;
3167 thetaDefault_ = ST::one()/(2*ST::one());
3168
3169 this->setStepperName("EDIRK 2 Stage Theta Method");
3170 this->setStepperType("EDIRK 2 Stage Theta Method");
3171 this->setTheta(theta);
3172 this->setupTableau();
3173 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3174 useEmbedded, zeroInitialGuess, stepperRKAppAction);
3175 }
3176
3177 void setTheta(Scalar theta)
3178 {
3179 TEUCHOS_TEST_FOR_EXCEPTION(
3180 theta == Teuchos::ScalarTraits<Scalar>::zero(), std::logic_error,
3181 "'theta' can not be zero, as it makes this stepper explicit. \n"
3182 "Try using the 'RK Forward Euler' stepper.\n");
3183 theta_ = theta;
3184 this->isInitialized_ = false;
3185 this->setupTableau();
3186 }
3187
3188 Scalar getTheta() const { return theta_; }
3189
3190 std::string getDescription() const
3191 {
3192 std::ostringstream Description;
3193 Description << this->getStepperType() << "\n"
3194 << "Computer Methods for ODEs and DAEs\n"
3195 << "U. M. Ascher and L. R. Petzold\n"
3196 << "p. 113\n"
3197 << "c = [ 0 1 ]'\n"
3198 << "A = [ 0 0 ]\n"
3199 << " [ 1-theta theta ]\n"
3200 << "b = [ 1-theta theta ]'";
3201 return Description.str();
3202 }
3203
3204 Teuchos::RCP<const Teuchos::ParameterList>
3206 {
3207 auto pl = this->getValidParametersBasicDIRK();
3208
3209 pl->template set<double>("theta", getTheta(),
3210 "Valid values are 0 < theta <= 1, where theta = 0 "
3211 "implies Forward Euler, theta = 1/2 implies trapezoidal "
3212 "method (default), and theta = 1 implies Backward Euler. "
3213 "For theta != 1/2, this method is first-order accurate, "
3214 "and with theta = 1/2, it is second-order accurate. "
3215 "This method is A-stable, but becomes L-stable with theta=1.");
3216
3217 return pl;
3218 }
3219
3220 void setUseFSAL(bool a) { this->useFSAL_ = a; this->isInitialized_ = false; }
3221
3222protected:
3223
3225 {
3226 typedef Teuchos::ScalarTraits<Scalar> ST;
3227 const Scalar one = ST::one();
3228 const Scalar zero = ST::zero();
3229
3230 int NumStages = 2;
3231 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3232 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3233 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3234
3235 // Fill A:
3236 A(0,0) = zero; A(0,1) = zero;
3237 A(1,0) = Teuchos::as<Scalar>( one - theta_ ); A(1,1) = theta_;
3238
3239 // Fill b:
3240 b(0) = Teuchos::as<Scalar>( one - theta_ );
3241 b(1) = theta_;
3242
3243 // Fill c:
3244 c(0) = zero;
3245 c(1) = one;
3246
3247 int order = 1;
3248 if ( std::abs((theta_-thetaDefault_)/theta_) < 1.0e-08 ) order = 2;
3249
3250 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
3251 this->getStepperType(),A,b,c,order,1,2));
3252 this->tableau_->setTVD(true);
3253 this->tableau_->setTVDCoeff(2.0);
3254 }
3255
3256 private:
3258 Scalar theta_;
3259};
3260
3261
3263// ------------------------------------------------------------------------
3264template<class Scalar>
3265Teuchos::RCP<StepperEDIRK_2StageTheta<Scalar> >
3267 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3268 Teuchos::RCP<Teuchos::ParameterList> pl)
3269{
3270 auto stepper = Teuchos::rcp(new StepperEDIRK_2StageTheta<Scalar>());
3271 stepper->setStepperDIRKValues(pl);
3272
3273 if (pl != Teuchos::null) {
3274 stepper->setTheta(pl->get<double>("theta", 0.5));
3275 }
3276
3277 if (model != Teuchos::null) {
3278 stepper->setModel(model);
3279 stepper->initialize();
3280 }
3281
3282 return stepper;
3283}
3284
3285
3286// ----------------------------------------------------------------------------
3307template<class Scalar>
3309 virtual public StepperDIRK<Scalar>
3310{
3311public:
3318 {
3319 this->setStepperName("RK Trapezoidal Rule");
3320 this->setStepperType("RK Trapezoidal Rule");
3321 this->setupTableau();
3322 this->setupDefault();
3323 this->setUseFSAL( true);
3324 this->setICConsistency( "Consistent");
3325 this->setICConsistencyCheck( false);
3326 }
3327
3329 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3330 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
3331 bool useFSAL,
3332 std::string ICConsistency,
3333 bool ICConsistencyCheck,
3334 bool useEmbedded,
3335 bool zeroInitialGuess,
3336 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3337 {
3338 this->setStepperName("RK Trapezoidal Rule");
3339 this->setStepperType("RK Trapezoidal Rule");
3340 this->setupTableau();
3341 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3342 useEmbedded, zeroInitialGuess, stepperRKAppAction);
3343 }
3344
3345 std::string getDescription() const
3346 {
3347 std::ostringstream Description;
3348 Description << this->getStepperType() << "\n"
3349 << "Also known as Crank-Nicolson Method.\n"
3350 << "c = [ 0 1 ]'\n"
3351 << "A = [ 0 0 ]\n"
3352 << " [ 1/2 1/2 ]\n"
3353 << "b = [ 1/2 1/2 ]'";
3354 return Description.str();
3355 }
3356
3357 void setUseFSAL(bool a) { this->useFSAL_ = a; this->isInitialized_ = false; }
3358
3359protected:
3360
3362 {
3363 typedef Teuchos::ScalarTraits<Scalar> ST;
3364 const Scalar one = ST::one();
3365 const Scalar zero = ST::zero();
3366 const Scalar onehalf = ST::one()/(2*ST::one());
3367
3368 int NumStages = 2;
3369 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3370 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3371 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3372
3373 // Fill A:
3374 A(0,0) = zero; A(0,1) = zero;
3375 A(1,0) = onehalf; A(1,1) = onehalf;
3376
3377 // Fill b:
3378 b(0) = onehalf;
3379 b(1) = onehalf;
3380
3381 // Fill c:
3382 c(0) = zero;
3383 c(1) = one;
3384
3385 int order = 2;
3386
3387 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
3388 this->getStepperType(),A,b,c,order,order,order));
3389 this->tableau_->setTVD(true);
3390 this->tableau_->setTVDCoeff(2.0);
3391 }
3392};
3393
3394
3396// ------------------------------------------------------------------------
3397template<class Scalar>
3398Teuchos::RCP<StepperEDIRK_TrapezoidalRule<Scalar> >
3400 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3401 Teuchos::RCP<Teuchos::ParameterList> pl)
3402{
3403 auto stepper = Teuchos::rcp(new StepperEDIRK_TrapezoidalRule<Scalar>());
3404
3405 // Test for aliases.
3406 if (pl != Teuchos::null) {
3407 auto stepperType =
3408 pl->get<std::string>("Stepper Type", stepper->getStepperType());
3409
3410 TEUCHOS_TEST_FOR_EXCEPTION(
3411 stepperType != stepper->getStepperType() &&
3412 stepperType != "RK Crank-Nicolson", std::logic_error,
3413 " ParameterList 'Stepper Type' (='" + stepperType +"')\n"
3414 " does not match type for this Stepper (='" + stepper->getStepperType() +
3415 "')\n or one of its aliases ('RK Crank-Nicolson').\n");
3416
3417 // Reset default StepperType.
3418 pl->set<std::string>("Stepper Type", stepper->getStepperType());
3419 }
3420
3421 stepper->setStepperDIRKValues(pl);
3422
3423 if (model != Teuchos::null) {
3424 stepper->setModel(model);
3425 stepper->initialize();
3426 }
3427
3428 return stepper;
3429}
3430
3431
3432// ----------------------------------------------------------------------------
3459template<class Scalar>
3461 virtual public StepperDIRK<Scalar>
3462{
3463public:
3470 {
3471 this->setStepperName("RK Implicit Midpoint");
3472 this->setStepperType("RK Implicit Midpoint");
3473 this->setupTableau();
3474 this->setupDefault();
3475 this->setUseFSAL( false);
3476 this->setICConsistency( "None");
3477 this->setICConsistencyCheck( false);
3478 }
3479
3481 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3482 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
3483 bool useFSAL,
3484 std::string ICConsistency,
3485 bool ICConsistencyCheck,
3486 bool useEmbedded,
3487 bool zeroInitialGuess,
3488 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3489 {
3490 this->setStepperName("RK Implicit Midpoint");
3491 this->setStepperType("RK Implicit Midpoint");
3492 this->setupTableau();
3493 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3494 useEmbedded, zeroInitialGuess, stepperRKAppAction);
3495 }
3496
3497 std::string getDescription() const
3498 {
3499 std::ostringstream Description;
3500 Description << this->getStepperType() << "\n"
3501 << "A-stable\n"
3502 << "Solving Ordinary Differential Equations II:\n"
3503 << "Stiff and Differential-Algebraic Problems,\n"
3504 << "2nd Revised Edition\n"
3505 << "E. Hairer and G. Wanner\n"
3506 << "Table 5.2, pg 72\n"
3507 << "Solving Ordinary Differential Equations I:\n"
3508 << "Nonstiff Problems, 2nd Revised Edition\n"
3509 << "E. Hairer, S. P. Norsett, and G. Wanner\n"
3510 << "Table 7.1, pg 205\n"
3511 << "c = [ 1/2 ]'\n"
3512 << "A = [ 1/2 ]\n"
3513 << "b = [ 1 ]'";
3514 return Description.str();
3515 }
3516
3517protected:
3518
3520 {
3521 typedef Teuchos::ScalarTraits<Scalar> ST;
3522 int NumStages = 1;
3523 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3524 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3525 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3526 const Scalar onehalf = ST::one()/(2*ST::one());
3527 const Scalar one = ST::one();
3528
3529 // Fill A:
3530 A(0,0) = onehalf;
3531
3532 // Fill b:
3533 b(0) = one;
3534
3535 // Fill c:
3536 c(0) = onehalf;
3537
3538 int order = 2;
3539
3540 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
3541 this->getStepperType(),A,b,c,order,order,order));
3542 this->tableau_->setTVD(true);
3543 this->tableau_->setTVDCoeff(2.0);
3544 }
3545};
3546
3547
3549// ------------------------------------------------------------------------
3550template<class Scalar>
3551Teuchos::RCP<StepperSDIRK_ImplicitMidpoint<Scalar> >
3553 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3554 Teuchos::RCP<Teuchos::ParameterList> pl)
3555{
3556 auto stepper = Teuchos::rcp(new StepperSDIRK_ImplicitMidpoint<Scalar>());
3557 stepper->setStepperDIRKValues(pl);
3558
3559 if (model != Teuchos::null) {
3560 stepper->setModel(model);
3561 stepper->initialize();
3562 }
3563
3564 return stepper;
3565}
3566
3567
3568// ----------------------------------------------------------------------------
3589template<class Scalar>
3591 virtual public StepperDIRK<Scalar>
3592{
3593 public:
3595 {
3596 this->setStepperName("SSPDIRK22");
3597 this->setStepperType("SSPDIRK22");
3598 this->setupTableau();
3599 this->setupDefault();
3600 this->setUseFSAL( false);
3601 this->setICConsistency( "None");
3602 this->setICConsistencyCheck( false);
3603 }
3604
3606 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3607 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
3608 bool useFSAL,
3609 std::string ICConsistency,
3610 bool ICConsistencyCheck,
3611 bool useEmbedded,
3612 bool zeroInitialGuess,
3613 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3614 {
3615 this->setStepperName("SSPDIRK22");
3616 this->setStepperType("SSPDIRK22");
3617 this->setupTableau();
3618 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3619 useEmbedded, zeroInitialGuess, stepperRKAppAction);
3620 }
3621
3622 std::string getDescription() const
3623 {
3624 std::ostringstream Description;
3625 Description << this->getStepperType() << "\n"
3626 << "Strong Stability Preserving Diagonally-Implicit RK (stage=2, order=2)\n"
3627 << "SSP-Coef = 4\n"
3628 << "c = [ 1/4 3/4 ]'\n"
3629 << "A = [ 1/4 ]\n"
3630 << " [ 1/2 1/4 ]\n"
3631 << "b = [ 1/2 1/2 ]\n" << std::endl;
3632 return Description.str();
3633 }
3634
3635protected:
3636
3638 {
3639 typedef Teuchos::ScalarTraits<Scalar> ST;
3640 using Teuchos::as;
3641 const int NumStages = 2;
3642 const int order = 2;
3643 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3644 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3645 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3646
3647 const Scalar one = ST::one();
3648 const Scalar zero = ST::zero();
3649 const Scalar onehalf = one/(2*one);
3650 const Scalar onefourth = one/(4*one);
3651
3652 // Fill A:
3653 A(0,0) = A(1,1) = onefourth;
3654 A(0,1) = zero;
3655 A(1,0) = onehalf;
3656
3657 // Fill b:
3658 b(0) = b(1) = onehalf;
3659
3660 // Fill c:
3661 c(0) = A(0,0);
3662 c(1) = A(1,0) + A(1,1);
3663
3664 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
3665 this->getStepperType(),A,b,c,order,order,order));
3666 this->tableau_->setTVD(true);
3667 this->tableau_->setTVDCoeff(4.0);
3668 }
3669};
3670
3671
3673// ------------------------------------------------------------------------
3674template<class Scalar>
3675Teuchos::RCP<StepperSDIRK_SSPDIRK22<Scalar> >
3677 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3678 Teuchos::RCP<Teuchos::ParameterList> pl)
3679{
3680 auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK22<Scalar>());
3681 stepper->setStepperDIRKValues(pl);
3682
3683 if (model != Teuchos::null) {
3684 stepper->setModel(model);
3685 stepper->initialize();
3686 }
3687
3688 return stepper;
3689}
3690
3691
3692// ----------------------------------------------------------------------------
3714template<class Scalar>
3716 virtual public StepperDIRK<Scalar>
3717{
3718 public:
3720 {
3721 this->setStepperName("SSPDIRK32");
3722 this->setStepperType("SSPDIRK32");
3723 this->setupTableau();
3724 this->setupDefault();
3725 this->setUseFSAL( false);
3726 this->setICConsistency( "None");
3727 this->setICConsistencyCheck( false);
3728 }
3729
3731 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3732 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
3733 bool useFSAL,
3734 std::string ICConsistency,
3735 bool ICConsistencyCheck,
3736 bool useEmbedded,
3737 bool zeroInitialGuess,
3738 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3739 {
3740 this->setStepperName("SSPDIRK32");
3741 this->setStepperType("SSPDIRK32");
3742 this->setupTableau();
3743 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3744 useEmbedded, zeroInitialGuess, stepperRKAppAction);
3745 }
3746
3747 std::string getDescription() const
3748 {
3749 std::ostringstream Description;
3750 Description << this->getStepperType() << "\n"
3751 << "Strong Stability Preserving Diagonally-Implicit RK (stage=3, order=2)\n"
3752 << "SSP-Coef = 6\n"
3753 << "c = [ 1/6 1/2 5/6 ]'\n"
3754 << "A = [ 1/6 ]\n"
3755 << " [ 1/3 1/6 ]\n"
3756 << " [ 1/3 1/3 1/6 ]\n"
3757 << "b = [ 1/3 1/3 1/3 ]\n" << std::endl;
3758 return Description.str();
3759 }
3760
3761protected:
3762
3764 {
3765
3766 typedef Teuchos::ScalarTraits<Scalar> ST;
3767 using Teuchos::as;
3768 const int NumStages = 3;
3769 const int order = 2;
3770 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3771 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3772 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3773
3774 const Scalar one = ST::one();
3775 const Scalar zero = ST::zero();
3776 const Scalar onethird = one/(3*one);
3777 const Scalar onesixth = one/(6*one);
3778
3779 // Fill A:
3780 A(0,0) = A(1,1) = A(2,2) = onesixth;
3781 A(1,0) = A(2,0) = A(2,1) = onethird;
3782 A(0,1) = A(0,2) = A(1,2) = zero;
3783
3784 // Fill b:
3785 b(0) = b(1) = b(2) = onethird;
3786
3787 // Fill c:
3788 c(0) = A(0,0);
3789 c(1) = A(1,0) + A(1,1);
3790 c(2) = A(2,0) + A(2,1) + A(2,2);
3791
3792 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
3793 this->getStepperType(),A,b,c,order,order,order));
3794 this->tableau_->setTVD(true);
3795 this->tableau_->setTVDCoeff(6.0);
3796 }
3797};
3798
3799
3801// ------------------------------------------------------------------------
3802template<class Scalar>
3803Teuchos::RCP<StepperSDIRK_SSPDIRK32<Scalar> >
3805 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3806 Teuchos::RCP<Teuchos::ParameterList> pl)
3807{
3808 auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK32<Scalar>());
3809 stepper->setStepperDIRKValues(pl);
3810
3811 if (model != Teuchos::null) {
3812 stepper->setModel(model);
3813 stepper->initialize();
3814 }
3815
3816 return stepper;
3817}
3818
3819
3820// ----------------------------------------------------------------------------
3841template<class Scalar>
3843 virtual public StepperDIRK<Scalar>
3844{
3845 public:
3847 {
3848 this->setStepperName("SSPDIRK23");
3849 this->setStepperType("SSPDIRK23");
3850 this->setupTableau();
3851 this->setupDefault();
3852 this->setUseFSAL( false);
3853 this->setICConsistency( "None");
3854 this->setICConsistencyCheck( false);
3855 }
3856
3858 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3859 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
3860 bool useFSAL,
3861 std::string ICConsistency,
3862 bool ICConsistencyCheck,
3863 bool useEmbedded,
3864 bool zeroInitialGuess,
3865 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3866 {
3867 this->setStepperName("SSPDIRK23");
3868 this->setStepperType("SSPDIRK23");
3869 this->setupTableau();
3870 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3871 useEmbedded, zeroInitialGuess, stepperRKAppAction);
3872 }
3873
3874 std::string getDescription() const
3875 {
3876 std::ostringstream Description;
3877 Description << this->getStepperType() << "\n"
3878 << "Strong Stability Preserving Diagonally-Implicit RK (stage=2, order=3)\n"
3879 << "SSP-Coef = 1 + sqrt( 3 )\n"
3880 << "c = [ 1/(3 + sqrt( 3 )) (1/6)(3 + sqrt( 3 )) ] '\n"
3881 << "A = [ 1/(3 + sqrt( 3 )) ] \n"
3882 << " [ 1/sqrt( 3 ) 1/(3 + sqrt( 3 )) ] \n"
3883 << "b = [ 1/2 1/2 ] \n" << std::endl;
3884 return Description.str();
3885 }
3886
3887protected:
3888
3890 {
3891
3892 typedef Teuchos::ScalarTraits<Scalar> ST;
3893 using Teuchos::as;
3894 const int NumStages = 2;
3895 const int order = 3;
3896 const Scalar sspcoef = 2.7321;
3897 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3898 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3899 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3900
3901 const Scalar one = ST::one();
3902 const Scalar zero = ST::zero();
3903 const Scalar onehalf = one/(2*one);
3904 const Scalar rootthree = ST::squareroot(3*one);
3905
3906 // Fill A:
3907 A(0,0) = A(1,1) = one/(3*one + rootthree);
3908 A(1,0) = one/rootthree;
3909 A(0,1) = zero;
3910
3911 // Fill b:
3912 b(0) = b(1) = onehalf;
3913
3914 // Fill c:
3915 c(0) = A(0,0);
3916 c(1) = A(1,0) + A(1,1);
3917
3918 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
3919 this->getStepperType(),A,b,c,order,order,order));
3920 this->tableau_->setTVD(true);
3921 this->tableau_->setTVDCoeff(sspcoef);
3922 }
3923};
3924
3925
3927// ------------------------------------------------------------------------
3928template<class Scalar>
3929Teuchos::RCP<StepperSDIRK_SSPDIRK23<Scalar> >
3931 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
3932 Teuchos::RCP<Teuchos::ParameterList> pl)
3933{
3934 auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK23<Scalar>());
3935 stepper->setStepperDIRKValues(pl);
3936
3937 if (model != Teuchos::null) {
3938 stepper->setModel(model);
3939 stepper->initialize();
3940 }
3941
3942 return stepper;
3943}
3944
3945
3946// ----------------------------------------------------------------------------
3968template<class Scalar>
3970 virtual public StepperDIRK<Scalar>
3971{
3972 public:
3974 {
3975 this->setStepperName("SSPDIRK33");
3976 this->setStepperType("SSPDIRK33");
3977 this->setupTableau();
3978 this->setupDefault();
3979 this->setUseFSAL( false);
3980 this->setICConsistency( "None");
3981 this->setICConsistencyCheck( false);
3982 }
3983
3985 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
3986 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
3987 bool useFSAL,
3988 std::string ICConsistency,
3989 bool ICConsistencyCheck,
3990 bool useEmbedded,
3991 bool zeroInitialGuess,
3992 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
3993 {
3994 this->setStepperName("SSPDIRK33");
3995 this->setStepperType("SSPDIRK33");
3996 this->setupTableau();
3997 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
3998 useEmbedded, zeroInitialGuess, stepperRKAppAction);
3999 }
4000
4001 std::string getDescription() const
4002 {
4003 std::ostringstream Description;
4004 Description << this->getStepperType() << "\n"
4005 << "Strong Stability Preserving Diagonally-Implicit RK (stage=3, order=3)\n"
4006 << "SSP-Coef = 2 + 2 sqrt(2)\n"
4007 << "c = [ 1/( 4 + 2 sqrt(2) 1/2 (1/4)(2 + sqrt(2) ] '\n"
4008 << "A = [ 1/( 4 + 2 sqrt(2) ] \n"
4009 << " [ 1/(2 sqrt(2) 1/( 4 + 2 sqrt(2) ] \n"
4010 << " [ 1/(2 sqrt(2) 1/(2 sqrt(2) 1/( 4 + 2 sqrt(2) ] \n"
4011 << "b = [ 1/3 1/3 1/3 ] \n"
4012 << std::endl;
4013 return Description.str();
4014 }
4015
4016protected:
4017
4019 {
4020
4021 typedef Teuchos::ScalarTraits<Scalar> ST;
4022 using Teuchos::as;
4023 const int NumStages = 3;
4024 const int order = 3;
4025 const Scalar sspcoef= 4.8284;
4026 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
4027 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
4028 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
4029
4030 const Scalar one = ST::one();
4031 const Scalar zero = ST::zero();
4032 const Scalar onethird = one/(3*one);
4033 const Scalar rootwo = ST::squareroot(2*one);
4034
4035 // Fill A:
4036 A(0,0) = A(1,1) = A(2,2) = one / (4*one + 2*rootwo);
4037 A(1,0) = A(2,0) = A(2,1) = one / (2*rootwo);
4038 A(0,1) = A(0,2) = A(1,2) = zero;
4039
4040 // Fill b:
4041 b(0) = b(1) = b(2) = onethird;
4042
4043 // Fill c:
4044 c(0) = A(0,0);
4045 c(1) = A(1,0) + A(1,1);
4046 c(2) = A(2,0) + A(2,1) + A(2,2);
4047
4048 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
4049 this->getStepperType(),A,b,c,order,order,order));
4050 this->tableau_->setTVD(true);
4051 this->tableau_->setTVDCoeff(sspcoef);
4052 }
4053};
4054
4055
4057// ------------------------------------------------------------------------
4058template<class Scalar>
4059Teuchos::RCP<StepperSDIRK_SSPDIRK33<Scalar> >
4061 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4062 Teuchos::RCP<Teuchos::ParameterList> pl)
4063{
4064 auto stepper = Teuchos::rcp(new StepperSDIRK_SSPDIRK33<Scalar>());
4065 stepper->setStepperDIRKValues(pl);
4066
4067 if (model != Teuchos::null) {
4068 stepper->setModel(model);
4069 stepper->initialize();
4070 }
4071
4072 return stepper;
4073}
4074
4075
4076// ----------------------------------------------------------------------------
4097template<class Scalar>
4099 virtual public StepperDIRK<Scalar>
4100{
4101public:
4108 {
4109 this->setStepperName("RK Implicit 1 Stage 1st order Radau IA");
4110 this->setStepperType("RK Implicit 1 Stage 1st order Radau IA");
4111 this->setupTableau();
4112 this->setupDefault();
4113 this->setUseFSAL( false);
4114 this->setICConsistency( "None");
4115 this->setICConsistencyCheck( false);
4116 }
4117
4119 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4120 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
4121 bool useFSAL,
4122 std::string ICConsistency,
4123 bool ICConsistencyCheck,
4124 bool useEmbedded,
4125 bool zeroInitialGuess,
4126 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4127 {
4128 this->setStepperName("RK Implicit 1 Stage 1st order Radau IA");
4129 this->setStepperType("RK Implicit 1 Stage 1st order Radau IA");
4130 this->setupTableau();
4131 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4132 useEmbedded, zeroInitialGuess, stepperRKAppAction);
4133 }
4134
4135 std::string getDescription() const
4136 {
4137 std::ostringstream Description;
4138 Description << this->getStepperType() << "\n"
4139 << "A-stable\n"
4140 << "Solving Ordinary Differential Equations II:\n"
4141 << "Stiff and Differential-Algebraic Problems,\n"
4142 << "2nd Revised Edition\n"
4143 << "E. Hairer and G. Wanner\n"
4144 << "Table 5.3, pg 73\n"
4145 << "c = [ 0 ]'\n"
4146 << "A = [ 1 ]\n"
4147 << "b = [ 1 ]'";
4148 return Description.str();
4149 }
4150
4151protected:
4152
4154 {
4155 typedef Teuchos::ScalarTraits<Scalar> ST;
4156 int NumStages = 1;
4157 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
4158 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
4159 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
4160 const Scalar one = ST::one();
4161 const Scalar zero = ST::zero();
4162 A(0,0) = one;
4163 b(0) = one;
4164 c(0) = zero;
4165 int order = 1;
4166
4167 auto emptyBStar = Teuchos::SerialDenseVector<int,Scalar>();
4168 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
4169 this->getStepperType(),A,b,c,order,order,order,emptyBStar,false));
4170 }
4171};
4172
4173
4175// ------------------------------------------------------------------------
4176template<class Scalar>
4177Teuchos::RCP<StepperDIRK_1Stage1stOrderRadauIA<Scalar> >
4179 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4180 Teuchos::RCP<Teuchos::ParameterList> pl)
4181{
4182 auto stepper = Teuchos::rcp(new StepperDIRK_1Stage1stOrderRadauIA<Scalar>());
4183 stepper->setStepperDIRKValues(pl);
4184
4185 if (model != Teuchos::null) {
4186 stepper->setModel(model);
4187 stepper->initialize();
4188 }
4189
4190 return stepper;
4191}
4192
4193
4194// ----------------------------------------------------------------------------
4217template<class Scalar>
4219 virtual public StepperDIRK<Scalar>
4220{
4221public:
4228 {
4229 this->setStepperName("RK Implicit 2 Stage 2nd order Lobatto IIIB");
4230 this->setStepperType("RK Implicit 2 Stage 2nd order Lobatto IIIB");
4231 this->setupTableau();
4232 this->setupDefault();
4233 this->setUseFSAL( false);
4234 this->setICConsistency( "None");
4235 this->setICConsistencyCheck( false);
4236 }
4237
4239 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4240 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
4241 bool useFSAL,
4242 std::string ICConsistency,
4243 bool ICConsistencyCheck,
4244 bool useEmbedded,
4245 bool zeroInitialGuess,
4246 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4247 {
4248 this->setStepperName("RK Implicit 2 Stage 2nd order Lobatto IIIB");
4249 this->setStepperType("RK Implicit 2 Stage 2nd order Lobatto IIIB");
4250 this->setupTableau();
4251 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4252 useEmbedded, zeroInitialGuess, stepperRKAppAction);
4253 }
4254
4255 std::string getDescription() const
4256 {
4257 std::ostringstream Description;
4258 Description << this->getStepperType() << "\n"
4259 << "A-stable\n"
4260 << "Solving Ordinary Differential Equations II:\n"
4261 << "Stiff and Differential-Algebraic Problems,\n"
4262 << "2nd Revised Edition\n"
4263 << "E. Hairer and G. Wanner\n"
4264 << "Table 5.9, pg 76\n"
4265 << "c = [ 0 1 ]'\n"
4266 << "A = [ 1/2 0 ]\n"
4267 << " [ 1/2 0 ]\n"
4268 << "b = [ 1/2 1/2 ]'";
4269 return Description.str();
4270 }
4271
4272protected:
4273
4275 {
4276 typedef Teuchos::ScalarTraits<Scalar> ST;
4277 using Teuchos::as;
4278 int NumStages = 2;
4279 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
4280 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
4281 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
4282 const Scalar zero = ST::zero();
4283 const Scalar one = ST::one();
4284
4285 // Fill A:
4286 A(0,0) = as<Scalar>( one/(2*one) );
4287 A(0,1) = zero;
4288 A(1,0) = as<Scalar>( one/(2*one) );
4289 A(1,1) = zero;
4290
4291 // Fill b:
4292 b(0) = as<Scalar>( one/(2*one) );
4293 b(1) = as<Scalar>( one/(2*one) );
4294
4295 // Fill c:
4296 c(0) = zero;
4297 c(1) = one;
4298 int order = 2;
4299
4300 auto emptyBStar = Teuchos::SerialDenseVector<int,Scalar>();
4301 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
4302 this->getStepperType(),A,b,c,order,order,order,emptyBStar,false));
4303 this->tableau_->setTVD(true);
4304 this->tableau_->setTVDCoeff(2.0);
4305 //TODO: fix this
4306 }
4307
4308};
4309
4310
4312// ------------------------------------------------------------------------
4313template<class Scalar>
4314Teuchos::RCP<StepperDIRK_2Stage2ndOrderLobattoIIIB<Scalar> >
4316 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4317 Teuchos::RCP<Teuchos::ParameterList> pl)
4318{
4319 auto stepper = Teuchos::rcp(new StepperDIRK_2Stage2ndOrderLobattoIIIB<Scalar>());
4320 stepper->setStepperDIRKValues(pl);
4321
4322 if (model != Teuchos::null) {
4323 stepper->setModel(model);
4324 stepper->initialize();
4325 }
4326
4327 return stepper;
4328}
4329
4330
4331// ----------------------------------------------------------------------------
4357template<class Scalar>
4359 virtual public StepperDIRK<Scalar>
4360{
4361public:
4368 {
4369 this->setStepperName("SDIRK 5 Stage 4th order");
4370 this->setStepperType("SDIRK 5 Stage 4th order");
4371 this->setupTableau();
4372 this->setupDefault();
4373 this->setUseFSAL( false);
4374 this->setICConsistency( "None");
4375 this->setICConsistencyCheck( false);
4376 }
4377
4379 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4380 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
4381 bool useFSAL,
4382 std::string ICConsistency,
4383 bool ICConsistencyCheck,
4384 bool useEmbedded,
4385 bool zeroInitialGuess,
4386 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4387 {
4388 this->setStepperName("SDIRK 5 Stage 4th order");
4389 this->setStepperType("SDIRK 5 Stage 4th order");
4390 this->setupTableau();
4391 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4392 useEmbedded, zeroInitialGuess, stepperRKAppAction);
4393 }
4394
4395 std::string getDescription() const
4396 {
4397 std::ostringstream Description;
4398 Description << this->getStepperType() << "\n"
4399 << "L-stable\n"
4400 << "Solving Ordinary Differential Equations II:\n"
4401 << "Stiff and Differential-Algebraic Problems,\n"
4402 << "2nd Revised Edition\n"
4403 << "E. Hairer and G. Wanner\n"
4404 << "pg100 \n"
4405 << "c = [ 1/4 3/4 11/20 1/2 1 ]'\n"
4406 << "A = [ 1/4 ]\n"
4407 << " [ 1/2 1/4 ]\n"
4408 << " [ 17/50 -1/25 1/4 ]\n"
4409 << " [ 371/1360 -137/2720 15/544 1/4 ]\n"
4410 << " [ 25/24 -49/48 125/16 -85/12 1/4 ]\n"
4411 << "b = [ 25/24 -49/48 125/16 -85/12 1/4 ]'";
4412 // << "b = [ 59/48 -17/96 225/32 -85/12 0 ]'";
4413 return Description.str();
4414 }
4415
4416protected:
4417
4419 {
4420 typedef Teuchos::ScalarTraits<Scalar> ST;
4421 using Teuchos::as;
4422 int NumStages = 5;
4423 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
4424 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
4425 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
4426 const Scalar zero = ST::zero();
4427 const Scalar one = ST::one();
4428 const Scalar onequarter = as<Scalar>( one/(4*one) );
4429
4430 // Fill A:
4431 A(0,0) = onequarter;
4432 A(0,1) = zero;
4433 A(0,2) = zero;
4434 A(0,3) = zero;
4435 A(0,4) = zero;
4436
4437 A(1,0) = as<Scalar>( one / (2*one) );
4438 A(1,1) = onequarter;
4439 A(1,2) = zero;
4440 A(1,3) = zero;
4441 A(1,4) = zero;
4442
4443 A(2,0) = as<Scalar>( 17*one/(50*one) );
4444 A(2,1) = as<Scalar>( -one/(25*one) );
4445 A(2,2) = onequarter;
4446 A(2,3) = zero;
4447 A(2,4) = zero;
4448
4449 A(3,0) = as<Scalar>( 371*one/(1360*one) );
4450 A(3,1) = as<Scalar>( -137*one/(2720*one) );
4451 A(3,2) = as<Scalar>( 15*one/(544*one) );
4452 A(3,3) = onequarter;
4453 A(3,4) = zero;
4454
4455 A(4,0) = as<Scalar>( 25*one/(24*one) );
4456 A(4,1) = as<Scalar>( -49*one/(48*one) );
4457 A(4,2) = as<Scalar>( 125*one/(16*one) );
4458 A(4,3) = as<Scalar>( -85*one/(12*one) );
4459 A(4,4) = onequarter;
4460
4461 // Fill b:
4462 b(0) = as<Scalar>( 25*one/(24*one) );
4463 b(1) = as<Scalar>( -49*one/(48*one) );
4464 b(2) = as<Scalar>( 125*one/(16*one) );
4465 b(3) = as<Scalar>( -85*one/(12*one) );
4466 b(4) = onequarter;
4467
4468 /*
4469 // Alternate version
4470 b(0) = as<Scalar>( 59*one/(48*one) );
4471 b(1) = as<Scalar>( -17*one/(96*one) );
4472 b(2) = as<Scalar>( 225*one/(32*one) );
4473 b(3) = as<Scalar>( -85*one/(12*one) );
4474 b(4) = zero;
4475 */
4476
4477 // Fill c:
4478 c(0) = onequarter;
4479 c(1) = as<Scalar>( 3*one/(4*one) );
4480 c(2) = as<Scalar>( 11*one/(20*one) );
4481 c(3) = as<Scalar>( one/(2*one) );
4482 c(4) = one;
4483
4484 int order = 4;
4485
4486 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
4487 this->getStepperType(),A,b,c,order,order,order));
4488 }
4489};
4490
4491
4493// ------------------------------------------------------------------------
4494template<class Scalar>
4495Teuchos::RCP<StepperSDIRK_5Stage4thOrder<Scalar> >
4497 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4498 Teuchos::RCP<Teuchos::ParameterList> pl)
4499{
4500 auto stepper = Teuchos::rcp(new StepperSDIRK_5Stage4thOrder<Scalar>());
4501 stepper->setStepperDIRKValues(pl);
4502
4503 if (model != Teuchos::null) {
4504 stepper->setModel(model);
4505 stepper->initialize();
4506 }
4507
4508 return stepper;
4509}
4510
4511
4512// ----------------------------------------------------------------------------
4537template<class Scalar>
4539 virtual public StepperDIRK<Scalar>
4540{
4541public:
4548 {
4549 this->setStepperName("SDIRK 3 Stage 4th order");
4550 this->setStepperType("SDIRK 3 Stage 4th order");
4551 this->setupTableau();
4552 this->setupDefault();
4553 this->setUseFSAL( false);
4554 this->setICConsistency( "None");
4555 this->setICConsistencyCheck( false);
4556 }
4557
4559 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4560 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
4561 bool useFSAL,
4562 std::string ICConsistency,
4563 bool ICConsistencyCheck,
4564 bool useEmbedded,
4565 bool zeroInitialGuess,
4566 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4567 {
4568 this->setStepperName("SDIRK 3 Stage 4th order");
4569 this->setStepperType("SDIRK 3 Stage 4th order");
4570 this->setupTableau();
4571 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4572 useEmbedded, zeroInitialGuess, stepperRKAppAction);
4573 }
4574
4575 std::string getDescription() const
4576 {
4577 std::ostringstream Description;
4578 Description << this->getStepperType() << "\n"
4579 << "A-stable\n"
4580 << "Solving Ordinary Differential Equations II:\n"
4581 << "Stiff and Differential-Algebraic Problems,\n"
4582 << "2nd Revised Edition\n"
4583 << "E. Hairer and G. Wanner\n"
4584 << "p. 100 \n"
4585 << "gamma = (1/sqrt(3))*cos(pi/18)+1/2\n"
4586 << "delta = 1/(6*(2*gamma-1)^2)\n"
4587 << "c = [ gamma 1/2 1-gamma ]'\n"
4588 << "A = [ gamma ]\n"
4589 << " [ 1/2-gamma gamma ]\n"
4590 << " [ 2*gamma 1-4*gamma gamma ]\n"
4591 << "b = [ delta 1-2*delta delta ]'";
4592 return Description.str();
4593 }
4594
4595protected:
4596
4598 {
4599 typedef Teuchos::ScalarTraits<Scalar> ST;
4600 using Teuchos::as;
4601 int NumStages = 3;
4602 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
4603 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
4604 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
4605 const Scalar zero = ST::zero();
4606 const Scalar one = ST::one();
4607 const Scalar pi = as<Scalar>(4*one)*std::atan(one);
4608 const Scalar gamma = as<Scalar>( one/ST::squareroot(3*one)*std::cos(pi/(18*one))+one/(2*one) );
4609 const Scalar delta = as<Scalar>( one/(6*one*std::pow(2*gamma-one,2*one)) );
4610
4611 // Fill A:
4612 A(0,0) = gamma;
4613 A(0,1) = zero;
4614 A(0,2) = zero;
4615
4616 A(1,0) = as<Scalar>( one/(2*one) - gamma );
4617 A(1,1) = gamma;
4618 A(1,2) = zero;
4619
4620 A(2,0) = as<Scalar>( 2*gamma );
4621 A(2,1) = as<Scalar>( one - 4*gamma );
4622 A(2,2) = gamma;
4623
4624 // Fill b:
4625 b(0) = delta;
4626 b(1) = as<Scalar>( one-2*delta );
4627 b(2) = delta;
4628
4629 // Fill c:
4630 c(0) = gamma;
4631 c(1) = as<Scalar>( one/(2*one) );
4632 c(2) = as<Scalar>( one - gamma );
4633
4634 int order = 4;
4635
4636 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
4637 this->getStepperType(),A,b,c,order,order,order));
4638 }
4639};
4640
4641
4643// ------------------------------------------------------------------------
4644template<class Scalar>
4645Teuchos::RCP<StepperSDIRK_3Stage4thOrder<Scalar> >
4647 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4648 Teuchos::RCP<Teuchos::ParameterList> pl)
4649{
4650 auto stepper = Teuchos::rcp(new StepperSDIRK_3Stage4thOrder<Scalar>());
4651 stepper->setStepperDIRKValues(pl);
4652
4653 if (model != Teuchos::null) {
4654 stepper->setModel(model);
4655 stepper->initialize();
4656 }
4657
4658 return stepper;
4659}
4660
4661
4662// ----------------------------------------------------------------------------
4688template<class Scalar>
4690 virtual public StepperDIRK<Scalar>
4691{
4692public:
4699 {
4700 this->setStepperName("SDIRK 5 Stage 5th order");
4701 this->setStepperType("SDIRK 5 Stage 5th order");
4702 this->setupTableau();
4703 this->setupDefault();
4704 this->setUseFSAL( false);
4705 this->setICConsistency( "None");
4706 this->setICConsistencyCheck( false);
4707 }
4708
4710 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4711 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
4712 bool useFSAL,
4713 std::string ICConsistency,
4714 bool ICConsistencyCheck,
4715 bool useEmbedded,
4716 bool zeroInitialGuess,
4717 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4718 {
4719 this->setStepperName("SDIRK 5 Stage 5th order");
4720 this->setStepperType("SDIRK 5 Stage 5th order");
4721 this->setupTableau();
4722 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4723 useEmbedded, zeroInitialGuess, stepperRKAppAction);
4724 }
4725
4726 std::string getDescription() const
4727 {
4728 std::ostringstream Description;
4729 Description << this->getStepperType() << "\n"
4730 << "Solving Ordinary Differential Equations II:\n"
4731 << "Stiff and Differential-Algebraic Problems,\n"
4732 << "2nd Revised Edition\n"
4733 << "E. Hairer and G. Wanner\n"
4734 << "pg101 \n"
4735 << "c = [ (6-sqrt(6))/10 ]\n"
4736 << " [ (6+9*sqrt(6))/35 ]\n"
4737 << " [ 1 ]\n"
4738 << " [ (4-sqrt(6))/10 ]\n"
4739 << " [ (4+sqrt(6))/10 ]\n"
4740 << "A = [ A1 A2 A3 A4 A5 ]\n"
4741 << " A1 = [ (6-sqrt(6))/10 ]\n"
4742 << " [ (-6+5*sqrt(6))/14 ]\n"
4743 << " [ (888+607*sqrt(6))/2850 ]\n"
4744 << " [ (3153-3082*sqrt(6))/14250 ]\n"
4745 << " [ (-32583+14638*sqrt(6))/71250 ]\n"
4746 << " A2 = [ 0 ]\n"
4747 << " [ (6-sqrt(6))/10 ]\n"
4748 << " [ (126-161*sqrt(6))/1425 ]\n"
4749 << " [ (3213+1148*sqrt(6))/28500 ]\n"
4750 << " [ (-17199+364*sqrt(6))/142500 ]\n"
4751 << " A3 = [ 0 ]\n"
4752 << " [ 0 ]\n"
4753 << " [ (6-sqrt(6))/10 ]\n"
4754 << " [ (-267+88*sqrt(6))/500 ]\n"
4755 << " [ (1329-544*sqrt(6))/2500 ]\n"
4756 << " A4 = [ 0 ]\n"
4757 << " [ 0 ]\n"
4758 << " [ 0 ]\n"
4759 << " [ (6-sqrt(6))/10 ]\n"
4760 << " [ (-96+131*sqrt(6))/625 ]\n"
4761 << " A5 = [ 0 ]\n"
4762 << " [ 0 ]\n"
4763 << " [ 0 ]\n"
4764 << " [ 0 ]\n"
4765 << " [ (6-sqrt(6))/10 ]\n"
4766 << "b = [ 0 ]\n"
4767 << " [ 0 ]\n"
4768 << " [ 1/9 ]\n"
4769 << " [ (16-sqrt(6))/36 ]\n"
4770 << " [ (16+sqrt(6))/36 ]'";
4771 return Description.str();
4772 }
4773
4774protected:
4775
4777 {
4778 typedef Teuchos::ScalarTraits<Scalar> ST;
4779 using Teuchos::as;
4780 int NumStages = 5;
4781 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
4782 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
4783 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
4784 const Scalar zero = ST::zero();
4785 const Scalar one = ST::one();
4786 const Scalar sqrt6 = ST::squareroot(as<Scalar>(6*one));
4787 const Scalar gamma = as<Scalar>( (6*one - sqrt6) / (10*one) ); // diagonal
4788
4789 // Fill A:
4790 A(0,0) = gamma;
4791 A(0,1) = zero;
4792 A(0,2) = zero;
4793 A(0,3) = zero;
4794 A(0,4) = zero;
4795
4796 A(1,0) = as<Scalar>( (-6*one+5*one*sqrt6)/(14*one) );
4797 A(1,1) = gamma;
4798 A(1,2) = zero;
4799 A(1,3) = zero;
4800 A(1,4) = zero;
4801
4802 A(2,0) = as<Scalar>( (888*one+607*one*sqrt6)/(2850*one) );
4803 A(2,1) = as<Scalar>( (126*one-161*one*sqrt6)/(1425*one) );
4804 A(2,2) = gamma;
4805 A(2,3) = zero;
4806 A(2,4) = zero;
4807
4808 A(3,0) = as<Scalar>( (3153*one-3082*one*sqrt6)/(14250*one) );
4809 A(3,1) = as<Scalar>( (3213*one+1148*one*sqrt6)/(28500*one) );
4810 A(3,2) = as<Scalar>( (-267*one+88*one*sqrt6)/(500*one) );
4811 A(3,3) = gamma;
4812 A(3,4) = zero;
4813
4814 A(4,0) = as<Scalar>( (-32583*one+14638*one*sqrt6)/(71250*one) );
4815 A(4,1) = as<Scalar>( (-17199*one+364*one*sqrt6)/(142500*one) );
4816 A(4,2) = as<Scalar>( (1329*one-544*one*sqrt6)/(2500*one) );
4817 A(4,3) = as<Scalar>( (-96*one+131*sqrt6)/(625*one) );
4818 A(4,4) = gamma;
4819
4820 // Fill b:
4821 b(0) = zero;
4822 b(1) = zero;
4823 b(2) = as<Scalar>( one/(9*one) );
4824 b(3) = as<Scalar>( (16*one-sqrt6)/(36*one) );
4825 b(4) = as<Scalar>( (16*one+sqrt6)/(36*one) );
4826
4827 // Fill c:
4828 c(0) = gamma;
4829 c(1) = as<Scalar>( (6*one+9*one*sqrt6)/(35*one) );
4830 c(2) = one;
4831 c(3) = as<Scalar>( (4*one-sqrt6)/(10*one) );
4832 c(4) = as<Scalar>( (4*one+sqrt6)/(10*one) );
4833
4834 int order = 5;
4835
4836 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
4837 this->getStepperType(),A,b,c,order,order,order));
4838 }
4839};
4840
4841
4843// ------------------------------------------------------------------------
4844template<class Scalar>
4845Teuchos::RCP<StepperSDIRK_5Stage5thOrder<Scalar> >
4847 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4848 Teuchos::RCP<Teuchos::ParameterList> pl)
4849{
4850 auto stepper = Teuchos::rcp(new StepperSDIRK_5Stage5thOrder<Scalar>());
4851 stepper->setStepperDIRKValues(pl);
4852
4853 if (model != Teuchos::null) {
4854 stepper->setModel(model);
4855 stepper->initialize();
4856 }
4857
4858 return stepper;
4859}
4860
4861
4862// ----------------------------------------------------------------------------
4881template<class Scalar>
4883 virtual public StepperDIRK<Scalar>
4884{
4885public:
4892 {
4893 this->setStepperName("SDIRK 2(1) Pair");
4894 this->setStepperType("SDIRK 2(1) Pair");
4895 this->setupTableau();
4896 this->setupDefault();
4897 this->setUseFSAL( false);
4898 this->setICConsistency( "None");
4899 this->setICConsistencyCheck( false);
4900 }
4901
4903 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
4904 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
4905 bool useFSAL,
4906 std::string ICConsistency,
4907 bool ICConsistencyCheck,
4908 bool useEmbedded,
4909 bool zeroInitialGuess,
4910 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
4911 {
4912 this->setStepperName("SDIRK 2(1) Pair");
4913 this->setStepperType("SDIRK 2(1) Pair");
4914 this->setupTableau();
4915 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
4916 useEmbedded, zeroInitialGuess, stepperRKAppAction);
4917 }
4918
4919 std::string getDescription() const
4920 {
4921 std::ostringstream Description;
4922 Description << this->getStepperType() << "\n"
4923 << "c = [ 1 0 ]'\n"
4924 << "A = [ 1 ]\n"
4925 << " [ -1 1 ]\n"
4926 << "b = [ 1/2 1/2 ]'\n"
4927 << "bstar = [ 1 0 ]'";
4928 return Description.str();
4929 }
4930
4931protected:
4932
4934 {
4935 typedef Teuchos::ScalarTraits<Scalar> ST;
4936 using Teuchos::as;
4937 int NumStages = 2;
4938 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
4939 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
4940 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
4941 Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);
4942
4943 const Scalar one = ST::one();
4944 const Scalar zero = ST::zero();
4945
4946 // Fill A:
4947 A(0,0) = one; A(0,1) = zero;
4948 A(1,0) = -one; A(1,1) = one;
4949
4950 // Fill b:
4951 b(0) = as<Scalar>(one/(2*one));
4952 b(1) = as<Scalar>(one/(2*one));
4953
4954 // Fill c:
4955 c(0) = one;
4956 c(1) = zero;
4957
4958 // Fill bstar
4959 bstar(0) = one;
4960 bstar(1) = zero;
4961 int order = 2;
4962
4963 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
4964 this->getStepperType(),A,b,c,order,order,order,bstar));
4965 }
4966};
4967
4968
4970// ------------------------------------------------------------------------
4971template<class Scalar>
4972Teuchos::RCP<StepperSDIRK_21Pair<Scalar> >
4974 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
4975 Teuchos::RCP<Teuchos::ParameterList> pl)
4976{
4977 auto stepper = Teuchos::rcp(new StepperSDIRK_21Pair<Scalar>());
4978 stepper->setStepperDIRKValues(pl);
4979
4980 if (model != Teuchos::null) {
4981 stepper->setModel(model);
4982 stepper->initialize();
4983 }
4984
4985 return stepper;
4986}
4987
4988
4989// ----------------------------------------------------------------------------
5022template<class Scalar>
5024 virtual public StepperDIRK<Scalar>
5025{
5026public:
5033 {
5034 this->setStepperName("General DIRK");
5035 this->setStepperType("General DIRK");
5036 this->setupTableau();
5037 this->setupDefault();
5038 this->setUseFSAL( false);
5039 this->setICConsistency( "None");
5040 this->setICConsistencyCheck( false);
5041 }
5042
5044 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
5045 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
5046 bool useFSAL,
5047 std::string ICConsistency,
5048 bool ICConsistencyCheck,
5049 bool useEmbedded,
5050 bool zeroInitialGuess,
5051 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
5052 const Teuchos::SerialDenseMatrix<int,Scalar>& A,
5053 const Teuchos::SerialDenseVector<int,Scalar>& b,
5054 const Teuchos::SerialDenseVector<int,Scalar>& c,
5055 const int order,
5056 const int orderMin,
5057 const int orderMax,
5058 const Teuchos::SerialDenseVector<int,Scalar>& bstar)
5059 {
5060 this->setStepperName("General DIRK");
5061 this->setStepperType("General DIRK");
5062 this->setTableau(A,b,c,order,orderMin,orderMax,bstar);
5063
5064 TEUCHOS_TEST_FOR_EXCEPTION(
5065 this->tableau_->isImplicit() != true, std::logic_error,
5066 "Error - General DIRK did not receive a DIRK Butcher Tableau!\n");
5067
5068 this->setup(appModel, solver, useFSAL, ICConsistency, ICConsistencyCheck,
5069 useEmbedded, zeroInitialGuess, stepperRKAppAction);
5070 }
5071
5072 std::string getDescription() const
5073 {
5074 std::stringstream Description;
5075 Description << this->getStepperType() << "\n"
5076 << "The format of the Butcher Tableau parameter list is\n"
5077 << " <Parameter name=\"A\" type=\"string\" value=\"# # # ;\n"
5078 << " # # # ;\n"
5079 << " # # #\"/>\n"
5080 << " <Parameter name=\"b\" type=\"string\" value=\"# # #\"/>\n"
5081 << " <Parameter name=\"c\" type=\"string\" value=\"# # #\"/>\n\n"
5082 << "Note the number of stages is implicit in the number of entries.\n"
5083 << "The number of stages must be consistent.\n"
5084 << "\n"
5085 << "Default tableau is 'SDIRK 2 Stage 2nd order':\n"
5086 << " Computer Methods for ODEs and DAEs\n"
5087 << " U. M. Ascher and L. R. Petzold\n"
5088 << " p. 106\n"
5089 << " gamma = (2-sqrt(2))/2\n"
5090 << " c = [ gamma 1 ]'\n"
5091 << " A = [ gamma 0 ]\n"
5092 << " [ 1-gamma gamma ]\n"
5093 << " b = [ 1-gamma gamma ]'";
5094 return Description.str();
5095 }
5096
5098 {
5099 if (this->tableau_ == Teuchos::null) {
5100 // Set tableau to the default if null, otherwise keep current tableau.
5101 auto stepper = Teuchos::rcp(new StepperSDIRK_2Stage2ndOrder<Scalar>());
5102 auto t = stepper->getTableau();
5103 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
5104 this->getStepperType(),
5105 t->A(),t->b(),t->c(),
5106 t->order(),t->orderMin(),t->orderMax(),
5107 t->bstar()));
5108 this->isInitialized_ = false;
5109 }
5110 }
5111
5112 void setTableau(const Teuchos::SerialDenseMatrix<int,Scalar>& A,
5113 const Teuchos::SerialDenseVector<int,Scalar>& b,
5114 const Teuchos::SerialDenseVector<int,Scalar>& c,
5115 const int order,
5116 const int orderMin,
5117 const int orderMax,
5118 const Teuchos::SerialDenseVector<int,Scalar>&
5119 bstar = Teuchos::SerialDenseVector<int,Scalar>())
5120 {
5121 this->tableau_ = Teuchos::rcp(new RKButcherTableau<Scalar>(
5122 this->getStepperType(),A,b,c,order,orderMin,orderMax,bstar));
5123 this->isInitialized_ = false;
5124 }
5125
5126 Teuchos::RCP<const Teuchos::ParameterList>
5128 {
5129 auto pl = this->getValidParametersBasicDIRK();
5130
5131 // Tableau ParameterList
5132 Teuchos::SerialDenseMatrix<int,Scalar> A = this->tableau_->A();
5133 Teuchos::SerialDenseVector<int,Scalar> b = this->tableau_->b();
5134 Teuchos::SerialDenseVector<int,Scalar> c = this->tableau_->c();
5135 Teuchos::SerialDenseVector<int,Scalar> bstar = this->tableau_->bstar();
5136
5137 Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
5138
5139 std::ostringstream Astream;
5140 Astream.precision(15);
5141 for(int i = 0; i < A.numRows(); i++) {
5142 for(int j = 0; j < A.numCols()-1; j++) {
5143 Astream << A(i,j) << " ";
5144 }
5145 Astream << A(i,A.numCols()-1);
5146 if ( i != A.numRows()-1 ) Astream << "; ";
5147 }
5148 tableauPL->set<std::string>("A", Astream.str());
5149
5150 std::ostringstream bstream;
5151 bstream.precision(15);
5152 for(int i = 0; i < b.length()-1; i++) {
5153 bstream << b(i) << " ";
5154 }
5155 bstream << b(b.length()-1);
5156 tableauPL->set<std::string>("b", bstream.str());
5157
5158 std::ostringstream cstream;
5159 cstream.precision(15);
5160 for(int i = 0; i < c.length()-1; i++) {
5161 cstream << c(i) << " ";
5162 }
5163 cstream << c(c.length()-1);
5164 tableauPL->set<std::string>("c", cstream.str());
5165
5166 tableauPL->set<int>("order", this->tableau_->order());
5167
5168 if ( bstar.length() == 0 ) {
5169 tableauPL->set("bstar", "");
5170 } else {
5171 std::ostringstream bstarstream;
5172 bstarstream.precision(15);
5173 for(int i = 0; i < bstar.length()-1; i++) {
5174 bstarstream << bstar(i) << " ";
5175 }
5176 bstarstream << bstar(bstar.length()-1);
5177 tableauPL->set<std::string>("bstar", bstarstream.str());
5178 }
5179
5180 pl->set("Tableau", *tableauPL);
5181
5182 return pl;
5183 }
5184};
5185
5186
5188// ------------------------------------------------------------------------
5189template<class Scalar>
5190Teuchos::RCP<StepperDIRK_General<Scalar> >
5192 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
5193 Teuchos::RCP<Teuchos::ParameterList> pl)
5194{
5195 auto stepper = Teuchos::rcp(new StepperDIRK_General<Scalar>());
5196 stepper->setStepperDIRKValues(pl);
5197
5198 if (pl != Teuchos::null) {
5199 if (pl->isParameter("Tableau")) {
5200 auto t = stepper->createTableau(pl);
5201 stepper->setTableau( t->A(),t->b(),t->c(),
5202 t->order(),t->orderMin(),t->orderMax(),
5203 t->bstar() );
5204 }
5205 }
5206 TEUCHOS_TEST_FOR_EXCEPTION(
5207 stepper->getTableau()->isDIRK() != true, std::logic_error,
5208 "Error - General DIRK did not receive a DIRK Butcher Tableau!\n");
5209
5210 if (model != Teuchos::null) {
5211 stepper->setModel(model);
5212 stepper->initialize();
5213 }
5214
5215 return stepper;
5216}
5217
5218
5219} // namespace Tempus
5220
5221
5222#endif // Tempus_StepperRKButcherTableau_hpp
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)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
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))
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)
Backward Euler Runge-Kutta Butcher Tableau.
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)
General Implicit 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)
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< const Teuchos::ParameterList > getValidParameters() const
Diagonally Implicit Runge-Kutta (DIRK) time stepper.
virtual void setupDefault()
Default setup for constructor.
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasicDIRK() const
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.
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)
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))
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
RK Trapezoidal Rule (A.K.A. RK Crank-Nicolson)
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)
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_3Stage3rdOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Explicit RK 3/8th Rule Butcher Tableau.
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)
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)
Runge-Kutta 4th order Butcher Tableau.
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)
RK Explicit 5 Stage 3rd order by Kinnmark and Gray.
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)
Explicit RK Bogacki-Shampine Butcher Tableau.
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)
Forward Euler Runge-Kutta Butcher Tableau.
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)
General Explicit Runge-Kutta Butcher Tableau.
StepperERK_General(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, 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, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction=Teuchos::null)
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 >())
virtual std::string getDefaultICConsistency() const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Explicit RK Merson Butcher Tableau.
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)
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)
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)
Strong Stability Preserving Explicit RK Butcher Tableau.
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)
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)
Explicit Runge-Kutta time stepper.
virtual void setup(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Setup for constructor.
virtual void setupDefault()
Default setup for constructor.
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasicERK() const
Teuchos::RCP< RKButcherTableau< Scalar > > tableau_
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)
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
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
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))
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)
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)
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)
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)
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)
Strong Stability Preserving Diagonally-Implicit RK Butcher Tableau.
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)
Strong Stability Preserving Diagonally-Implicit RK Butcher Tableau.
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)
Strong Stability Preserving Diagonally-Implicit RK Butcher Tableau.
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)
Strong Stability Preserving Diagonally-Implicit RK Butcher Tableau.
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)
bool isInitialized_
True if stepper's member data is initialized.
void setICConsistencyCheck(bool c)
void setStepperName(std::string s)
Set the stepper name.
bool useFSAL_
Use First-Same-As-Last (FSAL) principle.
std::string getStepperType() const
Get the stepper type. The stepper type is used as an identifier for the stepper, and can only be set ...
virtual void setUseFSAL(bool a)
void setStepperType(std::string s)
Set the stepper type.
void setICConsistency(std::string s)
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.
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.
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.
Teuchos::RCP< StepperSDIRK_SSPDIRK32< Scalar > > createStepperSDIRK_SSPDIRK32(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
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_Midpoint< Scalar > > createStepperERK_Midpoint(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
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< StepperERK_4Stage4thOrder< Scalar > > createStepperERK_4Stage4thOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< StepperERK_4Stage3rdOrderRunge< Scalar > > createStepperERK_4Stage3rdOrderRunge(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
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.
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.
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.
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< StepperERK_Ralston< Scalar > > createStepperERK_Ralston(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
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.
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< StepperDIRK_1StageTheta< Scalar > > createStepperDIRK_1StageTheta(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
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_2Stage3rdOrder< Scalar > > createStepperSDIRK_2Stage3rdOrder(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
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.
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.
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.
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< StepperSDIRK_ImplicitMidpoint< Scalar > > createStepperSDIRK_ImplicitMidpoint(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
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.
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.
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.
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.
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_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_3Stage3rdOrder< Scalar > > createStepperERK_3Stage3rdOrder(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.
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< StepperDIRK_General< Scalar > > createStepperDIRK_General(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
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.