Rythmos - Transient Integration for Differential Equations  Version of the Day
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 
52 namespace 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 
108 template<class Scalar>
109 class 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 = " << this->A() << std::endl;
204  out << "b = " << this->b() << std::endl;
205  out << "c = " << 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
239 template<class Scalar>
240 RCP<RKButcherTableauBase<Scalar> > rKButcherTableau()
241 {
242  return(rcp(new RKButcherTableauDefaultBase<Scalar>()));
243 }
244 
245 // Nonmember constructor
246 template<class Scalar>
247 RCP<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 
261 template<class Scalar>
262 class 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 
291 template<class Scalar>
292 class 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 
319 template<class Scalar>
320 class 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 
393 template<class Scalar>
394 class 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 
468 template<class Scalar>
469 class 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 
541 template<class Scalar>
542 class 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 one = ST::one();
568  Scalar onefifth = ST::one()/(5*ST::one());
569  Scalar onefourth = ST::one()/(4*ST::one());
570  Scalar onethird = ST::one()/(3*ST::one());
571  Scalar twothirds = 2*ST::one()/(3*ST::one());
572  Scalar threefourths = 3*ST::one()/(4*ST::one());
573  Scalar zero = ST::zero();
574 
575  // Fill A:
576  myA(0,0) = zero;
577  myA(0,1) = zero;
578  myA(0,2) = zero;
579  myA(0,3) = zero;
580  myA(0,4) = zero;
581 
582  myA(1,0) = onefifth;
583  myA(1,1) = zero;
584  myA(1,2) = zero;
585  myA(1,3) = zero;
586  myA(1,4) = zero;
587 
588  myA(2,0) = zero;
589  myA(2,1) = onefifth;
590  myA(2,2) = zero;
591  myA(2,3) = zero;
592  myA(2,4) = zero;
593 
594  myA(3,0) = zero;
595  myA(3,1) = zero;
596  myA(3,2) = onethird;
597  myA(3,3) = zero;
598  myA(3,4) = zero;
599 
600  myA(4,0) = zero;
601  myA(4,1) = zero;
602  myA(4,2) = zero;
603  myA(4,3) = twothirds;
604  myA(4,4) = zero;
605 
606  // Fill b:
607  myb(0) = onefourth;
608  myb(1) = zero;
609  myb(2) = zero;
610  myb(3) = zero;
611  myb(4) = threefourths;
612 
613  // Fill myc:
614  myc(0) = zero;
615  myc(1) = onefifth;
616  myc(2) = onefifth;
617  myc(3) = onethird;
618  myc(4) = twothirds;
619 
620  this->setMyDescription(myDescription.str());
621  this->setMy_A(myA);
622  this->setMy_b(myb);
623  this->setMy_c(myc);
624  this->setMy_order(3);
625  }
626 };
627 
628 
629 template<class Scalar>
630 class Explicit3Stage3rdOrder_RKBT :
631  virtual public RKButcherTableauDefaultBase<Scalar>
632 {
633  public:
634  Explicit3Stage3rdOrder_RKBT()
635  {
636 
637  std::ostringstream myDescription;
638  myDescription << Explicit3Stage3rdOrder_name() << "\n"
639  << "c = [ 0 1/2 1 ]'\n"
640  << "A = [ 0 ]\n"
641  << " [ 1/2 0 ]\n"
642  << " [ -1 2 0 ]\n"
643  << "b = [ 1/6 4/6 1/6 ]'" << std::endl;
644  typedef ScalarTraits<Scalar> ST;
645  Scalar one = ST::one();
646  Scalar two = as<Scalar>(2*ST::one());
647  Scalar zero = ST::zero();
648  Scalar onehalf = ST::one()/(2*ST::one());
649  Scalar onesixth = ST::one()/(6*ST::one());
650  Scalar foursixth = 4*ST::one()/(6*ST::one());
651 
652  int myNumStages = 3;
653  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
654  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
655  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
656 
657  // Fill myA:
658  myA(0,0) = zero;
659  myA(0,1) = zero;
660  myA(0,2) = zero;
661 
662  myA(1,0) = onehalf;
663  myA(1,1) = zero;
664  myA(1,2) = zero;
665 
666  myA(2,0) = -one;
667  myA(2,1) = two;
668  myA(2,2) = zero;
669 
670  // Fill myb:
671  myb(0) = onesixth;
672  myb(1) = foursixth;
673  myb(2) = onesixth;
674 
675  // fill b_c_
676  myc(0) = zero;
677  myc(1) = onehalf;
678  myc(2) = one;
679 
680  this->setMyDescription(myDescription.str());
681  this->setMy_A(myA);
682  this->setMy_b(myb);
683  this->setMy_c(myc);
684  this->setMy_order(3);
685  }
686 };
687 
688 template<class Scalar>
689 class Explicit2Stage2ndOrderTVD_RKBT :
690  virtual public RKButcherTableauDefaultBase<Scalar>
691 {
692  public:
693  Explicit2Stage2ndOrderTVD_RKBT()
694  {
695 
696  std::ostringstream myDescription;
697  myDescription << Explicit2Stage2ndOrderTVD_name() << "\n"
698  << "Sigal Gottlieb, David Ketcheson and Chi-Wang Shu\n"
699  << "`Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations'\n"
700  << "World Scientific, 2011\n"
701  << "pp. 15\n"
702  << "c = [ 0 1 ]'\n"
703  << "A = [ 0 ]\n"
704  << " [ 1 0 ]\n"
705  << "b = [ 1/2 1/2]'\n"
706  << "This is also written in the following set of updates.\n"
707  << "u1 = u^n + dt L(u^n)\n"
708  << "u^(n+1) = u^n/2 + u1/2 + dt L(u1)/2"
709  << std::endl;
710  typedef ScalarTraits<Scalar> ST;
711  Scalar one = ST::one();
712  Scalar zero = ST::zero();
713  Scalar onehalf = ST::one()/(2*ST::one());
714 
715  int myNumStages = 2;
716  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
717  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
718  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
719 
720  // Fill myA:
721  myA(0,0) = zero;
722  myA(0,1) = zero;
723 
724  myA(1,0) = one;
725  myA(1,1) = zero;
726 
727  // Fill myb:
728  myb(0) = onehalf;
729  myb(1) = onehalf;
730 
731  // fill b_c_
732  myc(0) = zero;
733  myc(1) = one;
734 
735  this->setMyDescription(myDescription.str());
736  this->setMy_A(myA);
737  this->setMy_b(myb);
738  this->setMy_c(myc);
739  this->setMy_order(2);
740  }
741 };
742 
743 template<class Scalar>
744 class Explicit3Stage3rdOrderTVD_RKBT :
745  virtual public RKButcherTableauDefaultBase<Scalar>
746 {
747  public:
748  Explicit3Stage3rdOrderTVD_RKBT()
749  {
750 
751  std::ostringstream myDescription;
752  myDescription << Explicit3Stage3rdOrderTVD_name() << "\n"
753  << "Sigal Gottlieb and Chi-Wang Shu\n"
754  << "`Total Variation Diminishing Runge-Kutta Schemes'\n"
755  << "Mathematics of Computation\n"
756  << "Volume 67, Number 221, January 1998, pp. 73-85\n"
757  << "c = [ 0 1 1/2 ]'\n"
758  << "A = [ 0 ]\n"
759  << " [ 1 0 ]\n"
760  << " [ 1/4 1/4 0 ]\n"
761  << "b = [ 1/6 1/6 4/6 ]'\n"
762  << "This is also written in the following set of updates.\n"
763  << "u1 = u^n + dt L(u^n)\n"
764  << "u2 = 3 u^n/4 + u1/4 + dt L(u1)/4\n"
765  << "u^(n+1) = u^n/3 + 2 u2/2 + 2 dt L(u2)/3"
766  << std::endl;
767  typedef ScalarTraits<Scalar> ST;
768  Scalar one = ST::one();
769  Scalar zero = ST::zero();
770  Scalar onehalf = ST::one()/(2*ST::one());
771  Scalar onefourth = ST::one()/(4*ST::one());
772  Scalar onesixth = ST::one()/(6*ST::one());
773  Scalar foursixth = 4*ST::one()/(6*ST::one());
774 
775  int myNumStages = 3;
776  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
777  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
778  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
779 
780  // Fill myA:
781  myA(0,0) = zero;
782  myA(0,1) = zero;
783  myA(0,2) = zero;
784 
785  myA(1,0) = one;
786  myA(1,1) = zero;
787  myA(1,2) = zero;
788 
789  myA(2,0) = onefourth;
790  myA(2,1) = onefourth;
791  myA(2,2) = zero;
792 
793  // Fill myb:
794  myb(0) = onesixth;
795  myb(1) = onesixth;
796  myb(2) = foursixth;
797 
798  // fill b_c_
799  myc(0) = zero;
800  myc(1) = one;
801  myc(2) = onehalf;
802 
803  this->setMyDescription(myDescription.str());
804  this->setMy_A(myA);
805  this->setMy_b(myb);
806  this->setMy_c(myc);
807  this->setMy_order(3);
808  }
809 };
810 
811 
812 template<class Scalar>
813 class Explicit3Stage3rdOrderHeun_RKBT :
814  virtual public RKButcherTableauDefaultBase<Scalar>
815 {
816  public:
817  Explicit3Stage3rdOrderHeun_RKBT()
818  {
819  std::ostringstream myDescription;
820  myDescription << Explicit3Stage3rdOrderHeun_name() << "\n"
821  << "Solving Ordinary Differential Equations I:\n"
822  << "Nonstiff Problems, 2nd Revised Edition\n"
823  << "E. Hairer, S.P. Norsett, G. Wanner\n"
824  << "Table 1.1, pg 135\n"
825  << "c = [ 0 1/3 2/3 ]'\n"
826  << "A = [ 0 ] \n"
827  << " [ 1/3 0 ]\n"
828  << " [ 0 2/3 0 ]\n"
829  << "b = [ 1/4 0 3/4 ]'" << std::endl;
830  typedef ScalarTraits<Scalar> ST;
831  Scalar one = ST::one();
832  Scalar zero = ST::zero();
833  Scalar onethird = one/(3*one);
834  Scalar twothirds = 2*one/(3*one);
835  Scalar onefourth = one/(4*one);
836  Scalar threefourths = 3*one/(4*one);
837 
838  int myNumStages = 3;
839  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
840  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
841  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
842 
843  // Fill myA:
844  myA(0,0) = zero;
845  myA(0,1) = zero;
846  myA(0,2) = zero;
847 
848  myA(1,0) = onethird;
849  myA(1,1) = zero;
850  myA(1,2) = zero;
851 
852  myA(2,0) = zero;
853  myA(2,1) = twothirds;
854  myA(2,2) = zero;
855 
856  // Fill myb:
857  myb(0) = onefourth;
858  myb(1) = zero;
859  myb(2) = threefourths;
860 
861  // fill b_c_
862  myc(0) = zero;
863  myc(1) = onethird;
864  myc(2) = twothirds;
865 
866  this->setMyDescription(myDescription.str());
867  this->setMy_A(myA);
868  this->setMy_b(myb);
869  this->setMy_c(myc);
870  this->setMy_order(3);
871  }
872 };
873 
874 
875 template<class Scalar>
876 class Explicit2Stage2ndOrderRunge_RKBT :
877  virtual public RKButcherTableauDefaultBase<Scalar>
878 {
879  public:
880  Explicit2Stage2ndOrderRunge_RKBT()
881  {
882  std::ostringstream myDescription;
883  myDescription << Explicit2Stage2ndOrderRunge_name() << "\n"
884  << "Also known as Explicit Midpoint\n"
885  << "Solving Ordinary Differential Equations I:\n"
886  << "Nonstiff Problems, 2nd Revised Edition\n"
887  << "E. Hairer, S.P. Norsett, G. Wanner\n"
888  << "Table 1.1, pg 135\n"
889  << "c = [ 0 1/2 ]'\n"
890  << "A = [ 0 ]\n"
891  << " [ 1/2 0 ]\n"
892  << "b = [ 0 1 ]'" << std::endl;
893  typedef ScalarTraits<Scalar> ST;
894  Scalar one = ST::one();
895  Scalar zero = ST::zero();
896  Scalar onehalf = ST::one()/(2*ST::one());
897 
898  int myNumStages = 2;
899  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
900  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
901  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
902 
903  // Fill myA:
904  myA(0,0) = zero;
905  myA(0,1) = zero;
906 
907  myA(1,0) = onehalf;
908  myA(1,1) = zero;
909 
910  // Fill myb:
911  myb(0) = zero;
912  myb(1) = one;
913 
914  // fill b_c_
915  myc(0) = zero;
916  myc(1) = onehalf;
917 
918  this->setMyDescription(myDescription.str());
919  this->setMy_A(myA);
920  this->setMy_b(myb);
921  this->setMy_c(myc);
922  this->setMy_order(2);
923  }
924 };
925 
926 
927 template<class Scalar>
928 class ExplicitTrapezoidal_RKBT :
929  virtual public RKButcherTableauDefaultBase<Scalar>
930 {
931  public:
932  ExplicitTrapezoidal_RKBT()
933  {
934  std::ostringstream myDescription;
935  myDescription << ExplicitTrapezoidal_name() << "\n"
936  << "c = [ 0 1 ]'\n"
937  << "A = [ 0 ]\n"
938  << " [ 1 0 ]\n"
939  << "b = [ 1/2 1/2 ]'" << std::endl;
940  typedef ScalarTraits<Scalar> ST;
941  Scalar one = ST::one();
942  Scalar zero = ST::zero();
943  Scalar onehalf = ST::one()/(2*ST::one());
944 
945  int myNumStages = 2;
946  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
947  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
948  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
949 
950  // Fill myA:
951  myA(0,0) = zero;
952  myA(0,1) = zero;
953 
954  myA(1,0) = one;
955  myA(1,1) = zero;
956 
957  // Fill myb:
958  myb(0) = onehalf;
959  myb(1) = onehalf;
960 
961  // fill b_c_
962  myc(0) = zero;
963  myc(1) = one;
964 
965  this->setMyDescription(myDescription.str());
966  this->setMy_A(myA);
967  this->setMy_b(myb);
968  this->setMy_c(myc);
969  this->setMy_order(2);
970  }
971 };
972 
973 
974 template<class Scalar>
975 class SDIRK2Stage2ndOrder_RKBT :
976  virtual public RKButcherTableauDefaultBase<Scalar>
977 {
978  public:
979  SDIRK2Stage2ndOrder_RKBT()
980  {
981  std::ostringstream myDescription;
982  myDescription << SDIRK2Stage2ndOrder_name() << "\n"
983  << "Computer Methods for ODEs and DAEs\n"
984  << "U. M. Ascher and L. R. Petzold\n"
985  << "p. 106\n"
986  << "gamma = (2+-sqrt(2))/2\n"
987  << "c = [ gamma 1 ]'\n"
988  << "A = [ gamma 0 ]\n"
989  << " [ 1-gamma gamma ]\n"
990  << "b = [ 1-gamma gamma ]'" << std::endl;
991 
992  this->setMyDescription(myDescription.str());
993  typedef ScalarTraits<Scalar> ST;
994  Scalar one = ST::one();
995  gamma_default_ = as<Scalar>( (2*one - ST::squareroot(2*one))/(2*one) );
996  gamma_ = gamma_default_;
997  this->setupData();
998 
999  RCP<ParameterList> validPL = Teuchos::parameterList();
1000  validPL->set("Description","",this->getMyDescription());
1001  validPL->set<double>("gamma",gamma_default_,
1002  "The default value is gamma = (2-sqrt(2))/2. "
1003  "This will produce an L-stable 2nd order method with the stage "
1004  "times within the timestep. Other values of gamma will still "
1005  "produce an L-stable scheme, but will only be 1st order accurate.");
1006  Teuchos::setupVerboseObjectSublist(&*validPL);
1007  this->setMyValidParameterList(validPL);
1008  }
1009  void setupData()
1010  {
1011  typedef ScalarTraits<Scalar> ST;
1012  int myNumStages = 2;
1013  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1014  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1015  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1016  Scalar one = ST::one();
1017  Scalar zero = ST::zero();
1018  myA(0,0) = gamma_;
1019  myA(0,1) = zero;
1020  myA(1,0) = as<Scalar>( one - gamma_ );
1021  myA(1,1) = gamma_;
1022  myb(0) = as<Scalar>( one - gamma_ );
1023  myb(1) = gamma_;
1024  myc(0) = gamma_;
1025  myc(1) = one;
1026 
1027  this->setMy_A(myA);
1028  this->setMy_b(myb);
1029  this->setMy_c(myc);
1030  this->setMy_order(2);
1031  }
1032  void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
1033  {
1034  TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1035  paramList->validateParameters(*this->getValidParameters());
1036  Teuchos::readVerboseObjectSublist(&*paramList,this);
1037  gamma_ = paramList->get<double>("gamma",gamma_default_);
1038  this->setupData();
1039  this->setMyParamList(paramList);
1040  }
1041  private:
1042  Scalar gamma_default_;
1043  Scalar gamma_;
1044 };
1045 
1046 
1047 // 04/07/09 tscoffe: I verified manually that the Convergence Testing passes
1048 // with gamma_default_ = -1.
1049 template<class Scalar>
1050 class SDIRK2Stage3rdOrder_RKBT :
1051  virtual public RKButcherTableauDefaultBase<Scalar>
1052 {
1053  public:
1054  SDIRK2Stage3rdOrder_RKBT()
1055  {
1056  std::ostringstream myDescription;
1057  myDescription << SDIRK2Stage3rdOrder_name() << "\n"
1058  << "Solving Ordinary Differential Equations I:\n"
1059  << "Nonstiff Problems, 2nd Revised Edition\n"
1060  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1061  << "Table 7.2, pg 207\n"
1062  << "gamma = (3+-sqrt(3))/6 -> 3rd order and A-stable\n"
1063  << "gamma = (2+-sqrt(2))/2 -> 2nd order and L-stable\n"
1064  << "c = [ gamma 1-gamma ]'\n"
1065  << "A = [ gamma 0 ]\n"
1066  << " [ 1-2*gamma gamma ]\n"
1067  << "b = [ 1/2 1/2 ]'" << std::endl;
1068 
1069  this->setMyDescription(myDescription.str());
1070  thirdOrderAStable_default_ = true;
1071  secondOrderLStable_default_ = false;
1072  thirdOrderAStable_ = true;
1073  secondOrderLStable_ = false;
1074  typedef ScalarTraits<Scalar> ST;
1075  Scalar one = ST::one();
1076  gamma_default_ = as<Scalar>( (3*one + ST::squareroot(3*one))/(6*one) );
1077  gamma_ = gamma_default_;
1078  this->setupData();
1079 
1080  RCP<ParameterList> validPL = Teuchos::parameterList();
1081  validPL->set("Description","",this->getMyDescription());
1082  validPL->set("3rd Order A-stable",thirdOrderAStable_default_,
1083  "If true, set gamma to gamma = (3+sqrt(3))/6 to obtain "
1084  "a 3rd order A-stable scheme. '3rd Order A-stable' and "
1085  "'2nd Order L-stable' can not both be true.");
1086  validPL->set("2nd Order L-stable",secondOrderLStable_default_,
1087  "If true, set gamma to gamma = (2+sqrt(2))/2 to obtain "
1088  "a 2nd order L-stable scheme. '3rd Order A-stable' and "
1089  "'2nd Order L-stable' can not both be true.");
1090  validPL->set("gamma",gamma_default_,
1091  "If both '3rd Order A-stable' and '2nd Order L-stable' "
1092  "are false, gamma will be used. The default value is the "
1093  "'3rd Order A-stable' gamma value, (3+sqrt(3))/6.");
1094 
1095  Teuchos::setupVerboseObjectSublist(&*validPL);
1096  this->setMyValidParameterList(validPL);
1097  }
1098  void setupData()
1099  {
1100  typedef ScalarTraits<Scalar> ST;
1101  int myNumStages = 2;
1102  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1103  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1104  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1105  Scalar one = ST::one();
1106  Scalar zero = ST::zero();
1107  Scalar gamma;
1108  if (thirdOrderAStable_)
1109  gamma = as<Scalar>( (3*one + ST::squareroot(3*one))/(6*one) );
1110  else if (secondOrderLStable_)
1111  gamma = as<Scalar>( (2*one + ST::squareroot(2*one))/(2*one) );
1112  else
1113  gamma = gamma_;
1114  myA(0,0) = gamma;
1115  myA(0,1) = zero;
1116  myA(1,0) = as<Scalar>( one - 2*gamma );
1117  myA(1,1) = gamma;
1118  myb(0) = as<Scalar>( one/(2*one) );
1119  myb(1) = as<Scalar>( one/(2*one) );
1120  myc(0) = gamma;
1121  myc(1) = as<Scalar>( one - gamma );
1122 
1123  this->setMy_A(myA);
1124  this->setMy_b(myb);
1125  this->setMy_c(myc);
1126  this->setMy_order(3);
1127  }
1128  void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
1129  {
1130  TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1131  paramList->validateParameters(*this->getValidParameters());
1132  Teuchos::readVerboseObjectSublist(&*paramList,this);
1133  thirdOrderAStable_ = paramList->get("3rd Order A-stable",
1134  thirdOrderAStable_default_);
1135  secondOrderLStable_ = paramList->get("2nd Order L-stable",
1136  secondOrderLStable_default_);
1137  TEUCHOS_TEST_FOR_EXCEPTION(
1138  thirdOrderAStable_ && secondOrderLStable_, std::logic_error,
1139  "'3rd Order A-stable' and '2nd Order L-stable' can not both be true.");
1140  gamma_ = paramList->get("gamma",gamma_default_);
1141 
1142  this->setupData();
1143  this->setMyParamList(paramList);
1144  }
1145  private:
1146  bool thirdOrderAStable_default_;
1147  bool thirdOrderAStable_;
1148  bool secondOrderLStable_default_;
1149  bool secondOrderLStable_;
1150  Scalar gamma_default_;
1151  Scalar gamma_;
1152 };
1153 
1154 
1155 template<class Scalar>
1156 class DIRK2Stage3rdOrder_RKBT :
1157  virtual public RKButcherTableauDefaultBase<Scalar>
1158 {
1159  public:
1160  DIRK2Stage3rdOrder_RKBT()
1161  {
1162 
1163  std::ostringstream myDescription;
1164  myDescription << DIRK2Stage3rdOrder_name() << "\n"
1165  << "Hammer & Hollingsworth method\n"
1166  << "Solving Ordinary Differential Equations I:\n"
1167  << "Nonstiff Problems, 2nd Revised Edition\n"
1168  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1169  << "Table 7.1, pg 205\n"
1170  << "c = [ 0 2/3 ]'\n"
1171  << "A = [ 0 0 ]\n"
1172  << " [ 1/3 1/3 ]\n"
1173  << "b = [ 1/4 3/4 ]'" << std::endl;
1174  typedef ScalarTraits<Scalar> ST;
1175  int myNumStages = 2;
1176  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1177  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1178  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1179  Scalar one = ST::one();
1180  Scalar zero = ST::zero();
1181  myA(0,0) = zero;
1182  myA(0,1) = zero;
1183  myA(1,0) = as<Scalar>( one/(3*one) );
1184  myA(1,1) = as<Scalar>( one/(3*one) );
1185  myb(0) = as<Scalar>( one/(4*one) );
1186  myb(1) = as<Scalar>( 3*one/(4*one) );
1187  myc(0) = zero;
1188  myc(1) = as<Scalar>( 2*one/(3*one) );
1189  this->setMyDescription(myDescription.str());
1190  this->setMy_A(myA);
1191  this->setMy_b(myb);
1192  this->setMy_c(myc);
1193  this->setMy_order(3);
1194  }
1195 };
1196 
1197 
1198 template<class Scalar>
1199 class Implicit3Stage6thOrderKuntzmannButcher_RKBT :
1200  virtual public RKButcherTableauDefaultBase<Scalar>
1201 {
1202  public:
1203  Implicit3Stage6thOrderKuntzmannButcher_RKBT()
1204  {
1205 
1206  std::ostringstream myDescription;
1207  myDescription << Implicit3Stage6thOrderKuntzmannButcher_name() << "\n"
1208  << "Kuntzmann & Butcher method\n"
1209  << "Solving Ordinary Differential Equations I:\n"
1210  << "Nonstiff Problems, 2nd Revised Edition\n"
1211  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1212  << "Table 7.4, pg 209\n"
1213  << "c = [ 1/2-sqrt(15)/10 1/2 1/2-sqrt(15)/10 ]'\n"
1214  << "A = [ 5/36 2/9-sqrt(15)/15 5/36-sqrt(15)/30 ]\n"
1215  << " [ 5/36+sqrt(15)/24 2/9 5/36-sqrt(15)/24 ]\n"
1216  << " [ 5/36+sqrt(15)/30 2/9+sqrt(15)/15 5/36 ]\n"
1217  << "b = [ 5/18 4/9 5/18 ]'" << std::endl;
1218  typedef ScalarTraits<Scalar> ST;
1219  int myNumStages = 3;
1220  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1221  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1222  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1223  Scalar one = ST::one();
1224  myA(0,0) = as<Scalar>( 5*one/(36*one) );
1225  myA(0,1) = as<Scalar>( 2*one/(9*one) - ST::squareroot(15*one)/(15*one) );
1226  myA(0,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(30*one) );
1227  myA(1,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(24*one) );
1228  myA(1,1) = as<Scalar>( 2*one/(9*one) );
1229  myA(1,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(24*one) );
1230  myA(2,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(30*one) );
1231  myA(2,1) = as<Scalar>( 2*one/(9*one) + ST::squareroot(15*one)/(15*one) );
1232  myA(2,2) = as<Scalar>( 5*one/(36*one) );
1233  myb(0) = as<Scalar>( 5*one/(18*one) );
1234  myb(1) = as<Scalar>( 4*one/(9*one) );
1235  myb(2) = as<Scalar>( 5*one/(18*one) );
1236  myc(0) = as<Scalar>( one/(2*one)-ST::squareroot(15*one)/(10*one) );
1237  myc(1) = as<Scalar>( one/(2*one) );
1238  myc(2) = as<Scalar>( one/(2*one)+ST::squareroot(15*one)/(10*one) );
1239  this->setMyDescription(myDescription.str());
1240  this->setMy_A(myA);
1241  this->setMy_b(myb);
1242  this->setMy_c(myc);
1243  this->setMy_order(6);
1244  }
1245 };
1246 
1247 
1248 template<class Scalar>
1249 class Implicit4Stage8thOrderKuntzmannButcher_RKBT :
1250  virtual public RKButcherTableauDefaultBase<Scalar>
1251 {
1252  public:
1253  Implicit4Stage8thOrderKuntzmannButcher_RKBT()
1254  {
1255 
1256  std::ostringstream myDescription;
1257  myDescription << Implicit4Stage8thOrderKuntzmannButcher_name() << "\n"
1258  << "Kuntzmann & Butcher method\n"
1259  << "Solving Ordinary Differential Equations I:\n"
1260  << "Nonstiff Problems, 2nd Revised Edition\n"
1261  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1262  << "Table 7.5, pg 209\n"
1263  << "c = [ 1/2-w2 1/2-w2p 1/2+w2p 1/2+w2 ]'\n"
1264  << "A = [ w1 w1p-w3+w4p w1p-w3-w4p w1-w5 ]\n"
1265  << " [ w1-w3p+w4 w1p w1p-w5p w1-w3p-w4 ]\n"
1266  << " [ w1+w3p+w4 w1p+w5p w1p w1+w3p-w4 ]\n"
1267  << " [ w1+w5 w1p+w3+w4p w1p+w3-w4p w1 ]\n"
1268  << "b = [ 2*w1 2*w1p 2*w1p 2*w1 ]'\n"
1269  << "w1 = 1/8-sqrt(30)/144\n"
1270  << "w2 = (1/2)*sqrt((15+2*sqrt(3))/35)\n"
1271  << "w3 = w2*(1/6+sqrt(30)/24)\n"
1272  << "w4 = w2*(1/21+5*sqrt(30)/168)\n"
1273  << "w5 = w2-2*w3\n"
1274  << "w1p = 1/8+sqrt(30)/144\n"
1275  << "w2p = (1/2)*sqrt((15-2*sqrt(3))/35)\n"
1276  << "w3p = w2*(1/6-sqrt(30)/24)\n"
1277  << "w4p = w2*(1/21-5*sqrt(30)/168)\n"
1278  << "w5p = w2p-2*w3p" << std::endl;
1279  typedef ScalarTraits<Scalar> ST;
1280  int myNumStages = 4;
1281  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1282  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1283  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1284  Scalar one = ST::one();
1285  Scalar onehalf = as<Scalar>( one/(2*one) );
1286  Scalar w1 = as<Scalar>( one/(8*one) - ST::squareroot(30*one)/(144*one) );
1287  Scalar w2 = as<Scalar>( (one/(2*one))*ST::squareroot((15*one+2*one*ST::squareroot(30*one))/(35*one)) );
1288  Scalar w3 = as<Scalar>( w2*(one/(6*one)+ST::squareroot(30*one)/(24*one)) );
1289  Scalar w4 = as<Scalar>( w2*(one/(21*one)+5*one*ST::squareroot(30*one)/(168*one)) );
1290  Scalar w5 = as<Scalar>( w2-2*w3 );
1291  Scalar w1p = as<Scalar>( one/(8*one) + ST::squareroot(30*one)/(144*one) );
1292  Scalar w2p = as<Scalar>( (one/(2*one))*ST::squareroot((15*one-2*one*ST::squareroot(30*one))/(35*one)) );
1293  Scalar w3p = as<Scalar>( w2p*(one/(6*one)-ST::squareroot(30*one)/(24*one)) );
1294  Scalar w4p = as<Scalar>( w2p*(one/(21*one)-5*one*ST::squareroot(30*one)/(168*one)) );
1295  Scalar w5p = as<Scalar>( w2p-2*w3p );
1296  myA(0,0) = w1;
1297  myA(0,1) = w1p-w3+w4p;
1298  myA(0,2) = w1p-w3-w4p;
1299  myA(0,3) = w1-w5;
1300  myA(1,0) = w1-w3p+w4;
1301  myA(1,1) = w1p;
1302  myA(1,2) = w1p-w5p;
1303  myA(1,3) = w1-w3p-w4;
1304  myA(2,0) = w1+w3p+w4;
1305  myA(2,1) = w1p+w5p;
1306  myA(2,2) = w1p;
1307  myA(2,3) = w1+w3p-w4;
1308  myA(3,0) = w1+w5;
1309  myA(3,1) = w1p+w3+w4p;
1310  myA(3,2) = w1p+w3-w4p;
1311  myA(3,3) = w1;
1312  myb(0) = 2*w1;
1313  myb(1) = 2*w1p;
1314  myb(2) = 2*w1p;
1315  myb(3) = 2*w1;
1316  myc(0) = onehalf - w2;
1317  myc(1) = onehalf - w2p;
1318  myc(2) = onehalf + w2p;
1319  myc(3) = onehalf + w2;
1320  this->setMyDescription(myDescription.str());
1321  this->setMy_A(myA);
1322  this->setMy_b(myb);
1323  this->setMy_c(myc);
1324  this->setMy_order(8);
1325  }
1326 };
1327 
1328 
1329 template<class Scalar>
1330 class Implicit2Stage4thOrderHammerHollingsworth_RKBT :
1331  virtual public RKButcherTableauDefaultBase<Scalar>
1332 {
1333  public:
1334  Implicit2Stage4thOrderHammerHollingsworth_RKBT()
1335  {
1336 
1337  std::ostringstream myDescription;
1338  myDescription << Implicit2Stage4thOrderHammerHollingsworth_name() << "\n"
1339  << "Hammer & Hollingsworth method\n"
1340  << "Solving Ordinary Differential Equations I:\n"
1341  << "Nonstiff Problems, 2nd Revised Edition\n"
1342  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1343  << "Table 7.3, pg 207\n"
1344  << "c = [ 1/2-sqrt(3)/6 1/2+sqrt(3)/6 ]'\n"
1345  << "A = [ 1/4 1/4-sqrt(3)/6 ]\n"
1346  << " [ 1/4+sqrt(3)/6 1/4 ]\n"
1347  << "b = [ 1/2 1/2 ]'" << std::endl;
1348  typedef ScalarTraits<Scalar> ST;
1349  int myNumStages = 2;
1350  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1351  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1352  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1353  Scalar one = ST::one();
1354  Scalar onequarter = as<Scalar>( one/(4*one) );
1355  Scalar onehalf = as<Scalar>( one/(2*one) );
1356  myA(0,0) = onequarter;
1357  myA(0,1) = as<Scalar>( onequarter-ST::squareroot(3*one)/(6*one) );
1358  myA(1,0) = as<Scalar>( onequarter+ST::squareroot(3*one)/(6*one) );
1359  myA(1,1) = onequarter;
1360  myb(0) = onehalf;
1361  myb(1) = onehalf;
1362  myc(0) = as<Scalar>( onehalf - ST::squareroot(3*one)/(6*one) );
1363  myc(1) = as<Scalar>( onehalf + ST::squareroot(3*one)/(6*one) );
1364  this->setMyDescription(myDescription.str());
1365  this->setMy_A(myA);
1366  this->setMy_b(myb);
1367  this->setMy_c(myc);
1368  this->setMy_order(4);
1369  }
1370 };
1371 
1372 
1373 template<class Scalar>
1374 class IRK1StageTheta_RKBT :
1375  virtual public RKButcherTableauDefaultBase<Scalar>
1376 {
1377  public:
1378  IRK1StageTheta_RKBT()
1379  {
1380 
1381  std::ostringstream myDescription;
1382  myDescription << IRK1StageTheta_name() << "\n"
1383  << "Non-standard finite-difference methods\n"
1384  << "in dynamical systems, P. Kama,\n"
1385  << "Dissertation, University of Pretoria, pg. 49.\n"
1386  << "Comment: Generalized Implicit Midpoint Method\n"
1387  << "c = [ theta ]'\n"
1388  << "A = [ theta ]\n"
1389  << "b = [ 1 ]'\n"
1390  << std::endl;
1391 
1392  this->setMyDescription(myDescription.str());
1393  typedef ScalarTraits<Scalar> ST;
1394  theta_default_ = ST::one()/(2*ST::one());
1395  theta_ = theta_default_;
1396  this->setupData();
1397 
1398  RCP<ParameterList> validPL = Teuchos::parameterList();
1399  validPL->set("Description","",this->getMyDescription());
1400  validPL->set<double>("theta",theta_default_,
1401  "Valid values are 0 <= theta <= 1, where theta = 0 "
1402  "implies Forward Euler, theta = 1/2 implies midpoint "
1403  "method, and theta = 1 implies Backward Euler. "
1404  "For theta != 1/2, this method is first-order accurate, "
1405  "and with theta = 1/2, it is second-order accurate. "
1406  "This method is A-stable, but becomes L-stable with theta=1.");
1407  Teuchos::setupVerboseObjectSublist(&*validPL);
1408  this->setMyValidParameterList(validPL);
1409  }
1410 
1411  void setupData()
1412  {
1413  typedef ScalarTraits<Scalar> ST;
1414  int myNumStages = 1;
1415  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1416  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1417  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1418  myA(0,0) = theta_;
1419  myb(0) = ST::one();
1420  myc(0) = theta_;
1421  this->setMy_A(myA);
1422  this->setMy_b(myb);
1423  this->setMy_c(myc);
1424  if (theta_ == theta_default_) this->setMy_order(2);
1425  else this->setMy_order(1);
1426  }
1427 
1428  void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
1429  {
1430  TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1431  paramList->validateParameters(*this->getValidParameters());
1432  Teuchos::readVerboseObjectSublist(&*paramList,this);
1433  theta_ = paramList->get<double>("theta",theta_default_);
1434  this->setupData();
1435  this->setMyParamList(paramList);
1436  }
1437  private:
1438  Scalar theta_default_;
1439  Scalar theta_;
1440 };
1441 
1442 
1443 template<class Scalar>
1444 class IRK2StageTheta_RKBT :
1445  virtual public RKButcherTableauDefaultBase<Scalar>
1446 {
1447  public:
1448  IRK2StageTheta_RKBT()
1449  {
1450  std::ostringstream myDescription;
1451  myDescription << IRK2StageTheta_name() << "\n"
1452  << "Computer Methods for ODEs and DAEs\n"
1453  << "U. M. Ascher and L. R. Petzold\n"
1454  << "p. 113\n"
1455  << "c = [ 0 1 ]'\n"
1456  << "A = [ 0 0 ]\n"
1457  << " [ 1-theta theta ]\n"
1458  << "b = [ 1-theta theta ]'\n"
1459  << std::endl;
1460 
1461  this->setMyDescription(myDescription.str());
1462  typedef ScalarTraits<Scalar> ST;
1463  theta_default_ = ST::one()/(2*ST::one());
1464  theta_ = theta_default_;
1465  this->setupData();
1466 
1467  RCP<ParameterList> validPL = Teuchos::parameterList();
1468  validPL->set("Description","",this->getMyDescription());
1469  validPL->set<double>("theta",theta_default_,
1470  "Valid values are 0 < theta <= 1, where theta = 0 "
1471  "implies Forward Euler, theta = 1/2 implies trapezoidal "
1472  "method, and theta = 1 implies Backward Euler. "
1473  "For theta != 1/2, this method is first-order accurate, "
1474  "and with theta = 1/2, it is second-order accurate. "
1475  "This method is A-stable, but becomes L-stable with theta=1.");
1476  Teuchos::setupVerboseObjectSublist(&*validPL);
1477  this->setMyValidParameterList(validPL);
1478  }
1479  void setupData()
1480  {
1481  typedef ScalarTraits<Scalar> ST;
1482  int myNumStages = 2;
1483  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1484  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1485  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1486  Scalar one = ST::one();
1487  Scalar zero = ST::zero();
1488  myA(0,0) = zero;
1489  myA(0,1) = zero;
1490  myA(1,0) = as<Scalar>( one - theta_ );
1491  myA(1,1) = theta_;
1492  myb(0) = as<Scalar>( one - theta_ );
1493  myb(1) = theta_;
1494  myc(0) = theta_;
1495  myc(1) = one;
1496 
1497  this->setMy_A(myA);
1498  this->setMy_b(myb);
1499  this->setMy_c(myc);
1500  TEUCHOS_TEST_FOR_EXCEPTION(
1501  theta_ == zero, std::logic_error,
1502  "'theta' can not be zero, as it makes this IRK stepper explicit.");
1503  if (theta_ == theta_default_) this->setMy_order(2);
1504  else this->setMy_order(1);
1505  }
1506  void setParameterList(RCP<Teuchos::ParameterList> const& paramList)
1507  {
1508  TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
1509  paramList->validateParameters(*this->getValidParameters());
1510  Teuchos::readVerboseObjectSublist(&*paramList,this);
1511  theta_ = paramList->get<double>("theta",theta_default_);
1512  this->setupData();
1513  this->setMyParamList(paramList);
1514  }
1515  private:
1516  Scalar theta_default_;
1517  Scalar theta_;
1518 };
1519 
1520 
1521 template<class Scalar>
1522 class Implicit1Stage2ndOrderGauss_RKBT :
1523  virtual public RKButcherTableauDefaultBase<Scalar>
1524 {
1525  public:
1526  Implicit1Stage2ndOrderGauss_RKBT()
1527  {
1528 
1529  std::ostringstream myDescription;
1530  myDescription << Implicit1Stage2ndOrderGauss_name() << "\n"
1531  << "A-stable\n"
1532  << "Solving Ordinary Differential Equations II:\n"
1533  << "Stiff and Differential-Algebraic Problems,\n"
1534  << "2nd Revised Edition\n"
1535  << "E. Hairer and G. Wanner\n"
1536  << "Table 5.2, pg 72\n"
1537  << "Also: Implicit midpoint rule\n"
1538  << "Solving Ordinary Differential Equations I:\n"
1539  << "Nonstiff Problems, 2nd Revised Edition\n"
1540  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1541  << "Table 7.1, pg 205\n"
1542  << "c = [ 1/2 ]'\n"
1543  << "A = [ 1/2 ]\n"
1544  << "b = [ 1 ]'" << std::endl;
1545  typedef ScalarTraits<Scalar> ST;
1546  int myNumStages = 1;
1547  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1548  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1549  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1550  Scalar onehalf = ST::one()/(2*ST::one());
1551  Scalar one = ST::one();
1552  myA(0,0) = onehalf;
1553  myb(0) = one;
1554  myc(0) = onehalf;
1555  this->setMyDescription(myDescription.str());
1556  this->setMy_A(myA);
1557  this->setMy_b(myb);
1558  this->setMy_c(myc);
1559  this->setMy_order(2);
1560  }
1561 };
1562 
1563 
1564 template<class Scalar>
1565 class Implicit2Stage4thOrderGauss_RKBT :
1566  virtual public RKButcherTableauDefaultBase<Scalar>
1567 {
1568  public:
1569  Implicit2Stage4thOrderGauss_RKBT()
1570  {
1571 
1572  std::ostringstream myDescription;
1573  myDescription << Implicit2Stage4thOrderGauss_name() << "\n"
1574  << "A-stable\n"
1575  << "Solving Ordinary Differential Equations II:\n"
1576  << "Stiff and Differential-Algebraic Problems,\n"
1577  << "2nd Revised Edition\n"
1578  << "E. Hairer and G. Wanner\n"
1579  << "Table 5.2, pg 72\n"
1580  << "c = [ 1/2-sqrt(3)/6 1/2+sqrt(3)/6 ]'\n"
1581  << "A = [ 1/4 1/4-sqrt(3)/6 ]\n"
1582  << " [ 1/4+sqrt(3)/6 1/4 ]\n"
1583  << "b = [ 1/2 1/2 ]'" << std::endl;
1584  typedef ScalarTraits<Scalar> ST;
1585  int myNumStages = 2;
1586  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1587  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1588  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1589  Scalar one = ST::one();
1590  Scalar onehalf = as<Scalar>(one/(2*one));
1591  Scalar three = as<Scalar>(3*one);
1592  Scalar six = as<Scalar>(6*one);
1593  Scalar onefourth = as<Scalar>(one/(4*one));
1594  Scalar alpha = ST::squareroot(three)/six;
1595 
1596  myA(0,0) = onefourth;
1597  myA(0,1) = onefourth-alpha;
1598  myA(1,0) = onefourth+alpha;
1599  myA(1,1) = onefourth;
1600  myb(0) = onehalf;
1601  myb(1) = onehalf;
1602  myc(0) = onehalf-alpha;
1603  myc(1) = onehalf+alpha;
1604  this->setMyDescription(myDescription.str());
1605  this->setMy_A(myA);
1606  this->setMy_b(myb);
1607  this->setMy_c(myc);
1608  this->setMy_order(4);
1609  }
1610 };
1611 
1612 
1613 template<class Scalar>
1614 class Implicit3Stage6thOrderGauss_RKBT :
1615  virtual public RKButcherTableauDefaultBase<Scalar>
1616 {
1617  public:
1618  Implicit3Stage6thOrderGauss_RKBT()
1619  {
1620 
1621  std::ostringstream myDescription;
1622  myDescription << Implicit3Stage6thOrderGauss_name() << "\n"
1623  << "A-stable\n"
1624  << "Solving Ordinary Differential Equations II:\n"
1625  << "Stiff and Differential-Algebraic Problems,\n"
1626  << "2nd Revised Edition\n"
1627  << "E. Hairer and G. Wanner\n"
1628  << "Table 5.2, pg 72\n"
1629  << "c = [ 1/2-sqrt(15)/10 1/2 1/2+sqrt(15)/10 ]'\n"
1630  << "A = [ 5/36 2/9-sqrt(15)/15 5/36-sqrt(15)/30 ]\n"
1631  << " [ 5/36+sqrt(15)/24 2/9 5/36-sqrt(15)/24 ]\n"
1632  << " [ 5/36+sqrt(15)/30 2/9+sqrt(15)/15 5/36 ]\n"
1633  << "b = [ 5/18 4/9 5/18 ]'" << std::endl;
1634  typedef ScalarTraits<Scalar> ST;
1635  int myNumStages = 3;
1636  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1637  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1638  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1639  Scalar one = ST::one();
1640  Scalar ten = as<Scalar>(10*one);
1641  Scalar fifteen = as<Scalar>(15*one);
1642  Scalar twentyfour = as<Scalar>(24*one);
1643  Scalar thirty = as<Scalar>(30*one);
1644  Scalar sqrt15over10 = as<Scalar>(ST::squareroot(fifteen)/ten);
1645  Scalar sqrt15over15 = as<Scalar>(ST::squareroot(fifteen)/fifteen);
1646  Scalar sqrt15over24 = as<Scalar>(ST::squareroot(fifteen)/twentyfour);
1647  Scalar sqrt15over30 = as<Scalar>(ST::squareroot(fifteen)/thirty);
1648 
1649  myA(0,0) = as<Scalar>(5*one/(36*one));
1650  myA(0,1) = as<Scalar>(2*one/(9*one))-sqrt15over15;
1651  myA(0,2) = as<Scalar>(5*one/(36*one))-sqrt15over30;
1652  myA(1,0) = as<Scalar>(5*one/(36*one))+sqrt15over24;
1653  myA(1,1) = as<Scalar>(2*one/(9*one));
1654  myA(1,2) = as<Scalar>(5*one/(36*one))-sqrt15over24;
1655  myA(2,0) = as<Scalar>(5*one/(36*one))+sqrt15over30;
1656  myA(2,1) = as<Scalar>(2*one/(9*one))+sqrt15over15;
1657  myA(2,2) = as<Scalar>(5*one/(36*one));
1658  myb(0) = as<Scalar>(5*one/(18*one));
1659  myb(1) = as<Scalar>(4*one/(9*one));
1660  myb(2) = as<Scalar>(5*one/(18*one));
1661  myc(0) = as<Scalar>(one/(2*one))-sqrt15over10;
1662  myc(1) = as<Scalar>(one/(2*one));
1663  myc(2) = as<Scalar>(one/(2*one))+sqrt15over10;
1664  this->setMyDescription(myDescription.str());
1665  this->setMy_A(myA);
1666  this->setMy_b(myb);
1667  this->setMy_c(myc);
1668  this->setMy_order(6);
1669  }
1670 };
1671 
1672 
1673 template<class Scalar>
1674 class Implicit1Stage1stOrderRadauA_RKBT :
1675  virtual public RKButcherTableauDefaultBase<Scalar>
1676 {
1677  public:
1678  Implicit1Stage1stOrderRadauA_RKBT()
1679  {
1680 
1681  std::ostringstream myDescription;
1682  myDescription << Implicit1Stage1stOrderRadauA_name() << "\n"
1683  << "A-stable\n"
1684  << "Solving Ordinary Differential Equations II:\n"
1685  << "Stiff and Differential-Algebraic Problems,\n"
1686  << "2nd Revised Edition\n"
1687  << "E. Hairer and G. Wanner\n"
1688  << "Table 5.3, pg 73\n"
1689  << "c = [ 0 ]'\n"
1690  << "A = [ 1 ]\n"
1691  << "b = [ 1 ]'" << std::endl;
1692  typedef ScalarTraits<Scalar> ST;
1693  int myNumStages = 1;
1694  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1695  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1696  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1697  Scalar one = ST::one();
1698  Scalar zero = ST::zero();
1699  myA(0,0) = one;
1700  myb(0) = one;
1701  myc(0) = zero;
1702  this->setMyDescription(myDescription.str());
1703  this->setMy_A(myA);
1704  this->setMy_b(myb);
1705  this->setMy_c(myc);
1706  this->setMy_order(1);
1707  }
1708 };
1709 
1710 
1711 template<class Scalar>
1712 class Implicit2Stage3rdOrderRadauA_RKBT :
1713  virtual public RKButcherTableauDefaultBase<Scalar>
1714 {
1715  public:
1716  Implicit2Stage3rdOrderRadauA_RKBT()
1717  {
1718 
1719  std::ostringstream myDescription;
1720  myDescription << Implicit2Stage3rdOrderRadauA_name() << "\n"
1721  << "A-stable\n"
1722  << "Solving Ordinary Differential Equations II:\n"
1723  << "Stiff and Differential-Algebraic Problems,\n"
1724  << "2nd Revised Edition\n"
1725  << "E. Hairer and G. Wanner\n"
1726  << "Table 5.3, pg 73\n"
1727  << "c = [ 0 2/3 ]'\n"
1728  << "A = [ 1/4 -1/4 ]\n"
1729  << " [ 1/4 5/12 ]\n"
1730  << "b = [ 1/4 3/4 ]'" << std::endl;
1731  typedef ScalarTraits<Scalar> ST;
1732  int myNumStages = 2;
1733  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1734  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1735  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1736  Scalar zero = ST::zero();
1737  Scalar one = ST::one();
1738  myA(0,0) = as<Scalar>(one/(4*one));
1739  myA(0,1) = as<Scalar>(-one/(4*one));
1740  myA(1,0) = as<Scalar>(one/(4*one));
1741  myA(1,1) = as<Scalar>(5*one/(12*one));
1742  myb(0) = as<Scalar>(one/(4*one));
1743  myb(1) = as<Scalar>(3*one/(4*one));
1744  myc(0) = zero;
1745  myc(1) = as<Scalar>(2*one/(3*one));
1746  this->setMyDescription(myDescription.str());
1747  this->setMy_A(myA);
1748  this->setMy_b(myb);
1749  this->setMy_c(myc);
1750  this->setMy_order(3);
1751  }
1752 };
1753 
1754 
1755 template<class Scalar>
1756 class Implicit3Stage5thOrderRadauA_RKBT :
1757  virtual public RKButcherTableauDefaultBase<Scalar>
1758 {
1759  public:
1760  Implicit3Stage5thOrderRadauA_RKBT()
1761  {
1762 
1763  std::ostringstream myDescription;
1764  myDescription << Implicit3Stage5thOrderRadauA_name() << "\n"
1765  << "A-stable\n"
1766  << "Solving Ordinary Differential Equations II:\n"
1767  << "Stiff and Differential-Algebraic Problems,\n"
1768  << "2nd Revised Edition\n"
1769  << "E. Hairer and G. Wanner\n"
1770  << "Table 5.4, pg 73\n"
1771  << "c = [ 0 (6-sqrt(6))/10 (6+sqrt(6))/10 ]'\n"
1772  << "A = [ 1/9 (-1-sqrt(6))/18 (-1+sqrt(6))/18 ]\n"
1773  << " [ 1/9 (88+7*sqrt(6))/360 (88-43*sqrt(6))/360 ]\n"
1774  << " [ 1/9 (88+43*sqrt(6))/360 (88-7*sqrt(6))/360 ]\n"
1775  << "b = [ 1/9 (16+sqrt(6))/36 (16-sqrt(6))/36 ]'" << std::endl;
1776  typedef ScalarTraits<Scalar> ST;
1777  int myNumStages = 3;
1778  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1779  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1780  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1781  Scalar zero = ST::zero();
1782  Scalar one = ST::one();
1783  myA(0,0) = as<Scalar>(one/(9*one));
1784  myA(0,1) = as<Scalar>( (-one-ST::squareroot(6*one))/(18*one) );
1785  myA(0,2) = as<Scalar>( (-one+ST::squareroot(6*one))/(18*one) );
1786  myA(1,0) = as<Scalar>(one/(9*one));
1787  myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
1788  myA(1,2) = as<Scalar>( (88*one-43*one*ST::squareroot(6*one))/(360*one) );
1789  myA(2,0) = as<Scalar>(one/(9*one));
1790  myA(2,1) = as<Scalar>( (88*one+43*one*ST::squareroot(6*one))/(360*one) );
1791  myA(2,2) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
1792  myb(0) = as<Scalar>(one/(9*one));
1793  myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1794  myb(2) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1795  myc(0) = zero;
1796  myc(1) = as<Scalar>( (6*one-ST::squareroot(6*one))/(10*one) );
1797  myc(2) = as<Scalar>( (6*one+ST::squareroot(6*one))/(10*one) );
1798  this->setMyDescription(myDescription.str());
1799  this->setMy_A(myA);
1800  this->setMy_b(myb);
1801  this->setMy_c(myc);
1802  this->setMy_order(5);
1803  }
1804 };
1805 
1806 
1807 template<class Scalar>
1808 class Implicit1Stage1stOrderRadauB_RKBT :
1809  virtual public RKButcherTableauDefaultBase<Scalar>
1810 {
1811  public:
1812  Implicit1Stage1stOrderRadauB_RKBT()
1813  {
1814 
1815  std::ostringstream myDescription;
1816  myDescription << Implicit1Stage1stOrderRadauB_name() << "\n"
1817  << "A-stable\n"
1818  << "Solving Ordinary Differential Equations II:\n"
1819  << "Stiff and Differential-Algebraic Problems,\n"
1820  << "2nd Revised Edition\n"
1821  << "E. Hairer and G. Wanner\n"
1822  << "Table 5.5, pg 74\n"
1823  << "c = [ 1 ]'\n"
1824  << "A = [ 1 ]\n"
1825  << "b = [ 1 ]'" << std::endl;
1826  typedef ScalarTraits<Scalar> ST;
1827  int myNumStages = 1;
1828  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1829  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1830  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1831  Scalar one = ST::one();
1832  myA(0,0) = one;
1833  myb(0) = one;
1834  myc(0) = one;
1835  this->setMyDescription(myDescription.str());
1836  this->setMy_A(myA);
1837  this->setMy_b(myb);
1838  this->setMy_c(myc);
1839  this->setMy_order(1);
1840  }
1841 };
1842 
1843 
1844 template<class Scalar>
1845 class Implicit2Stage3rdOrderRadauB_RKBT :
1846  virtual public RKButcherTableauDefaultBase<Scalar>
1847 {
1848  public:
1849  Implicit2Stage3rdOrderRadauB_RKBT()
1850  {
1851 
1852  std::ostringstream myDescription;
1853  myDescription << Implicit2Stage3rdOrderRadauB_name() << "\n"
1854  << "A-stable\n"
1855  << "Solving Ordinary Differential Equations II:\n"
1856  << "Stiff and Differential-Algebraic Problems,\n"
1857  << "2nd Revised Edition\n"
1858  << "E. Hairer and G. Wanner\n"
1859  << "Table 5.5, pg 74\n"
1860  << "c = [ 1/3 1 ]'\n"
1861  << "A = [ 5/12 -1/12 ]\n"
1862  << " [ 3/4 1/4 ]\n"
1863  << "b = [ 3/4 1/4 ]'" << std::endl;
1864  typedef ScalarTraits<Scalar> ST;
1865  int myNumStages = 2;
1866  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1867  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1868  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1869  Scalar one = ST::one();
1870  myA(0,0) = as<Scalar>( 5*one/(12*one) );
1871  myA(0,1) = as<Scalar>( -one/(12*one) );
1872  myA(1,0) = as<Scalar>( 3*one/(4*one) );
1873  myA(1,1) = as<Scalar>( one/(4*one) );
1874  myb(0) = as<Scalar>( 3*one/(4*one) );
1875  myb(1) = as<Scalar>( one/(4*one) );
1876  myc(0) = as<Scalar>( one/(3*one) );
1877  myc(1) = one;
1878  this->setMyDescription(myDescription.str());
1879  this->setMy_A(myA);
1880  this->setMy_b(myb);
1881  this->setMy_c(myc);
1882  this->setMy_order(3);
1883  }
1884 };
1885 
1886 
1887 template<class Scalar>
1888 class Implicit3Stage5thOrderRadauB_RKBT :
1889  virtual public RKButcherTableauDefaultBase<Scalar>
1890 {
1891  public:
1892  Implicit3Stage5thOrderRadauB_RKBT()
1893  {
1894 
1895  std::ostringstream myDescription;
1896  myDescription << Implicit3Stage5thOrderRadauB_name() << "\n"
1897  << "A-stable\n"
1898  << "Solving Ordinary Differential Equations II:\n"
1899  << "Stiff and Differential-Algebraic Problems,\n"
1900  << "2nd Revised Edition\n"
1901  << "E. Hairer and G. Wanner\n"
1902  << "Table 5.6, pg 74\n"
1903  << "c = [ (4-sqrt(6))/10 (4+sqrt(6))/10 1 ]'\n"
1904  << "A = [ A1 A2 A3 ]\n"
1905  << " A1 = [ (88-7*sqrt(6))/360 ]\n"
1906  << " [ (296+169*sqrt(6))/1800 ]\n"
1907  << " [ (16-sqrt(6))/36 ]\n"
1908  << " A2 = [ (296-169*sqrt(6))/1800 ]\n"
1909  << " [ (88+7*sqrt(6))/360 ]\n"
1910  << " [ (16+sqrt(6))/36 ]\n"
1911  << " A3 = [ (-2+3*sqrt(6))/225 ]\n"
1912  << " [ (-2-3*sqrt(6))/225 ]\n"
1913  << " [ 1/9 ]\n"
1914  << "b = [ (16-sqrt(6))/36 (16+sqrt(6))/36 1/9 ]'"
1915  << std::endl;
1916  typedef ScalarTraits<Scalar> ST;
1917  int myNumStages = 3;
1918  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1919  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1920  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1921  Scalar one = ST::one();
1922  myA(0,0) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
1923  myA(0,1) = as<Scalar>( (296*one-169*one*ST::squareroot(6*one))/(1800*one) );
1924  myA(0,2) = as<Scalar>( (-2*one+3*one*ST::squareroot(6*one))/(225*one) );
1925  myA(1,0) = as<Scalar>( (296*one+169*one*ST::squareroot(6*one))/(1800*one) );
1926  myA(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
1927  myA(1,2) = as<Scalar>( (-2*one-3*one*ST::squareroot(6*one))/(225*one) );
1928  myA(2,0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1929  myA(2,1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1930  myA(2,2) = as<Scalar>( one/(9*one) );
1931  myb(0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
1932  myb(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
1933  myb(2) = as<Scalar>( one/(9*one) );
1934  myc(0) = as<Scalar>( (4*one-ST::squareroot(6*one))/(10*one) );
1935  myc(1) = as<Scalar>( (4*one+ST::squareroot(6*one))/(10*one) );
1936  myc(2) = one;
1937  this->setMyDescription(myDescription.str());
1938  this->setMy_A(myA);
1939  this->setMy_b(myb);
1940  this->setMy_c(myc);
1941  this->setMy_order(5);
1942  }
1943 };
1944 
1945 
1946 template<class Scalar>
1947 class Implicit2Stage2ndOrderLobattoA_RKBT :
1948  virtual public RKButcherTableauDefaultBase<Scalar>
1949 {
1950  public:
1951  Implicit2Stage2ndOrderLobattoA_RKBT()
1952  {
1953 
1954  std::ostringstream myDescription;
1955  myDescription << Implicit2Stage2ndOrderLobattoA_name() << "\n"
1956  << "A-stable\n"
1957  << "Solving Ordinary Differential Equations II:\n"
1958  << "Stiff and Differential-Algebraic Problems,\n"
1959  << "2nd Revised Edition\n"
1960  << "E. Hairer and G. Wanner\n"
1961  << "Table 5.7, pg 75\n"
1962  << "c = [ 0 1 ]'\n"
1963  << "A = [ 0 0 ]\n"
1964  << " [ 1/2 1/2 ]\n"
1965  << "b = [ 1/2 1/2 ]'" << std::endl;
1966  typedef ScalarTraits<Scalar> ST;
1967  int myNumStages = 2;
1968  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
1969  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
1970  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
1971  Scalar zero = ST::zero();
1972  Scalar one = ST::one();
1973  myA(0,0) = zero;
1974  myA(0,1) = zero;
1975  myA(1,0) = as<Scalar>( one/(2*one) );
1976  myA(1,1) = as<Scalar>( one/(2*one) );
1977  myb(0) = as<Scalar>( one/(2*one) );
1978  myb(1) = as<Scalar>( one/(2*one) );
1979  myc(0) = zero;
1980  myc(1) = one;
1981  this->setMyDescription(myDescription.str());
1982  this->setMy_A(myA);
1983  this->setMy_b(myb);
1984  this->setMy_c(myc);
1985  this->setMy_order(2);
1986  }
1987 };
1988 
1989 
1990 template<class Scalar>
1991 class Implicit3Stage4thOrderLobattoA_RKBT :
1992  virtual public RKButcherTableauDefaultBase<Scalar>
1993 {
1994  public:
1995  Implicit3Stage4thOrderLobattoA_RKBT()
1996  {
1997 
1998  std::ostringstream myDescription;
1999  myDescription << Implicit3Stage4thOrderLobattoA_name() << "\n"
2000  << "A-stable\n"
2001  << "Solving Ordinary Differential Equations II:\n"
2002  << "Stiff and Differential-Algebraic Problems,\n"
2003  << "2nd Revised Edition\n"
2004  << "E. Hairer and G. Wanner\n"
2005  << "Table 5.7, pg 75\n"
2006  << "c = [ 0 1/2 1 ]'\n"
2007  << "A = [ 0 0 0 ]\n"
2008  << " [ 5/24 1/3 -1/24 ]\n"
2009  << " [ 1/6 2/3 1/6 ]\n"
2010  << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
2011  typedef ScalarTraits<Scalar> ST;
2012  int myNumStages = 3;
2013  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2014  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2015  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2016  Scalar zero = ST::zero();
2017  Scalar one = ST::one();
2018  myA(0,0) = zero;
2019  myA(0,1) = zero;
2020  myA(0,2) = zero;
2021  myA(1,0) = as<Scalar>( (5*one)/(24*one) );
2022  myA(1,1) = as<Scalar>( (one)/(3*one) );
2023  myA(1,2) = as<Scalar>( (-one)/(24*one) );
2024  myA(2,0) = as<Scalar>( (one)/(6*one) );
2025  myA(2,1) = as<Scalar>( (2*one)/(3*one) );
2026  myA(2,2) = as<Scalar>( (1*one)/(6*one) );
2027  myb(0) = as<Scalar>( (one)/(6*one) );
2028  myb(1) = as<Scalar>( (2*one)/(3*one) );
2029  myb(2) = as<Scalar>( (1*one)/(6*one) );
2030  myc(0) = zero;
2031  myc(1) = as<Scalar>( one/(2*one) );
2032  myc(2) = one;
2033  this->setMyDescription(myDescription.str());
2034  this->setMy_A(myA);
2035  this->setMy_b(myb);
2036  this->setMy_c(myc);
2037  this->setMy_order(4);
2038  }
2039 };
2040 
2041 
2042 template<class Scalar>
2043 class Implicit4Stage6thOrderLobattoA_RKBT :
2044  virtual public RKButcherTableauDefaultBase<Scalar>
2045 {
2046  public:
2047  Implicit4Stage6thOrderLobattoA_RKBT()
2048  {
2049 
2050  using Teuchos::as;
2051  typedef Teuchos::ScalarTraits<Scalar> ST;
2052 
2053  std::ostringstream myDescription;
2054  myDescription << Implicit4Stage6thOrderLobattoA_name() << "\n"
2055  << "A-stable\n"
2056  << "Solving Ordinary Differential Equations II:\n"
2057  << "Stiff and Differential-Algebraic Problems,\n"
2058  << "2nd Revised Edition\n"
2059  << "E. Hairer and G. Wanner\n"
2060  << "Table 5.8, pg 75\n"
2061  << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
2062  << "A = [ A1 A2 A3 A4 ]\n"
2063  << " A1 = [ 0 ]\n"
2064  << " [ (11+sqrt(5)/120 ]\n"
2065  << " [ (11-sqrt(5)/120 ]\n"
2066  << " [ 1/12 ]\n"
2067  << " A2 = [ 0 ]\n"
2068  << " [ (25-sqrt(5))/120 ]\n"
2069  << " [ (25+13*sqrt(5))/120 ]\n"
2070  << " [ 5/12 ]\n"
2071  << " A3 = [ 0 ]\n"
2072  << " [ (25-13*sqrt(5))/120 ]\n"
2073  << " [ (25+sqrt(5))/120 ]\n"
2074  << " [ 5/12 ]\n"
2075  << " A4 = [ 0 ]\n"
2076  << " [ (-1+sqrt(5))/120 ]\n"
2077  << " [ (-1-sqrt(5))/120 ]\n"
2078  << " [ 1/12 ]\n"
2079  << "b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
2080  typedef ScalarTraits<Scalar> ST;
2081  int myNumStages = 4;
2082  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2083  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2084  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2085  Scalar zero = ST::zero();
2086  Scalar one = ST::one();
2087  myA(0,0) = zero;
2088  myA(0,1) = zero;
2089  myA(0,2) = zero;
2090  myA(0,3) = zero;
2091  myA(1,0) = as<Scalar>( (11*one+ST::squareroot(5*one))/(120*one) );
2092  myA(1,1) = as<Scalar>( (25*one-ST::squareroot(5*one))/(120*one) );
2093  myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5*one))/(120*one) );
2094  myA(1,3) = as<Scalar>( (-one+ST::squareroot(5*one))/(120*one) );
2095  myA(2,0) = as<Scalar>( (11*one-ST::squareroot(5*one))/(120*one) );
2096  myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5*one))/(120*one) );
2097  myA(2,2) = as<Scalar>( (25*one+ST::squareroot(5*one))/(120*one) );
2098  myA(2,3) = as<Scalar>( (-one-ST::squareroot(5*one))/(120*one) );
2099  myA(3,0) = as<Scalar>( one/(12*one) );
2100  myA(3,1) = as<Scalar>( 5*one/(12*one) );
2101  myA(3,2) = as<Scalar>( 5*one/(12*one) );
2102  myA(3,3) = as<Scalar>( one/(12*one) );
2103  myb(0) = as<Scalar>( one/(12*one) );
2104  myb(1) = as<Scalar>( 5*one/(12*one) );
2105  myb(2) = as<Scalar>( 5*one/(12*one) );
2106  myb(3) = as<Scalar>( one/(12*one) );
2107  myc(0) = zero;
2108  myc(1) = as<Scalar>( (5*one-ST::squareroot(5))/(10*one) );
2109  myc(2) = as<Scalar>( (5*one+ST::squareroot(5))/(10*one) );
2110  myc(3) = one;
2111  this->setMyDescription(myDescription.str());
2112  this->setMy_A(myA);
2113  this->setMy_b(myb);
2114  this->setMy_c(myc);
2115  this->setMy_order(6);
2116  }
2117 };
2118 
2119 
2120 template<class Scalar>
2121 class Implicit2Stage2ndOrderLobattoB_RKBT :
2122  virtual public RKButcherTableauDefaultBase<Scalar>
2123 {
2124  public:
2125  Implicit2Stage2ndOrderLobattoB_RKBT()
2126  {
2127 
2128  std::ostringstream myDescription;
2129  myDescription << Implicit2Stage2ndOrderLobattoB_name() << "\n"
2130  << "A-stable\n"
2131  << "Solving Ordinary Differential Equations II:\n"
2132  << "Stiff and Differential-Algebraic Problems,\n"
2133  << "2nd Revised Edition\n"
2134  << "E. Hairer and G. Wanner\n"
2135  << "Table 5.9, pg 76\n"
2136  << "c = [ 0 1 ]'\n"
2137  << "A = [ 1/2 0 ]\n"
2138  << " [ 1/2 0 ]\n"
2139  << "b = [ 1/2 1/2 ]'" << std::endl;
2140  typedef ScalarTraits<Scalar> ST;
2141  int myNumStages = 2;
2142  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2143  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2144  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2145  Scalar zero = ST::zero();
2146  Scalar one = ST::one();
2147  myA(0,0) = as<Scalar>( one/(2*one) );
2148  myA(0,1) = zero;
2149  myA(1,0) = as<Scalar>( one/(2*one) );
2150  myA(1,1) = zero;
2151  myb(0) = as<Scalar>( one/(2*one) );
2152  myb(1) = as<Scalar>( one/(2*one) );
2153  myc(0) = zero;
2154  myc(1) = one;
2155  this->setMyDescription(myDescription.str());
2156  this->setMy_A(myA);
2157  this->setMy_b(myb);
2158  this->setMy_c(myc);
2159  this->setMy_order(2);
2160  }
2161 };
2162 
2163 
2164 template<class Scalar>
2165 class Implicit3Stage4thOrderLobattoB_RKBT :
2166  virtual public RKButcherTableauDefaultBase<Scalar>
2167 {
2168  public:
2169  Implicit3Stage4thOrderLobattoB_RKBT()
2170  {
2171 
2172  std::ostringstream myDescription;
2173  myDescription << Implicit3Stage4thOrderLobattoB_name() << "\n"
2174  << "A-stable\n"
2175  << "Solving Ordinary Differential Equations II:\n"
2176  << "Stiff and Differential-Algebraic Problems,\n"
2177  << "2nd Revised Edition\n"
2178  << "E. Hairer and G. Wanner\n"
2179  << "Table 5.9, pg 76\n"
2180  << "c = [ 0 1/2 1 ]'\n"
2181  << "A = [ 1/6 -1/6 0 ]\n"
2182  << " [ 1/6 1/3 0 ]\n"
2183  << " [ 1/6 5/6 0 ]\n"
2184  << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
2185  typedef ScalarTraits<Scalar> ST;
2186  int myNumStages = 3;
2187  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2188  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2189  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2190  Scalar zero = ST::zero();
2191  Scalar one = ST::one();
2192  myA(0,0) = as<Scalar>( one/(6*one) );
2193  myA(0,1) = as<Scalar>( -one/(6*one) );
2194  myA(0,2) = zero;
2195  myA(1,0) = as<Scalar>( one/(6*one) );
2196  myA(1,1) = as<Scalar>( one/(3*one) );
2197  myA(1,2) = zero;
2198  myA(2,0) = as<Scalar>( one/(6*one) );
2199  myA(2,1) = as<Scalar>( 5*one/(6*one) );
2200  myA(2,2) = zero;
2201  myb(0) = as<Scalar>( one/(6*one) );
2202  myb(1) = as<Scalar>( 2*one/(3*one) );
2203  myb(2) = as<Scalar>( one/(6*one) );
2204  myc(0) = zero;
2205  myc(1) = as<Scalar>( one/(2*one) );
2206  myc(2) = one;
2207  this->setMyDescription(myDescription.str());
2208  this->setMy_A(myA);
2209  this->setMy_b(myb);
2210  this->setMy_c(myc);
2211  this->setMy_order(4);
2212  }
2213 };
2214 
2215 
2216 template<class Scalar>
2217 class Implicit4Stage6thOrderLobattoB_RKBT :
2218  virtual public RKButcherTableauDefaultBase<Scalar>
2219 {
2220  public:
2221  Implicit4Stage6thOrderLobattoB_RKBT()
2222  {
2223 
2224  std::ostringstream myDescription;
2225  myDescription << Implicit4Stage6thOrderLobattoB_name() << "\n"
2226  << "A-stable\n"
2227  << "Solving Ordinary Differential Equations II:\n"
2228  << "Stiff and Differential-Algebraic Problems,\n"
2229  << "2nd Revised Edition\n"
2230  << "E. Hairer and G. Wanner\n"
2231  << "Table 5.10, pg 76\n"
2232  << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
2233  << "A = [ 1/12 (-1-sqrt(5))/24 (-1+sqrt(5))/24 0 ]\n"
2234  << " [ 1/12 (25+sqrt(5))/120 (25-13*sqrt(5))/120 0 ]\n"
2235  << " [ 1/12 (25+13*sqrt(5))/120 (25-sqrt(5))/120 0 ]\n"
2236  << " [ 1/12 (11-sqrt(5))/24 (11+sqrt(5))/24 0 ]\n"
2237  << "b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
2238  typedef ScalarTraits<Scalar> ST;
2239  int myNumStages = 4;
2240  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2241  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2242  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2243  Scalar zero = ST::zero();
2244  Scalar one = ST::one();
2245  myA(0,0) = as<Scalar>( one/(12*one) );
2246  myA(0,1) = as<Scalar>( (-one-ST::squareroot(5))/(24*one) );
2247  myA(0,2) = as<Scalar>( (-one+ST::squareroot(5))/(24*one) );
2248  myA(0,3) = zero;
2249  myA(1,0) = as<Scalar>( one/(12*one) );
2250  myA(1,1) = as<Scalar>( (25*one+ST::squareroot(5))/(120*one) );
2251  myA(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5))/(120*one) );
2252  myA(1,3) = zero;
2253  myA(2,0) = as<Scalar>( one/(12*one) );
2254  myA(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5))/(120*one) );
2255  myA(2,2) = as<Scalar>( (25*one-ST::squareroot(5))/(120*one) );
2256  myA(2,3) = zero;
2257  myA(3,0) = as<Scalar>( one/(12*one) );
2258  myA(3,1) = as<Scalar>( (11*one-ST::squareroot(5*one))/(24*one) );
2259  myA(3,2) = as<Scalar>( (11*one+ST::squareroot(5*one))/(24*one) );
2260  myA(3,3) = zero;
2261  myb(0) = as<Scalar>( one/(12*one) );
2262  myb(1) = as<Scalar>( 5*one/(12*one) );
2263  myb(2) = as<Scalar>( 5*one/(12*one) );
2264  myb(3) = as<Scalar>( one/(12*one) );
2265  myc(0) = zero;
2266  myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
2267  myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
2268  myc(3) = one;
2269  this->setMyDescription(myDescription.str());
2270  this->setMy_A(myA);
2271  this->setMy_b(myb);
2272  this->setMy_c(myc);
2273  this->setMy_order(6);
2274  }
2275 };
2276 
2277 
2278 template<class Scalar>
2279 class Implicit2Stage2ndOrderLobattoC_RKBT :
2280  virtual public RKButcherTableauDefaultBase<Scalar>
2281 {
2282  public:
2283  Implicit2Stage2ndOrderLobattoC_RKBT()
2284  {
2285 
2286  std::ostringstream myDescription;
2287  myDescription << Implicit2Stage2ndOrderLobattoC_name() << "\n"
2288  << "A-stable\n"
2289  << "Solving Ordinary Differential Equations II:\n"
2290  << "Stiff and Differential-Algebraic Problems,\n"
2291  << "2nd Revised Edition\n"
2292  << "E. Hairer and G. Wanner\n"
2293  << "Table 5.11, pg 76\n"
2294  << "c = [ 0 1 ]'\n"
2295  << "A = [ 1/2 -1/2 ]\n"
2296  << " [ 1/2 1/2 ]\n"
2297  << "b = [ 1/2 1/2 ]'" << std::endl;
2298  typedef ScalarTraits<Scalar> ST;
2299  int myNumStages = 2;
2300  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2301  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2302  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2303  Scalar zero = ST::zero();
2304  Scalar one = ST::one();
2305  myA(0,0) = as<Scalar>( one/(2*one) );
2306  myA(0,1) = as<Scalar>( -one/(2*one) );
2307  myA(1,0) = as<Scalar>( one/(2*one) );
2308  myA(1,1) = as<Scalar>( one/(2*one) );
2309  myb(0) = as<Scalar>( one/(2*one) );
2310  myb(1) = as<Scalar>( one/(2*one) );
2311  myc(0) = zero;
2312  myc(1) = one;
2313  this->setMyDescription(myDescription.str());
2314  this->setMy_A(myA);
2315  this->setMy_b(myb);
2316  this->setMy_c(myc);
2317  this->setMy_order(2);
2318  }
2319 };
2320 
2321 
2322 template<class Scalar>
2323 class Implicit3Stage4thOrderLobattoC_RKBT :
2324  virtual public RKButcherTableauDefaultBase<Scalar>
2325 {
2326  public:
2327  Implicit3Stage4thOrderLobattoC_RKBT()
2328  {
2329 
2330  std::ostringstream myDescription;
2331  myDescription << Implicit3Stage4thOrderLobattoC_name() << "\n"
2332  << "A-stable\n"
2333  << "Solving Ordinary Differential Equations II:\n"
2334  << "Stiff and Differential-Algebraic Problems,\n"
2335  << "2nd Revised Edition\n"
2336  << "E. Hairer and G. Wanner\n"
2337  << "Table 5.11, pg 76\n"
2338  << "c = [ 0 1/2 1 ]'\n"
2339  << "A = [ 1/6 -1/3 1/6 ]\n"
2340  << " [ 1/6 5/12 -1/12 ]\n"
2341  << " [ 1/6 2/3 1/6 ]\n"
2342  << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
2343  typedef ScalarTraits<Scalar> ST;
2344  int myNumStages = 3;
2345  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2346  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2347  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2348  Scalar zero = ST::zero();
2349  Scalar one = ST::one();
2350  myA(0,0) = as<Scalar>( one/(6*one) );
2351  myA(0,1) = as<Scalar>( -one/(3*one) );
2352  myA(0,2) = as<Scalar>( one/(6*one) );
2353  myA(1,0) = as<Scalar>( one/(6*one) );
2354  myA(1,1) = as<Scalar>( 5*one/(12*one) );
2355  myA(1,2) = as<Scalar>( -one/(12*one) );
2356  myA(2,0) = as<Scalar>( one/(6*one) );
2357  myA(2,1) = as<Scalar>( 2*one/(3*one) );
2358  myA(2,2) = as<Scalar>( one/(6*one) );
2359  myb(0) = as<Scalar>( one/(6*one) );
2360  myb(1) = as<Scalar>( 2*one/(3*one) );
2361  myb(2) = as<Scalar>( one/(6*one) );
2362  myc(0) = zero;
2363  myc(1) = as<Scalar>( one/(2*one) );
2364  myc(2) = one;
2365  this->setMyDescription(myDescription.str());
2366  this->setMy_A(myA);
2367  this->setMy_b(myb);
2368  this->setMy_c(myc);
2369  this->setMy_order(4);
2370  }
2371 };
2372 
2373 
2374 template<class Scalar>
2375 class Implicit4Stage6thOrderLobattoC_RKBT :
2376  virtual public RKButcherTableauDefaultBase<Scalar>
2377 {
2378  public:
2379  Implicit4Stage6thOrderLobattoC_RKBT()
2380  {
2381 
2382  std::ostringstream myDescription;
2383  myDescription << Implicit4Stage6thOrderLobattoC_name() << "\n"
2384  << "A-stable\n"
2385  << "Solving Ordinary Differential Equations II:\n"
2386  << "Stiff and Differential-Algebraic Problems,\n"
2387  << "2nd Revised Edition\n"
2388  << "E. Hairer and G. Wanner\n"
2389  << "Table 5.12, pg 76\n"
2390  << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
2391  << "A = [ 1/12 -sqrt(5)/12 sqrt(5)/12 -1/12 ]\n"
2392  << " [ 1/12 1/4 (10-7*sqrt(5))/60 sqrt(5)/60 ]\n"
2393  << " [ 1/12 (10+7*sqrt(5))/60 1/4 -sqrt(5)/60 ]\n"
2394  << " [ 1/12 5/12 5/12 1/12 ]\n"
2395  << "b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
2396  typedef ScalarTraits<Scalar> ST;
2397  int myNumStages = 4;
2398  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2399  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2400  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2401  Scalar zero = ST::zero();
2402  Scalar one = ST::one();
2403  myA(0,0) = as<Scalar>( one/(12*one) );
2404  myA(0,1) = as<Scalar>( -ST::squareroot(5*one)/(12*one) );
2405  myA(0,2) = as<Scalar>( ST::squareroot(5*one)/(12*one) );
2406  myA(0,3) = as<Scalar>( -one/(12*one) );
2407  myA(1,0) = as<Scalar>( one/(12*one) );
2408  myA(1,1) = as<Scalar>( one/(4*one) );
2409  myA(1,2) = as<Scalar>( (10*one-7*one*ST::squareroot(5*one))/(60*one) );
2410  myA(1,3) = as<Scalar>( ST::squareroot(5*one)/(60*one) );
2411  myA(2,0) = as<Scalar>( one/(12*one) );
2412  myA(2,1) = as<Scalar>( (10*one+7*one*ST::squareroot(5*one))/(60*one) );
2413  myA(2,2) = as<Scalar>( one/(4*one) );
2414  myA(2,3) = as<Scalar>( -ST::squareroot(5*one)/(60*one) );
2415  myA(3,0) = as<Scalar>( one/(12*one) );
2416  myA(3,1) = as<Scalar>( 5*one/(12*one) );
2417  myA(3,2) = as<Scalar>( 5*one/(12*one) );
2418  myA(3,3) = as<Scalar>( one/(12*one) );
2419  myb(0) = as<Scalar>( one/(12*one) );
2420  myb(1) = as<Scalar>( 5*one/(12*one) );
2421  myb(2) = as<Scalar>( 5*one/(12*one) );
2422  myb(3) = as<Scalar>( one/(12*one) );
2423  myc(0) = zero;
2424  myc(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
2425  myc(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
2426  myc(3) = one;
2427  this->setMyDescription(myDescription.str());
2428  this->setMy_A(myA);
2429  this->setMy_b(myb);
2430  this->setMy_c(myc);
2431  this->setMy_order(6);
2432  }
2433 };
2434 
2435 
2436 
2437 template<class Scalar>
2438 class SDIRK5Stage5thOrder_RKBT :
2439  virtual public RKButcherTableauDefaultBase<Scalar>
2440 {
2441  public:
2442  SDIRK5Stage5thOrder_RKBT()
2443  {
2444 
2445  std::ostringstream myDescription;
2446  myDescription << SDIRK5Stage5thOrder_name() << "\n"
2447  << "A-stable\n"
2448  << "Solving Ordinary Differential Equations II:\n"
2449  << "Stiff and Differential-Algebraic Problems,\n"
2450  << "2nd Revised Edition\n"
2451  << "E. Hairer and G. Wanner\n"
2452  << "pg101 \n"
2453  << "c = [ (6-sqrt(6))/10 ]\n"
2454  << " [ (6+9*sqrt(6))/35 ]\n"
2455  << " [ 1 ]\n"
2456  << " [ (4-sqrt(6))/10 ]\n"
2457  << " [ (4+sqrt(6))/10 ]\n"
2458  << "A = [ A1 A2 A3 A4 A5 ]\n"
2459  << " A1 = [ (6-sqrt(6))/10 ]\n"
2460  << " [ (-6+5*sqrt(6))/14 ]\n"
2461  << " [ (888+607*sqrt(6))/2850 ]\n"
2462  << " [ (3153-3082*sqrt(6))/14250 ]\n"
2463  << " [ (-32583+14638*sqrt(6))/71250 ]\n"
2464  << " A2 = [ 0 ]\n"
2465  << " [ (6-sqrt(6))/10 ]\n"
2466  << " [ (126-161*sqrt(6))/1425 ]\n"
2467  << " [ (3213+1148*sqrt(6))/28500 ]\n"
2468  << " [ (-17199+364*sqrt(6))/142500 ]\n"
2469  << " A3 = [ 0 ]\n"
2470  << " [ 0 ]\n"
2471  << " [ (6-sqrt(6))/10 ]\n"
2472  << " [ (-267+88*sqrt(6))/500 ]\n"
2473  << " [ (1329-544*sqrt(6))/2500 ]\n"
2474  << " A4 = [ 0 ]\n"
2475  << " [ 0 ]\n"
2476  << " [ 0 ]\n"
2477  << " [ (6-sqrt(6))/10 ]\n"
2478  << " [ (-96+131*sqrt(6))/625 ]\n"
2479  << " A5 = [ 0 ]\n"
2480  << " [ 0 ]\n"
2481  << " [ 0 ]\n"
2482  << " [ 0 ]\n"
2483  << " [ (6-sqrt(6))/10 ]\n"
2484  << "b = [ 0 ]\n"
2485  << " [ 0 ]\n"
2486  << " [ 1/9 ]\n"
2487  << " [ (16-sqrt(6))/36 ]\n"
2488  << " [ (16+sqrt(6))/36 ]" << std::endl;
2489  typedef ScalarTraits<Scalar> ST;
2490  int myNumStages = 5;
2491  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2492  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2493  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2494  Scalar zero = ST::zero();
2495  Scalar one = ST::one();
2496  Scalar sqrt6 = ST::squareroot(as<Scalar>(6*one));
2497  Scalar gamma = as<Scalar>( (6*one - sqrt6) / (10*one) ); // diagonal
2498  myA(0,0) = gamma;
2499  myA(0,1) = zero;
2500  myA(0,2) = zero;
2501  myA(0,3) = zero;
2502  myA(0,4) = zero;
2503 
2504  myA(1,0) = as<Scalar>( (-6*one+5*one*sqrt6)/(14*one) );
2505  myA(1,1) = gamma;
2506  myA(1,2) = zero;
2507  myA(1,3) = zero;
2508  myA(1,4) = zero;
2509 
2510  myA(2,0) = as<Scalar>( (888*one+607*one*sqrt6)/(2850*one) );
2511  myA(2,1) = as<Scalar>( (126*one-161*one*sqrt6)/(1425*one) );
2512  myA(2,2) = gamma;
2513  myA(2,3) = zero;
2514  myA(2,4) = zero;
2515 
2516  myA(3,0) = as<Scalar>( (3153*one-3082*one*sqrt6)/(14250*one) );
2517  myA(3,1) = as<Scalar>( (3213*one+1148*one*sqrt6)/(28500*one) );
2518  myA(3,2) = as<Scalar>( (-267*one+88*one*sqrt6)/(500*one) );
2519  myA(3,3) = gamma;
2520  myA(3,4) = zero;
2521 
2522  myA(4,0) = as<Scalar>( (-32583*one+14638*one*sqrt6)/(71250*one) );
2523  myA(4,1) = as<Scalar>( (-17199*one+364*one*sqrt6)/(142500*one) );
2524  myA(4,2) = as<Scalar>( (1329*one-544*one*sqrt6)/(2500*one) );
2525  myA(4,3) = as<Scalar>( (-96*one+131*sqrt6)/(625*one) );
2526  myA(4,4) = gamma;
2527 
2528  myb(0) = zero;
2529  myb(1) = zero;
2530  myb(2) = as<Scalar>( one/(9*one) );
2531  myb(3) = as<Scalar>( (16*one-sqrt6)/(36*one) );
2532  myb(4) = as<Scalar>( (16*one+sqrt6)/(36*one) );
2533 
2534  myc(0) = gamma;
2535  myc(1) = as<Scalar>( (6*one+9*one*sqrt6)/(35*one) );
2536  myc(2) = one;
2537  myc(3) = as<Scalar>( (4*one-sqrt6)/(10*one) );
2538  myc(4) = as<Scalar>( (4*one+sqrt6)/(10*one) );
2539 
2540  this->setMyDescription(myDescription.str());
2541  this->setMy_A(myA);
2542  this->setMy_b(myb);
2543  this->setMy_c(myc);
2544  this->setMy_order(5);
2545  }
2546 };
2547 
2548 
2549 template<class Scalar>
2550 class SDIRK5Stage4thOrder_RKBT :
2551  virtual public RKButcherTableauDefaultBase<Scalar>
2552 {
2553  public:
2554  SDIRK5Stage4thOrder_RKBT()
2555  {
2556 
2557  std::ostringstream myDescription;
2558  myDescription << SDIRK5Stage4thOrder_name() << "\n"
2559  << "L-stable\n"
2560  << "Solving Ordinary Differential Equations II:\n"
2561  << "Stiff and Differential-Algebraic Problems,\n"
2562  << "2nd Revised Edition\n"
2563  << "E. Hairer and G. Wanner\n"
2564  << "pg100 \n"
2565  << "c = [ 1/4 3/4 11/20 1/2 1 ]'\n"
2566  << "A = [ 1/4 ]\n"
2567  << " [ 1/2 1/4 ]\n"
2568  << " [ 17/50 -1/25 1/4 ]\n"
2569  << " [ 371/1360 -137/2720 15/544 1/4 ]\n"
2570  << " [ 25/24 -49/48 125/16 -85/12 1/4 ]\n"
2571  << "b = [ 25/24 -49/48 125/16 -85/12 1/4 ]'\n"
2572  << "b' = [ 59/48 -17/96 225/32 -85/12 0 ]'" << std::endl;
2573  typedef ScalarTraits<Scalar> ST;
2574  int myNumStages = 5;
2575  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2576  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2577  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2578  Scalar zero = ST::zero();
2579  Scalar one = ST::one();
2580  Scalar onequarter = as<Scalar>( one/(4*one) );
2581  myA(0,0) = onequarter;
2582  myA(0,1) = zero;
2583  myA(0,2) = zero;
2584  myA(0,3) = zero;
2585  myA(0,4) = zero;
2586 
2587  myA(1,0) = as<Scalar>( one / (2*one) );
2588  myA(1,1) = onequarter;
2589  myA(1,2) = zero;
2590  myA(1,3) = zero;
2591  myA(1,4) = zero;
2592 
2593  myA(2,0) = as<Scalar>( 17*one/(50*one) );
2594  myA(2,1) = as<Scalar>( -one/(25*one) );
2595  myA(2,2) = onequarter;
2596  myA(2,3) = zero;
2597  myA(2,4) = zero;
2598 
2599  myA(3,0) = as<Scalar>( 371*one/(1360*one) );
2600  myA(3,1) = as<Scalar>( -137*one/(2720*one) );
2601  myA(3,2) = as<Scalar>( 15*one/(544*one) );
2602  myA(3,3) = onequarter;
2603  myA(3,4) = zero;
2604 
2605  myA(4,0) = as<Scalar>( 25*one/(24*one) );
2606  myA(4,1) = as<Scalar>( -49*one/(48*one) );
2607  myA(4,2) = as<Scalar>( 125*one/(16*one) );
2608  myA(4,3) = as<Scalar>( -85*one/(12*one) );
2609  myA(4,4) = onequarter;
2610 
2611  myb(0) = as<Scalar>( 25*one/(24*one) );
2612  myb(1) = as<Scalar>( -49*one/(48*one) );
2613  myb(2) = as<Scalar>( 125*one/(16*one) );
2614  myb(3) = as<Scalar>( -85*one/(12*one) );
2615  myb(4) = onequarter;
2616 
2617  /*
2618  // Alternate version
2619  myb(0) = as<Scalar>( 59*one/(48*one) );
2620  myb(1) = as<Scalar>( -17*one/(96*one) );
2621  myb(2) = as<Scalar>( 225*one/(32*one) );
2622  myb(3) = as<Scalar>( -85*one/(12*one) );
2623  myb(4) = zero;
2624  */
2625  myc(0) = onequarter;
2626  myc(1) = as<Scalar>( 3*one/(4*one) );
2627  myc(2) = as<Scalar>( 11*one/(20*one) );
2628  myc(3) = as<Scalar>( one/(2*one) );
2629  myc(4) = one;
2630 
2631  this->setMyDescription(myDescription.str());
2632  this->setMy_A(myA);
2633  this->setMy_b(myb);
2634  this->setMy_c(myc);
2635  this->setMy_order(4);
2636  }
2637 };
2638 
2639 
2640 template<class Scalar>
2641 class SDIRK3Stage4thOrder_RKBT :
2642  virtual public RKButcherTableauDefaultBase<Scalar>
2643 {
2644  public:
2645  SDIRK3Stage4thOrder_RKBT()
2646  {
2647 
2648  std::ostringstream myDescription;
2649  myDescription << SDIRK3Stage4thOrder_name() << "\n"
2650  << "A-stable\n"
2651  << "Solving Ordinary Differential Equations II:\n"
2652  << "Stiff and Differential-Algebraic Problems,\n"
2653  << "2nd Revised Edition\n"
2654  << "E. Hairer and G. Wanner\n"
2655  << "pg100 \n"
2656  << "gamma = (1/sqrt(3))*cos(pi/18)+1/2\n"
2657  << "delta = 1/(6*(2*gamma-1)^2)\n"
2658  << "c = [ gamma 1/2 1-gamma ]'\n"
2659  << "A = [ gamma ]\n"
2660  << " [ 1/2-gamma gamma ]\n"
2661  << " [ 2*gamma 1-4*gamma gamma ]\n"
2662  << "b = [ delta 1-2*delta delta ]'" << std::endl;
2663  typedef ScalarTraits<Scalar> ST;
2664  int myNumStages = 3;
2665  Teuchos::SerialDenseMatrix<int,Scalar> myA(myNumStages,myNumStages);
2666  Teuchos::SerialDenseVector<int,Scalar> myb(myNumStages);
2667  Teuchos::SerialDenseVector<int,Scalar> myc(myNumStages);
2668  Scalar zero = ST::zero();
2669  Scalar one = ST::one();
2670  Scalar pi = as<Scalar>(4*one)*std::atan(one);
2671  Scalar gamma = as<Scalar>( one/ST::squareroot(3*one)*std::cos(pi/(18*one))+one/(2*one) );
2672  Scalar delta = as<Scalar>( one/(6*one*std::pow(2*gamma-one,2*one)) );
2673  myA(0,0) = gamma;
2674  myA(0,1) = zero;
2675  myA(0,2) = zero;
2676 
2677  myA(1,0) = as<Scalar>( one/(2*one) - gamma );
2678  myA(1,1) = gamma;
2679  myA(1,2) = zero;
2680 
2681  myA(2,0) = as<Scalar>( 2*gamma );
2682  myA(2,1) = as<Scalar>( one - 4*gamma );
2683  myA(2,2) = gamma;
2684 
2685  myb(0) = delta;
2686  myb(1) = as<Scalar>( one-2*delta );
2687  myb(2) = delta;
2688 
2689  myc(0) = gamma;
2690  myc(1) = as<Scalar>( one/(2*one) );
2691  myc(2) = as<Scalar>( one - gamma );
2692 
2693  this->setMyDescription(myDescription.str());
2694  this->setMy_A(myA);
2695  this->setMy_b(myb);
2696  this->setMy_c(myc);
2697  this->setMy_order(4);
2698  }
2699 };
2700 
2701 
2702 } // namespace Rythmos
2703 
2704 
2705 #endif // RYTHMOS_RK_BUTCHER_TABLEAU_HPP