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