Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_RKButcherTableau.hpp
1//@HEADER
2// ***********************************************************************
3//
4// Rythmos Package
5// Copyright (2006) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25//
26// ***********************************************************************
27//@HEADER
28
29
30#ifndef RYTHMOS_RK_BUTCHER_TABLEAU_HPP
31#define RYTHMOS_RK_BUTCHER_TABLEAU_HPP
32
33// disable clang warnings
34#ifdef __clang__
35#pragma clang system_header
36#endif
37
38#include "Rythmos_Types.hpp"
39#include "Rythmos_RKButcherTableauBase.hpp"
40
41#include "Teuchos_Assert.hpp"
42#include "Teuchos_as.hpp"
43#include "Teuchos_StandardParameterEntryValidators.hpp"
44#include "Teuchos_Describable.hpp"
45#include "Teuchos_VerboseObject.hpp"
46#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
47#include "Teuchos_ParameterListAcceptor.hpp"
48#include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
49
50#include "Thyra_ProductVectorBase.hpp"
51
52namespace Rythmos {
53
54 using Teuchos::as;
55
56 inline const std::string RKBT_ForwardEuler_name() { return "Forward Euler"; } // done
57 inline const std::string RKBT_BackwardEuler_name() { return "Backward Euler"; } // done
58 inline const std::string Explicit4Stage_name() { return "Explicit 4 Stage"; } // done
59 inline const std::string Explicit3_8Rule_name() { return "Explicit 3/8 Rule"; } // done
60
61 inline const std::string ExplicitTrapezoidal_name() { return "Explicit Trapezoidal"; } // done
62 inline const std::string Explicit2Stage2ndOrderRunge_name() { return "Explicit 2 Stage 2nd order by Runge"; } // done
63 inline const std::string Explicit2Stage2ndOrderTVD_name() { return "Explicit 2 Stage 2nd order TVD"; } // done
64 inline const std::string Explicit3Stage3rdOrderHeun_name() { return "Explicit 3 Stage 3rd order by Heun"; } // done
65 inline const std::string Explicit3Stage3rdOrder_name() { return "Explicit 3 Stage 3rd order"; } // done
66 inline const std::string Explicit3Stage3rdOrderTVD_name() { return "Explicit 3 Stage 3rd order TVD"; } // done
67 inline const std::string Explicit4Stage3rdOrderRunge_name() { return "Explicit 4 Stage 3rd order by Runge"; } // done
68 inline const std::string Explicit5Stage3rdOrderKandG_name() { return "Explicit 5 Stage 3rd order by Kinnmark and Gray"; } // done
69
70 inline const std::string IRK1StageTheta_name() { return "IRK 1 Stage Theta Method"; } // done
71 inline const std::string IRK2StageTheta_name() { return "IRK 2 Stage Theta Method"; } // done
72 inline const std::string Implicit1Stage2ndOrderGauss_name() { return "Implicit 1 Stage 2nd order Gauss"; } // done
73 inline const std::string Implicit2Stage4thOrderGauss_name() { return "Implicit 2 Stage 4th order Gauss"; } // done
74 inline const std::string Implicit3Stage6thOrderGauss_name() { return "Implicit 3 Stage 6th order Gauss"; } // done
75
76 inline const std::string Implicit1Stage1stOrderRadauA_name() { return "Implicit 1 Stage 1st order Radau left"; } // done
77 inline const std::string Implicit2Stage3rdOrderRadauA_name() { return "Implicit 2 Stage 3rd order Radau left"; } // done
78 inline const std::string Implicit3Stage5thOrderRadauA_name() { return "Implicit 3 Stage 5th order Radau left"; } // done
79
80 inline const std::string Implicit1Stage1stOrderRadauB_name() { return "Implicit 1 Stage 1st order Radau right"; } // done
81 inline const std::string Implicit2Stage3rdOrderRadauB_name() { return "Implicit 2 Stage 3rd order Radau right"; } // done
82 inline const std::string Implicit3Stage5thOrderRadauB_name() { return "Implicit 3 Stage 5th order Radau right"; } // done
83
84 inline const std::string Implicit2Stage2ndOrderLobattoA_name() { return "Implicit 2 Stage 2nd order Lobatto A"; } // done
85 inline const std::string Implicit3Stage4thOrderLobattoA_name() { return "Implicit 3 Stage 4th order Lobatto A"; } // done
86 inline const std::string Implicit4Stage6thOrderLobattoA_name() { return "Implicit 4 Stage 6th order Lobatto A"; } // done
87
88 inline const std::string Implicit2Stage2ndOrderLobattoB_name() { return "Implicit 2 Stage 2nd order Lobatto B"; } // done
89 inline const std::string Implicit3Stage4thOrderLobattoB_name() { return "Implicit 3 Stage 4th order Lobatto B"; } // done
90 inline const std::string Implicit4Stage6thOrderLobattoB_name() { return "Implicit 4 Stage 6th order Lobatto B"; } // done
91
92 inline const std::string Implicit2Stage2ndOrderLobattoC_name() { return "Implicit 2 Stage 2nd order Lobatto C"; } // done
93 inline const std::string Implicit3Stage4thOrderLobattoC_name() { return "Implicit 3 Stage 4th order Lobatto C"; } // done
94 inline const std::string Implicit4Stage6thOrderLobattoC_name() { return "Implicit 4 Stage 6th order Lobatto C"; } // done
95
96 inline const std::string Implicit2Stage4thOrderHammerHollingsworth_name() { return "Implicit 2 Stage 4th Order Hammer & Hollingsworth"; } // done
97 inline const std::string Implicit3Stage6thOrderKuntzmannButcher_name() { return "Implicit 3 Stage 6th Order Kuntzmann & Butcher"; } // done
98 inline const std::string Implicit4Stage8thOrderKuntzmannButcher_name() { return "Implicit 4 Stage 8th Order Kuntzmann & Butcher"; } // done
99
100 inline const std::string DIRK2Stage3rdOrder_name() { return "Diagonal IRK 2 Stage 3rd order"; } // done
101
102 inline const std::string SDIRK2Stage2ndOrder_name() { return "Singly Diagonal IRK 2 Stage 2nd order"; } // done
103 inline const std::string SDIRK2Stage3rdOrder_name() { return "Singly Diagonal IRK 2 Stage 3rd order"; } // done
104 inline const std::string SDIRK5Stage5thOrder_name() { return "Singly Diagonal IRK 5 Stage 5th order"; } // done
105 inline const std::string SDIRK5Stage4thOrder_name() { return "Singly Diagonal IRK 5 Stage 4th order"; } // done
106 inline const std::string SDIRK3Stage4thOrder_name() { return "Singly Diagonal IRK 3 Stage 4th order"; } // done
107
108template<class Scalar>
109class RKButcherTableauDefaultBase :
110 virtual public RKButcherTableauBase<Scalar>,
111 virtual public Teuchos::ParameterListAcceptorDefaultBase
112{
113 public:
115 virtual int numStages() const { return A_.numRows(); }
117 virtual const Teuchos::SerialDenseMatrix<int,Scalar>& A() const { return A_; }
119 virtual const Teuchos::SerialDenseVector<int,Scalar>& b() const { return b_; }
121 virtual const Teuchos::SerialDenseVector<int,Scalar>& bhat() const { return bhat_ ; }
123 virtual const Teuchos::SerialDenseVector<int,Scalar>& c() const { return c_; }
125 virtual int order() const { return order_; }
127 virtual bool isEmbeddedMethod() const { return isEmbedded_; } // returns whether the stepper is Embedded or not (Sidafa)
129 virtual void setDescription(std::string longDescription) { longDescription_ = longDescription; }
130
132 virtual void initialize(
133 const Teuchos::SerialDenseMatrix<int,Scalar>& A_in,
134 const Teuchos::SerialDenseVector<int,Scalar>& b_in,
135 const Teuchos::SerialDenseVector<int,Scalar>& c_in,
136 const int order_in,
137 const std::string& longDescription_in,
138 bool isEmbedded = false, /* (default) tell the stepper the RK is an embedded method */
139 const Teuchos::SerialDenseVector<int,Scalar>& bhat_in = Teuchos::SerialDenseVector<int,Scalar>() /* (default) */
140 )
141 {
142 const int numStages_in = A_in.numRows();
143 TEUCHOS_ASSERT_EQUALITY( A_in.numRows(), numStages_in );
144 TEUCHOS_ASSERT_EQUALITY( A_in.numCols(), numStages_in );
145 TEUCHOS_ASSERT_EQUALITY( b_in.length(), numStages_in );
146 TEUCHOS_ASSERT_EQUALITY( c_in.length(), numStages_in );
147 TEUCHOS_ASSERT( order_in > 0 );
148 A_ = A_in;
149 b_ = b_in;
150 c_ = c_in;
151 order_ = order_in;
152 longDescription_ = longDescription_in;
153
154 /* Sidafa */
155 if (isEmbedded) {
156 TEUCHOS_ASSERT_EQUALITY( bhat_in.length(), numStages_in );
157 bhat_ = bhat_in;
158 isEmbedded_ = true;
159 }
160 }
161
162 /* \brief Redefined from Teuchos::ParameterListAcceptorDefaultBase */
164
166 virtual void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
167 {
168 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
169 paramList->validateParameters(*this->getValidParameters());
170 Teuchos::readVerboseObjectSublist(&*paramList,this);
171 setMyParamList(paramList);
172 }
173
175 virtual RCP<const Teuchos::ParameterList> getValidParameters() const
176 {
177 if (is_null(validPL_)) {
178 validPL_ = Teuchos::parameterList();
179 validPL_->set("Description","",this->getMyDescription());
180 Teuchos::setupVerboseObjectSublist(&*validPL_);
181 }
182 return validPL_;
183 }
184
186
187 /* \brief Redefined from Teuchos::Describable */
189
191 virtual std::string description() const { return "Rythmos::RKButcherTableauDefaultBase"; }
192
194 virtual void describe(
195 Teuchos::FancyOStream &out,
196 const Teuchos::EVerbosityLevel verbLevel
197 ) const
198 {
199 if (verbLevel != Teuchos::VERB_NONE) {
200 out << this->description() << std::endl;
201 out << this->getMyDescription() << std::endl;
202 out << "number of Stages = " << this->numStages() << std::endl;
203 out << "A = " << printMat(this->A()) << std::endl;
204 out << "b = " << printMat(this->b()) << std::endl;
205 out << "c = " << printMat(this->c()) << std::endl;
206 out << "order = " << this->order() << std::endl;
207 }
208 }
209
211
212 protected:
213 void setMyDescription(std::string longDescription) { longDescription_ = longDescription; }
214 const std::string& getMyDescription() const { return longDescription_; }
215
216 void setMy_A(const Teuchos::SerialDenseMatrix<int,Scalar>& new_A) { A_ = new_A; }
217 void setMy_b(const Teuchos::SerialDenseVector<int,Scalar>& new_b) { b_ = new_b; }
218 void setMy_c(const Teuchos::SerialDenseVector<int,Scalar>& new_c) { c_ = new_c; }
219 void setMy_order(const int& new_order) { order_ = new_order; }
220
221 void setMyValidParameterList( const RCP<ParameterList> validPL ) { validPL_ = validPL; }
222 RCP<ParameterList> getMyNonconstValidParameterList() { return validPL_; }
223
224 private:
225 Teuchos::SerialDenseMatrix<int,Scalar> A_;
226 Teuchos::SerialDenseVector<int,Scalar> b_;
227 Teuchos::SerialDenseVector<int,Scalar> c_;
228 int order_;
229 std::string longDescription_;
230 mutable RCP<ParameterList> validPL_;
231
232 /* Sidafa - Embedded method parameters */
233 Teuchos::SerialDenseVector<int,Scalar> bhat_;
234 bool isEmbedded_ = false;
235};
236
237
238// Nonmember constructor
239template<class Scalar>
240RCP<RKButcherTableauBase<Scalar> > rKButcherTableau()
241{
242 return(rcp(new RKButcherTableauDefaultBase<Scalar>()));
243}
244
245// Nonmember constructor
246template<class Scalar>
247RCP<RKButcherTableauBase<Scalar> > rKButcherTableau(
248 const Teuchos::SerialDenseMatrix<int,Scalar>& A_in,
249 const Teuchos::SerialDenseVector<int,Scalar>& b_in,
250 const Teuchos::SerialDenseVector<int,Scalar>& c_in,
251 int order_in,
252 const std::string& description_in = ""
253 )
254{
255 RCP<RKButcherTableauDefaultBase<Scalar> > rkbt = rcp(new RKButcherTableauDefaultBase<Scalar>());
256 rkbt->initialize(A_in,b_in,c_in,order_in,description_in);
257 return(rkbt);
258}
259
260
261template<class Scalar>
262class BackwardEuler_RKBT :
263 virtual public RKButcherTableauDefaultBase<Scalar>
264{
265 public:
266 BackwardEuler_RKBT()
267 {
268 std::ostringstream myDescription;
269 myDescription << RKBT_BackwardEuler_name() << "\n"
270 << "c = [ 1 ]'\n"
271 << "A = [ 1 ]\n"
272 << "b = [ 1 ]'" << std::endl;
273 typedef ScalarTraits<Scalar> ST;
274 Teuchos::SerialDenseMatrix<int,Scalar> myA(1,1);
275 myA(0,0) = ST::one();
276 Teuchos::SerialDenseVector<int,Scalar> myb(1);
277 myb(0) = ST::one();
278 Teuchos::SerialDenseVector<int,Scalar> myc(1);
279 myc(0) = ST::one();
280
281 this->setMyDescription(myDescription.str());
282 this->setMy_A(myA);
283 this->setMy_b(myb);
284 this->setMy_c(myc);
285 this->setMy_order(1);
286 }
287};
288
289
290
291template<class Scalar>
292class ForwardEuler_RKBT :
293 virtual public RKButcherTableauDefaultBase<Scalar>
294{
295 public:
296
297 ForwardEuler_RKBT()
298 {
299 std::ostringstream myDescription;
300 myDescription << RKBT_ForwardEuler_name() << "\n"
301 << "c = [ 0 ]'\n"
302 << "A = [ 0 ]\n"
303 << "b = [ 1 ]'" << std::endl;
304 typedef ScalarTraits<Scalar> ST;
305 Teuchos::SerialDenseMatrix<int,Scalar> myA(1,1);
306 Teuchos::SerialDenseVector<int,Scalar> myb(1);
307 myb(0) = ST::one();
308 Teuchos::SerialDenseVector<int,Scalar> myc(1);
309
310 this->setMyDescription(myDescription.str());
311 this->setMy_A(myA);
312 this->setMy_b(myb);
313 this->setMy_c(myc);
314 this->setMy_order(1);
315 }
316};
317
318
319template<class Scalar>
320class Explicit4Stage4thOrder_RKBT :
321 virtual public RKButcherTableauDefaultBase<Scalar>
322{
323 public:
324 Explicit4Stage4thOrder_RKBT()
325 {
326 std::ostringstream myDescription;
327 myDescription << Explicit4Stage_name() << "\n"
328 << "\"The\" Runge-Kutta Method (explicit):\n"
329 << "Solving Ordinary Differential Equations I:\n"
330 << "Nonstiff Problems, 2nd Revised Edition\n"
331 << "E. Hairer, S.P. Norsett, G. Wanner\n"
332 << "Table 1.2, pg 138\n"
333 << "c = [ 0 1/2 1/2 1 ]'\n"
334 << "A = [ 0 ] \n"
335 << " [ 1/2 0 ]\n"
336 << " [ 0 1/2 0 ]\n"
337 << " [ 0 0 1 0 ]\n"
338 << "b = [ 1/6 1/3 1/3 1/6 ]'" << std::endl;
339 typedef ScalarTraits<Scalar> ST;
340 Scalar one = ST::one();
341 Scalar zero = ST::zero();
342 Scalar onehalf = ST::one()/(2*ST::one());
343 Scalar onesixth = ST::one()/(6*ST::one());
344 Scalar onethird = ST::one()/(3*ST::one());
345
346 int myNumStages = 4;
347 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
348 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
349 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
350
351 // Fill A:
352 myA(0,0) = zero;
353 myA(0,1) = zero;
354 myA(0,2) = zero;
355 myA(0,3) = zero;
356
357 myA(1,0) = onehalf;
358 myA(1,1) = zero;
359 myA(1,2) = zero;
360 myA(1,3) = zero;
361
362 myA(2,0) = zero;
363 myA(2,1) = onehalf;
364 myA(2,2) = zero;
365 myA(2,3) = zero;
366
367 myA(3,0) = zero;
368 myA(3,1) = zero;
369 myA(3,2) = one;
370 myA(3,3) = zero;
371
372 // Fill myb:
373 myb(0) = onesixth;
374 myb(1) = onethird;
375 myb(2) = onethird;
376 myb(3) = onesixth;
377
378 // fill b_c_
379 myc(0) = zero;
380 myc(1) = onehalf;
381 myc(2) = onehalf;
382 myc(3) = one;
383
384 this->setMyDescription(myDescription.str());
385 this->setMy_A(myA);
386 this->setMy_b(myb);
387 this->setMy_c(myc);
388 this->setMy_order(4);
389 }
390};
391
392
393template<class Scalar>
394class Explicit3_8Rule_RKBT :
395 virtual public RKButcherTableauDefaultBase<Scalar>
396{
397 public:
398 Explicit3_8Rule_RKBT()
399 {
400
401 std::ostringstream myDescription;
402 myDescription << Explicit3_8Rule_name() << "\n"
403 << "Solving Ordinary Differential Equations I:\n"
404 << "Nonstiff Problems, 2nd Revised Edition\n"
405 << "E. Hairer, S.P. Norsett, G. Wanner\n"
406 << "Table 1.2, pg 138\n"
407 << "c = [ 0 1/3 2/3 1 ]'\n"
408 << "A = [ 0 ]\n"
409 << " [ 1/3 0 ]\n"
410 << " [-1/3 1 0 ]\n"
411 << " [ 1 -1 1 0 ]\n"
412 << "b = [ 1/8 3/8 3/8 1/8 ]'" << std::endl;
413 typedef ScalarTraits<Scalar> ST;
414 int myNumStages = 4;
415 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
416 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
417 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
418
419 Scalar one = ST::one();
420 Scalar zero = ST::zero();
421 Scalar one_third = as<Scalar>(ST::one()/(3*ST::one()));
422 Scalar two_third = as<Scalar>(2*ST::one()/(3*ST::one()));
423 Scalar one_eighth = as<Scalar>(ST::one()/(8*ST::one()));
424 Scalar three_eighth = as<Scalar>(3*ST::one()/(8*ST::one()));
425
426 // Fill myA:
427 myA(0,0) = zero;
428 myA(0,1) = zero;
429 myA(0,2) = zero;
430 myA(0,3) = zero;
431
432 myA(1,0) = one_third;
433 myA(1,1) = zero;
434 myA(1,2) = zero;
435 myA(1,3) = zero;
436
437 myA(2,0) = as<Scalar>(-one_third);
438 myA(2,1) = one;
439 myA(2,2) = zero;
440 myA(2,3) = zero;
441
442 myA(3,0) = one;
443 myA(3,1) = as<Scalar>(-one);
444 myA(3,2) = one;
445 myA(3,3) = zero;
446
447 // Fill myb:
448 myb(0) = one_eighth;
449 myb(1) = three_eighth;
450 myb(2) = three_eighth;
451 myb(3) = one_eighth;
452
453 // Fill myc:
454 myc(0) = zero;
455 myc(1) = one_third;
456 myc(2) = two_third;
457 myc(3) = one;
458
459 this->setMyDescription(myDescription.str());
460 this->setMy_A(myA);
461 this->setMy_b(myb);
462 this->setMy_c(myc);
463 this->setMy_order(4);
464 }
465};
466
467
468template<class Scalar>
469class Explicit4Stage3rdOrderRunge_RKBT :
470 virtual public RKButcherTableauDefaultBase<Scalar>
471{
472 public:
473 Explicit4Stage3rdOrderRunge_RKBT()
474 {
475
476 std::ostringstream myDescription;
477 myDescription << Explicit4Stage3rdOrderRunge_name() << "\n"
478 << "Solving Ordinary Differential Equations I:\n"
479 << "Nonstiff Problems, 2nd Revised Edition\n"
480 << "E. Hairer, S.P. Norsett, G. Wanner\n"
481 << "Table 1.1, pg 135\n"
482 << "c = [ 0 1/2 1 1 ]'\n"
483 << "A = [ 0 ]\n"
484 << " [ 1/2 0 ]\n"
485 << " [ 0 1 0 ]\n"
486 << " [ 0 0 1 0 ]\n"
487 << "b = [ 1/6 2/3 0 1/6 ]'" << std::endl;
488 typedef ScalarTraits<Scalar> ST;
489 int myNumStages = 4;
490 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
491 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
492 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
493
494 Scalar one = ST::one();
495 Scalar onehalf = ST::one()/(2*ST::one());
496 Scalar onesixth = ST::one()/(6*ST::one());
497 Scalar twothirds = 2*ST::one()/(3*ST::one());
498 Scalar zero = ST::zero();
499
500 // Fill A:
501 myA(0,0) = zero;
502 myA(0,1) = zero;
503 myA(0,2) = zero;
504 myA(0,3) = zero;
505
506 myA(1,0) = onehalf;
507 myA(1,1) = zero;
508 myA(1,2) = zero;
509 myA(1,3) = zero;
510
511 myA(2,0) = zero;
512 myA(2,1) = one;
513 myA(2,2) = zero;
514 myA(2,3) = zero;
515
516 myA(3,0) = zero;
517 myA(3,1) = zero;
518 myA(3,2) = one;
519 myA(3,3) = zero;
520
521 // Fill b:
522 myb(0) = onesixth;
523 myb(1) = twothirds;
524 myb(2) = zero;
525 myb(3) = onesixth;
526
527 // Fill myc:
528 myc(0) = zero;
529 myc(1) = onehalf;
530 myc(2) = one;
531 myc(3) = one;
532
533 this->setMyDescription(myDescription.str());
534 this->setMy_A(myA);
535 this->setMy_b(myb);
536 this->setMy_c(myc);
537 this->setMy_order(3);
538 }
539};
540
541template<class Scalar>
542class Explicit5Stage3rdOrderKandG_RKBT :
543 virtual public RKButcherTableauDefaultBase<Scalar>
544{
545 public:
546 Explicit5Stage3rdOrderKandG_RKBT()
547 {
548
549 std::ostringstream myDescription;
550 myDescription << Explicit5Stage3rdOrderKandG_name() << "\n"
551 << "Kinnmark & Gray 5 stage, 3rd order scheme \n"
552 << "Modified by P. Ullrich. From the prim_advance_mod.F90 \n"
553 << "routine in the HOMME atmosphere model code.\n"
554 << "c = [ 0 1/5 1/5 1/3 2/3 ]'\n"
555 << "A = [ 0 ]\n"
556 << " [ 1/5 0 ]\n"
557 << " [ 0 1/5 0 ]\n"
558 << " [ 0 0 1/3 0 ]\n"
559 << " [ 0 0 0 2/3 0 ]\n"
560 << "b = [ 1/4 0 0 0 3/4 ]'" << std::endl;
561 typedef ScalarTraits<Scalar> ST;
562 int myNumStages = 5;
563 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
564 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
565 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
566
567 Scalar onefifth = ST::one()/(5*ST::one());
568 Scalar onefourth = ST::one()/(4*ST::one());
569 Scalar onethird = ST::one()/(3*ST::one());
570 Scalar twothirds = 2*ST::one()/(3*ST::one());
571 Scalar threefourths = 3*ST::one()/(4*ST::one());
572 Scalar zero = ST::zero();
573
574 // Fill A:
575 myA(0,0) = zero;
576 myA(0,1) = zero;
577 myA(0,2) = zero;
578 myA(0,3) = zero;
579 myA(0,4) = zero;
580
581 myA(1,0) = onefifth;
582 myA(1,1) = zero;
583 myA(1,2) = zero;
584 myA(1,3) = zero;
585 myA(1,4) = zero;
586
587 myA(2,0) = zero;
588 myA(2,1) = onefifth;
589 myA(2,2) = zero;
590 myA(2,3) = zero;
591 myA(2,4) = zero;
592
593 myA(3,0) = zero;
594 myA(3,1) = zero;
595 myA(3,2) = onethird;
596 myA(3,3) = zero;
597 myA(3,4) = zero;
598
599 myA(4,0) = zero;
600 myA(4,1) = zero;
601 myA(4,2) = zero;
602 myA(4,3) = twothirds;
603 myA(4,4) = zero;
604
605 // Fill b:
606 myb(0) = onefourth;
607 myb(1) = zero;
608 myb(2) = zero;
609 myb(3) = zero;
610 myb(4) = threefourths;
611
612 // Fill myc:
613 myc(0) = zero;
614 myc(1) = onefifth;
615 myc(2) = onefifth;
616 myc(3) = onethird;
617 myc(4) = twothirds;
618
619 this->setMyDescription(myDescription.str());
620 this->setMy_A(myA);
621 this->setMy_b(myb);
622 this->setMy_c(myc);
623 this->setMy_order(3);
624 }
625};
626
627
628template<class Scalar>
629class Explicit3Stage3rdOrder_RKBT :
630 virtual public RKButcherTableauDefaultBase<Scalar>
631{
632 public:
633 Explicit3Stage3rdOrder_RKBT()
634 {
635
636 std::ostringstream myDescription;
637 myDescription << Explicit3Stage3rdOrder_name() << "\n"
638 << "c = [ 0 1/2 1 ]'\n"
639 << "A = [ 0 ]\n"
640 << " [ 1/2 0 ]\n"
641 << " [ -1 2 0 ]\n"
642 << "b = [ 1/6 4/6 1/6 ]'" << std::endl;
643 typedef ScalarTraits<Scalar> ST;
644 Scalar one = ST::one();
645 Scalar two = as<Scalar>(2*ST::one());
646 Scalar zero = ST::zero();
647 Scalar onehalf = ST::one()/(2*ST::one());
648 Scalar onesixth = ST::one()/(6*ST::one());
649 Scalar foursixth = 4*ST::one()/(6*ST::one());
650
651 int myNumStages = 3;
652 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
653 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
654 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
655
656 // Fill myA:
657 myA(0,0) = zero;
658 myA(0,1) = zero;
659 myA(0,2) = zero;
660
661 myA(1,0) = onehalf;
662 myA(1,1) = zero;
663 myA(1,2) = zero;
664
665 myA(2,0) = -one;
666 myA(2,1) = two;
667 myA(2,2) = zero;
668
669 // Fill myb:
670 myb(0) = onesixth;
671 myb(1) = foursixth;
672 myb(2) = onesixth;
673
674 // fill b_c_
675 myc(0) = zero;
676 myc(1) = onehalf;
677 myc(2) = one;
678
679 this->setMyDescription(myDescription.str());
680 this->setMy_A(myA);
681 this->setMy_b(myb);
682 this->setMy_c(myc);
683 this->setMy_order(3);
684 }
685};
686
687template<class Scalar>
688class Explicit2Stage2ndOrderTVD_RKBT :
689 virtual public RKButcherTableauDefaultBase<Scalar>
690{
691 public:
692 Explicit2Stage2ndOrderTVD_RKBT()
693 {
694
695 std::ostringstream myDescription;
696 myDescription << Explicit2Stage2ndOrderTVD_name() << "\n"
697 << "Sigal Gottlieb, David Ketcheson and Chi-Wang Shu\n"
698 << "`Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations'\n"
699 << "World Scientific, 2011\n"
700 << "pp. 15\n"
701 << "c = [ 0 1 ]'\n"
702 << "A = [ 0 ]\n"
703 << " [ 1 0 ]\n"
704 << "b = [ 1/2 1/2]'\n"
705 << "This is also written in the following set of updates.\n"
706 << "u1 = u^n + dt L(u^n)\n"
707 << "u^(n+1) = u^n/2 + u1/2 + dt L(u1)/2"
708 << std::endl;
709 typedef ScalarTraits<Scalar> ST;
710 Scalar one = ST::one();
711 Scalar zero = ST::zero();
712 Scalar onehalf = ST::one()/(2*ST::one());
713
714 int myNumStages = 2;
715 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
716 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
717 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
718
719 // Fill myA:
720 myA(0,0) = zero;
721 myA(0,1) = zero;
722
723 myA(1,0) = one;
724 myA(1,1) = zero;
725
726 // Fill myb:
727 myb(0) = onehalf;
728 myb(1) = onehalf;
729
730 // fill b_c_
731 myc(0) = zero;
732 myc(1) = one;
733
734 this->setMyDescription(myDescription.str());
735 this->setMy_A(myA);
736 this->setMy_b(myb);
737 this->setMy_c(myc);
738 this->setMy_order(2);
739 }
740};
741
742template<class Scalar>
743class Explicit3Stage3rdOrderTVD_RKBT :
744 virtual public RKButcherTableauDefaultBase<Scalar>
745{
746 public:
747 Explicit3Stage3rdOrderTVD_RKBT()
748 {
749
750 std::ostringstream myDescription;
751 myDescription << Explicit3Stage3rdOrderTVD_name() << "\n"
752 << "Sigal Gottlieb and Chi-Wang Shu\n"
753 << "`Total Variation Diminishing Runge-Kutta Schemes'\n"
754 << "Mathematics of Computation\n"
755 << "Volume 67, Number 221, January 1998, pp. 73-85\n"
756 << "c = [ 0 1 1/2 ]'\n"
757 << "A = [ 0 ]\n"
758 << " [ 1 0 ]\n"
759 << " [ 1/4 1/4 0 ]\n"
760 << "b = [ 1/6 1/6 4/6 ]'\n"
761 << "This is also written in the following set of updates.\n"
762 << "u1 = u^n + dt L(u^n)\n"
763 << "u2 = 3 u^n/4 + u1/4 + dt L(u1)/4\n"
764 << "u^(n+1) = u^n/3 + 2 u2/2 + 2 dt L(u2)/3"
765 << std::endl;
766 typedef ScalarTraits<Scalar> ST;
767 Scalar one = ST::one();
768 Scalar zero = ST::zero();
769 Scalar onehalf = ST::one()/(2*ST::one());
770 Scalar onefourth = ST::one()/(4*ST::one());
771 Scalar onesixth = ST::one()/(6*ST::one());
772 Scalar foursixth = 4*ST::one()/(6*ST::one());
773
774 int myNumStages = 3;
775 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
776 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
777 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
778
779 // Fill myA:
780 myA(0,0) = zero;
781 myA(0,1) = zero;
782 myA(0,2) = zero;
783
784 myA(1,0) = one;
785 myA(1,1) = zero;
786 myA(1,2) = zero;
787
788 myA(2,0) = onefourth;
789 myA(2,1) = onefourth;
790 myA(2,2) = zero;
791
792 // Fill myb:
793 myb(0) = onesixth;
794 myb(1) = onesixth;
795 myb(2) = foursixth;
796
797 // fill b_c_
798 myc(0) = zero;
799 myc(1) = one;
800 myc(2) = onehalf;
801
802 this->setMyDescription(myDescription.str());
803 this->setMy_A(myA);
804 this->setMy_b(myb);
805 this->setMy_c(myc);
806 this->setMy_order(3);
807 }
808};
809
810
811template<class Scalar>
812class Explicit3Stage3rdOrderHeun_RKBT :
813 virtual public RKButcherTableauDefaultBase<Scalar>
814{
815 public:
816 Explicit3Stage3rdOrderHeun_RKBT()
817 {
818 std::ostringstream myDescription;
819 myDescription << Explicit3Stage3rdOrderHeun_name() << "\n"
820 << "Solving Ordinary Differential Equations I:\n"
821 << "Nonstiff Problems, 2nd Revised Edition\n"
822 << "E. Hairer, S.P. Norsett, G. Wanner\n"
823 << "Table 1.1, pg 135\n"
824 << "c = [ 0 1/3 2/3 ]'\n"
825 << "A = [ 0 ] \n"
826 << " [ 1/3 0 ]\n"
827 << " [ 0 2/3 0 ]\n"
828 << "b = [ 1/4 0 3/4 ]'" << std::endl;
829 typedef ScalarTraits<Scalar> ST;
830 Scalar one = ST::one();
831 Scalar zero = ST::zero();
832 Scalar onethird = one/(3*one);
833 Scalar twothirds = 2*one/(3*one);
834 Scalar onefourth = one/(4*one);
835 Scalar threefourths = 3*one/(4*one);
836
837 int myNumStages = 3;
838 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
839 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
840 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
841
842 // Fill myA:
843 myA(0,0) = zero;
844 myA(0,1) = zero;
845 myA(0,2) = zero;
846
847 myA(1,0) = onethird;
848 myA(1,1) = zero;
849 myA(1,2) = zero;
850
851 myA(2,0) = zero;
852 myA(2,1) = twothirds;
853 myA(2,2) = zero;
854
855 // Fill myb:
856 myb(0) = onefourth;
857 myb(1) = zero;
858 myb(2) = threefourths;
859
860 // fill b_c_
861 myc(0) = zero;
862 myc(1) = onethird;
863 myc(2) = twothirds;
864
865 this->setMyDescription(myDescription.str());
866 this->setMy_A(myA);
867 this->setMy_b(myb);
868 this->setMy_c(myc);
869 this->setMy_order(3);
870 }
871};
872
873
874template<class Scalar>
875class Explicit2Stage2ndOrderRunge_RKBT :
876 virtual public RKButcherTableauDefaultBase<Scalar>
877{
878 public:
879 Explicit2Stage2ndOrderRunge_RKBT()
880 {
881 std::ostringstream myDescription;
882 myDescription << Explicit2Stage2ndOrderRunge_name() << "\n"
883 << "Also known as Explicit Midpoint\n"
884 << "Solving Ordinary Differential Equations I:\n"
885 << "Nonstiff Problems, 2nd Revised Edition\n"
886 << "E. Hairer, S.P. Norsett, G. Wanner\n"
887 << "Table 1.1, pg 135\n"
888 << "c = [ 0 1/2 ]'\n"
889 << "A = [ 0 ]\n"
890 << " [ 1/2 0 ]\n"
891 << "b = [ 0 1 ]'" << std::endl;
892 typedef ScalarTraits<Scalar> ST;
893 Scalar one = ST::one();
894 Scalar zero = ST::zero();
895 Scalar onehalf = ST::one()/(2*ST::one());
896
897 int myNumStages = 2;
898 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
899 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
900 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
901
902 // Fill myA:
903 myA(0,0) = zero;
904 myA(0,1) = zero;
905
906 myA(1,0) = onehalf;
907 myA(1,1) = zero;
908
909 // Fill myb:
910 myb(0) = zero;
911 myb(1) = one;
912
913 // fill b_c_
914 myc(0) = zero;
915 myc(1) = onehalf;
916
917 this->setMyDescription(myDescription.str());
918 this->setMy_A(myA);
919 this->setMy_b(myb);
920 this->setMy_c(myc);
921 this->setMy_order(2);
922 }
923};
924
925
926template<class Scalar>
927class ExplicitTrapezoidal_RKBT :
928 virtual public RKButcherTableauDefaultBase<Scalar>
929{
930 public:
931 ExplicitTrapezoidal_RKBT()
932 {
933 std::ostringstream myDescription;
934 myDescription << ExplicitTrapezoidal_name() << "\n"
935 << "c = [ 0 1 ]'\n"
936 << "A = [ 0 ]\n"
937 << " [ 1 0 ]\n"
938 << "b = [ 1/2 1/2 ]'" << std::endl;
939 typedef ScalarTraits<Scalar> ST;
940 Scalar one = ST::one();
941 Scalar zero = ST::zero();
942 Scalar onehalf = ST::one()/(2*ST::one());
943
944 int myNumStages = 2;
945 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
946 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
947 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
948
949 // Fill myA:
950 myA(0,0) = zero;
951 myA(0,1) = zero;
952
953 myA(1,0) = one;
954 myA(1,1) = zero;
955
956 // Fill myb:
957 myb(0) = onehalf;
958 myb(1) = onehalf;
959
960 // fill b_c_
961 myc(0) = zero;
962 myc(1) = one;
963
964 this->setMyDescription(myDescription.str());
965 this->setMy_A(myA);
966 this->setMy_b(myb);
967 this->setMy_c(myc);
968 this->setMy_order(2);
969 }
970};
971
972
973template<class Scalar>
974class SDIRK2Stage2ndOrder_RKBT :
975 virtual public RKButcherTableauDefaultBase<Scalar>
976{
977 public:
978 SDIRK2Stage2ndOrder_RKBT()
979 {
980 std::ostringstream myDescription;
981 myDescription << SDIRK2Stage2ndOrder_name() << "\n"
982 << "Computer Methods for ODEs and DAEs\n"
983 << "U. M. Ascher and L. R. Petzold\n"
984 << "p. 106\n"
985 << "gamma = (2+-sqrt(2))/2\n"
986 << "c = [ gamma 1 ]'\n"
987 << "A = [ gamma 0 ]\n"
988 << " [ 1-gamma gamma ]\n"
989 << "b = [ 1-gamma gamma ]'" << std::endl;
990
991 this->setMyDescription(myDescription.str());
992 typedef ScalarTraits<Scalar> ST;
993 Scalar one = ST::one();
994 gamma_default_ = as<Scalar>( (2*one - ST::squareroot(2*one))/(2*one) );
995 gamma_ = gamma_default_;
996 this->setupData();
997
998 RCP<ParameterList> validPL = Teuchos::parameterList();
999 validPL->set("Description","",this->getMyDescription());
1000 validPL->set<double>("gamma",gamma_default_,
1001 "The default value is gamma = (2-sqrt(2))/2. "
1002 "This will produce an L-stable 2nd order method with the stage "
1003 "times within the timestep. Other values of gamma will still "
1004 "produce an L-stable scheme, but will only be 1st order accurate.");
1005 Teuchos::setupVerboseObjectSublist(&*validPL);
1006 this->setMyValidParameterList(validPL);
1007 }
1008 void setupData()
1009 {
1010 typedef ScalarTraits<Scalar> ST;
1011 int myNumStages = 2;
1012 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1013 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1014 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1015 Scalar one = ST::one();
1016 Scalar zero = ST::zero();
1017 myA(0,0) = gamma_;
1018 myA(0,1) = zero;
1019 myA(1,0) = as<Scalar>( one - gamma_ );
1020 myA(1,1) = gamma_;
1021 myb(0) = as<Scalar>( one - gamma_ );
1022 myb(1) = gamma_;
1023 myc(0) = gamma_;
1024 myc(1) = one;
1025
1026 this->setMy_A(myA);
1027 this->setMy_b(myb);
1028 this->setMy_c(myc);
1029 this->setMy_order(2);
1030 }
1031 void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
1032 {
1033 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1034 paramList->validateParameters(*this->getValidParameters());
1035 Teuchos::readVerboseObjectSublist(&*paramList,this);
1036 gamma_ = paramList->get<double>("gamma",gamma_default_);
1037 this->setupData();
1038 this->setMyParamList(paramList);
1039 }
1040 private:
1041 Scalar gamma_default_;
1042 Scalar gamma_;
1043};
1044
1045
1046// 04/07/09 tscoffe: I verified manually that the Convergence Testing passes
1047// with gamma_default_ = -1.
1048template<class Scalar>
1049class SDIRK2Stage3rdOrder_RKBT :
1050 virtual public RKButcherTableauDefaultBase<Scalar>
1051{
1052 public:
1053 SDIRK2Stage3rdOrder_RKBT()
1054 {
1055 std::ostringstream myDescription;
1056 myDescription << SDIRK2Stage3rdOrder_name() << "\n"
1057 << "Solving Ordinary Differential Equations I:\n"
1058 << "Nonstiff Problems, 2nd Revised Edition\n"
1059 << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1060 << "Table 7.2, pg 207\n"
1061 << "gamma = (3+-sqrt(3))/6 -> 3rd order and A-stable\n"
1062 << "gamma = (2+-sqrt(2))/2 -> 2nd order and L-stable\n"
1063 << "c = [ gamma 1-gamma ]'\n"
1064 << "A = [ gamma 0 ]\n"
1065 << " [ 1-2*gamma gamma ]\n"
1066 << "b = [ 1/2 1/2 ]'" << std::endl;
1067
1068 this->setMyDescription(myDescription.str());
1069 thirdOrderAStable_default_ = true;
1070 secondOrderLStable_default_ = false;
1071 thirdOrderAStable_ = true;
1072 secondOrderLStable_ = false;
1073 typedef ScalarTraits<Scalar> ST;
1074 Scalar one = ST::one();
1075 gamma_default_ = as<Scalar>( (3*one + ST::squareroot(3*one))/(6*one) );
1076 gamma_ = gamma_default_;
1077 this->setupData();
1078
1079 RCP<ParameterList> validPL = Teuchos::parameterList();
1080 validPL->set("Description","",this->getMyDescription());
1081 validPL->set("3rd Order A-stable",thirdOrderAStable_default_,
1082 "If true, set gamma to gamma = (3+sqrt(3))/6 to obtain "
1083 "a 3rd order A-stable scheme. '3rd Order A-stable' and "
1084 "'2nd Order L-stable' can not both be true.");
1085 validPL->set("2nd Order L-stable",secondOrderLStable_default_,
1086 "If true, set gamma to gamma = (2+sqrt(2))/2 to obtain "
1087 "a 2nd order L-stable scheme. '3rd Order A-stable' and "
1088 "'2nd Order L-stable' can not both be true.");
1089 validPL->set("gamma",gamma_default_,
1090 "If both '3rd Order A-stable' and '2nd Order L-stable' "
1091 "are false, gamma will be used. The default value is the "
1092 "'3rd Order A-stable' gamma value, (3+sqrt(3))/6.");
1093
1094 Teuchos::setupVerboseObjectSublist(&*validPL);
1095 this->setMyValidParameterList(validPL);
1096 }
1097 void setupData()
1098 {
1099 typedef ScalarTraits<Scalar> ST;
1100 int myNumStages = 2;
1101 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1102 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1103 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1104 Scalar one = ST::one();
1105 Scalar zero = ST::zero();
1106 Scalar gamma;
1107 if (thirdOrderAStable_)
1108 gamma = as<Scalar>( (3*one + ST::squareroot(3*one))/(6*one) );
1109 else if (secondOrderLStable_)
1110 gamma = as<Scalar>( (2*one + ST::squareroot(2*one))/(2*one) );
1111 else
1112 gamma = gamma_;
1113 myA(0,0) = gamma;
1114 myA(0,1) = zero;
1115 myA(1,0) = as<Scalar>( one - 2*gamma );
1116 myA(1,1) = gamma;
1117 myb(0) = as<Scalar>( one/(2*one) );
1118 myb(1) = as<Scalar>( one/(2*one) );
1119 myc(0) = gamma;
1120 myc(1) = as<Scalar>( one - gamma );
1121
1122 this->setMy_A(myA);
1123 this->setMy_b(myb);
1124 this->setMy_c(myc);
1125 this->setMy_order(3);
1126 }
1127 void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
1128 {
1129 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1130 paramList->validateParameters(*this->getValidParameters());
1131 Teuchos::readVerboseObjectSublist(&*paramList,this);
1132 thirdOrderAStable_ = paramList->get("3rd Order A-stable",
1133 thirdOrderAStable_default_);
1134 secondOrderLStable_ = paramList->get("2nd Order L-stable",
1135 secondOrderLStable_default_);
1136 TEUCHOS_TEST_FOR_EXCEPTION(
1137 thirdOrderAStable_ && secondOrderLStable_, std::logic_error,
1138 "'3rd Order A-stable' and '2nd Order L-stable' can not both be true.");
1139 gamma_ = paramList->get("gamma",gamma_default_);
1140
1141 this->setupData();
1142 this->setMyParamList(paramList);
1143 }
1144 private:
1145 bool thirdOrderAStable_default_;
1146 bool thirdOrderAStable_;
1147 bool secondOrderLStable_default_;
1148 bool secondOrderLStable_;
1149 Scalar gamma_default_;
1150 Scalar gamma_;
1151};
1152
1153
1154template<class Scalar>
1155class DIRK2Stage3rdOrder_RKBT :
1156 virtual public RKButcherTableauDefaultBase<Scalar>
1157{
1158 public:
1159 DIRK2Stage3rdOrder_RKBT()
1160 {
1161
1162 std::ostringstream myDescription;
1163 myDescription << DIRK2Stage3rdOrder_name() << "\n"
1164 << "Hammer & Hollingsworth method\n"
1165 << "Solving Ordinary Differential Equations I:\n"
1166 << "Nonstiff Problems, 2nd Revised Edition\n"
1167 << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1168 << "Table 7.1, pg 205\n"
1169 << "c = [ 0 2/3 ]'\n"
1170 << "A = [ 0 0 ]\n"
1171 << " [ 1/3 1/3 ]\n"
1172 << "b = [ 1/4 3/4 ]'" << std::endl;
1173 typedef ScalarTraits<Scalar> ST;
1174 int myNumStages = 2;
1175 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1176 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1177 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1178 Scalar one = ST::one();
1179 Scalar zero = ST::zero();
1180 myA(0,0) = zero;
1181 myA(0,1) = zero;
1182 myA(1,0) = as<Scalar>( one/(3*one) );
1183 myA(1,1) = as<Scalar>( one/(3*one) );
1184 myb(0) = as<Scalar>( one/(4*one) );
1185 myb(1) = as<Scalar>( 3*one/(4*one) );
1186 myc(0) = zero;
1187 myc(1) = as<Scalar>( 2*one/(3*one) );
1188 this->setMyDescription(myDescription.str());
1189 this->setMy_A(myA);
1190 this->setMy_b(myb);
1191 this->setMy_c(myc);
1192 this->setMy_order(3);
1193 }
1194};
1195
1196
1197template<class Scalar>
1198class Implicit3Stage6thOrderKuntzmannButcher_RKBT :
1199 virtual public RKButcherTableauDefaultBase<Scalar>
1200{
1201 public:
1202 Implicit3Stage6thOrderKuntzmannButcher_RKBT()
1203 {
1204
1205 std::ostringstream myDescription;
1206 myDescription << Implicit3Stage6thOrderKuntzmannButcher_name() << "\n"
1207 << "Kuntzmann & Butcher method\n"
1208 << "Solving Ordinary Differential Equations I:\n"
1209 << "Nonstiff Problems, 2nd Revised Edition\n"
1210 << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1211 << "Table 7.4, pg 209\n"
1212 << "c = [ 1/2-sqrt(15)/10 1/2 1/2-sqrt(15)/10 ]'\n"
1213 << "A = [ 5/36 2/9-sqrt(15)/15 5/36-sqrt(15)/30 ]\n"
1214 << " [ 5/36+sqrt(15)/24 2/9 5/36-sqrt(15)/24 ]\n"
1215 << " [ 5/36+sqrt(15)/30 2/9+sqrt(15)/15 5/36 ]\n"
1216 << "b = [ 5/18 4/9 5/18 ]'" << std::endl;
1217 typedef ScalarTraits<Scalar> ST;
1218 int myNumStages = 3;
1219 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1220 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1221 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1222 Scalar one = ST::one();
1223 myA(0,0) = as<Scalar>( 5*one/(36*one) );
1224 myA(0,1) = as<Scalar>( 2*one/(9*one) - ST::squareroot(15*one)/(15*one) );
1225 myA(0,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(30*one) );
1226 myA(1,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(24*one) );
1227 myA(1,1) = as<Scalar>( 2*one/(9*one) );
1228 myA(1,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(24*one) );
1229 myA(2,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(30*one) );
1230 myA(2,1) = as<Scalar>( 2*one/(9*one) + ST::squareroot(15*one)/(15*one) );
1231 myA(2,2) = as<Scalar>( 5*one/(36*one) );
1232 myb(0) = as<Scalar>( 5*one/(18*one) );
1233 myb(1) = as<Scalar>( 4*one/(9*one) );
1234 myb(2) = as<Scalar>( 5*one/(18*one) );
1235 myc(0) = as<Scalar>( one/(2*one)-ST::squareroot(15*one)/(10*one) );
1236 myc(1) = as<Scalar>( one/(2*one) );
1237 myc(2) = as<Scalar>( one/(2*one)+ST::squareroot(15*one)/(10*one) );
1238 this->setMyDescription(myDescription.str());
1239 this->setMy_A(myA);
1240 this->setMy_b(myb);
1241 this->setMy_c(myc);
1242 this->setMy_order(6);
1243 }
1244};
1245
1246
1247template<class Scalar>
1248class Implicit4Stage8thOrderKuntzmannButcher_RKBT :
1249 virtual public RKButcherTableauDefaultBase<Scalar>
1250{
1251 public:
1252 Implicit4Stage8thOrderKuntzmannButcher_RKBT()
1253 {
1254
1255 std::ostringstream myDescription;
1256 myDescription << Implicit4Stage8thOrderKuntzmannButcher_name() << "\n"
1257 << "Kuntzmann & Butcher method\n"
1258 << "Solving Ordinary Differential Equations I:\n"
1259 << "Nonstiff Problems, 2nd Revised Edition\n"
1260 << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1261 << "Table 7.5, pg 209\n"
1262 << "c = [ 1/2-w2 1/2-w2p 1/2+w2p 1/2+w2 ]'\n"
1263 << "A = [ w1 w1p-w3+w4p w1p-w3-w4p w1-w5 ]\n"
1264 << " [ w1-w3p+w4 w1p w1p-w5p w1-w3p-w4 ]\n"
1265 << " [ w1+w3p+w4 w1p+w5p w1p w1+w3p-w4 ]\n"
1266 << " [ w1+w5 w1p+w3+w4p w1p+w3-w4p w1 ]\n"
1267 << "b = [ 2*w1 2*w1p 2*w1p 2*w1 ]'\n"
1268 << "w1 = 1/8-sqrt(30)/144\n"
1269 << "w2 = (1/2)*sqrt((15+2*sqrt(3))/35)\n"
1270 << "w3 = w2*(1/6+sqrt(30)/24)\n"
1271 << "w4 = w2*(1/21+5*sqrt(30)/168)\n"
1272 << "w5 = w2-2*w3\n"
1273 << "w1p = 1/8+sqrt(30)/144\n"
1274 << "w2p = (1/2)*sqrt((15-2*sqrt(3))/35)\n"
1275 << "w3p = w2*(1/6-sqrt(30)/24)\n"
1276 << "w4p = w2*(1/21-5*sqrt(30)/168)\n"
1277 << "w5p = w2p-2*w3p" << std::endl;
1278 typedef ScalarTraits<Scalar> ST;
1279 int myNumStages = 4;
1280 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1281 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1282 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1283 Scalar one = ST::one();
1284 Scalar onehalf = as<Scalar>( one/(2*one) );
1285 Scalar w1 = as<Scalar>( one/(8*one) - ST::squareroot(30*one)/(144*one) );
1286 Scalar w2 = as<Scalar>( (one/(2*one))*ST::squareroot((15*one+2*one*ST::squareroot(30*one))/(35*one)) );
1287 Scalar w3 = as<Scalar>( w2*(one/(6*one)+ST::squareroot(30*one)/(24*one)) );
1288 Scalar w4 = as<Scalar>( w2*(one/(21*one)+5*one*ST::squareroot(30*one)/(168*one)) );
1289 Scalar w5 = as<Scalar>( w2-2*w3 );
1290 Scalar w1p = as<Scalar>( one/(8*one) + ST::squareroot(30*one)/(144*one) );
1291 Scalar w2p = as<Scalar>( (one/(2*one))*ST::squareroot((15*one-2*one*ST::squareroot(30*one))/(35*one)) );
1292 Scalar w3p = as<Scalar>( w2p*(one/(6*one)-ST::squareroot(30*one)/(24*one)) );
1293 Scalar w4p = as<Scalar>( w2p*(one/(21*one)-5*one*ST::squareroot(30*one)/(168*one)) );
1294 Scalar w5p = as<Scalar>( w2p-2*w3p );
1295 myA(0,0) = w1;
1296 myA(0,1) = w1p-w3+w4p;
1297 myA(0,2) = w1p-w3-w4p;
1298 myA(0,3) = w1-w5;
1299 myA(1,0) = w1-w3p+w4;
1300 myA(1,1) = w1p;
1301 myA(1,2) = w1p-w5p;
1302 myA(1,3) = w1-w3p-w4;
1303 myA(2,0) = w1+w3p+w4;
1304 myA(2,1) = w1p+w5p;
1305 myA(2,2) = w1p;
1306 myA(2,3) = w1+w3p-w4;
1307 myA(3,0) = w1+w5;
1308 myA(3,1) = w1p+w3+w4p;
1309 myA(3,2) = w1p+w3-w4p;
1310 myA(3,3) = w1;
1311 myb(0) = 2*w1;
1312 myb(1) = 2*w1p;
1313 myb(2) = 2*w1p;
1314 myb(3) = 2*w1;
1315 myc(0) = onehalf - w2;
1316 myc(1) = onehalf - w2p;
1317 myc(2) = onehalf + w2p;
1318 myc(3) = onehalf + w2;
1319 this->setMyDescription(myDescription.str());
1320 this->setMy_A(myA);
1321 this->setMy_b(myb);
1322 this->setMy_c(myc);
1323 this->setMy_order(8);
1324 }
1325};
1326
1327
1328template<class Scalar>
1329class Implicit2Stage4thOrderHammerHollingsworth_RKBT :
1330 virtual public RKButcherTableauDefaultBase<Scalar>
1331{
1332 public:
1333 Implicit2Stage4thOrderHammerHollingsworth_RKBT()
1334 {
1335
1336 std::ostringstream myDescription;
1337 myDescription << Implicit2Stage4thOrderHammerHollingsworth_name() << "\n"
1338 << "Hammer & Hollingsworth method\n"
1339 << "Solving Ordinary Differential Equations I:\n"
1340 << "Nonstiff Problems, 2nd Revised Edition\n"
1341 << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1342 << "Table 7.3, pg 207\n"
1343 << "c = [ 1/2-sqrt(3)/6 1/2+sqrt(3)/6 ]'\n"
1344 << "A = [ 1/4 1/4-sqrt(3)/6 ]\n"
1345 << " [ 1/4+sqrt(3)/6 1/4 ]\n"
1346 << "b = [ 1/2 1/2 ]'" << std::endl;
1347 typedef ScalarTraits<Scalar> ST;
1348 int myNumStages = 2;
1349 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1350 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1351 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1352 Scalar one = ST::one();
1353 Scalar onequarter = as<Scalar>( one/(4*one) );
1354 Scalar onehalf = as<Scalar>( one/(2*one) );
1355 myA(0,0) = onequarter;
1356 myA(0,1) = as<Scalar>( onequarter-ST::squareroot(3*one)/(6*one) );
1357 myA(1,0) = as<Scalar>( onequarter+ST::squareroot(3*one)/(6*one) );
1358 myA(1,1) = onequarter;
1359 myb(0) = onehalf;
1360 myb(1) = onehalf;
1361 myc(0) = as<Scalar>( onehalf - ST::squareroot(3*one)/(6*one) );
1362 myc(1) = as<Scalar>( onehalf + ST::squareroot(3*one)/(6*one) );
1363 this->setMyDescription(myDescription.str());
1364 this->setMy_A(myA);
1365 this->setMy_b(myb);
1366 this->setMy_c(myc);
1367 this->setMy_order(4);
1368 }
1369};
1370
1371
1372template<class Scalar>
1373class IRK1StageTheta_RKBT :
1374 virtual public RKButcherTableauDefaultBase<Scalar>
1375{
1376 public:
1377 IRK1StageTheta_RKBT()
1378 {
1379
1380 std::ostringstream myDescription;
1381 myDescription << IRK1StageTheta_name() << "\n"
1382 << "Non-standard finite-difference methods\n"
1383 << "in dynamical systems, P. Kama,\n"
1384 << "Dissertation, University of Pretoria, pg. 49.\n"
1385 << "Comment: Generalized Implicit Midpoint Method\n"
1386 << "c = [ theta ]'\n"
1387 << "A = [ theta ]\n"
1388 << "b = [ 1 ]'\n"
1389 << std::endl;
1390
1391 this->setMyDescription(myDescription.str());
1392 typedef ScalarTraits<Scalar> ST;
1393 theta_default_ = ST::one()/(2*ST::one());
1394 theta_ = theta_default_;
1395 this->setupData();
1396
1397 RCP<ParameterList> validPL = Teuchos::parameterList();
1398 validPL->set("Description","",this->getMyDescription());
1399 validPL->set<double>("theta",theta_default_,
1400 "Valid values are 0 <= theta <= 1, where theta = 0 "
1401 "implies Forward Euler, theta = 1/2 implies midpoint "
1402 "method, and theta = 1 implies Backward Euler. "
1403 "For theta != 1/2, this method is first-order accurate, "
1404 "and with theta = 1/2, it is second-order accurate. "
1405 "This method is A-stable, but becomes L-stable with theta=1.");
1406 Teuchos::setupVerboseObjectSublist(&*validPL);
1407 this->setMyValidParameterList(validPL);
1408 }
1409
1410 void setupData()
1411 {
1412 typedef ScalarTraits<Scalar> ST;
1413 int myNumStages = 1;
1414 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1415 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1416 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1417 myA(0,0) = theta_;
1418 myb(0) = ST::one();
1419 myc(0) = theta_;
1420 this->setMy_A(myA);
1421 this->setMy_b(myb);
1422 this->setMy_c(myc);
1423 if (theta_ == theta_default_) this->setMy_order(2);
1424 else this->setMy_order(1);
1425 }
1426
1427 void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
1428 {
1429 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1430 paramList->validateParameters(*this->getValidParameters());
1431 Teuchos::readVerboseObjectSublist(&*paramList,this);
1432 theta_ = paramList->get<double>("theta",theta_default_);
1433 this->setupData();
1434 this->setMyParamList(paramList);
1435 }
1436 private:
1437 Scalar theta_default_;
1438 Scalar theta_;
1439};
1440
1441
1442template<class Scalar>
1443class IRK2StageTheta_RKBT :
1444 virtual public RKButcherTableauDefaultBase<Scalar>
1445{
1446 public:
1447 IRK2StageTheta_RKBT()
1448 {
1449 std::ostringstream myDescription;
1450 myDescription << IRK2StageTheta_name() << "\n"
1451 << "Computer Methods for ODEs and DAEs\n"
1452 << "U. M. Ascher and L. R. Petzold\n"
1453 << "p. 113\n"
1454 << "c = [ 0 1 ]'\n"
1455 << "A = [ 0 0 ]\n"
1456 << " [ 1-theta theta ]\n"
1457 << "b = [ 1-theta theta ]'\n"
1458 << std::endl;
1459
1460 this->setMyDescription(myDescription.str());
1461 typedef ScalarTraits<Scalar> ST;
1462 theta_default_ = ST::one()/(2*ST::one());
1463 theta_ = theta_default_;
1464 this->setupData();
1465
1466 RCP<ParameterList> validPL = Teuchos::parameterList();
1467 validPL->set("Description","",this->getMyDescription());
1468 validPL->set<double>("theta",theta_default_,
1469 "Valid values are 0 < theta <= 1, where theta = 0 "
1470 "implies Forward Euler, theta = 1/2 implies trapezoidal "
1471 "method, and theta = 1 implies Backward Euler. "
1472 "For theta != 1/2, this method is first-order accurate, "
1473 "and with theta = 1/2, it is second-order accurate. "
1474 "This method is A-stable, but becomes L-stable with theta=1.");
1475 Teuchos::setupVerboseObjectSublist(&*validPL);
1476 this->setMyValidParameterList(validPL);
1477 }
1478 void setupData()
1479 {
1480 typedef ScalarTraits<Scalar> ST;
1481 int myNumStages = 2;
1482 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1483 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1484 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1485 Scalar one = ST::one();
1486 Scalar zero = ST::zero();
1487 myA(0,0) = zero;
1488 myA(0,1) = zero;
1489 myA(1,0) = as<Scalar>( one - theta_ );
1490 myA(1,1) = theta_;
1491 myb(0) = as<Scalar>( one - theta_ );
1492 myb(1) = theta_;
1493 myc(0) = theta_;
1494 myc(1) = one;
1495
1496 this->setMy_A(myA);
1497 this->setMy_b(myb);
1498 this->setMy_c(myc);
1499 TEUCHOS_TEST_FOR_EXCEPTION(
1500 theta_ == zero, std::logic_error,
1501 "'theta' can not be zero, as it makes this IRK stepper explicit.");
1502 if (theta_ == theta_default_) this->setMy_order(2);
1503 else this->setMy_order(1);
1504 }
1505 void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
1506 {
1507 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1508 paramList->validateParameters(*this->getValidParameters());
1509 Teuchos::readVerboseObjectSublist(&*paramList,this);
1510 theta_ = paramList->get<double>("theta",theta_default_);
1511 this->setupData();
1512 this->setMyParamList(paramList);
1513 }
1514 private:
1515 Scalar theta_default_;
1516 Scalar theta_;
1517};
1518
1519
1520template<class Scalar>
1521class Implicit1Stage2ndOrderGauss_RKBT :
1522 virtual public RKButcherTableauDefaultBase<Scalar>
1523{
1524 public:
1525 Implicit1Stage2ndOrderGauss_RKBT()
1526 {
1527
1528 std::ostringstream myDescription;
1529 myDescription << Implicit1Stage2ndOrderGauss_name() << "\n"
1530 << "A-stable\n"
1531 << "Solving Ordinary Differential Equations II:\n"
1532 << "Stiff and Differential-Algebraic Problems,\n"
1533 << "2nd Revised Edition\n"
1534 << "E. Hairer and G. Wanner\n"
1535 << "Table 5.2, pg 72\n"
1536 << "Also: Implicit midpoint rule\n"
1537 << "Solving Ordinary Differential Equations I:\n"
1538 << "Nonstiff Problems, 2nd Revised Edition\n"
1539 << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1540 << "Table 7.1, pg 205\n"
1541 << "c = [ 1/2 ]'\n"
1542 << "A = [ 1/2 ]\n"
1543 << "b = [ 1 ]'" << std::endl;
1544 typedef ScalarTraits<Scalar> ST;
1545 int myNumStages = 1;
1546 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1547 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1548 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1549 Scalar onehalf = ST::one()/(2*ST::one());
1550 Scalar one = ST::one();
1551 myA(0,0) = onehalf;
1552 myb(0) = one;
1553 myc(0) = onehalf;
1554 this->setMyDescription(myDescription.str());
1555 this->setMy_A(myA);
1556 this->setMy_b(myb);
1557 this->setMy_c(myc);
1558 this->setMy_order(2);
1559 }
1560};
1561
1562
1563template<class Scalar>
1564class Implicit2Stage4thOrderGauss_RKBT :
1565 virtual public RKButcherTableauDefaultBase<Scalar>
1566{
1567 public:
1568 Implicit2Stage4thOrderGauss_RKBT()
1569 {
1570
1571 std::ostringstream myDescription;
1572 myDescription << Implicit2Stage4thOrderGauss_name() << "\n"
1573 << "A-stable\n"
1574 << "Solving Ordinary Differential Equations II:\n"
1575 << "Stiff and Differential-Algebraic Problems,\n"
1576 << "2nd Revised Edition\n"
1577 << "E. Hairer and G. Wanner\n"
1578 << "Table 5.2, pg 72\n"
1579 << "c = [ 1/2-sqrt(3)/6 1/2+sqrt(3)/6 ]'\n"
1580 << "A = [ 1/4 1/4-sqrt(3)/6 ]\n"
1581 << " [ 1/4+sqrt(3)/6 1/4 ]\n"
1582 << "b = [ 1/2 1/2 ]'" << std::endl;
1583 typedef ScalarTraits<Scalar> ST;
1584 int myNumStages = 2;
1585 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1586 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1587 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1588 Scalar one = ST::one();
1589 Scalar onehalf = as<Scalar>(one/(2*one));
1590 Scalar three = as<Scalar>(3*one);
1591 Scalar six = as<Scalar>(6*one);
1592 Scalar onefourth = as<Scalar>(one/(4*one));
1593 Scalar alpha = ST::squareroot(three)/six;
1594
1595 myA(0,0) = onefourth;
1596 myA(0,1) = onefourth-alpha;
1597 myA(1,0) = onefourth+alpha;
1598 myA(1,1) = onefourth;
1599 myb(0) = onehalf;
1600 myb(1) = onehalf;
1601 myc(0) = onehalf-alpha;
1602 myc(1) = onehalf+alpha;
1603 this->setMyDescription(myDescription.str());
1604 this->setMy_A(myA);
1605 this->setMy_b(myb);
1606 this->setMy_c(myc);
1607 this->setMy_order(4);
1608 }
1609};
1610
1611
1612template<class Scalar>
1613class Implicit3Stage6thOrderGauss_RKBT :
1614 virtual public RKButcherTableauDefaultBase<Scalar>
1615{
1616 public:
1617 Implicit3Stage6thOrderGauss_RKBT()
1618 {
1619
1620 std::ostringstream myDescription;
1621 myDescription << Implicit3Stage6thOrderGauss_name() << "\n"
1622 << "A-stable\n"
1623 << "Solving Ordinary Differential Equations II:\n"
1624 << "Stiff and Differential-Algebraic Problems,\n"
1625 << "2nd Revised Edition\n"
1626 << "E. Hairer and G. Wanner\n"
1627 << "Table 5.2, pg 72\n"
1628 << "c = [ 1/2-sqrt(15)/10 1/2 1/2+sqrt(15)/10 ]'\n"
1629 << "A = [ 5/36 2/9-sqrt(15)/15 5/36-sqrt(15)/30 ]\n"
1630 << " [ 5/36+sqrt(15)/24 2/9 5/36-sqrt(15)/24 ]\n"
1631 << " [ 5/36+sqrt(15)/30 2/9+sqrt(15)/15 5/36 ]\n"
1632 << "b = [ 5/18 4/9 5/18 ]'" << std::endl;
1633 typedef ScalarTraits<Scalar> ST;
1634 int myNumStages = 3;
1635 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1636 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1637 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1638 Scalar one = ST::one();
1639 Scalar ten = as<Scalar>(10*one);
1640 Scalar fifteen = as<Scalar>(15*one);
1641 Scalar twentyfour = as<Scalar>(24*one);
1642 Scalar thirty = as<Scalar>(30*one);
1643 Scalar sqrt15over10 = as<Scalar>(ST::squareroot(fifteen)/ten);
1644 Scalar sqrt15over15 = as<Scalar>(ST::squareroot(fifteen)/fifteen);
1645 Scalar sqrt15over24 = as<Scalar>(ST::squareroot(fifteen)/twentyfour);
1646 Scalar sqrt15over30 = as<Scalar>(ST::squareroot(fifteen)/thirty);
1647
1648 myA(0,0) = as<Scalar>(5*one/(36*one));
1649 myA(0,1) = as<Scalar>(2*one/(9*one))-sqrt15over15;
1650 myA(0,2) = as<Scalar>(5*one/(36*one))-sqrt15over30;
1651 myA(1,0) = as<Scalar>(5*one/(36*one))+sqrt15over24;
1652 myA(1,1) = as<Scalar>(2*one/(9*one));
1653 myA(1,2) = as<Scalar>(5*one/(36*one))-sqrt15over24;
1654 myA(2,0) = as<Scalar>(5*one/(36*one))+sqrt15over30;
1655 myA(2,1) = as<Scalar>(2*one/(9*one))+sqrt15over15;
1656 myA(2,2) = as<Scalar>(5*one/(36*one));
1657 myb(0) = as<Scalar>(5*one/(18*one));
1658 myb(1) = as<Scalar>(4*one/(9*one));
1659 myb(2) = as<Scalar>(5*one/(18*one));
1660 myc(0) = as<Scalar>(one/(2*one))-sqrt15over10;
1661 myc(1) = as<Scalar>(one/(2*one));
1662 myc(2) = as<Scalar>(one/(2*one))+sqrt15over10;
1663 this->setMyDescription(myDescription.str());
1664 this->setMy_A(myA);
1665 this->setMy_b(myb);
1666 this->setMy_c(myc);
1667 this->setMy_order(6);
1668 }
1669};
1670
1671
1672template<class Scalar>
1673class Implicit1Stage1stOrderRadauA_RKBT :
1674 virtual public RKButcherTableauDefaultBase<Scalar>
1675{
1676 public:
1677 Implicit1Stage1stOrderRadauA_RKBT()
1678 {
1679
1680 std::ostringstream myDescription;
1681 myDescription << Implicit1Stage1stOrderRadauA_name() << "\n"
1682 << "A-stable\n"
1683 << "Solving Ordinary Differential Equations II:\n"
1684 << "Stiff and Differential-Algebraic Problems,\n"
1685 << "2nd Revised Edition\n"
1686 << "E. Hairer and G. Wanner\n"
1687 << "Table 5.3, pg 73\n"
1688 << "c = [ 0 ]'\n"
1689 << "A = [ 1 ]\n"
1690 << "b = [ 1 ]'" << std::endl;
1691 typedef ScalarTraits<Scalar> ST;
1692 int myNumStages = 1;
1693 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1694 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1695 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1696 Scalar one = ST::one();
1697 Scalar zero = ST::zero();
1698 myA(0,0) = one;
1699 myb(0) = one;
1700 myc(0) = zero;
1701 this->setMyDescription(myDescription.str());
1702 this->setMy_A(myA);
1703 this->setMy_b(myb);
1704 this->setMy_c(myc);
1705 this->setMy_order(1);
1706 }
1707};
1708
1709
1710template<class Scalar>
1711class Implicit2Stage3rdOrderRadauA_RKBT :
1712 virtual public RKButcherTableauDefaultBase<Scalar>
1713{
1714 public:
1715 Implicit2Stage3rdOrderRadauA_RKBT()
1716 {
1717
1718 std::ostringstream myDescription;
1719 myDescription << Implicit2Stage3rdOrderRadauA_name() << "\n"
1720 << "A-stable\n"
1721 << "Solving Ordinary Differential Equations II:\n"
1722 << "Stiff and Differential-Algebraic Problems,\n"
1723 << "2nd Revised Edition\n"
1724 << "E. Hairer and G. Wanner\n"
1725 << "Table 5.3, pg 73\n"
1726 << "c = [ 0 2/3 ]'\n"
1727 << "A = [ 1/4 -1/4 ]\n"
1728 << " [ 1/4 5/12 ]\n"
1729 << "b = [ 1/4 3/4 ]'" << std::endl;
1730 typedef ScalarTraits<Scalar> ST;
1731 int myNumStages = 2;
1732 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1733 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1734 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1735 Scalar zero = ST::zero();
1736 Scalar one = ST::one();
1737 myA(0,0) = as<Scalar>(one/(4*one));
1738 myA(0,1) = as<Scalar>(-one/(4*one));
1739 myA(1,0) = as<Scalar>(one/(4*one));
1740 myA(1,1) = as<Scalar>(5*one/(12*one));
1741 myb(0) = as<Scalar>(one/(4*one));
1742 myb(1) = as<Scalar>(3*one/(4*one));
1743 myc(0) = zero;
1744 myc(1) = as<Scalar>(2*one/(3*one));
1745 this->setMyDescription(myDescription.str());
1746 this->setMy_A(myA);
1747 this->setMy_b(myb);
1748 this->setMy_c(myc);
1749 this->setMy_order(3);
1750 }
1751};
1752
1753
1754template<class Scalar>
1755class Implicit3Stage5thOrderRadauA_RKBT :
1756 virtual public RKButcherTableauDefaultBase<Scalar>
1757{
1758 public:
1759 Implicit3Stage5thOrderRadauA_RKBT()
1760 {
1761
1762 std::ostringstream myDescription;
1763 myDescription << Implicit3Stage5thOrderRadauA_name() << "\n"
1764 << "A-stable\n"
1765 << "Solving Ordinary Differential Equations II:\n"
1766 << "Stiff and Differential-Algebraic Problems,\n"
1767 << "2nd Revised Edition\n"
1768 << "E. Hairer and G. Wanner\n"
1769 << "Table 5.4, pg 73\n"
1770 << "c = [ 0 (6-sqrt(6))/10 (6+sqrt(6))/10 ]'\n"
1771 << "A = [ 1/9 (-1-sqrt(6))/18 (-1+sqrt(6))/18 ]\n"
1772 << " [ 1/9 (88+7*sqrt(6))/360 (88-43*sqrt(6))/360 ]\n"
1773 << " [ 1/9 (88+43*sqrt(6))/360 (88-7*sqrt(6))/360 ]\n"
1774 << "b = [ 1/9 (16+sqrt(6))/36 (16-sqrt(6))/36 ]'" << std::endl;
1775 typedef ScalarTraits<Scalar> ST;
1776 int myNumStages = 3;
1777 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1778 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1779 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1780 Scalar zero = ST::zero();
1781 Scalar one = ST::one();
1782 myA(0,0) = as<Scalar>(one/(9*one));
1783 myA(0,1) = as<Scalar>( (-one-ST::squareroot(6*one))/(18*one) );
1784 myA(0,2) = as<Scalar>( (-one+ST::squareroot(6*one))/(18*one) );
1785 myA(1,0) = as<Scalar>(one/(9*one));
1786 myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
1787 myA(1,2) = as<Scalar>( (88*one-43*one*ST::squareroot(6*one))/(360*one) );
1788 myA(2,0) = as<Scalar>(one/(9*one));
1789 myA(2,1) = as<Scalar>( (88*one+43*one*ST::squareroot(6*one))/(360*one) );
1790 myA(2,2) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
1791 myb(0) = as<Scalar>(one/(9*one));
1792 myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1793 myb(2) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1794 myc(0) = zero;
1795 myc(1) = as<Scalar>( (6*one-ST::squareroot(6*one))/(10*one) );
1796 myc(2) = as<Scalar>( (6*one+ST::squareroot(6*one))/(10*one) );
1797 this->setMyDescription(myDescription.str());
1798 this->setMy_A(myA);
1799 this->setMy_b(myb);
1800 this->setMy_c(myc);
1801 this->setMy_order(5);
1802 }
1803};
1804
1805
1806template<class Scalar>
1807class Implicit1Stage1stOrderRadauB_RKBT :
1808 virtual public RKButcherTableauDefaultBase<Scalar>
1809{
1810 public:
1811 Implicit1Stage1stOrderRadauB_RKBT()
1812 {
1813
1814 std::ostringstream myDescription;
1815 myDescription << Implicit1Stage1stOrderRadauB_name() << "\n"
1816 << "A-stable\n"
1817 << "Solving Ordinary Differential Equations II:\n"
1818 << "Stiff and Differential-Algebraic Problems,\n"
1819 << "2nd Revised Edition\n"
1820 << "E. Hairer and G. Wanner\n"
1821 << "Table 5.5, pg 74\n"
1822 << "c = [ 1 ]'\n"
1823 << "A = [ 1 ]\n"
1824 << "b = [ 1 ]'" << std::endl;
1825 typedef ScalarTraits<Scalar> ST;
1826 int myNumStages = 1;
1827 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1828 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1829 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1830 Scalar one = ST::one();
1831 myA(0,0) = one;
1832 myb(0) = one;
1833 myc(0) = one;
1834 this->setMyDescription(myDescription.str());
1835 this->setMy_A(myA);
1836 this->setMy_b(myb);
1837 this->setMy_c(myc);
1838 this->setMy_order(1);
1839 }
1840};
1841
1842
1843template<class Scalar>
1844class Implicit2Stage3rdOrderRadauB_RKBT :
1845 virtual public RKButcherTableauDefaultBase<Scalar>
1846{
1847 public:
1848 Implicit2Stage3rdOrderRadauB_RKBT()
1849 {
1850
1851 std::ostringstream myDescription;
1852 myDescription << Implicit2Stage3rdOrderRadauB_name() << "\n"
1853 << "A-stable\n"
1854 << "Solving Ordinary Differential Equations II:\n"
1855 << "Stiff and Differential-Algebraic Problems,\n"
1856 << "2nd Revised Edition\n"
1857 << "E. Hairer and G. Wanner\n"
1858 << "Table 5.5, pg 74\n"
1859 << "c = [ 1/3 1 ]'\n"
1860 << "A = [ 5/12 -1/12 ]\n"
1861 << " [ 3/4 1/4 ]\n"
1862 << "b = [ 3/4 1/4 ]'" << std::endl;
1863 typedef ScalarTraits<Scalar> ST;
1864 int myNumStages = 2;
1865 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1866 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1867 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1868 Scalar one = ST::one();
1869 myA(0,0) = as<Scalar>( 5*one/(12*one) );
1870 myA(0,1) = as<Scalar>( -one/(12*one) );
1871 myA(1,0) = as<Scalar>( 3*one/(4*one) );
1872 myA(1,1) = as<Scalar>( one/(4*one) );
1873 myb(0) = as<Scalar>( 3*one/(4*one) );
1874 myb(1) = as<Scalar>( one/(4*one) );
1875 myc(0) = as<Scalar>( one/(3*one) );
1876 myc(1) = one;
1877 this->setMyDescription(myDescription.str());
1878 this->setMy_A(myA);
1879 this->setMy_b(myb);
1880 this->setMy_c(myc);
1881 this->setMy_order(3);
1882 }
1883};
1884
1885
1886template<class Scalar>
1887class Implicit3Stage5thOrderRadauB_RKBT :
1888 virtual public RKButcherTableauDefaultBase<Scalar>
1889{
1890 public:
1891 Implicit3Stage5thOrderRadauB_RKBT()
1892 {
1893
1894 std::ostringstream myDescription;
1895 myDescription << Implicit3Stage5thOrderRadauB_name() << "\n"
1896 << "A-stable\n"
1897 << "Solving Ordinary Differential Equations II:\n"
1898 << "Stiff and Differential-Algebraic Problems,\n"
1899 << "2nd Revised Edition\n"
1900 << "E. Hairer and G. Wanner\n"
1901 << "Table 5.6, pg 74\n"
1902 << "c = [ (4-sqrt(6))/10 (4+sqrt(6))/10 1 ]'\n"
1903 << "A = [ A1 A2 A3 ]\n"
1904 << " A1 = [ (88-7*sqrt(6))/360 ]\n"
1905 << " [ (296+169*sqrt(6))/1800 ]\n"
1906 << " [ (16-sqrt(6))/36 ]\n"
1907 << " A2 = [ (296-169*sqrt(6))/1800 ]\n"
1908 << " [ (88+7*sqrt(6))/360 ]\n"
1909 << " [ (16+sqrt(6))/36 ]\n"
1910 << " A3 = [ (-2+3*sqrt(6))/225 ]\n"
1911 << " [ (-2-3*sqrt(6))/225 ]\n"
1912 << " [ 1/9 ]\n"
1913 << "b = [ (16-sqrt(6))/36 (16+sqrt(6))/36 1/9 ]'"
1914 << std::endl;
1915 typedef ScalarTraits<Scalar> ST;
1916 int myNumStages = 3;
1917 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1918 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1919 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1920 Scalar one = ST::one();
1921 myA(0,0) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
1922 myA(0,1) = as<Scalar>( (296*one-169*one*ST::squareroot(6*one))/(1800*one) );
1923 myA(0,2) = as<Scalar>( (-2*one+3*one*ST::squareroot(6*one))/(225*one) );
1924 myA(1,0) = as<Scalar>( (296*one+169*one*ST::squareroot(6*one))/(1800*one) );
1925 myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
1926 myA(1,2) = as<Scalar>( (-2*one-3*one*ST::squareroot(6*one))/(225*one) );
1927 myA(2,0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1928 myA(2,1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1929 myA(2,2) = as<Scalar>( one/(9*one) );
1930 myb(0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1931 myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1932 myb(2) = as<Scalar>( one/(9*one) );
1933 myc(0) = as<Scalar>( (4*one-ST::squareroot(6*one))/(10*one) );
1934 myc(1) = as<Scalar>( (4*one+ST::squareroot(6*one))/(10*one) );
1935 myc(2) = one;
1936 this->setMyDescription(myDescription.str());
1937 this->setMy_A(myA);
1938 this->setMy_b(myb);
1939 this->setMy_c(myc);
1940 this->setMy_order(5);
1941 }
1942};
1943
1944
1945template<class Scalar>
1946class Implicit2Stage2ndOrderLobattoA_RKBT :
1947 virtual public RKButcherTableauDefaultBase<Scalar>
1948{
1949 public:
1950 Implicit2Stage2ndOrderLobattoA_RKBT()
1951 {
1952
1953 std::ostringstream myDescription;
1954 myDescription << Implicit2Stage2ndOrderLobattoA_name() << "\n"
1955 << "A-stable\n"
1956 << "Solving Ordinary Differential Equations II:\n"
1957 << "Stiff and Differential-Algebraic Problems,\n"
1958 << "2nd Revised Edition\n"
1959 << "E. Hairer and G. Wanner\n"
1960 << "Table 5.7, pg 75\n"
1961 << "c = [ 0 1 ]'\n"
1962 << "A = [ 0 0 ]\n"
1963 << " [ 1/2 1/2 ]\n"
1964 << "b = [ 1/2 1/2 ]'" << std::endl;
1965 typedef ScalarTraits<Scalar> ST;
1966 int myNumStages = 2;
1967 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1968 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1969 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1970 Scalar zero = ST::zero();
1971 Scalar one = ST::one();
1972 myA(0,0) = zero;
1973 myA(0,1) = zero;
1974 myA(1,0) = as<Scalar>( one/(2*one) );
1975 myA(1,1) = as<Scalar>( one/(2*one) );
1976 myb(0) = as<Scalar>( one/(2*one) );
1977 myb(1) = as<Scalar>( one/(2*one) );
1978 myc(0) = zero;
1979 myc(1) = one;
1980 this->setMyDescription(myDescription.str());
1981 this->setMy_A(myA);
1982 this->setMy_b(myb);
1983 this->setMy_c(myc);
1984 this->setMy_order(2);
1985 }
1986};
1987
1988
1989template<class Scalar>
1990class Implicit3Stage4thOrderLobattoA_RKBT :
1991 virtual public RKButcherTableauDefaultBase<Scalar>
1992{
1993 public:
1994 Implicit3Stage4thOrderLobattoA_RKBT()
1995 {
1996
1997 std::ostringstream myDescription;
1998 myDescription << Implicit3Stage4thOrderLobattoA_name() << "\n"
1999 << "A-stable\n"
2000 << "Solving Ordinary Differential Equations II:\n"
2001 << "Stiff and Differential-Algebraic Problems,\n"
2002 << "2nd Revised Edition\n"
2003 << "E. Hairer and G. Wanner\n"
2004 << "Table 5.7, pg 75\n"
2005 << "c = [ 0 1/2 1 ]'\n"
2006 << "A = [ 0 0 0 ]\n"
2007 << " [ 5/24 1/3 -1/24 ]\n"
2008 << " [ 1/6 2/3 1/6 ]\n"
2009 << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
2010 typedef ScalarTraits<Scalar> ST;
2011 int myNumStages = 3;
2012 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2013 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2014 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2015 Scalar zero = ST::zero();
2016 Scalar one = ST::one();
2017 myA(0,0) = zero;
2018 myA(0,1) = zero;
2019 myA(0,2) = zero;
2020 myA(1,0) = as<Scalar>( (5*one)/(24*one) );
2021 myA(1,1) = as<Scalar>( (one)/(3*one) );
2022 myA(1,2) = as<Scalar>( (-one)/(24*one) );
2023 myA(2,0) = as<Scalar>( (one)/(6*one) );
2024 myA(2,1) = as<Scalar>( (2*one)/(3*one) );
2025 myA(2,2) = as<Scalar>( (1*one)/(6*one) );
2026 myb(0) = as<Scalar>( (one)/(6*one) );
2027 myb(1) = as<Scalar>( (2*one)/(3*one) );
2028 myb(2) = as<Scalar>( (1*one)/(6*one) );
2029 myc(0) = zero;
2030 myc(1) = as<Scalar>( one/(2*one) );
2031 myc(2) = one;
2032 this->setMyDescription(myDescription.str());
2033 this->setMy_A(myA);
2034 this->setMy_b(myb);
2035 this->setMy_c(myc);
2036 this->setMy_order(4);
2037 }
2038};
2039
2040
2041template<class Scalar>
2042class Implicit4Stage6thOrderLobattoA_RKBT :
2043 virtual public RKButcherTableauDefaultBase<Scalar>
2044{
2045 public:
2046 Implicit4Stage6thOrderLobattoA_RKBT()
2047 {
2048
2049 using Teuchos::as;
2050 typedef Teuchos::ScalarTraits<Scalar> ST;
2051
2052 std::ostringstream myDescription;
2053 myDescription << Implicit4Stage6thOrderLobattoA_name() << "\n"
2054 << "A-stable\n"
2055 << "Solving Ordinary Differential Equations II:\n"
2056 << "Stiff and Differential-Algebraic Problems,\n"
2057 << "2nd Revised Edition\n"
2058 << "E. Hairer and G. Wanner\n"
2059 << "Table 5.8, pg 75\n"
2060 << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
2061 << "A = [ A1 A2 A3 A4 ]\n"
2062 << " A1 = [ 0 ]\n"
2063 << " [ (11+sqrt(5)/120 ]\n"
2064 << " [ (11-sqrt(5)/120 ]\n"
2065 << " [ 1/12 ]\n"
2066 << " A2 = [ 0 ]\n"
2067 << " [ (25-sqrt(5))/120 ]\n"
2068 << " [ (25+13*sqrt(5))/120 ]\n"
2069 << " [ 5/12 ]\n"
2070 << " A3 = [ 0 ]\n"
2071 << " [ (25-13*sqrt(5))/120 ]\n"
2072 << " [ (25+sqrt(5))/120 ]\n"
2073 << " [ 5/12 ]\n"
2074 << " A4 = [ 0 ]\n"
2075 << " [ (-1+sqrt(5))/120 ]\n"
2076 << " [ (-1-sqrt(5))/120 ]\n"
2077 << " [ 1/12 ]\n"
2078 << "b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
2079 typedef ScalarTraits<Scalar> ST;
2080 int myNumStages = 4;
2081 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2082 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2083 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2084 Scalar zero = ST::zero();
2085 Scalar one = ST::one();
2086 myA(0,0) = zero;
2087 myA(0,1) = zero;
2088 myA(0,2) = zero;
2089 myA(0,3) = zero;
2090 myA(1,0) = as<Scalar>( (11*one+ST::squareroot(5*one))/(120*one) );
2091 myA(1,1) = as<Scalar>( (25*one-ST::squareroot(5*one))/(120*one) );
2092 myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5*one))/(120*one) );
2093 myA(1,3) = as<Scalar>( (-one+ST::squareroot(5*one))/(120*one) );
2094 myA(2,0) = as<Scalar>( (11*one-ST::squareroot(5*one))/(120*one) );
2095 myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5*one))/(120*one) );
2096 myA(2,2) = as<Scalar>( (25*one+ST::squareroot(5*one))/(120*one) );
2097 myA(2,3) = as<Scalar>( (-one-ST::squareroot(5*one))/(120*one) );
2098 myA(3,0) = as<Scalar>( one/(12*one) );
2099 myA(3,1) = as<Scalar>( 5*one/(12*one) );
2100 myA(3,2) = as<Scalar>( 5*one/(12*one) );
2101 myA(3,3) = as<Scalar>( one/(12*one) );
2102 myb(0) = as<Scalar>( one/(12*one) );
2103 myb(1) = as<Scalar>( 5*one/(12*one) );
2104 myb(2) = as<Scalar>( 5*one/(12*one) );
2105 myb(3) = as<Scalar>( one/(12*one) );
2106 myc(0) = zero;
2107 myc(1) = as<Scalar>( (5*one-ST::squareroot(5))/(10*one) );
2108 myc(2) = as<Scalar>( (5*one+ST::squareroot(5))/(10*one) );
2109 myc(3) = one;
2110 this->setMyDescription(myDescription.str());
2111 this->setMy_A(myA);
2112 this->setMy_b(myb);
2113 this->setMy_c(myc);
2114 this->setMy_order(6);
2115 }
2116};
2117
2118
2119template<class Scalar>
2120class Implicit2Stage2ndOrderLobattoB_RKBT :
2121 virtual public RKButcherTableauDefaultBase<Scalar>
2122{
2123 public:
2124 Implicit2Stage2ndOrderLobattoB_RKBT()
2125 {
2126
2127 std::ostringstream myDescription;
2128 myDescription << Implicit2Stage2ndOrderLobattoB_name() << "\n"
2129 << "A-stable\n"
2130 << "Solving Ordinary Differential Equations II:\n"
2131 << "Stiff and Differential-Algebraic Problems,\n"
2132 << "2nd Revised Edition\n"
2133 << "E. Hairer and G. Wanner\n"
2134 << "Table 5.9, pg 76\n"
2135 << "c = [ 0 1 ]'\n"
2136 << "A = [ 1/2 0 ]\n"
2137 << " [ 1/2 0 ]\n"
2138 << "b = [ 1/2 1/2 ]'" << std::endl;
2139 typedef ScalarTraits<Scalar> ST;
2140 int myNumStages = 2;
2141 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2142 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2143 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2144 Scalar zero = ST::zero();
2145 Scalar one = ST::one();
2146 myA(0,0) = as<Scalar>( one/(2*one) );
2147 myA(0,1) = zero;
2148 myA(1,0) = as<Scalar>( one/(2*one) );
2149 myA(1,1) = zero;
2150 myb(0) = as<Scalar>( one/(2*one) );
2151 myb(1) = as<Scalar>( one/(2*one) );
2152 myc(0) = zero;
2153 myc(1) = one;
2154 this->setMyDescription(myDescription.str());
2155 this->setMy_A(myA);
2156 this->setMy_b(myb);
2157 this->setMy_c(myc);
2158 this->setMy_order(2);
2159 }
2160};
2161
2162
2163template<class Scalar>
2164class Implicit3Stage4thOrderLobattoB_RKBT :
2165 virtual public RKButcherTableauDefaultBase<Scalar>
2166{
2167 public:
2168 Implicit3Stage4thOrderLobattoB_RKBT()
2169 {
2170
2171 std::ostringstream myDescription;
2172 myDescription << Implicit3Stage4thOrderLobattoB_name() << "\n"
2173 << "A-stable\n"
2174 << "Solving Ordinary Differential Equations II:\n"
2175 << "Stiff and Differential-Algebraic Problems,\n"
2176 << "2nd Revised Edition\n"
2177 << "E. Hairer and G. Wanner\n"
2178 << "Table 5.9, pg 76\n"
2179 << "c = [ 0 1/2 1 ]'\n"
2180 << "A = [ 1/6 -1/6 0 ]\n"
2181 << " [ 1/6 1/3 0 ]\n"
2182 << " [ 1/6 5/6 0 ]\n"
2183 << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
2184 typedef ScalarTraits<Scalar> ST;
2185 int myNumStages = 3;
2186 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2187 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2188 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2189 Scalar zero = ST::zero();
2190 Scalar one = ST::one();
2191 myA(0,0) = as<Scalar>( one/(6*one) );
2192 myA(0,1) = as<Scalar>( -one/(6*one) );
2193 myA(0,2) = zero;
2194 myA(1,0) = as<Scalar>( one/(6*one) );
2195 myA(1,1) = as<Scalar>( one/(3*one) );
2196 myA(1,2) = zero;
2197 myA(2,0) = as<Scalar>( one/(6*one) );
2198 myA(2,1) = as<Scalar>( 5*one/(6*one) );
2199 myA(2,2) = zero;
2200 myb(0) = as<Scalar>( one/(6*one) );
2201 myb(1) = as<Scalar>( 2*one/(3*one) );
2202 myb(2) = as<Scalar>( one/(6*one) );
2203 myc(0) = zero;
2204 myc(1) = as<Scalar>( one/(2*one) );
2205 myc(2) = one;
2206 this->setMyDescription(myDescription.str());
2207 this->setMy_A(myA);
2208 this->setMy_b(myb);
2209 this->setMy_c(myc);
2210 this->setMy_order(4);
2211 }
2212};
2213
2214
2215template<class Scalar>
2216class Implicit4Stage6thOrderLobattoB_RKBT :
2217 virtual public RKButcherTableauDefaultBase<Scalar>
2218{
2219 public:
2220 Implicit4Stage6thOrderLobattoB_RKBT()
2221 {
2222
2223 std::ostringstream myDescription;
2224 myDescription << Implicit4Stage6thOrderLobattoB_name() << "\n"
2225 << "A-stable\n"
2226 << "Solving Ordinary Differential Equations II:\n"
2227 << "Stiff and Differential-Algebraic Problems,\n"
2228 << "2nd Revised Edition\n"
2229 << "E. Hairer and G. Wanner\n"
2230 << "Table 5.10, pg 76\n"
2231 << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
2232 << "A = [ 1/12 (-1-sqrt(5))/24 (-1+sqrt(5))/24 0 ]\n"
2233 << " [ 1/12 (25+sqrt(5))/120 (25-13*sqrt(5))/120 0 ]\n"
2234 << " [ 1/12 (25+13*sqrt(5))/120 (25-sqrt(5))/120 0 ]\n"
2235 << " [ 1/12 (11-sqrt(5))/24 (11+sqrt(5))/24 0 ]\n"
2236 << "b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
2237 typedef ScalarTraits<Scalar> ST;
2238 int myNumStages = 4;
2239 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2240 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2241 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2242 Scalar zero = ST::zero();
2243 Scalar one = ST::one();
2244 myA(0,0) = as<Scalar>( one/(12*one) );
2245 myA(0,1) = as<Scalar>( (-one-ST::squareroot(5))/(24*one) );
2246 myA(0,2) = as<Scalar>( (-one+ST::squareroot(5))/(24*one) );
2247 myA(0,3) = zero;
2248 myA(1,0) = as<Scalar>( one/(12*one) );
2249 myA(1,1) = as<Scalar>( (25*one+ST::squareroot(5))/(120*one) );
2250 myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5))/(120*one) );
2251 myA(1,3) = zero;
2252 myA(2,0) = as<Scalar>( one/(12*one) );
2253 myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5))/(120*one) );
2254 myA(2,2) = as<Scalar>( (25*one-ST::squareroot(5))/(120*one) );
2255 myA(2,3) = zero;
2256 myA(3,0) = as<Scalar>( one/(12*one) );
2257 myA(3,1) = as<Scalar>( (11*one-ST::squareroot(5*one))/(24*one) );
2258 myA(3,2) = as<Scalar>( (11*one+ST::squareroot(5*one))/(24*one) );
2259 myA(3,3) = zero;
2260 myb(0) = as<Scalar>( one/(12*one) );
2261 myb(1) = as<Scalar>( 5*one/(12*one) );
2262 myb(2) = as<Scalar>( 5*one/(12*one) );
2263 myb(3) = as<Scalar>( one/(12*one) );
2264 myc(0) = zero;
2265 myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
2266 myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
2267 myc(3) = one;
2268 this->setMyDescription(myDescription.str());
2269 this->setMy_A(myA);
2270 this->setMy_b(myb);
2271 this->setMy_c(myc);
2272 this->setMy_order(6);
2273 }
2274};
2275
2276
2277template<class Scalar>
2278class Implicit2Stage2ndOrderLobattoC_RKBT :
2279 virtual public RKButcherTableauDefaultBase<Scalar>
2280{
2281 public:
2282 Implicit2Stage2ndOrderLobattoC_RKBT()
2283 {
2284
2285 std::ostringstream myDescription;
2286 myDescription << Implicit2Stage2ndOrderLobattoC_name() << "\n"
2287 << "A-stable\n"
2288 << "Solving Ordinary Differential Equations II:\n"
2289 << "Stiff and Differential-Algebraic Problems,\n"
2290 << "2nd Revised Edition\n"
2291 << "E. Hairer and G. Wanner\n"
2292 << "Table 5.11, pg 76\n"
2293 << "c = [ 0 1 ]'\n"
2294 << "A = [ 1/2 -1/2 ]\n"
2295 << " [ 1/2 1/2 ]\n"
2296 << "b = [ 1/2 1/2 ]'" << std::endl;
2297 typedef ScalarTraits<Scalar> ST;
2298 int myNumStages = 2;
2299 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2300 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2301 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2302 Scalar zero = ST::zero();
2303 Scalar one = ST::one();
2304 myA(0,0) = as<Scalar>( one/(2*one) );
2305 myA(0,1) = as<Scalar>( -one/(2*one) );
2306 myA(1,0) = as<Scalar>( one/(2*one) );
2307 myA(1,1) = as<Scalar>( one/(2*one) );
2308 myb(0) = as<Scalar>( one/(2*one) );
2309 myb(1) = as<Scalar>( one/(2*one) );
2310 myc(0) = zero;
2311 myc(1) = one;
2312 this->setMyDescription(myDescription.str());
2313 this->setMy_A(myA);
2314 this->setMy_b(myb);
2315 this->setMy_c(myc);
2316 this->setMy_order(2);
2317 }
2318};
2319
2320
2321template<class Scalar>
2322class Implicit3Stage4thOrderLobattoC_RKBT :
2323 virtual public RKButcherTableauDefaultBase<Scalar>
2324{
2325 public:
2326 Implicit3Stage4thOrderLobattoC_RKBT()
2327 {
2328
2329 std::ostringstream myDescription;
2330 myDescription << Implicit3Stage4thOrderLobattoC_name() << "\n"
2331 << "A-stable\n"
2332 << "Solving Ordinary Differential Equations II:\n"
2333 << "Stiff and Differential-Algebraic Problems,\n"
2334 << "2nd Revised Edition\n"
2335 << "E. Hairer and G. Wanner\n"
2336 << "Table 5.11, pg 76\n"
2337 << "c = [ 0 1/2 1 ]'\n"
2338 << "A = [ 1/6 -1/3 1/6 ]\n"
2339 << " [ 1/6 5/12 -1/12 ]\n"
2340 << " [ 1/6 2/3 1/6 ]\n"
2341 << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
2342 typedef ScalarTraits<Scalar> ST;
2343 int myNumStages = 3;
2344 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2345 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2346 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2347 Scalar zero = ST::zero();
2348 Scalar one = ST::one();
2349 myA(0,0) = as<Scalar>( one/(6*one) );
2350 myA(0,1) = as<Scalar>( -one/(3*one) );
2351 myA(0,2) = as<Scalar>( one/(6*one) );
2352 myA(1,0) = as<Scalar>( one/(6*one) );
2353 myA(1,1) = as<Scalar>( 5*one/(12*one) );
2354 myA(1,2) = as<Scalar>( -one/(12*one) );
2355 myA(2,0) = as<Scalar>( one/(6*one) );
2356 myA(2,1) = as<Scalar>( 2*one/(3*one) );
2357 myA(2,2) = as<Scalar>( one/(6*one) );
2358 myb(0) = as<Scalar>( one/(6*one) );
2359 myb(1) = as<Scalar>( 2*one/(3*one) );
2360 myb(2) = as<Scalar>( one/(6*one) );
2361 myc(0) = zero;
2362 myc(1) = as<Scalar>( one/(2*one) );
2363 myc(2) = one;
2364 this->setMyDescription(myDescription.str());
2365 this->setMy_A(myA);
2366 this->setMy_b(myb);
2367 this->setMy_c(myc);
2368 this->setMy_order(4);
2369 }
2370};
2371
2372
2373template<class Scalar>
2374class Implicit4Stage6thOrderLobattoC_RKBT :
2375 virtual public RKButcherTableauDefaultBase<Scalar>
2376{
2377 public:
2378 Implicit4Stage6thOrderLobattoC_RKBT()
2379 {
2380
2381 std::ostringstream myDescription;
2382 myDescription << Implicit4Stage6thOrderLobattoC_name() << "\n"
2383 << "A-stable\n"
2384 << "Solving Ordinary Differential Equations II:\n"
2385 << "Stiff and Differential-Algebraic Problems,\n"
2386 << "2nd Revised Edition\n"
2387 << "E. Hairer and G. Wanner\n"
2388 << "Table 5.12, pg 76\n"
2389 << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
2390 << "A = [ 1/12 -sqrt(5)/12 sqrt(5)/12 -1/12 ]\n"
2391 << " [ 1/12 1/4 (10-7*sqrt(5))/60 sqrt(5)/60 ]\n"
2392 << " [ 1/12 (10+7*sqrt(5))/60 1/4 -sqrt(5)/60 ]\n"
2393 << " [ 1/12 5/12 5/12 1/12 ]\n"
2394 << "b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
2395 typedef ScalarTraits<Scalar> ST;
2396 int myNumStages = 4;
2397 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2398 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2399 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2400 Scalar zero = ST::zero();
2401 Scalar one = ST::one();
2402 myA(0,0) = as<Scalar>( one/(12*one) );
2403 myA(0,1) = as<Scalar>( -ST::squareroot(5*one)/(12*one) );
2404 myA(0,2) = as<Scalar>( ST::squareroot(5*one)/(12*one) );
2405 myA(0,3) = as<Scalar>( -one/(12*one) );
2406 myA(1,0) = as<Scalar>( one/(12*one) );
2407 myA(1,1) = as<Scalar>( one/(4*one) );
2408 myA(1,2) = as<Scalar>( (10*one-7*one*ST::squareroot(5*one))/(60*one) );
2409 myA(1,3) = as<Scalar>( ST::squareroot(5*one)/(60*one) );
2410 myA(2,0) = as<Scalar>( one/(12*one) );
2411 myA(2,1) = as<Scalar>( (10*one+7*one*ST::squareroot(5*one))/(60*one) );
2412 myA(2,2) = as<Scalar>( one/(4*one) );
2413 myA(2,3) = as<Scalar>( -ST::squareroot(5*one)/(60*one) );
2414 myA(3,0) = as<Scalar>( one/(12*one) );
2415 myA(3,1) = as<Scalar>( 5*one/(12*one) );
2416 myA(3,2) = as<Scalar>( 5*one/(12*one) );
2417 myA(3,3) = as<Scalar>( one/(12*one) );
2418 myb(0) = as<Scalar>( one/(12*one) );
2419 myb(1) = as<Scalar>( 5*one/(12*one) );
2420 myb(2) = as<Scalar>( 5*one/(12*one) );
2421 myb(3) = as<Scalar>( one/(12*one) );
2422 myc(0) = zero;
2423 myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
2424 myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
2425 myc(3) = one;
2426 this->setMyDescription(myDescription.str());
2427 this->setMy_A(myA);
2428 this->setMy_b(myb);
2429 this->setMy_c(myc);
2430 this->setMy_order(6);
2431 }
2432};
2433
2434
2435
2436template<class Scalar>
2437class SDIRK5Stage5thOrder_RKBT :
2438 virtual public RKButcherTableauDefaultBase<Scalar>
2439{
2440 public:
2441 SDIRK5Stage5thOrder_RKBT()
2442 {
2443
2444 std::ostringstream myDescription;
2445 myDescription << SDIRK5Stage5thOrder_name() << "\n"
2446 << "A-stable\n"
2447 << "Solving Ordinary Differential Equations II:\n"
2448 << "Stiff and Differential-Algebraic Problems,\n"
2449 << "2nd Revised Edition\n"
2450 << "E. Hairer and G. Wanner\n"
2451 << "pg101 \n"
2452 << "c = [ (6-sqrt(6))/10 ]\n"
2453 << " [ (6+9*sqrt(6))/35 ]\n"
2454 << " [ 1 ]\n"
2455 << " [ (4-sqrt(6))/10 ]\n"
2456 << " [ (4+sqrt(6))/10 ]\n"
2457 << "A = [ A1 A2 A3 A4 A5 ]\n"
2458 << " A1 = [ (6-sqrt(6))/10 ]\n"
2459 << " [ (-6+5*sqrt(6))/14 ]\n"
2460 << " [ (888+607*sqrt(6))/2850 ]\n"
2461 << " [ (3153-3082*sqrt(6))/14250 ]\n"
2462 << " [ (-32583+14638*sqrt(6))/71250 ]\n"
2463 << " A2 = [ 0 ]\n"
2464 << " [ (6-sqrt(6))/10 ]\n"
2465 << " [ (126-161*sqrt(6))/1425 ]\n"
2466 << " [ (3213+1148*sqrt(6))/28500 ]\n"
2467 << " [ (-17199+364*sqrt(6))/142500 ]\n"
2468 << " A3 = [ 0 ]\n"
2469 << " [ 0 ]\n"
2470 << " [ (6-sqrt(6))/10 ]\n"
2471 << " [ (-267+88*sqrt(6))/500 ]\n"
2472 << " [ (1329-544*sqrt(6))/2500 ]\n"
2473 << " A4 = [ 0 ]\n"
2474 << " [ 0 ]\n"
2475 << " [ 0 ]\n"
2476 << " [ (6-sqrt(6))/10 ]\n"
2477 << " [ (-96+131*sqrt(6))/625 ]\n"
2478 << " A5 = [ 0 ]\n"
2479 << " [ 0 ]\n"
2480 << " [ 0 ]\n"
2481 << " [ 0 ]\n"
2482 << " [ (6-sqrt(6))/10 ]\n"
2483 << "b = [ 0 ]\n"
2484 << " [ 0 ]\n"
2485 << " [ 1/9 ]\n"
2486 << " [ (16-sqrt(6))/36 ]\n"
2487 << " [ (16+sqrt(6))/36 ]" << std::endl;
2488 typedef ScalarTraits<Scalar> ST;
2489 int myNumStages = 5;
2490 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2491 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2492 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2493 Scalar zero = ST::zero();
2494 Scalar one = ST::one();
2495 Scalar sqrt6 = ST::squareroot(as<Scalar>(6*one));
2496 Scalar gamma = as<Scalar>( (6*one - sqrt6) / (10*one) ); // diagonal
2497 myA(0,0) = gamma;
2498 myA(0,1) = zero;
2499 myA(0,2) = zero;
2500 myA(0,3) = zero;
2501 myA(0,4) = zero;
2502
2503 myA(1,0) = as<Scalar>( (-6*one+5*one*sqrt6)/(14*one) );
2504 myA(1,1) = gamma;
2505 myA(1,2) = zero;
2506 myA(1,3) = zero;
2507 myA(1,4) = zero;
2508
2509 myA(2,0) = as<Scalar>( (888*one+607*one*sqrt6)/(2850*one) );
2510 myA(2,1) = as<Scalar>( (126*one-161*one*sqrt6)/(1425*one) );
2511 myA(2,2) = gamma;
2512 myA(2,3) = zero;
2513 myA(2,4) = zero;
2514
2515 myA(3,0) = as<Scalar>( (3153*one-3082*one*sqrt6)/(14250*one) );
2516 myA(3,1) = as<Scalar>( (3213*one+1148*one*sqrt6)/(28500*one) );
2517 myA(3,2) = as<Scalar>( (-267*one+88*one*sqrt6)/(500*one) );
2518 myA(3,3) = gamma;
2519 myA(3,4) = zero;
2520
2521 myA(4,0) = as<Scalar>( (-32583*one+14638*one*sqrt6)/(71250*one) );
2522 myA(4,1) = as<Scalar>( (-17199*one+364*one*sqrt6)/(142500*one) );
2523 myA(4,2) = as<Scalar>( (1329*one-544*one*sqrt6)/(2500*one) );
2524 myA(4,3) = as<Scalar>( (-96*one+131*sqrt6)/(625*one) );
2525 myA(4,4) = gamma;
2526
2527 myb(0) = zero;
2528 myb(1) = zero;
2529 myb(2) = as<Scalar>( one/(9*one) );
2530 myb(3) = as<Scalar>( (16*one-sqrt6)/(36*one) );
2531 myb(4) = as<Scalar>( (16*one+sqrt6)/(36*one) );
2532
2533 myc(0) = gamma;
2534 myc(1) = as<Scalar>( (6*one+9*one*sqrt6)/(35*one) );
2535 myc(2) = one;
2536 myc(3) = as<Scalar>( (4*one-sqrt6)/(10*one) );
2537 myc(4) = as<Scalar>( (4*one+sqrt6)/(10*one) );
2538
2539 this->setMyDescription(myDescription.str());
2540 this->setMy_A(myA);
2541 this->setMy_b(myb);
2542 this->setMy_c(myc);
2543 this->setMy_order(5);
2544 }
2545};
2546
2547
2548template<class Scalar>
2549class SDIRK5Stage4thOrder_RKBT :
2550 virtual public RKButcherTableauDefaultBase<Scalar>
2551{
2552 public:
2553 SDIRK5Stage4thOrder_RKBT()
2554 {
2555
2556 std::ostringstream myDescription;
2557 myDescription << SDIRK5Stage4thOrder_name() << "\n"
2558 << "L-stable\n"
2559 << "Solving Ordinary Differential Equations II:\n"
2560 << "Stiff and Differential-Algebraic Problems,\n"
2561 << "2nd Revised Edition\n"
2562 << "E. Hairer and G. Wanner\n"
2563 << "pg100 \n"
2564 << "c = [ 1/4 3/4 11/20 1/2 1 ]'\n"
2565 << "A = [ 1/4 ]\n"
2566 << " [ 1/2 1/4 ]\n"
2567 << " [ 17/50 -1/25 1/4 ]\n"
2568 << " [ 371/1360 -137/2720 15/544 1/4 ]\n"
2569 << " [ 25/24 -49/48 125/16 -85/12 1/4 ]\n"
2570 << "b = [ 25/24 -49/48 125/16 -85/12 1/4 ]'\n"
2571 << "b' = [ 59/48 -17/96 225/32 -85/12 0 ]'" << std::endl;
2572 typedef ScalarTraits<Scalar> ST;
2573 int myNumStages = 5;
2574 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2575 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2576 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2577 Scalar zero = ST::zero();
2578 Scalar one = ST::one();
2579 Scalar onequarter = as<Scalar>( one/(4*one) );
2580 myA(0,0) = onequarter;
2581 myA(0,1) = zero;
2582 myA(0,2) = zero;
2583 myA(0,3) = zero;
2584 myA(0,4) = zero;
2585
2586 myA(1,0) = as<Scalar>( one / (2*one) );
2587 myA(1,1) = onequarter;
2588 myA(1,2) = zero;
2589 myA(1,3) = zero;
2590 myA(1,4) = zero;
2591
2592 myA(2,0) = as<Scalar>( 17*one/(50*one) );
2593 myA(2,1) = as<Scalar>( -one/(25*one) );
2594 myA(2,2) = onequarter;
2595 myA(2,3) = zero;
2596 myA(2,4) = zero;
2597
2598 myA(3,0) = as<Scalar>( 371*one/(1360*one) );
2599 myA(3,1) = as<Scalar>( -137*one/(2720*one) );
2600 myA(3,2) = as<Scalar>( 15*one/(544*one) );
2601 myA(3,3) = onequarter;
2602 myA(3,4) = zero;
2603
2604 myA(4,0) = as<Scalar>( 25*one/(24*one) );
2605 myA(4,1) = as<Scalar>( -49*one/(48*one) );
2606 myA(4,2) = as<Scalar>( 125*one/(16*one) );
2607 myA(4,3) = as<Scalar>( -85*one/(12*one) );
2608 myA(4,4) = onequarter;
2609
2610 myb(0) = as<Scalar>( 25*one/(24*one) );
2611 myb(1) = as<Scalar>( -49*one/(48*one) );
2612 myb(2) = as<Scalar>( 125*one/(16*one) );
2613 myb(3) = as<Scalar>( -85*one/(12*one) );
2614 myb(4) = onequarter;
2615
2616 /*
2617 // Alternate version
2618 myb(0) = as<Scalar>( 59*one/(48*one) );
2619 myb(1) = as<Scalar>( -17*one/(96*one) );
2620 myb(2) = as<Scalar>( 225*one/(32*one) );
2621 myb(3) = as<Scalar>( -85*one/(12*one) );
2622 myb(4) = zero;
2623 */
2624 myc(0) = onequarter;
2625 myc(1) = as<Scalar>( 3*one/(4*one) );
2626 myc(2) = as<Scalar>( 11*one/(20*one) );
2627 myc(3) = as<Scalar>( one/(2*one) );
2628 myc(4) = one;
2629
2630 this->setMyDescription(myDescription.str());
2631 this->setMy_A(myA);
2632 this->setMy_b(myb);
2633 this->setMy_c(myc);
2634 this->setMy_order(4);
2635 }
2636};
2637
2638
2639template<class Scalar>
2640class SDIRK3Stage4thOrder_RKBT :
2641 virtual public RKButcherTableauDefaultBase<Scalar>
2642{
2643 public:
2644 SDIRK3Stage4thOrder_RKBT()
2645 {
2646
2647 std::ostringstream myDescription;
2648 myDescription << SDIRK3Stage4thOrder_name() << "\n"
2649 << "A-stable\n"
2650 << "Solving Ordinary Differential Equations II:\n"
2651 << "Stiff and Differential-Algebraic Problems,\n"
2652 << "2nd Revised Edition\n"
2653 << "E. Hairer and G. Wanner\n"
2654 << "pg100 \n"
2655 << "gamma = (1/sqrt(3))*cos(pi/18)+1/2\n"
2656 << "delta = 1/(6*(2*gamma-1)^2)\n"
2657 << "c = [ gamma 1/2 1-gamma ]'\n"
2658 << "A = [ gamma ]\n"
2659 << " [ 1/2-gamma gamma ]\n"
2660 << " [ 2*gamma 1-4*gamma gamma ]\n"
2661 << "b = [ delta 1-2*delta delta ]'" << std::endl;
2662 typedef ScalarTraits<Scalar> ST;
2663 int myNumStages = 3;
2664 Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2665 Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2666 Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2667 Scalar zero = ST::zero();
2668 Scalar one = ST::one();
2669 Scalar pi = as<Scalar>(4*one)*std::atan(one);
2670 Scalar gamma = as<Scalar>( one/ST::squareroot(3*one)*std::cos(pi/(18*one))+one/(2*one) );
2671 Scalar delta = as<Scalar>( one/(6*one*std::pow(2*gamma-one,2*one)) );
2672 myA(0,0) = gamma;
2673 myA(0,1) = zero;
2674 myA(0,2) = zero;
2675
2676 myA(1,0) = as<Scalar>( one/(2*one) - gamma );
2677 myA(1,1) = gamma;
2678 myA(1,2) = zero;
2679
2680 myA(2,0) = as<Scalar>( 2*gamma );
2681 myA(2,1) = as<Scalar>( one - 4*gamma );
2682 myA(2,2) = gamma;
2683
2684 myb(0) = delta;
2685 myb(1) = as<Scalar>( one-2*delta );
2686 myb(2) = delta;
2687
2688 myc(0) = gamma;
2689 myc(1) = as<Scalar>( one/(2*one) );
2690 myc(2) = as<Scalar>( one - gamma );
2691
2692 this->setMyDescription(myDescription.str());
2693 this->setMy_A(myA);
2694 this->setMy_b(myb);
2695 this->setMy_c(myc);
2696 this->setMy_order(4);
2697 }
2698};
2699
2700
2701} // namespace Rythmos
2702
2703
2704#endif // RYTHMOS_RK_BUTCHER_TABLEAU_HPP