Tempus  Version of the Day
Time Integration
Tempus_RKButcherTableau.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ****************************************************************************
3 // Tempus: Copyright (2017) Sandia Corporation
4 //
5 // Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6 // ****************************************************************************
7 // @HEADER
8 
9 #ifndef Tempus_RKButcherTableau_hpp
10 #define Tempus_RKButcherTableau_hpp
11 
12 // disable clang warnings
13 #ifdef __clang__
14 #pragma clang system_header
15 #endif
16 
17 #include "Tempus_String_Utilities.hpp"
18 
19 #include "Teuchos_Assert.hpp"
20 #include "Teuchos_as.hpp"
21 #include "Teuchos_Describable.hpp"
22 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
23 #include "Teuchos_VerboseObject.hpp"
24 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
25 #include "Teuchos_SerialDenseMatrix.hpp"
26 #include "Teuchos_SerialDenseVector.hpp"
27 #include "Thyra_MultiVectorStdOps.hpp"
28 
29 
30 //#include "Thyra_ProductVectorBase.hpp"
31 
32 namespace Tempus {
33 
34 
35 /** \brief Runge-Kutta methods.
36  *
37  * This base class specifies the Butcher tableau which defines the
38  * Runge-Kutta (RK) method. Both explicit and implicit RK methods
39  * can be specied here, and of arbitrary number of stages and orders.
40  * Embedded methods are also supported.
41  *
42  * Since this is a generic RK class, no low-storage methods are
43  * incorporated here, however any RK method with a Butcher tableau
44  * can be created with the base class.
45  *
46  * There are over 40 derived RK methods that have been implemented,
47  * ranging from first order and eight order, and from single stage
48  * to 5 stages.
49  *
50  * This class was taken and modified from Rythmos' RKButcherTableau class.
51  */
52 template<class Scalar>
53 class RKButcherTableau :
54  virtual public Teuchos::Describable,
55  virtual public Teuchos::ParameterListAcceptorDefaultBase,
56  virtual public Teuchos::VerboseObject<RKButcherTableau<Scalar> >
57 {
58  public:
59  /** \brief Return the number of stages */
60  virtual std::size_t numStages() const { return A_.numRows(); }
61  /** \brief Return the matrix coefficients */
62  virtual const Teuchos::SerialDenseMatrix<int,Scalar>& A() const
63  { return A_; }
64  /** \brief Return the vector of quadrature weights */
65  virtual const Teuchos::SerialDenseVector<int,Scalar>& b() const
66  { return b_; }
67  /** \brief Return the vector of quadrature weights for embedded methods */
68  virtual const Teuchos::SerialDenseVector<int,Scalar>& bstar() const
69  { return bstar_ ; }
70  /** \brief Return the vector of stage positions */
71  virtual const Teuchos::SerialDenseVector<int,Scalar>& c() const
72  { return c_; }
73  /** \brief Return the order */
74  virtual int order() const { return order_; }
75  /** \brief Return the minimum order */
76  virtual int orderMin() const { return orderMin_; }
77  /** \brief Return the maximum order */
78  virtual int orderMax() const { return orderMax_; }
79  /** \brief Return true if the RK method is implicit */
80  virtual bool isImplicit() const { return isImplicit_; }
81  /** \brief Return true if the RK method is Diagonally Implicit */
82  virtual bool isDIRK() const { return isDIRK_; }
83  /** \brief Return true if the RK method has embedded capabilities */
84  virtual bool isEmbedded() const { return isEmbedded_; }
85 
86  virtual void initialize(
87  const Teuchos::SerialDenseMatrix<int,Scalar>& A,
88  const Teuchos::SerialDenseVector<int,Scalar>& b,
89  const Teuchos::SerialDenseVector<int,Scalar>& c,
90  const int order,
91  const std::string& longDescription,
92  bool isEmbedded = false,
93  const Teuchos::SerialDenseVector<int,Scalar>&
94  bstar = Teuchos::SerialDenseVector<int,Scalar>())
95  {
96  this->initialize(A,b,c,order,order,order,
97  longDescription,isEmbedded,bstar);
98  }
99 
100  virtual void initialize(
101  const Teuchos::SerialDenseMatrix<int,Scalar>& A,
102  const Teuchos::SerialDenseVector<int,Scalar>& b,
103  const Teuchos::SerialDenseVector<int,Scalar>& c,
104  const int order,
105  const int orderMin,
106  const int orderMax,
107  const std::string& longDescription,
108  bool isEmbedded = false,
109  const Teuchos::SerialDenseVector<int,Scalar>&
110  bstar = Teuchos::SerialDenseVector<int,Scalar>())
111  {
112  const int numStages = A.numRows();
113  TEUCHOS_ASSERT_EQUALITY( A.numCols(), numStages );
114  TEUCHOS_ASSERT_EQUALITY( b.length(), numStages );
115  TEUCHOS_ASSERT_EQUALITY( c.length(), numStages );
116  TEUCHOS_ASSERT( order > 0 );
117  A_ = A;
118  b_ = b;
119  c_ = c;
120  order_ = order;
121  orderMin_ = orderMin;
122  orderMax_ = orderMax;
123  this->set_isImplicit();
124  this->set_isDIRK();
125  longDescription_ = longDescription;
126 
127  if (isEmbedded) {
128  TEUCHOS_ASSERT_EQUALITY( bstar.length(), numStages );
129  bstar_ = bstar;
130  isEmbedded_ = true;
131  }
132  }
133 
134  /* \brief Redefined from Teuchos::ParameterListAcceptorDefaultBase */
135  //@{
136  virtual void setParameterList(
137  Teuchos::RCP<Teuchos::ParameterList> const& pList)
138  {
139  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
140  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
141  else pl = pList;
142  pl->validateParameters(*this->getValidParameters());
143  TEUCHOS_TEST_FOR_EXCEPTION(
144  pl->get<std::string>("Stepper Type") != this->description()
145  ,std::runtime_error,
146  " Stepper Type != \""+this->description()+"\"\n"
147  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
148  this->setMyParamList(pl);
149  this->rkbtPL_ = pl;
150  }
151 
152  virtual Teuchos::RCP<const Teuchos::ParameterList>
153  getValidParameters() const
154  {
155  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
156  pl->setName("Default Stepper - " + this->description());
157  pl->set<std::string>("Description", this->getDescription());
158  pl->set<std::string>("Stepper Type", this->description());
159  pl->set<bool>("Use Embedded", false);
160  return pl;
161  }
162  //@}
163 
164  /* \brief Redefined from Teuchos::Describable */
165  //@{
166  virtual std::string description() const = 0;
167 
168  virtual void describe( Teuchos::FancyOStream &out,
169  const Teuchos::EVerbosityLevel verbLevel) const
170  {
171  if (verbLevel != Teuchos::VERB_NONE) {
172  out << this->description() << std::endl;
173  out << this->getDescription() << std::endl;
174  out << "number of Stages = " << this->numStages() << std::endl;
175  out << "A = " << this->A() << std::endl;
176  out << "b = " << this->b() << std::endl;
177  out << "c = " << this->c() << std::endl;
178  out << "bstar = " << this->bstar() << std::endl;
179  out << "order = " << this->order() << std::endl;
180  out << "orderMin = " << this->orderMin() << std::endl;
181  out << "orderMax = " << this->orderMax() << std::endl;
182  out << "isImplicit = " << this->isImplicit() << std::endl;
183  out << "isDIRK = " << this->isDIRK() << std::endl;
184  out << "isEmbedded = " << this->isEmbedded() << std::endl;
185  }
186  }
187  //@}
188 
189  protected:
190  virtual void setDescription(std::string longD) { longDescription_ = longD; }
191  const std::string& getDescription() const { return longDescription_; }
192 
193  void set_A(const Teuchos::SerialDenseMatrix<int,Scalar>& A) { A_ = A; }
194  void set_b(const Teuchos::SerialDenseVector<int,Scalar>& b) { b_ = b; }
195  void set_c(const Teuchos::SerialDenseVector<int,Scalar>& c) { c_ = c; }
196  void set_order(const int& order) { order_ = order; }
197  void set_orderMin(const int& order) { orderMin_ = order; }
198  void set_orderMax(const int& order) { orderMax_ = order; }
199  void set_isImplicit() {
200  isImplicit_ = false;
201  for (size_t i = 0; i < this->numStages(); i++)
202  for (size_t j = i; j < this->numStages(); j++)
203  if (A_(i,j) != 0.0) isImplicit_ = true;
204  }
205  /// DIRK is defined as if a_ij = 0 for j>i and a_ii != 0 for at least one i.
206  void set_isDIRK() {
207  isDIRK_ = true;
208  bool nonZero = false;
209  for (size_t i = 0; i < this->numStages(); i++) {
210  if (A_(i,i) != 0.0) nonZero = true;
211  for (size_t j = i+1; j < this->numStages(); j++)
212  if (A_(i,j) != 0.0) isDIRK_ = false;
213  }
214  if (nonZero == false) isDIRK_ = false;
215  }
216 
217  private:
218  Teuchos::SerialDenseMatrix<int,Scalar> A_;
219  Teuchos::SerialDenseVector<int,Scalar> b_;
220  Teuchos::SerialDenseVector<int,Scalar> c_;
221  int order_;
222  int orderMin_;
223  int orderMax_;
224  bool isImplicit_;
225  bool isDIRK_;
226  std::string longDescription_;
227 
228  bool isEmbedded_ = false;
229  Teuchos::SerialDenseVector<int,Scalar> bstar_;
230 
231  protected:
232  Teuchos::RCP<Teuchos::ParameterList> rkbtPL_;
233 };
234 
235 // ----------------------------------------------------------------------------
236 // Nonmember constructor
237 template<class Scalar>
238 Teuchos::RCP<RKButcherTableau<Scalar> > rKButcherTableau()
239 {
240  return(rcp(new RKButcherTableau<Scalar>()));
241 }
242 
243 // Nonmember constructor
244 template<class Scalar>
245 Teuchos::RCP<RKButcherTableau<Scalar> > rKButcherTableau(
246  const Teuchos::SerialDenseMatrix<int,Scalar>& A,
247  const Teuchos::SerialDenseVector<int,Scalar>& b,
248  const Teuchos::SerialDenseVector<int,Scalar>& c,
249  int order,
250  const std::string& description = ""
251  )
252 {
253  Teuchos::RCP<RKButcherTableau<Scalar> > rkbt =
254  rcp(new RKButcherTableau<Scalar>());
255  rkbt->initialize(A,b,c,order,order,order,description);
256  return(rkbt);
257 }
258 
259 
260 // Nonmember constructor
261 template<class Scalar>
262 Teuchos::RCP<RKButcherTableau<Scalar> > rKButcherTableau(
263  const Teuchos::SerialDenseMatrix<int,Scalar>& A,
264  const Teuchos::SerialDenseVector<int,Scalar>& b,
265  const Teuchos::SerialDenseVector<int,Scalar>& c,
266  int order,
267  int orderMin,
268  int orderMax,
269  const std::string& description = ""
270  )
271 {
272  Teuchos::RCP<RKButcherTableau<Scalar> > rkbt =
273  rcp(new RKButcherTableau<Scalar>());
274  rkbt->initialize(A,b,c,order,orderMin,orderMax,description);
275  return(rkbt);
276 }
277 
278 template<class Scalar>
279 class General_RKButcherTableau :
280  virtual public RKButcherTableau<Scalar>
281 {
282  protected:
283  void parseGeneralPL(Teuchos::RCP<Teuchos::ParameterList> const& pList)
284  {
285  using Teuchos::as;
286 
287  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
288  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
289  else pl = pList;
290  // Can not validate because optional parameters (e.g., bstar).
291  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
292  TEUCHOS_TEST_FOR_EXCEPTION(
293  pl->get<std::string>("Stepper Type") != this->description()
294  ,std::runtime_error,
295  " Stepper Type != \""+this->description()+"\"\n"
296  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
297 
298  Teuchos::RCP<Teuchos::ParameterList> tableauPL = sublist(pl,"Tableau",true);
299  std::size_t numStages = 0;
300  int order = tableauPL->get<int>("order");
301  Teuchos::SerialDenseMatrix<int,Scalar> A;
302  Teuchos::SerialDenseVector<int,Scalar> b;
303  Teuchos::SerialDenseVector<int,Scalar> c;
304  Teuchos::SerialDenseVector<int,Scalar> bstar;
305 
306  // read in the A matrix
307  {
308  std::vector<std::string> A_row_tokens;
309  Tempus::StringTokenizer(A_row_tokens, tableauPL->get<std::string>("A"),
310  ";",true);
311 
312  // this is the only place where numStages is set
313  numStages = A_row_tokens.size();
314 
315  // allocate the matrix
316  A.shape(as<int>(numStages),as<int>(numStages));
317 
318  // fill the rows
319  for(std::size_t row=0;row<numStages;row++) {
320  // parse the row (tokenize on space)
321  std::vector<std::string> tokens;
322  Tempus::StringTokenizer(tokens,A_row_tokens[row]," ",true);
323 
324  std::vector<double> values;
325  Tempus::TokensToDoubles(values,tokens);
326 
327  TEUCHOS_TEST_FOR_EXCEPTION(values.size()!=numStages,std::runtime_error,
328  "Error parsing A matrix, wrong number of stages in row "
329  << row << "\n" + this->description());
330 
331  for(std::size_t col=0;col<numStages;col++)
332  A(row,col) = values[col];
333  }
334  }
335 
336  // size b and c vectors
337  b.size(as<int>(numStages));
338  c.size(as<int>(numStages));
339 
340  // read in the b vector
341  {
342  std::vector<std::string> tokens;
343  Tempus::StringTokenizer(tokens,tableauPL->get<std::string>("b")," ",true);
344  std::vector<double> values;
345  Tempus::TokensToDoubles(values,tokens);
346 
347  TEUCHOS_TEST_FOR_EXCEPTION(values.size()!=numStages,std::runtime_error,
348  "Error parsing b vector, wrong number of stages.\n"
349  + this->description());
350 
351  for(std::size_t i=0;i<numStages;i++)
352  b(i) = values[i];
353  }
354 
355  // read in the c vector
356  {
357  std::vector<std::string> tokens;
358  Tempus::StringTokenizer(tokens,tableauPL->get<std::string>("c")," ",true);
359  std::vector<double> values;
360  Tempus::TokensToDoubles(values,tokens);
361 
362  TEUCHOS_TEST_FOR_EXCEPTION(values.size()!=numStages,std::runtime_error,
363  "Error parsing c vector, wrong number of stages.\n"
364  + this->description());
365 
366  for(std::size_t i=0;i<numStages;i++)
367  c(i) = values[i];
368  }
369 
370  if (tableauPL->isParameter("bstar") and
371  tableauPL->get<std::string>("bstar") != "1.0") {
372  bstar.size(as<int>(numStages));
373  // read in the bstar vector
374  {
375  std::vector<std::string> tokens;
376  Tempus::StringTokenizer(
377  tokens, tableauPL->get<std::string>("bstar"), " ", true);
378  std::vector<double> values;
379  Tempus::TokensToDoubles(values,tokens);
380 
381  TEUCHOS_TEST_FOR_EXCEPTION(values.size()!=numStages,std::runtime_error,
382  "Error parsing bstar vector, wrong number of stages.\n"
383  + this->description());
384 
385  for(std::size_t i=0;i<numStages;i++)
386  bstar(i) = values[i];
387  }
388  this->initialize(A,b,c,order,this->getDescription(), true, bstar);
389  } else {
390  this->initialize(A,b,c,order,this->getDescription());
391  }
392 
393  this->setMyParamList(pl);
394  this->rkbtPL_ = pl;
395  }
396 
397 };
398 
399 // ----------------------------------------------------------------------------
400 /** \brief General Explicit Runge-Kutta Butcher Tableau
401  *
402  * The format of the Butcher Tableau parameter list is
403  \verbatim
404  <Parameter name="A" type="string" value="# # # ;
405  # # # ;
406  # # #">
407  <Parameter name="b" type="string" value="# # #">
408  <Parameter name="c" type="string" value="# # #">
409  \endverbatim
410  * Note the number of stages is implicit in the number of entries.
411  * The number of stages must be consistent.
412  *
413  * Default tableau is RK4 (order=4):
414  * \f[
415  * \begin{array}{c|c}
416  * c & A \\ \hline
417  * & b^T
418  * \end{array}
419  * \;\;\;\;\mbox{ where }\;\;\;\;
420  * \begin{array}{c|cccc} 0 & 0 & & & \\
421  * 1/2 & 1/2 & 0 & & \\
422  * 1/2 & 0 & 1/2 & 0 & \\
423  * 1 & 0 & 0 & 1 & 0 \\ \hline
424  * & 1/6 & 1/3 & 1/3 & 1/6 \end{array}
425  * \f]
426  */
427 template<class Scalar>
428 class GeneralExplicit_RKBT :
429  virtual public General_RKButcherTableau<Scalar>
430 {
431  public:
432  GeneralExplicit_RKBT()
433  {
434  std::stringstream Description;
435  Description << this->description() << "\n"
436  << "The format of the Butcher Tableau parameter list is\n"
437  << " <Parameter name=\"A\" type=\"string\" value=\"# # # ;\n"
438  << " # # # ;\n"
439  << " # # #\"/>\n"
440  << " <Parameter name=\"b\" type=\"string\" value=\"# # #\"/>\n"
441  << " <Parameter name=\"c\" type=\"string\" value=\"# # #\"/>\n\n"
442  << "Note the number of stages is implicit in the number of entries.\n"
443  << "The number of stages must be consistent.\n"
444  << "\n"
445  << "Default tableau is RK4 (order=4):\n"
446  << "c = [ 0 1/2 1/2 1 ]'\n"
447  << "A = [ 0 ]\n"
448  << " [ 1/2 0 ]\n"
449  << " [ 0 1/2 0 ]\n"
450  << " [ 0 0 1 0 ]\n"
451  << "b = [ 1/6 1/3 1/3 1/6 ]'" << std::endl;
452 
453  this->setDescription(Description.str());
454  this->setParameterList(Teuchos::null);
455  }
456 
457  virtual std::string description() const { return "General ERK"; }
458 
459  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
460  {
461  this->parseGeneralPL(pList);
462  TEUCHOS_TEST_FOR_EXCEPTION(this->isImplicit() == true, std::logic_error,
463  "Error - General ERK received an implicit Butcher Tableau!\n");
464  }
465 
466  Teuchos::RCP<const Teuchos::ParameterList>
467  getValidParameters() const
468  {
469  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
470  pl->setName("Default Stepper - " + this->description());
471  pl->set<std::string>("Description", this->getDescription());
472  pl->set<std::string>("Stepper Type", this->description());
473  pl->set<bool>("Use Embedded", false);
474 
475  // Tableau ParameterList
476  Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
477  tableauPL->set<std::string>("A",
478  "0.0 0.0 0.0 0.0; 0.5 0.0 0.0 0.0; 0.0 0.5 0.0 0.0; 0.0 0.0 1.0 0.0");
479  tableauPL->set<std::string>("b",
480  "0.166666666666667 0.333333333333333 0.333333333333333 0.166666666666667");
481  //tableauPL->set<std::string>("bstar", "1.0");
482  tableauPL->set<std::string>("c", "0.0 0.5 0.5 1.0");
483  tableauPL->set<int>("order", 4);
484  pl->set("Tableau", *tableauPL);
485 
486  return pl;
487  }
488 };
489 
490 
491 // ----------------------------------------------------------------------------
492 /** \brief Backward Euler Runge-Kutta Butcher Tableau
493  *
494  * The tableau for Backward Euler (order=1) is
495  * \f[
496  * \begin{array}{c|c}
497  * c & A \\ \hline
498  * & b^T
499  * \end{array}
500  * \;\;\;\;\mbox{ where }\;\;\;\;
501  * \begin{array}{c|c} 1 & 1 \\ \hline
502  * & 1 \end{array}
503  * \f]
504  */
505 template<class Scalar>
506 class BackwardEuler_RKBT :
507  virtual public RKButcherTableau<Scalar>
508 {
509  public:
510  BackwardEuler_RKBT()
511  {
512  std::ostringstream Description;
513  Description << this->description() << "\n"
514  << "c = [ 1 ]'\n"
515  << "A = [ 1 ]\n"
516  << "b = [ 1 ]'" << std::endl;
517 
518  this->setDescription(Description.str());
519  this->setParameterList(Teuchos::null);
520  }
521 
522  virtual std::string description() const { return "RK Backward Euler"; }
523 
524  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
525  {
526  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
527  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
528  else pl = pList;
529  // Can not validate because optional parameters (e.g., Solver Name).
530  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
531  TEUCHOS_TEST_FOR_EXCEPTION(
532  pl->get<std::string>("Stepper Type") != this->description()
533  ,std::runtime_error,
534  " Stepper Type != \""+this->description()+"\"\n"
535  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
536 
537  typedef Teuchos::ScalarTraits<Scalar> ST;
538  Teuchos::SerialDenseMatrix<int,Scalar> A(1,1);
539  A(0,0) = ST::one();
540  Teuchos::SerialDenseVector<int,Scalar> b(1);
541  b(0) = ST::one();
542  Teuchos::SerialDenseVector<int,Scalar> c(1);
543  c(0) = ST::one();
544  int order = 1;
545 
546  this->initialize(A,b,c,order,this->getDescription());
547  this->setMyParamList(pl);
548  this->rkbtPL_ = pl;
549  }
550 
551  Teuchos::RCP<const Teuchos::ParameterList>
552  getValidParameters() const
553  {
554  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
555  pl->setName("Default Stepper - " + this->description());
556  pl->set<std::string>("Description", this->getDescription());
557  pl->set<std::string>("Stepper Type", this->description());
558  pl->set<bool>("Use Embedded", false);
559  pl->set<std::string>("Solver Name", "",
560  "Name of ParameterList containing the solver specifications.");
561 
562  return pl;
563  }
564 };
565 
566 
567 // ----------------------------------------------------------------------------
568 /** \brief Forward Euler Runge-Kutta Butcher Tableau
569  *
570  * The tableau for Forward Euler (order=1) is
571  * \f[
572  * \begin{array}{c|c}
573  * c & A \\ \hline
574  * & b^T
575  * \end{array}
576  * \;\;\;\;\mbox{ where }\;\;\;\;
577  * \begin{array}{c|c}
578  * c & A \\ \hline
579  * & b^T
580  * \end{array}
581  * \;\;\;\;\mbox{ where }\;\;\;\;
582  * \begin{array}{c|c} 0 & 0 \\ \hline
583  * & 1 \end{array}
584  * \f]
585  */
586 template<class Scalar>
587 class ForwardEuler_RKBT :
588  virtual public RKButcherTableau<Scalar>
589 {
590  public:
591  ForwardEuler_RKBT()
592  {
593  std::ostringstream Description;
594  Description << this->description() << "\n"
595  << "c = [ 0 ]'\n"
596  << "A = [ 0 ]\n"
597  << "b = [ 1 ]'" << std::endl;
598  typedef Teuchos::ScalarTraits<Scalar> ST;
599  Teuchos::SerialDenseMatrix<int,Scalar> A(1,1);
600  Teuchos::SerialDenseVector<int,Scalar> b(1);
601  b(0) = ST::one();
602  Teuchos::SerialDenseVector<int,Scalar> c(1);
603  int order = 1;
604 
605  this->initialize(A,b,c,order,Description.str());
606  }
607  virtual std::string description() const { return "RK Forward Euler"; }
608 };
609 
610 
611 // ----------------------------------------------------------------------------
612 /** \brief Runge-Kutta 4th order Butcher Tableau
613  *
614  * The tableau for RK4 (order=4) is
615  * \f[
616  * \begin{array}{c|c}
617  * c & A \\ \hline
618  * & b^T
619  * \end{array}
620  * \;\;\;\;\mbox{ where }\;\;\;\;
621  * \begin{array}{c|cccc} 0 & 0 & & & \\
622  * 1/2 & 1/2 & 0 & & \\
623  * 1/2 & 0 & 1/2 & 0 & \\
624  * 1 & 0 & 0 & 1 & 0 \\ \hline
625  * & 1/6 & 1/3 & 1/3 & 1/6 \end{array}
626  * \f]
627  */
628 template<class Scalar>
629 class Explicit4Stage4thOrder_RKBT :
630  virtual public RKButcherTableau<Scalar>
631 {
632  public:
633  Explicit4Stage4thOrder_RKBT()
634  {
635  std::ostringstream Description;
636  Description << this->description() << "\n"
637  << "\"The\" Runge-Kutta Method (explicit):\n"
638  << "Solving Ordinary Differential Equations I:\n"
639  << "Nonstiff Problems, 2nd Revised Edition\n"
640  << "E. Hairer, S.P. Norsett, G. Wanner\n"
641  << "Table 1.2, pg 138\n"
642  << "c = [ 0 1/2 1/2 1 ]'\n"
643  << "A = [ 0 ] \n"
644  << " [ 1/2 0 ]\n"
645  << " [ 0 1/2 0 ]\n"
646  << " [ 0 0 1 0 ]\n"
647  << "b = [ 1/6 1/3 1/3 1/6 ]'" << std::endl;
648  typedef Teuchos::ScalarTraits<Scalar> ST;
649  const Scalar one = ST::one();
650  const Scalar zero = ST::zero();
651  const Scalar onehalf = one/(2*one);
652  const Scalar onesixth = one/(6*one);
653  const Scalar onethird = one/(3*one);
654 
655  int NumStages = 4;
656  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
657  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
658  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
659 
660  // Fill A:
661  A(0,0) = zero;
662  A(0,1) = zero;
663  A(0,2) = zero;
664  A(0,3) = zero;
665 
666  A(1,0) = onehalf;
667  A(1,1) = zero;
668  A(1,2) = zero;
669  A(1,3) = zero;
670 
671  A(2,0) = zero;
672  A(2,1) = onehalf;
673  A(2,2) = zero;
674  A(2,3) = zero;
675 
676  A(3,0) = zero;
677  A(3,1) = zero;
678  A(3,2) = one;
679  A(3,3) = zero;
680 
681  // Fill b:
682  b(0) = onesixth;
683  b(1) = onethird;
684  b(2) = onethird;
685  b(3) = onesixth;
686 
687  // fill b_c_
688  c(0) = zero;
689  c(1) = onehalf;
690  c(2) = onehalf;
691  c(3) = one;
692 
693  int order = 4;
694 
695  this->initialize(A,b,c,order,Description.str());
696  }
697  virtual std::string description() const { return "RK Explicit 4 Stage"; }
698 };
699 
700 
701 // ----------------------------------------------------------------------------
702 /** \brief Explicit RK Bogacki-Shampine Butcher Tableau
703  *
704  * The tableau (order=3(2)) is
705  * \f[
706  * \begin{array}{c|c}
707  * c & A \\ \hline
708  * & b^T \\ \hline
709  * & \hat{b}^T
710  * \end{array}
711  * \;\;\;\;\mbox{ where }\;\;\;\;
712  * \begin{array}{c|cccc} 0 & 0 & & & \\
713  * 1/3 & 1/2 & 0 & & \\
714  * 2/3 & 0 & 3/4 & 0 & \\
715  * 1 & 2/9 & 1/3 & 4/9 & 0 \\ \hline
716  * & 2/9 & 1/3 & 4/9 & 0 \\
717  * & 7/24 & 1/4 & 1/3 & 1/8 \end{array}
718  * \f]
719  * Reference: P. Bogacki and L.F. Shampine.
720  * A 3(2) pair of Runge–Kutta formulas.
721  * Applied Mathematics Letters, 2(4):321 – 325, 1989.
722  *
723  */
724 template<class Scalar>
725 class ExplicitBogackiShampine32_RKBT :
726  virtual public RKButcherTableau<Scalar>
727 {
728  public:
729  ExplicitBogackiShampine32_RKBT()
730  {
731  std::ostringstream Description;
732  Description << this->description() << "\n"
733  << "P. Bogacki and L.F. Shampine.\n"
734  << "A 3(2) pair of Runge–Kutta formulas.\n"
735  << "Applied Mathematics Letters, 2(4):321 – 325, 1989.\n"
736  << "c = [ 0 1/3 2/3 1 ]'\n"
737  << "A = [ 0 ]\n"
738  << " [ 1/2 0 ]\n"
739  << " [ 0 3/4 0 ]\n"
740  << " [ 2/9 1/3 4/9 0 ]\n"
741  << "b = [ 2/9 1/3 4/9 0 ]\n"
742  << "bstar = [ 7/24 1/4 1/3 1/8 ]\n" << std::endl;
743  typedef Teuchos::ScalarTraits<Scalar> ST;
744  using Teuchos::as;
745  int NumStages = 4;
746  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
747  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
748  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
749  Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);
750 
751  const Scalar one = ST::one();
752  const Scalar zero = ST::zero();
753 
754  // Fill A:
755  A(0,0) = zero;
756  A(0,1) = zero;
757  A(0,2) = zero;
758  A(0,3) = zero;
759 
760  A(1,0) = as<Scalar>(one/(2*one));
761  A(1,1) = zero;
762  A(1,2) = zero;
763  A(1,3) = zero;
764 
765  A(2,0) = zero;
766  A(2,1) = as<Scalar>(3*one/(4*one));
767  A(2,2) = zero;
768  A(2,3) = zero;
769 
770  A(3,0) = as<Scalar>(2*one/(9*one));
771  A(3,1) = as<Scalar>(1*one/(3*one));
772  A(3,2) = as<Scalar>(4*one/(9*one));
773  A(3,3) = zero;
774 
775  // Fill b:
776  b(0) = A(3,0);
777  b(1) = A(3,1);
778  b(2) = A(3,2);
779  b(3) = A(3,3);
780 
781  // Fill c:
782  c(0) = zero;
783  c(1) = as<Scalar>(1*one/(3*one));
784  c(2) = as<Scalar>(2*one/(3*one));
785  c(3) = one;
786 
787  // Fill bstar
788  bstar(0) = as<Scalar>(7.0*one/(24*one));
789  bstar(1) = as<Scalar>(1*one/(4*one));
790  bstar(2) = as<Scalar>(1*one/(3*one));
791  bstar(3) = as<Scalar>(1*one/(8*one));
792  int order = 3;
793 
794  this->initialize(A,b,c,order,Description.str(),true,bstar);
795  }
796  virtual std::string description() const { return "Bogacki-Shampine 3(2) Pair"; }
797 };
798 
799 
800 // ----------------------------------------------------------------------------
801 /** \brief Explicit RK Merson Butcher Tableau
802  *
803  * The tableau (order=4(5)) is
804  * \f[
805  * \begin{array}{c|c}
806  * c & A \\ \hline
807  * & b^T \\ \hline
808  * & \hat{b}^T
809  * \end{array}
810  * \;\;\;\;\mbox{ where }\;\;\;\;
811  * \begin{array}{c|cccc} 0 & 0 & & & & \\
812  * 1/3 & 1/3 & 0 & & & \\
813  * 1/3 & 1/6 & 1/6 & 0 & & \\
814  * 1/2 & 1/8 & 0 & 3/8 & & \\
815  * 1 & 1/2 & 0 & -3/2 & 2 & \\ \hline
816  * & 1/6 & 0 & 0 & 2/3 & 1/6 \\
817  * & 1/10 & 0 & 3/10 & 2/5 & 1/5 \end{array}
818  * \f]
819  * Reference: E. Hairer, S.P. Norsett, G. Wanner,
820  * "Solving Ordinary Differential Equations I:
821  * Nonstiff Problems", 2nd Revised Edition,
822  * Table 4.1, pg 167.
823  *
824  */
825 template<class Scalar>
826 class ExplicitMerson45_RKBT :
827  virtual public RKButcherTableau<Scalar>
828 {
829  public:
830  ExplicitMerson45_RKBT()
831  {
832  std::ostringstream Description;
833  Description << this->description() << "\n"
834  << "Solving Ordinary Differential Equations I:\n"
835  << "Nonstiff Problems, 2nd Revised Edition\n"
836  << "E. Hairer, S.P. Norsett, G. Wanner\n"
837  << "Table 4.1, pg 167\n"
838  << "c = [ 0 1/3 1/3 1/2 1 ]'\n"
839  << "A = [ 0 ]\n"
840  << " [ 1/3 0 ]\n"
841  << " [ 1/6 1/6 0 ]\n"
842  << " [ 1/8 0 3/8 0 ]\n"
843  << " [ 1/2 0 -3/2 2 0 ]\n"
844  << "b = [ 1/6 0 0 2/3 1/6 ]\n"
845  << "bstar = [ 1/10 0 3/10 2/5 1/5 ]\n" << std::endl;
846  typedef Teuchos::ScalarTraits<Scalar> ST;
847  using Teuchos::as;
848  int NumStages = 5;
849  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages, true);
850  Teuchos::SerialDenseVector<int,Scalar> b(NumStages, true);
851  Teuchos::SerialDenseVector<int,Scalar> c(NumStages, true);
852  Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages, true);
853 
854  const Scalar one = ST::one();
855  const Scalar zero = ST::zero();
856 
857  // Fill A:
858  A(1,0) = as<Scalar>(one/(3*one));;
859 
860  A(2,0) = as<Scalar>(one/(6*one));;
861  A(2,1) = as<Scalar>(one/(6*one));;
862 
863  A(3,0) = as<Scalar>(one/(8*one));;
864  A(3,2) = as<Scalar>(3*one/(8*one));;
865 
866  A(4,0) = as<Scalar>(one/(2*one));;
867  A(4,2) = as<Scalar>(-3*one/(2*one));;
868  A(4,3) = 2*one;
869 
870  // Fill b:
871  b(0) = as<Scalar>(one/(6*one));
872  b(3) = as<Scalar>(2*one/(3*one));
873  b(4) = as<Scalar>(one/(6*one));
874 
875  // Fill c:
876  c(0) = zero;
877  c(1) = as<Scalar>(1*one/(3*one));
878  c(2) = as<Scalar>(1*one/(3*one));
879  c(3) = as<Scalar>(1*one/(2*one));
880  c(4) = one;
881 
882  // Fill bstar
883  bstar(0) = as<Scalar>(1*one/(10*one));
884  bstar(2) = as<Scalar>(3*one/(10*one));
885  bstar(3) = as<Scalar>(2*one/(5*one));
886  bstar(4) = as<Scalar>(1*one/(5*one));
887  int order = 4;
888 
889  this->initialize(A,b,c,order,Description.str(),true,bstar);
890  }
891  virtual std::string description() const { return "Merson 4(5) Pair"; }
892 };
893 
894 // ----------------------------------------------------------------------------
895 /** \brief Explicit RK 3/8th Rule Butcher Tableau
896  *
897  * The tableau (order=4) is
898  * \f[
899  * \begin{array}{c|c}
900  * c & A \\ \hline
901  * & b^T
902  * \end{array}
903  * \;\;\;\;\mbox{ where }\;\;\;\;
904  * \begin{array}{c|cccc} 0 & 0 & & & \\
905  * 1/3 & 1/3 & 0 & & \\
906  * 2/3 &-1/3 & 1 & 0 & \\
907  * 1 & 1 & -1 & 1 & 0 \\ \hline
908  * & 1/8 & 3/8 & 3/8 & 1/8 \end{array}
909  * \f]
910  * Reference: E. Hairer, S.P. Norsett, G. Wanner,
911  * "Solving Ordinary Differential Equations I:
912  * Nonstiff Problems", 2nd Revised Edition,
913  * Table 1.2, pg 138.
914  */
915 template<class Scalar>
916 class Explicit3_8Rule_RKBT :
917  virtual public RKButcherTableau<Scalar>
918 {
919  public:
920  Explicit3_8Rule_RKBT()
921  {
922  std::ostringstream Description;
923  Description << this->description() << "\n"
924  << "Solving Ordinary Differential Equations I:\n"
925  << "Nonstiff Problems, 2nd Revised Edition\n"
926  << "E. Hairer, S.P. Norsett, G. Wanner\n"
927  << "Table 1.2, pg 138\n"
928  << "c = [ 0 1/3 2/3 1 ]'\n"
929  << "A = [ 0 ]\n"
930  << " [ 1/3 0 ]\n"
931  << " [-1/3 1 0 ]\n"
932  << " [ 1 -1 1 0 ]\n"
933  << "b = [ 1/8 3/8 3/8 1/8 ]'" << std::endl;
934  typedef Teuchos::ScalarTraits<Scalar> ST;
935  using Teuchos::as;
936  int NumStages = 4;
937  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
938  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
939  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
940 
941  const Scalar one = ST::one();
942  const Scalar zero = ST::zero();
943  const Scalar one_third = as<Scalar>(one/(3*one));
944  const Scalar two_third = as<Scalar>(2*one/(3*one));
945  const Scalar one_eighth = as<Scalar>(one/(8*one));
946  const Scalar three_eighth = as<Scalar>(3*one/(8*one));
947 
948  // Fill A:
949  A(0,0) = zero;
950  A(0,1) = zero;
951  A(0,2) = zero;
952  A(0,3) = zero;
953 
954  A(1,0) = one_third;
955  A(1,1) = zero;
956  A(1,2) = zero;
957  A(1,3) = zero;
958 
959  A(2,0) = as<Scalar>(-one_third);
960  A(2,1) = one;
961  A(2,2) = zero;
962  A(2,3) = zero;
963 
964  A(3,0) = one;
965  A(3,1) = as<Scalar>(-one);
966  A(3,2) = one;
967  A(3,3) = zero;
968 
969  // Fill b:
970  b(0) = one_eighth;
971  b(1) = three_eighth;
972  b(2) = three_eighth;
973  b(3) = one_eighth;
974 
975  // Fill c:
976  c(0) = zero;
977  c(1) = one_third;
978  c(2) = two_third;
979  c(3) = one;
980 
981  int order = 4;
982 
983  this->initialize(A,b,c,order,Description.str());
984  }
985  virtual std::string description() const { return "RK Explicit 3/8 Rule"; }
986 };
987 
988 
989 // ----------------------------------------------------------------------------
990 /** \brief RK Explicit 4 Stage 3rd order by Runge
991  *
992  * The tableau (order=3) is
993  * \f[
994  * \begin{array}{c|c}
995  * c & A \\ \hline
996  * & b^T
997  * \end{array}
998  * \;\;\;\;\mbox{ where }\;\;\;\;
999  * \begin{array}{c|cccc} 0 & 0 & & & \\
1000  * 1/2 & 1/2 & 0 & & \\
1001  * 1 & 0 & 1 & 0 & \\
1002  * 1 & 0 & 0 & 1 & 0 \\ \hline
1003  * & 1/6 & 2/3 & 0 & 1/6 \end{array}
1004  * \f]
1005  * Reference: E. Hairer, S.P. Norsett, G. Wanner,
1006  * "Solving Ordinary Differential Equations I:
1007  * Nonstiff Problems", 2nd Revised Edition,
1008  * Table 1.1, pg 135.
1009  */
1010 template<class Scalar>
1011 class Explicit4Stage3rdOrderRunge_RKBT :
1012  virtual public RKButcherTableau<Scalar>
1013 {
1014  public:
1015  Explicit4Stage3rdOrderRunge_RKBT()
1016  {
1017  std::ostringstream Description;
1018  Description << this->description() << "\n"
1019  << "Solving Ordinary Differential Equations I:\n"
1020  << "Nonstiff Problems, 2nd Revised Edition\n"
1021  << "E. Hairer, S.P. Norsett, G. Wanner\n"
1022  << "Table 1.1, pg 135\n"
1023  << "c = [ 0 1/2 1 1 ]'\n"
1024  << "A = [ 0 ]\n"
1025  << " [ 1/2 0 ]\n"
1026  << " [ 0 1 0 ]\n"
1027  << " [ 0 0 1 0 ]\n"
1028  << "b = [ 1/6 2/3 0 1/6 ]'" << std::endl;
1029  typedef Teuchos::ScalarTraits<Scalar> ST;
1030  int NumStages = 4;
1031  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1032  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1033  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1034 
1035  const Scalar one = ST::one();
1036  const Scalar onehalf = one/(2*one);
1037  const Scalar onesixth = one/(6*one);
1038  const Scalar twothirds = 2*one/(3*one);
1039  const Scalar zero = ST::zero();
1040 
1041  // Fill A:
1042  A(0,0) = zero;
1043  A(0,1) = zero;
1044  A(0,2) = zero;
1045  A(0,3) = zero;
1046 
1047  A(1,0) = onehalf;
1048  A(1,1) = zero;
1049  A(1,2) = zero;
1050  A(1,3) = zero;
1051 
1052  A(2,0) = zero;
1053  A(2,1) = one;
1054  A(2,2) = zero;
1055  A(2,3) = zero;
1056 
1057  A(3,0) = zero;
1058  A(3,1) = zero;
1059  A(3,2) = one;
1060  A(3,3) = zero;
1061 
1062  // Fill b:
1063  b(0) = onesixth;
1064  b(1) = twothirds;
1065  b(2) = zero;
1066  b(3) = onesixth;
1067 
1068  // Fill c:
1069  c(0) = zero;
1070  c(1) = onehalf;
1071  c(2) = one;
1072  c(3) = one;
1073 
1074  int order = 3;
1075 
1076  this->initialize(A,b,c,order,Description.str());
1077  }
1078  virtual std::string description() const
1079  { return "RK Explicit 4 Stage 3rd order by Runge"; }
1080 };
1081 
1082 
1083 // ----------------------------------------------------------------------------
1084 /** \brief RK Explicit 5 Stage 3rd order by Kinnmark and Gray
1085  *
1086  * The tableau (order=3) is
1087  * \f[
1088  * \begin{array}{c|c}
1089  * c & A \\ \hline
1090  * & b^T
1091  * \end{array}
1092  * \;\;\;\;\mbox{ where }\;\;\;\;
1093  * \begin{array}{c|ccccc} 0 & 0 & & & & \\
1094  * 1/5 & 1/5 & 0 & & & \\
1095  * 1/5 & 0 & 1/5 & 0 & & \\
1096  * 1/3 & 0 & 0 & 1/3 & 0 & \\
1097  * 2/3 & 0 & 0 & 0 & 2/3 & 0 \\ \hline
1098  * & 1/4 & 0 & 0 & 0 & 3/4 \end{array}
1099  * \f]
1100  * Reference: Modified by P. Ullrich. From the prim_advance_mod.F90
1101  * routine in the HOMME atmosphere model code.
1102  */
1103 template<class Scalar>
1104 class Explicit5Stage3rdOrderKandG_RKBT :
1105  virtual public RKButcherTableau<Scalar>
1106 {
1107  public:
1108  Explicit5Stage3rdOrderKandG_RKBT()
1109  {
1110  std::ostringstream Description;
1111  Description << this->description() << "\n"
1112  << "Kinnmark & Gray 5 stage, 3rd order scheme \n"
1113  << "Modified by P. Ullrich. From the prim_advance_mod.F90 \n"
1114  << "routine in the HOMME atmosphere model code.\n"
1115  << "c = [ 0 1/5 1/5 1/3 2/3 ]'\n"
1116  << "A = [ 0 ]\n"
1117  << " [ 1/5 0 ]\n"
1118  << " [ 0 1/5 0 ]\n"
1119  << " [ 0 0 1/3 0 ]\n"
1120  << " [ 0 0 0 2/3 0 ]\n"
1121  << "b = [ 1/4 0 0 0 3/4 ]'" << std::endl;
1122  typedef Teuchos::ScalarTraits<Scalar> ST;
1123  int NumStages = 5;
1124  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1125  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1126  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1127 
1128  const Scalar one = ST::one();
1129  const Scalar onefifth = one/(5*one);
1130  const Scalar onefourth = one/(4*one);
1131  const Scalar onethird = one/(3*one);
1132  const Scalar twothirds = 2*one/(3*one);
1133  const Scalar threefourths = 3*one/(4*one);
1134  const Scalar zero = ST::zero();
1135 
1136  // Fill A:
1137  A(0,0) = zero;
1138  A(0,1) = zero;
1139  A(0,2) = zero;
1140  A(0,3) = zero;
1141  A(0,4) = zero;
1142 
1143  A(1,0) = onefifth;
1144  A(1,1) = zero;
1145  A(1,2) = zero;
1146  A(1,3) = zero;
1147  A(1,4) = zero;
1148 
1149  A(2,0) = zero;
1150  A(2,1) = onefifth;
1151  A(2,2) = zero;
1152  A(2,3) = zero;
1153  A(2,4) = zero;
1154 
1155  A(3,0) = zero;
1156  A(3,1) = zero;
1157  A(3,2) = onethird;
1158  A(3,3) = zero;
1159  A(3,4) = zero;
1160 
1161  A(4,0) = zero;
1162  A(4,1) = zero;
1163  A(4,2) = zero;
1164  A(4,3) = twothirds;
1165  A(4,4) = zero;
1166 
1167  // Fill b:
1168  b(0) = onefourth;
1169  b(1) = zero;
1170  b(2) = zero;
1171  b(3) = zero;
1172  b(4) = threefourths;
1173 
1174  // Fill c:
1175  c(0) = zero;
1176  c(1) = onefifth;
1177  c(2) = onefifth;
1178  c(3) = onethird;
1179  c(4) = twothirds;
1180 
1181  int order = 3;
1182 
1183  this->initialize(A,b,c,order,Description.str());
1184  }
1185  virtual std::string description() const
1186  { return "RK Explicit 5 Stage 3rd order by Kinnmark and Gray"; }
1187 };
1188 
1189 
1190 // ----------------------------------------------------------------------------
1191 /** \brief RK Explicit 3 Stage 3rd order
1192  *
1193  * The tableau (order=3) is
1194  * \f[
1195  * \begin{array}{c|c}
1196  * c & A \\ \hline
1197  * & b^T
1198  * \end{array}
1199  * \;\;\;\;\mbox{ where }\;\;\;\;
1200  * \begin{array}{c|ccc} 0 & 0 & & \\
1201  * 1/2 & 1/2 & 0 & \\
1202  * 1 & -1 & 2 & 0 \\ \hline
1203  * & 1/6 & 4/6 & 1/6 \end{array}
1204  * \f]
1205  */
1206 template<class Scalar>
1207 class Explicit3Stage3rdOrder_RKBT :
1208  virtual public RKButcherTableau<Scalar>
1209 {
1210  public:
1211  Explicit3Stage3rdOrder_RKBT()
1212  {
1213  std::ostringstream Description;
1214  Description << this->description() << "\n"
1215  << "c = [ 0 1/2 1 ]'\n"
1216  << "A = [ 0 ]\n"
1217  << " [ 1/2 0 ]\n"
1218  << " [ -1 2 0 ]\n"
1219  << "b = [ 1/6 4/6 1/6 ]'" << std::endl;
1220  typedef Teuchos::ScalarTraits<Scalar> ST;
1221  const Scalar one = ST::one();
1222  const Scalar two = Teuchos::as<Scalar>(2*one);
1223  const Scalar zero = ST::zero();
1224  const Scalar onehalf = one/(2*one);
1225  const Scalar onesixth = one/(6*one);
1226  const Scalar foursixth = 4*one/(6*one);
1227 
1228  int NumStages = 3;
1229  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1230  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1231  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1232 
1233  // Fill A:
1234  A(0,0) = zero;
1235  A(0,1) = zero;
1236  A(0,2) = zero;
1237 
1238  A(1,0) = onehalf;
1239  A(1,1) = zero;
1240  A(1,2) = zero;
1241 
1242  A(2,0) = -one;
1243  A(2,1) = two;
1244  A(2,2) = zero;
1245 
1246  // Fill b:
1247  b(0) = onesixth;
1248  b(1) = foursixth;
1249  b(2) = onesixth;
1250 
1251  // fill b_c_
1252  c(0) = zero;
1253  c(1) = onehalf;
1254  c(2) = one;
1255 
1256  int order = 3;
1257 
1258  this->initialize(A,b,c,order,Description.str());
1259  }
1260  virtual std::string description() const
1261  { return "RK Explicit 3 Stage 3rd order"; }
1262 };
1263 
1264 
1265 // ----------------------------------------------------------------------------
1266 /** \brief RK Explicit 3 Stage 3rd order TVD
1267  *
1268  * The tableau (order=3) is
1269  * \f[
1270  * \begin{array}{c|c}
1271  * c & A \\ \hline
1272  * & b^T
1273  * \end{array}
1274  * \;\;\;\;\mbox{ where }\;\;\;\;
1275  * \begin{array}{c|ccc} 0 & 0 & & \\
1276  * 1 & 1 & 0 & \\
1277  * 1/2 & 1/4 & 1/4 & 0 \\ \hline
1278  * & 1/6 & 1/6 & 4/6 \end{array}
1279  * \f]
1280  * Reference: Sigal Gottlieb and Chi-Wang Shu,
1281  * 'Total Variation Diminishing Runge-Kutta Schemes',
1282  * Mathematics of Computation,
1283  * Volume 67, Number 221, January 1998, pp. 73-85.
1284  *
1285  * This is also written in the following set of updates.
1286  \verbatim
1287  u1 = u^n + dt L(u^n)
1288  u2 = 3 u^n/4 + u1/4 + dt L(u1)/4
1289  u^(n+1) = u^n/3 + 2 u2/2 + 2 dt L(u2)/3
1290  \endverbatim
1291  */
1292 template<class Scalar>
1293 class Explicit3Stage3rdOrderTVD_RKBT :
1294  virtual public RKButcherTableau<Scalar>
1295 {
1296  public:
1297  Explicit3Stage3rdOrderTVD_RKBT()
1298  {
1299  std::ostringstream Description;
1300  Description << this->description() << "\n"
1301  << "Sigal Gottlieb and Chi-Wang Shu\n"
1302  << "`Total Variation Diminishing Runge-Kutta Schemes'\n"
1303  << "Mathematics of Computation\n"
1304  << "Volume 67, Number 221, January 1998, pp. 73-85\n"
1305  << "c = [ 0 1 1/2 ]'\n"
1306  << "A = [ 0 ]\n"
1307  << " [ 1 0 ]\n"
1308  << " [ 1/4 1/4 0 ]\n"
1309  << "b = [ 1/6 1/6 4/6 ]'\n"
1310  << "This is also written in the following set of updates.\n"
1311  << "u1 = u^n + dt L(u^n)\n"
1312  << "u2 = 3 u^n/4 + u1/4 + dt L(u1)/4\n"
1313  << "u^(n+1) = u^n/3 + 2 u2/2 + 2 dt L(u2)/3"
1314  << std::endl;
1315  typedef Teuchos::ScalarTraits<Scalar> ST;
1316  const Scalar one = ST::one();
1317  const Scalar zero = ST::zero();
1318  const Scalar onehalf = one/(2*one);
1319  const Scalar onefourth = one/(4*one);
1320  const Scalar onesixth = one/(6*one);
1321  const Scalar foursixth = 4*one/(6*one);
1322 
1323  int NumStages = 3;
1324  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1325  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1326  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1327 
1328  // Fill A:
1329  A(0,0) = zero;
1330  A(0,1) = zero;
1331  A(0,2) = zero;
1332 
1333  A(1,0) = one;
1334  A(1,1) = zero;
1335  A(1,2) = zero;
1336 
1337  A(2,0) = onefourth;
1338  A(2,1) = onefourth;
1339  A(2,2) = zero;
1340 
1341  // Fill b:
1342  b(0) = onesixth;
1343  b(1) = onesixth;
1344  b(2) = foursixth;
1345 
1346  // fill b_c_
1347  c(0) = zero;
1348  c(1) = one;
1349  c(2) = onehalf;
1350 
1351  int order = 3;
1352 
1353  this->initialize(A,b,c,order,Description.str());
1354  }
1355  virtual std::string description() const
1356  { return "RK Explicit 3 Stage 3rd order TVD"; }
1357 };
1358 
1359 
1360 // ----------------------------------------------------------------------------
1361 /** \brief RK Explicit 3 Stage 3rd order by Heun
1362  *
1363  * The tableau (order=3) is
1364  * \f[
1365  * \begin{array}{c|c}
1366  * c & A \\ \hline
1367  * & b^T
1368  * \end{array}
1369  * \;\;\;\;\mbox{ where }\;\;\;\;
1370  * \begin{array}{c|ccc} 0 & 0 & & \\
1371  * 1/3 & 1/3 & 0 & \\
1372  * 2/3 & 0 & 2/3 & 0 \\ \hline
1373  * & 1/4 & 0 & 3/4 \end{array}
1374  * \f]
1375  * Reference: E. Hairer, S.P. Norsett, G. Wanner,
1376  * "Solving Ordinary Differential Equations I:
1377  * Nonstiff Problems", 2nd Revised Edition,
1378  * Table 1.1, pg 135.
1379  */
1380 template<class Scalar>
1381 class Explicit3Stage3rdOrderHeun_RKBT :
1382  virtual public RKButcherTableau<Scalar>
1383 {
1384  public:
1385  Explicit3Stage3rdOrderHeun_RKBT()
1386  {
1387  std::ostringstream Description;
1388  Description << this->description() << "\n"
1389  << "Solving Ordinary Differential Equations I:\n"
1390  << "Nonstiff Problems, 2nd Revised Edition\n"
1391  << "E. Hairer, S.P. Norsett, G. Wanner\n"
1392  << "Table 1.1, pg 135\n"
1393  << "c = [ 0 1/3 2/3 ]'\n"
1394  << "A = [ 0 ] \n"
1395  << " [ 1/3 0 ]\n"
1396  << " [ 0 2/3 0 ]\n"
1397  << "b = [ 1/4 0 3/4 ]'" << std::endl;
1398  typedef Teuchos::ScalarTraits<Scalar> ST;
1399  const Scalar one = ST::one();
1400  const Scalar zero = ST::zero();
1401  const Scalar onethird = one/(3*one);
1402  const Scalar twothirds = 2*one/(3*one);
1403  const Scalar onefourth = one/(4*one);
1404  const Scalar threefourths = 3*one/(4*one);
1405 
1406  int NumStages = 3;
1407  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1408  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1409  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1410 
1411  // Fill A:
1412  A(0,0) = zero;
1413  A(0,1) = zero;
1414  A(0,2) = zero;
1415 
1416  A(1,0) = onethird;
1417  A(1,1) = zero;
1418  A(1,2) = zero;
1419 
1420  A(2,0) = zero;
1421  A(2,1) = twothirds;
1422  A(2,2) = zero;
1423 
1424  // Fill b:
1425  b(0) = onefourth;
1426  b(1) = zero;
1427  b(2) = threefourths;
1428 
1429  // fill b_c_
1430  c(0) = zero;
1431  c(1) = onethird;
1432  c(2) = twothirds;
1433 
1434  int order = 3;
1435 
1436  this->initialize(A,b,c,order,Description.str());
1437  }
1438  virtual std::string description() const
1439  { return "RK Explicit 3 Stage 3rd order by Heun"; }
1440 };
1441 
1442 
1443 // ----------------------------------------------------------------------------
1444 /** \brief RK Explicit 2 Stage 2nd order by Runge
1445  *
1446  * The tableau (order=2) (also known as Explicit Midpoint) is
1447  * \f[
1448  * \begin{array}{c|c}
1449  * c & A \\ \hline
1450  * & b^T
1451  * \end{array}
1452  * \;\;\;\;\mbox{ where }\;\;\;\;
1453  * \begin{array}{c|cc} 0 & 0 & \\
1454  * 1/2 & 1/2 & 0 \\ \hline
1455  * & 0 & 1 \end{array}
1456  * \f]
1457  * Reference: E. Hairer, S.P. Norsett, G. Wanner,
1458  * "Solving Ordinary Differential Equations I:
1459  * Nonstiff Problems", 2nd Revised Edition,
1460  * Table 1.1, pg 135.
1461  */
1462 template<class Scalar>
1463 class Explicit2Stage2ndOrderRunge_RKBT :
1464  virtual public RKButcherTableau<Scalar>
1465 {
1466  public:
1467  Explicit2Stage2ndOrderRunge_RKBT()
1468  {
1469  std::ostringstream Description;
1470  Description << this->description() << "\n"
1471  << "Also known as Explicit Midpoint\n"
1472  << "Solving Ordinary Differential Equations I:\n"
1473  << "Nonstiff Problems, 2nd Revised Edition\n"
1474  << "E. Hairer, S.P. Norsett, G. Wanner\n"
1475  << "Table 1.1, pg 135\n"
1476  << "c = [ 0 1/2 ]'\n"
1477  << "A = [ 0 ]\n"
1478  << " [ 1/2 0 ]\n"
1479  << "b = [ 0 1 ]'" << std::endl;
1480  typedef Teuchos::ScalarTraits<Scalar> ST;
1481  const Scalar one = ST::one();
1482  const Scalar zero = ST::zero();
1483  const Scalar onehalf = one/(2*one);
1484 
1485  int NumStages = 2;
1486  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1487  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1488  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1489 
1490  // Fill A:
1491  A(0,0) = zero;
1492  A(0,1) = zero;
1493 
1494  A(1,0) = onehalf;
1495  A(1,1) = zero;
1496 
1497  // Fill b:
1498  b(0) = zero;
1499  b(1) = one;
1500 
1501  // fill b_c_
1502  c(0) = zero;
1503  c(1) = onehalf;
1504 
1505  int order = 2;
1506 
1507  this->initialize(A,b,c,order,Description.str());
1508  }
1509  virtual std::string description() const
1510  { return "RK Explicit 2 Stage 2nd order by Runge"; }
1511 };
1512 
1513 
1514 // ----------------------------------------------------------------------------
1515 /** \brief RK Explicit Trapezoidal
1516  *
1517  * The tableau (order=2) (also known as Explicit Midpoint) is
1518  * \f[
1519  * \begin{array}{c|c}
1520  * c & A \\ \hline
1521  * & b^T
1522  * \end{array}
1523  * \;\;\;\;\mbox{ where }\;\;\;\;
1524  * \begin{array}{c|cc} 0 & 0 & \\
1525  * 1 & 1 & 0 \\ \hline
1526  * & 1/2 & 1/2 \end{array}
1527  * \f]
1528  */
1529 template<class Scalar>
1530 class ExplicitTrapezoidal_RKBT :
1531  virtual public RKButcherTableau<Scalar>
1532 {
1533  public:
1534  ExplicitTrapezoidal_RKBT()
1535  {
1536  std::ostringstream Description;
1537  Description << this->description() << "\n"
1538  << "c = [ 0 1 ]'\n"
1539  << "A = [ 0 ]\n"
1540  << " [ 1 0 ]\n"
1541  << "b = [ 1/2 1/2 ]'" << std::endl;
1542  typedef Teuchos::ScalarTraits<Scalar> ST;
1543  const Scalar one = ST::one();
1544  const Scalar zero = ST::zero();
1545  const Scalar onehalf = one/(2*one);
1546 
1547  int NumStages = 2;
1548  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1549  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1550  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1551 
1552  // Fill A:
1553  A(0,0) = zero;
1554  A(0,1) = zero;
1555 
1556  A(1,0) = one;
1557  A(1,1) = zero;
1558 
1559  // Fill b:
1560  b(0) = onehalf;
1561  b(1) = onehalf;
1562 
1563  // fill b_c_
1564  c(0) = zero;
1565  c(1) = one;
1566 
1567  int order = 2;
1568 
1569  this->initialize(A,b,c,order,Description.str());
1570  }
1571  virtual std::string description() const { return "RK Explicit Trapezoidal"; }
1572 };
1573 
1574 
1575 // ----------------------------------------------------------------------------
1576 /** \brief General Implicit Runge-Kutta Butcher Tableau
1577  *
1578  * The format of the Butcher Tableau parameter list is
1579  \verbatim
1580  <Parameter name="A" type="string" value="# # # ;
1581  # # # ;
1582  # # #">
1583  <Parameter name="b" type="string" value="# # #">
1584  <Parameter name="c" type="string" value="# # #">
1585  \endverbatim
1586  * Note the number of stages is implicit in the number of entries.
1587  * The number of stages must be consistent.
1588  *
1589  * Default tableau is "SDIRK 2 Stage 2nd order":
1590  * \f[
1591  * \begin{array}{c|c}
1592  * c & A \\ \hline
1593  * & b^T
1594  * \end{array}
1595  * \;\;\;\;\mbox{ where }\;\;\;\;
1596  * \begin{array}{c|cc} \gamma & \gamma & \\
1597  * 1 & 1-\gamma & \gamma \\ \hline
1598  * & 1-\gamma & \gamma \end{array}
1599  * \f]
1600  * where \f$\gamma = (2\pm \sqrt{2})/2\f$. This will produce an
1601  * L-stable 2nd order method.
1602  *
1603  * Reference: U. M. Ascher and L. R. Petzold,
1604  * Computer Methods for ODEs and DAEs, p. 106.
1605  */
1606 template<class Scalar>
1607 class GeneralDIRK_RKBT :
1608  virtual public General_RKButcherTableau<Scalar>
1609 {
1610  public:
1611  GeneralDIRK_RKBT()
1612  {
1613  std::stringstream Description;
1614  Description << this->description() << "\n"
1615  << "The format of the Butcher Tableau parameter list is\n"
1616  << " <Parameter name=\"A\" type=\"string\" value=\"# # # ;\n"
1617  << " # # # ;\n"
1618  << " # # #\"/>\n"
1619  << " <Parameter name=\"b\" type=\"string\" value=\"# # #\"/>\n"
1620  << " <Parameter name=\"c\" type=\"string\" value=\"# # #\"/>\n\n"
1621  << "Note the number of stages is implicit in the number of entries.\n"
1622  << "The number of stages must be consistent.\n"
1623  << "\n"
1624  << "Default tableau is 'SDIRK 2 Stage 2nd order':\n"
1625  << " Computer Methods for ODEs and DAEs\n"
1626  << " U. M. Ascher and L. R. Petzold\n"
1627  << " p. 106\n"
1628  << " gamma = (2+-sqrt(2))/2\n"
1629  << " c = [ gamma 1 ]'\n"
1630  << " A = [ gamma 0 ]\n"
1631  << " [ 1-gamma gamma ]\n"
1632  << " b = [ 1-gamma gamma ]'" << std::endl;
1633 
1634  this->setDescription(Description.str());
1635  this->setParameterList(Teuchos::null);
1636  }
1637 
1638  virtual std::string description() const { return "General DIRK"; }
1639 
1640  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
1641  {
1642  this->parseGeneralPL(pList);
1643  TEUCHOS_TEST_FOR_EXCEPTION(this->isImplicit() != true, std::logic_error,
1644  "Error - General DIRK did not receive a DIRK Butcher Tableau!\n");
1645  }
1646 
1647  Teuchos::RCP<const Teuchos::ParameterList>
1648  getValidParameters() const
1649  {
1650  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1651  pl->setName("Default Stepper - " + this->description());
1652  pl->set<std::string>("Description", this->getDescription());
1653  pl->set<std::string>("Stepper Type", this->description());
1654  pl->set<bool>("Use Embedded", false);
1655 
1656  // Tableau ParameterList
1657  typedef Teuchos::ScalarTraits<Scalar> ST;
1658  std::string gamma = std::to_string((2.0 - ST::squareroot(2.0))/(2.0));
1659  std::string one_gamma = std::to_string(1.0-(2.0-ST::squareroot(2.0))/(2.0));
1660  Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
1661  tableauPL->set<std::string>("A", gamma + " 0.0; " + one_gamma + " "+gamma);
1662  tableauPL->set<std::string>("b", one_gamma + " " + gamma);
1663  tableauPL->set<std::string>("c", gamma + " 1.0");
1664  tableauPL->set<int>("order", 2);
1665  pl->set("Tableau", *tableauPL);
1666 
1667  return pl;
1668  }
1669 };
1670 
1671 
1672 // ----------------------------------------------------------------------------
1673 /** \brief SDIRK 1 Stage 1st order
1674  *
1675  * The tableau (order=1) (also known as Backward Euler) is
1676  * \f[
1677  * \begin{array}{c|c}
1678  * c & A \\ \hline
1679  * & b^T
1680  * \end{array}
1681  * \;\;\;\;\mbox{ where }\;\;\;\;
1682  * \begin{array}{c|c} 1 & 1 \\ \hline
1683  * & 1 \end{array}
1684  * \f]
1685  */
1686 template<class Scalar>
1687 class SDIRK1Stage1stOrder_RKBT :
1688  virtual public RKButcherTableau<Scalar>
1689 {
1690  public:
1691  SDIRK1Stage1stOrder_RKBT()
1692  {
1693  std::stringstream Description;
1694  Description << this->description() << "\n"
1695  << "c = [ 1 ]'\n"
1696  << "A = [ 1 ]\n"
1697  << "b = [ 1 ]'" << std::endl;
1698 
1699  this->setDescription(Description.str());
1700  this->setParameterList(Teuchos::null);
1701  }
1702 
1703  virtual std::string description() const { return "SDIRK 1 Stage 1st order"; }
1704 
1705  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
1706  {
1707  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1708  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
1709  else pl = pList;
1710  // Can not validate because optional parameters (e.g., Solver Name).
1711  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
1712  TEUCHOS_TEST_FOR_EXCEPTION(
1713  pl->get<std::string>("Stepper Type") != this->description()
1714  ,std::runtime_error,
1715  " Stepper Type != \""+this->description()+"\"\n"
1716  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
1717 
1718  typedef Teuchos::ScalarTraits<Scalar> ST;
1719  Teuchos::SerialDenseMatrix<int,Scalar> A(1,1);
1720  A(0,0) = ST::one();
1721  Teuchos::SerialDenseVector<int,Scalar> b(1);
1722  b(0) = ST::one();
1723  Teuchos::SerialDenseVector<int,Scalar> c(1);
1724  c(0) = ST::one();
1725  int order = 1;
1726 
1727  this->initialize(A,b,c,order,this->getDescription());
1728  this->setMyParamList(pl);
1729  this->rkbtPL_ = pl;
1730  }
1731 
1732  Teuchos::RCP<const Teuchos::ParameterList>
1733  getValidParameters() const
1734  {
1735  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1736  pl->setName("Default Stepper - " + this->description());
1737  pl->set<std::string>("Description", this->getDescription());
1738  pl->set<std::string>("Stepper Type", this->description());
1739  pl->set<bool>("Use Embedded", false);
1740  pl->set<std::string>("Solver Name", "",
1741  "Name of ParameterList containing the solver specifications.");
1742 
1743  return pl;
1744  }
1745 };
1746 
1747 
1748 // ----------------------------------------------------------------------------
1749 /** \brief SDIRK 2 Stage 2nd order
1750  *
1751  * The tableau (order=1 or 2) is
1752  * \f[
1753  * \begin{array}{c|c}
1754  * c & A \\ \hline
1755  * & b^T
1756  * \end{array}
1757  * \;\;\;\;\mbox{ where }\;\;\;\;
1758  * \begin{array}{c|cc} \gamma & \gamma & \\
1759  * 1 & 1-\gamma & \gamma \\ \hline
1760  * & 1-\gamma & \gamma \end{array}
1761  * \f]
1762  * The default value is \f$\gamma = (2\pm \sqrt{2})/2\f$.
1763  * This will produce an L-stable 2nd order method with the stage
1764  * times within the timestep. Other values of gamma will still
1765  * produce an L-stable scheme, but will only be 1st order accurate.
1766  * L-stability is guaranteed because \f$A_{sj} = b_j\f$.
1767  *
1768  * Reference: U. M. Ascher and L. R. Petzold,
1769  * Computer Methods for ODEs and DAEs, p. 106.
1770  */
1771 template<class Scalar>
1772 class SDIRK2Stage2ndOrder_RKBT :
1773  virtual public RKButcherTableau<Scalar>
1774 {
1775  public:
1776  SDIRK2Stage2ndOrder_RKBT()
1777  {
1778  std::ostringstream Description;
1779  Description << this->description() << "\n"
1780  << "Computer Methods for ODEs and DAEs\n"
1781  << "U. M. Ascher and L. R. Petzold\n"
1782  << "p. 106\n"
1783  << "gamma = (2+-sqrt(2))/2\n"
1784  << "c = [ gamma 1 ]'\n"
1785  << "A = [ gamma 0 ]\n"
1786  << " [ 1-gamma gamma ]\n"
1787  << "b = [ 1-gamma gamma ]'" << std::endl;
1788 
1789  typedef Teuchos::ScalarTraits<Scalar> ST;
1790  const Scalar one = ST::one();
1791  gamma_default_ = Teuchos::as<Scalar>((2*one-ST::squareroot(2*one))/(2*one));
1792  gamma_ = gamma_default_;
1793 
1794  this->setDescription(Description.str());
1795  this->setParameterList(Teuchos::null);
1796  }
1797 
1798  virtual std::string description() const { return "SDIRK 2 Stage 2nd order"; }
1799 
1800  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
1801  {
1802  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1803  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
1804  else pl = pList;
1805  // Can not validate because optional parameters (e.g., Solver Name).
1806  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
1807  TEUCHOS_TEST_FOR_EXCEPTION(
1808  pl->get<std::string>("Stepper Type") != this->description()
1809  ,std::runtime_error,
1810  " Stepper Type != \""+this->description()+"\"\n"
1811  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
1812 
1813  gamma_ = pl->get<double>("gamma", gamma_default_);
1814 
1815  typedef Teuchos::ScalarTraits<Scalar> ST;
1816  int NumStages = 2;
1817  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1818  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1819  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1820  const Scalar one = ST::one();
1821  const Scalar zero = ST::zero();
1822  A(0,0) = gamma_;
1823  A(0,1) = zero;
1824  A(1,0) = Teuchos::as<Scalar>( one - gamma_ );
1825  A(1,1) = gamma_;
1826  b(0) = Teuchos::as<Scalar>( one - gamma_ );
1827  b(1) = gamma_;
1828  c(0) = gamma_;
1829  c(1) = one;
1830 
1831  int order = 1;
1832  if ( std::abs((gamma_-gamma_default_)/gamma_) < 1.0e-08 ) order = 2;
1833 
1834  this->initialize(A,b,c,order,1,2,this->getDescription());
1835  this->setMyParamList(pl);
1836  this->rkbtPL_ = pl;
1837  }
1838 
1839  Teuchos::RCP<const Teuchos::ParameterList>
1840  getValidParameters() const
1841  {
1842  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1843  pl->setName("Default Stepper - " + this->description());
1844  pl->set<std::string>("Description", this->getDescription());
1845  pl->set<std::string>("Stepper Type", this->description());
1846  pl->set<bool>("Use Embedded", false);
1847  pl->set("Solver Name", "",
1848  "Name of ParameterList containing the solver specifications.");
1849  pl->set<double>("gamma",gamma_default_,
1850  "The default value is gamma = (2-sqrt(2))/2. "
1851  "This will produce an L-stable 2nd order method with the stage "
1852  "times within the timestep. Other values of gamma will still "
1853  "produce an L-stable scheme, but will only be 1st order accurate.");
1854 
1855  return pl;
1856  }
1857 
1858  private:
1859  Scalar gamma_default_;
1860  Scalar gamma_;
1861 };
1862 
1863 
1864 // ----------------------------------------------------------------------------
1865 /** \brief SDIRK 2 Stage 3rd order
1866  *
1867  * The tableau (order=2 or 3) is
1868  * \f[
1869  * \begin{array}{c|c}
1870  * c & A \\ \hline
1871  * & b^T
1872  * \end{array}
1873  * \;\;\;\;\mbox{ where }\;\;\;\;
1874  * \begin{array}{c|cc} \gamma & \gamma & \\
1875  * 1-\gamma & 1-2\gamma & \gamma \\ \hline
1876  * & 1/2 & 1/2 \end{array}
1877  * \f]
1878  * \f[
1879  * \gamma = \left\{ \begin{array}{cc}
1880  * (2\pm \sqrt{2})/2 & \mbox{then 2nd order and L-stable} \\
1881  * (3\pm \sqrt{3})/6 & \mbox{then 3rd order and A-stable}
1882  * \end{array} \right.
1883  * \f]
1884  * The default value is \f$\gamma = (3\pm \sqrt{3})/6\f$.
1885  *
1886  * Reference: E. Hairer, S. P. Norsett, and G. Wanner,
1887  * Solving Ordinary Differential Equations I:
1888  * Nonstiff Problems, 2nd Revised Edition,
1889  * Table 7.2, pg 207.
1890  */
1891 template<class Scalar>
1892 class SDIRK2Stage3rdOrder_RKBT :
1893  virtual public RKButcherTableau<Scalar>
1894 {
1895  public:
1896  SDIRK2Stage3rdOrder_RKBT()
1897  {
1898  std::ostringstream Description;
1899  Description << this->description() << "\n"
1900  << "Solving Ordinary Differential Equations I:\n"
1901  << "Nonstiff Problems, 2nd Revised Edition\n"
1902  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
1903  << "Table 7.2, pg 207\n"
1904  << "gamma = (3+-sqrt(3))/6 -> 3rd order and A-stable\n"
1905  << "gamma = (2+-sqrt(2))/2 -> 2nd order and L-stable\n"
1906  << "c = [ gamma 1-gamma ]'\n"
1907  << "A = [ gamma 0 ]\n"
1908  << " [ 1-2*gamma gamma ]\n"
1909  << "b = [ 1/2 1/2 ]'" << std::endl;
1910 
1911  thirdOrderAStable_default_ = true;
1912  secondOrderLStable_default_ = false;
1913  thirdOrderAStable_ = true;
1914  secondOrderLStable_ = false;
1915  typedef Teuchos::ScalarTraits<Scalar> ST;
1916  const Scalar one = ST::one();
1917  gamma_default_ =
1918  Teuchos::as<Scalar>( (3*one + ST::squareroot(3*one))/(6*one) );
1919  gamma_ = gamma_default_;
1920 
1921  this->setDescription(Description.str());
1922  this->setParameterList(Teuchos::null);
1923  }
1924 
1925  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
1926  {
1927  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1928  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
1929  else pl = pList;
1930  // Can not validate because optional parameters (e.g., bstar).
1931  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
1932  TEUCHOS_TEST_FOR_EXCEPTION(
1933  pl->get<std::string>("Stepper Type") != this->description()
1934  ,std::runtime_error,
1935  " Stepper Type != \""+this->description()+"\"\n"
1936  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
1937 
1938  thirdOrderAStable_ = pl->get<bool>("3rd Order A-stable",false);
1939  secondOrderLStable_ = pl->get<bool>("2nd Order L-stable",false);
1940  TEUCHOS_TEST_FOR_EXCEPTION(
1941  thirdOrderAStable_ && secondOrderLStable_, std::logic_error,
1942  "'3rd Order A-stable' and '2nd Order L-stable' can not both be true.");
1943  gamma_ = pl->get<double>("gamma", gamma_default_);
1944 
1945  typedef Teuchos::ScalarTraits<Scalar> ST;
1946  using Teuchos::as;
1947  int NumStages = 2;
1948  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
1949  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
1950  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
1951  const Scalar one = ST::one();
1952  const Scalar zero = ST::zero();
1953  const Scalar gammaLStable =
1954  as<Scalar>( (2*one + ST::squareroot(2*one))/(2*one) );
1955  if (thirdOrderAStable_)
1956  gamma_ = gamma_default_;
1957  else if (secondOrderLStable_)
1958  gamma_ = gammaLStable;
1959  A(0,0) = gamma_;
1960  A(0,1) = zero;
1961  A(1,0) = as<Scalar>( one - 2*gamma_ );
1962  A(1,1) = gamma_;
1963  b(0) = as<Scalar>( one/(2*one) );
1964  b(1) = as<Scalar>( one/(2*one) );
1965  c(0) = gamma_;
1966  c(1) = as<Scalar>( one - gamma_ );
1967 
1968  int order = 2;
1969  if ( std::abs((gamma_-gamma_default_)/gamma_) < 1.0e-08 ) {
1970  order = 3;
1971  thirdOrderAStable_ = true;
1972  secondOrderLStable_ = false;
1973  } else if ( std::abs((gamma_-gammaLStable)/gamma_) < 1.0e-08 ) {
1974  thirdOrderAStable_ = false;
1975  secondOrderLStable_ = true;
1976  } else {
1977  thirdOrderAStable_ = false;
1978  secondOrderLStable_ = false;
1979  }
1980 
1981  this->initialize(A,b,c,order,2,3,this->getDescription());
1982  this->setMyParamList(pl);
1983  this->rkbtPL_ = pl;
1984  }
1985 
1986  Teuchos::RCP<const Teuchos::ParameterList>
1987  getValidParameters() const
1988  {
1989  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1990  pl->setName("Default Stepper - " + this->description());
1991  pl->set<std::string>("Description", this->getDescription());
1992  pl->set<std::string>("Stepper Type", this->description());
1993  pl->set<bool>("Use Embedded", false);
1994  pl->set("Solver Name", "",
1995  "Name of ParameterList containing the solver specifications.");
1996  pl->set<bool>("3rd Order A-stable",thirdOrderAStable_default_,
1997  "If true, set gamma to gamma = (3+sqrt(3))/6 to obtain "
1998  "a 3rd order A-stable scheme. '3rd Order A-stable' and "
1999  "'2nd Order L-stable' can not both be true.");
2000  pl->set<bool>("2nd Order L-stable",secondOrderLStable_default_,
2001  "If true, set gamma to gamma = (2+sqrt(2))/2 to obtain "
2002  "a 2nd order L-stable scheme. '3rd Order A-stable' and "
2003  "'2nd Order L-stable' can not both be true.");
2004  pl->set<double>("gamma",gamma_default_,
2005  "If both '3rd Order A-stable' and '2nd Order L-stable' "
2006  "are false, gamma will be used. The default value is the "
2007  "'3rd Order A-stable' gamma value, (3+sqrt(3))/6.");
2008 
2009  return pl;
2010  }
2011 
2012  virtual std::string description() const { return "SDIRK 2 Stage 3rd order"; }
2013 
2014  private:
2015  bool thirdOrderAStable_default_;
2016  bool thirdOrderAStable_;
2017  bool secondOrderLStable_default_;
2018  bool secondOrderLStable_;
2019  Scalar gamma_default_;
2020  Scalar gamma_;
2021 };
2022 
2023 
2024 // ----------------------------------------------------------------------------
2025 /** \brief EDIRK 2 Stage 3rd order
2026  *
2027  * The tableau (order=3) is
2028  * \f[
2029  * \begin{array}{c|c}
2030  * c & A \\ \hline
2031  * & b^T
2032  * \end{array}
2033  * \;\;\;\;\mbox{ where }\;\;\;\;
2034  * \begin{array}{c|cc} 0 & 0 & \\
2035  * 2/3 & 1/3 & 1/3 \\ \hline
2036  * & 1/4 & 3/4 \end{array}
2037  * \f]
2038  * Reference: E. Hairer, S. P. Norsett, and G. Wanner,
2039  * Solving Ordinary Differential Equations I:
2040  * Nonstiff Problems, 2nd Revised Edition,
2041  * Table 7.1, pg 205.
2042  */
2043 template<class Scalar>
2044 class EDIRK2Stage3rdOrder_RKBT :
2045  virtual public RKButcherTableau<Scalar>
2046 {
2047  public:
2048  EDIRK2Stage3rdOrder_RKBT()
2049  {
2050  std::ostringstream Description;
2051  Description << this->description() << "\n"
2052  << "Hammer & Hollingsworth method\n"
2053  << "Solving Ordinary Differential Equations I:\n"
2054  << "Nonstiff Problems, 2nd Revised Edition\n"
2055  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
2056  << "Table 7.1, pg 205\n"
2057  << "c = [ 0 2/3 ]'\n"
2058  << "A = [ 0 0 ]\n"
2059  << " [ 1/3 1/3 ]\n"
2060  << "b = [ 1/4 3/4 ]'" << std::endl;
2061 
2062  this->setDescription(Description.str());
2063  this->setParameterList(Teuchos::null);
2064  }
2065 
2066  virtual std::string description() const { return "EDIRK 2 Stage 3rd order"; }
2067 
2068  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
2069  {
2070  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
2071  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
2072  else pl = pList;
2073  // Can not validate because optional parameters (e.g., Solver Name).
2074  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
2075  TEUCHOS_TEST_FOR_EXCEPTION(
2076  pl->get<std::string>("Stepper Type") != this->description()
2077  ,std::runtime_error,
2078  " Stepper Type != \""+this->description()+"\"\n"
2079  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
2080  this->setDescription(pl->get<std::string>("Description",""));
2081 
2082  typedef Teuchos::ScalarTraits<Scalar> ST;
2083  using Teuchos::as;
2084  int NumStages = 2;
2085  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2086  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2087  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2088  const Scalar one = ST::one();
2089  const Scalar zero = ST::zero();
2090  A(0,0) = zero;
2091  A(0,1) = zero;
2092  A(1,0) = as<Scalar>( one/(3*one) );
2093  A(1,1) = as<Scalar>( one/(3*one) );
2094  b(0) = as<Scalar>( one/(4*one) );
2095  b(1) = as<Scalar>( 3*one/(4*one) );
2096  c(0) = zero;
2097  c(1) = as<Scalar>( 2*one/(3*one) );
2098  int order = 3;
2099 
2100  this->initialize(A,b,c,order,this->getDescription());
2101  this->setMyParamList(pl);
2102  this->rkbtPL_ = pl;
2103  }
2104 
2105  Teuchos::RCP<const Teuchos::ParameterList>
2106  getValidParameters() const
2107  {
2108  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
2109  pl->setName("Default Stepper - " + this->description());
2110  pl->set<std::string>("Description", this->getDescription());
2111  pl->set<std::string>("Stepper Type", this->description());
2112  pl->set<bool>("Use Embedded", false);
2113  pl->set("Solver Name", "",
2114  "Name of ParameterList containing the solver specifications.");
2115 
2116  return pl;
2117  }
2118 
2119 };
2120 
2121 
2122 // ----------------------------------------------------------------------------
2123 template<class Scalar>
2124 class Implicit3Stage6thOrderKuntzmannButcher_RKBT :
2125  virtual public RKButcherTableau<Scalar>
2126 {
2127  public:
2128  Implicit3Stage6thOrderKuntzmannButcher_RKBT()
2129  {
2130  std::ostringstream Description;
2131  Description << this->description() << "\n"
2132  << "Kuntzmann & Butcher method\n"
2133  << "Solving Ordinary Differential Equations I:\n"
2134  << "Nonstiff Problems, 2nd Revised Edition\n"
2135  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
2136  << "Table 7.4, pg 209\n"
2137  << "c = [ 1/2-sqrt(15)/10 1/2 1/2+sqrt(15)/10 ]'\n"
2138  << "A = [ 5/36 2/9-sqrt(15)/15 5/36-sqrt(15)/30 ]\n"
2139  << " [ 5/36+sqrt(15)/24 2/9 5/36-sqrt(15)/24 ]\n"
2140  << " [ 5/36+sqrt(15)/30 2/9+sqrt(15)/15 5/36 ]\n"
2141  << "b = [ 5/18 4/9 5/18 ]'"
2142  << std::endl;
2143  typedef Teuchos::ScalarTraits<Scalar> ST;
2144  using Teuchos::as;
2145  int NumStages = 3;
2146  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2147  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2148  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2149  const Scalar one = ST::one();
2150  A(0,0) = as<Scalar>( 5*one/(36*one) );
2151  A(0,1) = as<Scalar>( 2*one/(9*one) - ST::squareroot(15*one)/(15*one) );
2152  A(0,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(30*one) );
2153  A(1,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(24*one) );
2154  A(1,1) = as<Scalar>( 2*one/(9*one) );
2155  A(1,2) = as<Scalar>( 5*one/(36*one) - ST::squareroot(15*one)/(24*one) );
2156  A(2,0) = as<Scalar>( 5*one/(36*one) + ST::squareroot(15*one)/(30*one) );
2157  A(2,1) = as<Scalar>( 2*one/(9*one) + ST::squareroot(15*one)/(15*one) );
2158  A(2,2) = as<Scalar>( 5*one/(36*one) );
2159  b(0) = as<Scalar>( 5*one/(18*one) );
2160  b(1) = as<Scalar>( 4*one/(9*one) );
2161  b(2) = as<Scalar>( 5*one/(18*one) );
2162  c(0) = as<Scalar>( one/(2*one)-ST::squareroot(15*one)/(10*one) );
2163  c(1) = as<Scalar>( one/(2*one) );
2164  c(2) = as<Scalar>( one/(2*one)+ST::squareroot(15*one)/(10*one) );
2165  int order = 6;
2166 
2167  this->initialize(A,b,c,order,Description.str());
2168  }
2169  virtual std::string description() const
2170  { return "RK Implicit 3 Stage 6th Order Kuntzmann & Butcher"; }
2171 };
2172 
2173 
2174 // ----------------------------------------------------------------------------
2175 template<class Scalar>
2176 class Implicit4Stage8thOrderKuntzmannButcher_RKBT :
2177  virtual public RKButcherTableau<Scalar>
2178 {
2179  public:
2180  Implicit4Stage8thOrderKuntzmannButcher_RKBT()
2181  {
2182  std::ostringstream Description;
2183  Description << this->description() << "\n"
2184  << "Kuntzmann & Butcher method\n"
2185  << "Solving Ordinary Differential Equations I:\n"
2186  << "Nonstiff Problems, 2nd Revised Edition\n"
2187  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
2188  << "Table 7.5, pg 209\n"
2189  << "c = [ 1/2-w2 1/2-w2p 1/2+w2p 1/2+w2 ]'\n"
2190  << "A = [ w1 w1p-w3+w4p w1p-w3-w4p w1-w5 ]\n"
2191  << " [ w1-w3p+w4 w1p w1p-w5p w1-w3p-w4 ]\n"
2192  << " [ w1+w3p+w4 w1p+w5p w1p w1+w3p-w4 ]\n"
2193  << " [ w1+w5 w1p+w3+w4p w1p+w3-w4p w1 ]\n"
2194  << "b = [ 2*w1 2*w1p 2*w1p 2*w1 ]'\n"
2195  << "w1 = 1/8-sqrt(30)/144\n"
2196  << "w2 = (1/2)*sqrt((15+2*sqrt(3))/35)\n"
2197  << "w3 = w2*(1/6+sqrt(30)/24)\n"
2198  << "w4 = w2*(1/21+5*sqrt(30)/168)\n"
2199  << "w5 = w2-2*w3\n"
2200  << "w1p = 1/8+sqrt(30)/144\n"
2201  << "w2p = (1/2)*sqrt((15-2*sqrt(3))/35)\n"
2202  << "w3p = w2*(1/6-sqrt(30)/24)\n"
2203  << "w4p = w2*(1/21-5*sqrt(30)/168)\n"
2204  << "w5p = w2p-2*w3p" << std::endl;
2205  typedef Teuchos::ScalarTraits<Scalar> ST;
2206  using Teuchos::as;
2207  int NumStages = 4;
2208  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2209  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2210  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2211  const Scalar one = ST::one();
2212  const Scalar onehalf = as<Scalar>( one/(2*one) );
2213  const Scalar w1 = as<Scalar>( one/(8*one) - ST::squareroot(30*one)/(144*one) );
2214  const Scalar w2 = as<Scalar>( (one/(2*one))*ST::squareroot((15*one+2*one*ST::squareroot(30*one))/(35*one)) );
2215  const Scalar w3 = as<Scalar>( w2*(one/(6*one)+ST::squareroot(30*one)/(24*one)) );
2216  const Scalar w4 = as<Scalar>( w2*(one/(21*one)+5*one*ST::squareroot(30*one)/(168*one)) );
2217  const Scalar w5 = as<Scalar>( w2-2*w3 );
2218  const Scalar w1p = as<Scalar>( one/(8*one) + ST::squareroot(30*one)/(144*one) );
2219  const Scalar w2p = as<Scalar>( (one/(2*one))*ST::squareroot((15*one-2*one*ST::squareroot(30*one))/(35*one)) );
2220  const Scalar w3p = as<Scalar>( w2p*(one/(6*one)-ST::squareroot(30*one)/(24*one)) );
2221  const Scalar w4p = as<Scalar>( w2p*(one/(21*one)-5*one*ST::squareroot(30*one)/(168*one)) );
2222  const Scalar w5p = as<Scalar>( w2p-2*w3p );
2223  A(0,0) = w1;
2224  A(0,1) = w1p-w3+w4p;
2225  A(0,2) = w1p-w3-w4p;
2226  A(0,3) = w1-w5;
2227  A(1,0) = w1-w3p+w4;
2228  A(1,1) = w1p;
2229  A(1,2) = w1p-w5p;
2230  A(1,3) = w1-w3p-w4;
2231  A(2,0) = w1+w3p+w4;
2232  A(2,1) = w1p+w5p;
2233  A(2,2) = w1p;
2234  A(2,3) = w1+w3p-w4;
2235  A(3,0) = w1+w5;
2236  A(3,1) = w1p+w3+w4p;
2237  A(3,2) = w1p+w3-w4p;
2238  A(3,3) = w1;
2239  b(0) = 2*w1;
2240  b(1) = 2*w1p;
2241  b(2) = 2*w1p;
2242  b(3) = 2*w1;
2243  c(0) = onehalf - w2;
2244  c(1) = onehalf - w2p;
2245  c(2) = onehalf + w2p;
2246  c(3) = onehalf + w2;
2247  int order = 8;
2248 
2249  this->initialize(A,b,c,order,Description.str());
2250  }
2251  virtual std::string description() const
2252  { return "RK Implicit 4 Stage 8th Order Kuntzmann & Butcher"; }
2253 };
2254 
2255 
2256 // ----------------------------------------------------------------------------
2257 template<class Scalar>
2258 class Implicit2Stage4thOrderHammerHollingsworth_RKBT :
2259  virtual public RKButcherTableau<Scalar>
2260 {
2261  public:
2262  Implicit2Stage4thOrderHammerHollingsworth_RKBT()
2263  {
2264  std::ostringstream Description;
2265  Description << this->description() << "\n"
2266  << "Hammer & Hollingsworth method\n"
2267  << "Solving Ordinary Differential Equations I:\n"
2268  << "Nonstiff Problems, 2nd Revised Edition\n"
2269  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
2270  << "Table 7.3, pg 207\n"
2271  << "c = [ 1/2-sqrt(3)/6 1/2+sqrt(3)/6 ]'\n"
2272  << "A = [ 1/4 1/4-sqrt(3)/6 ]\n"
2273  << " [ 1/4+sqrt(3)/6 1/4 ]\n"
2274  << "b = [ 1/2 1/2 ]'" << std::endl;
2275  typedef Teuchos::ScalarTraits<Scalar> ST;
2276  using Teuchos::as;
2277  int NumStages = 2;
2278  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2279  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2280  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2281  const Scalar one = ST::one();
2282  const Scalar onequarter = as<Scalar>( one/(4*one) );
2283  const Scalar onehalf = as<Scalar>( one/(2*one) );
2284  A(0,0) = onequarter;
2285  A(0,1) = as<Scalar>( onequarter-ST::squareroot(3*one)/(6*one) );
2286  A(1,0) = as<Scalar>( onequarter+ST::squareroot(3*one)/(6*one) );
2287  A(1,1) = onequarter;
2288  b(0) = onehalf;
2289  b(1) = onehalf;
2290  c(0) = as<Scalar>( onehalf - ST::squareroot(3*one)/(6*one) );
2291  c(1) = as<Scalar>( onehalf + ST::squareroot(3*one)/(6*one) );
2292  int order = 4;
2293 
2294  this->initialize(A,b,c,order,Description.str());
2295  }
2296  virtual std::string description() const
2297  { return "RK Implicit 2 Stage 4th Order Hammer & Hollingsworth"; }
2298 };
2299 
2300 
2301 // ----------------------------------------------------------------------------
2302 template<class Scalar>
2303 class IRK1StageTheta_RKBT :
2304  virtual public RKButcherTableau<Scalar>
2305 {
2306  public:
2307  IRK1StageTheta_RKBT()
2308  {
2309  std::ostringstream Description;
2310  Description << this->description() << "\n"
2311  << "Non-standard finite-difference methods\n"
2312  << "in dynamical systems, P. Kama,\n"
2313  << "Dissertation, University of Pretoria, pg. 49.\n"
2314  << "Comment: Generalized Implicit Midpoint Method\n"
2315  << "c = [ theta ]'\n"
2316  << "A = [ theta ]\n"
2317  << "b = [ 1 ]'\n"
2318  << std::endl;
2319 
2320  typedef Teuchos::ScalarTraits<Scalar> ST;
2321  theta_default_ = ST::one()/(2*ST::one());
2322  theta_ = theta_default_;
2323 
2324  this->setDescription(Description.str());
2325  this->setParameterList(Teuchos::null);
2326  }
2327 
2328  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
2329  {
2330  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
2331  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
2332  else pl = pList;
2333  // Can not validate because optional parameters (e.g., Solver Name).
2334  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
2335  TEUCHOS_TEST_FOR_EXCEPTION(
2336  pl->get<std::string>("Stepper Type") != this->description() and
2337  pl->get<std::string>("Stepper Type") != "Implicit Midpoint"
2338  ,std::runtime_error,
2339  " Stepper Type != \""+this->description()+"\"\n"
2340  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
2341  theta_ = pl->get<double>("theta",theta_default_);
2342 
2343  typedef Teuchos::ScalarTraits<Scalar> ST;
2344  int NumStages = 1;
2345  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2346  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2347  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2348  A(0,0) = theta_;
2349  b(0) = ST::one();
2350  c(0) = theta_;
2351 
2352  int order = 1;
2353  if (theta_ == theta_default_) order = 2;
2354 
2355  this->initialize(A, b, c, order, 1, 2, this->getDescription());
2356  this->setMyParamList(pl);
2357  this->rkbtPL_ = pl;
2358  }
2359 
2360  virtual std::string description() const { return "IRK 1 Stage Theta Method"; }
2361 
2362  Teuchos::RCP<const Teuchos::ParameterList>
2363  getValidParameters() const
2364  {
2365  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
2366  pl->setName("Default Stepper - " + this->description());
2367  pl->set<std::string>("Description", this->getDescription());
2368  pl->set<std::string>("Stepper Type", this->description());
2369  pl->set<bool>("Use Embedded", false);
2370  pl->set<std::string>("Solver Name", "",
2371  "Name of ParameterList containing the solver specifications.");
2372  pl->set<double>("theta",theta_default_,
2373  "Valid values are 0 <= theta <= 1, where theta = 0 "
2374  "implies Forward Euler, theta = 1/2 implies implicit midpoint "
2375  "method (default), and theta = 1 implies Backward Euler. "
2376  "For theta != 1/2, this method is first-order accurate, "
2377  "and with theta = 1/2, it is second-order accurate. "
2378  "This method is A-stable, but becomes L-stable with theta=1.");
2379 
2380  return pl;
2381  }
2382 
2383  private:
2384  Scalar theta_default_;
2385  Scalar theta_;
2386 };
2387 
2388 
2389 // ----------------------------------------------------------------------------
2390 template<class Scalar>
2391 class EDIRK2StageTheta_RKBT :
2392  virtual public RKButcherTableau<Scalar>
2393 {
2394  public:
2395  EDIRK2StageTheta_RKBT()
2396  {
2397  std::ostringstream Description;
2398  Description << this->description() << "\n"
2399  << "Computer Methods for ODEs and DAEs\n"
2400  << "U. M. Ascher and L. R. Petzold\n"
2401  << "p. 113\n"
2402  << "c = [ 0 1 ]'\n"
2403  << "A = [ 0 0 ]\n"
2404  << " [ 1-theta theta ]\n"
2405  << "b = [ 1-theta theta ]'\n"
2406  << std::endl;
2407 
2408  typedef Teuchos::ScalarTraits<Scalar> ST;
2409  theta_default_ = ST::one()/(2*ST::one());
2410  theta_ = theta_default_;
2411 
2412  this->setDescription(Description.str());
2413  this->setParameterList(Teuchos::null);
2414  }
2415 
2416  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
2417  {
2418  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
2419  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
2420  else pl = pList;
2421  // Can not validate because optional parameters (e.g., Solver Name).
2422  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
2423  TEUCHOS_TEST_FOR_EXCEPTION(
2424  pl->get<std::string>("Stepper Type") != this->description()
2425  ,std::runtime_error,
2426  " Stepper Type != \""+this->description()+"\"\n"
2427  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
2428  theta_ = pl->get<double>("theta",theta_default_);
2429 
2430  typedef Teuchos::ScalarTraits<Scalar> ST;
2431  const Scalar one = ST::one();
2432  const Scalar zero = ST::zero();
2433  TEUCHOS_TEST_FOR_EXCEPTION(
2434  theta_ == zero, std::logic_error,
2435  "'theta' can not be zero, as it makes this IRK stepper explicit.");
2436 
2437  int NumStages = 2;
2438  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2439  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2440  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2441  A(0,0) = zero;
2442  A(0,1) = zero;
2443  A(1,0) = Teuchos::as<Scalar>( one - theta_ );
2444  A(1,1) = theta_;
2445  b(0) = Teuchos::as<Scalar>( one - theta_ );
2446  b(1) = theta_;
2447  c(0) = theta_;
2448  c(1) = one;
2449 
2450  int order = 1;
2451  if (theta_ == theta_default_) order = 2;
2452 
2453  this->initialize(A, b, c, order, 1, 2, this->getDescription());
2454  this->setMyParamList(pl);
2455  this->rkbtPL_ = pl;
2456  }
2457 
2458  Teuchos::RCP<const Teuchos::ParameterList>
2459  getValidParameters() const
2460  {
2461  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
2462  pl->setName("Default Stepper - " + this->description());
2463  pl->set<std::string>("Description", this->getDescription());
2464  pl->set<std::string>("Stepper Type", this->description());
2465  pl->set<bool>("Use Embedded", false);
2466  pl->set<std::string>("Solver Name", "",
2467  "Name of ParameterList containing the solver specifications.");
2468  pl->set<double>("theta",theta_default_,
2469  "Valid values are 0 < theta <= 1, where theta = 0 "
2470  "implies Forward Euler, theta = 1/2 implies trapezoidal "
2471  "method (default), and theta = 1 implies Backward Euler. "
2472  "For theta != 1/2, this method is first-order accurate, "
2473  "and with theta = 1/2, it is second-order accurate. "
2474  "This method is A-stable, but becomes L-stable with theta=1.");
2475 
2476  return pl;
2477  }
2478 
2479  virtual std::string description() const {return "EDIRK 2 Stage Theta Method";}
2480 
2481  private:
2482  Scalar theta_default_;
2483  Scalar theta_;
2484 };
2485 
2486 
2487 // ----------------------------------------------------------------------------
2488 template<class Scalar>
2489 class Implicit1Stage2ndOrderGauss_RKBT :
2490  virtual public RKButcherTableau<Scalar>
2491 {
2492  public:
2493  Implicit1Stage2ndOrderGauss_RKBT()
2494  {
2495  std::ostringstream Description;
2496  Description << this->description() << "\n"
2497  << "A-stable\n"
2498  << "Solving Ordinary Differential Equations II:\n"
2499  << "Stiff and Differential-Algebraic Problems,\n"
2500  << "2nd Revised Edition\n"
2501  << "E. Hairer and G. Wanner\n"
2502  << "Table 5.2, pg 72\n"
2503  << "Also: Implicit midpoint rule\n"
2504  << "Solving Ordinary Differential Equations I:\n"
2505  << "Nonstiff Problems, 2nd Revised Edition\n"
2506  << "E. Hairer, S. P. Norsett, and G. Wanner\n"
2507  << "Table 7.1, pg 205\n"
2508  << "c = [ 1/2 ]'\n"
2509  << "A = [ 1/2 ]\n"
2510  << "b = [ 1 ]'" << std::endl;
2511  typedef Teuchos::ScalarTraits<Scalar> ST;
2512  int NumStages = 1;
2513  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2514  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2515  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2516  const Scalar onehalf = ST::one()/(2*ST::one());
2517  const Scalar one = ST::one();
2518  A(0,0) = onehalf;
2519  b(0) = one;
2520  c(0) = onehalf;
2521  int order = 2;
2522 
2523  this->initialize(A,b,c,order,Description.str());
2524  }
2525  virtual std::string description() const
2526  { return "RK Implicit 1 Stage 2nd order Gauss"; }
2527 };
2528 
2529 
2530 // ----------------------------------------------------------------------------
2531 template<class Scalar>
2532 class Implicit2Stage4thOrderGauss_RKBT :
2533  virtual public RKButcherTableau<Scalar>
2534 {
2535  public:
2536  Implicit2Stage4thOrderGauss_RKBT()
2537  {
2538  std::ostringstream Description;
2539  Description << this->description() << "\n"
2540  << "A-stable\n"
2541  << "Solving Ordinary Differential Equations II:\n"
2542  << "Stiff and Differential-Algebraic Problems,\n"
2543  << "2nd Revised Edition\n"
2544  << "E. Hairer and G. Wanner\n"
2545  << "Table 5.2, pg 72\n"
2546  << "c = [ 1/2-sqrt(3)/6 1/2+sqrt(3)/6 ]'\n"
2547  << "A = [ 1/4 1/4-sqrt(3)/6 ]\n"
2548  << " [ 1/4+sqrt(3)/6 1/4 ]\n"
2549  << "b = [ 1/2 1/2 ]'" << std::endl;
2550  typedef Teuchos::ScalarTraits<Scalar> ST;
2551  using Teuchos::as;
2552  int NumStages = 2;
2553  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2554  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2555  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2556  const Scalar one = ST::one();
2557  const Scalar onehalf = as<Scalar>(one/(2*one));
2558  const Scalar three = as<Scalar>(3*one);
2559  const Scalar six = as<Scalar>(6*one);
2560  const Scalar onefourth = as<Scalar>(one/(4*one));
2561  const Scalar alpha = ST::squareroot(three)/six;
2562 
2563  A(0,0) = onefourth;
2564  A(0,1) = onefourth-alpha;
2565  A(1,0) = onefourth+alpha;
2566  A(1,1) = onefourth;
2567  b(0) = onehalf;
2568  b(1) = onehalf;
2569  c(0) = onehalf-alpha;
2570  c(1) = onehalf+alpha;
2571  int order = 4;
2572 
2573  this->initialize(A,b,c,order,Description.str());
2574  }
2575  virtual std::string description() const
2576  { return "RK Implicit 2 Stage 4th order Gauss"; }
2577 };
2578 
2579 
2580 // ----------------------------------------------------------------------------
2581 template<class Scalar>
2582 class Implicit3Stage6thOrderGauss_RKBT :
2583  virtual public RKButcherTableau<Scalar>
2584 {
2585  public:
2586  Implicit3Stage6thOrderGauss_RKBT()
2587  {
2588  std::ostringstream Description;
2589  Description << this->description() << "\n"
2590  << "A-stable\n"
2591  << "Solving Ordinary Differential Equations II:\n"
2592  << "Stiff and Differential-Algebraic Problems,\n"
2593  << "2nd Revised Edition\n"
2594  << "E. Hairer and G. Wanner\n"
2595  << "Table 5.2, pg 72\n"
2596  << "c = [ 1/2-sqrt(15)/10 1/2 1/2+sqrt(15)/10 ]'\n"
2597  << "A = [ 5/36 2/9-sqrt(15)/15 5/36-sqrt(15)/30 ]\n"
2598  << " [ 5/36+sqrt(15)/24 2/9 5/36-sqrt(15)/24 ]\n"
2599  << " [ 5/36+sqrt(15)/30 2/9+sqrt(15)/15 5/36 ]\n"
2600  << "b = [ 5/18 4/9 5/18 ]'"
2601  << std::endl;
2602  typedef Teuchos::ScalarTraits<Scalar> ST;
2603  using Teuchos::as;
2604  int NumStages = 3;
2605  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2606  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2607  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2608  const Scalar one = ST::one();
2609  const Scalar ten = as<Scalar>(10*one);
2610  const Scalar fifteen = as<Scalar>(15*one);
2611  const Scalar twentyfour = as<Scalar>(24*one);
2612  const Scalar thirty = as<Scalar>(30*one);
2613  const Scalar sqrt15over10 = as<Scalar>(ST::squareroot(fifteen)/ten);
2614  const Scalar sqrt15over15 = as<Scalar>(ST::squareroot(fifteen)/fifteen);
2615  const Scalar sqrt15over24 = as<Scalar>(ST::squareroot(fifteen)/twentyfour);
2616  const Scalar sqrt15over30 = as<Scalar>(ST::squareroot(fifteen)/thirty);
2617 
2618  A(0,0) = as<Scalar>(5*one/(36*one));
2619  A(0,1) = as<Scalar>(2*one/(9*one))-sqrt15over15;
2620  A(0,2) = as<Scalar>(5*one/(36*one))-sqrt15over30;
2621  A(1,0) = as<Scalar>(5*one/(36*one))+sqrt15over24;
2622  A(1,1) = as<Scalar>(2*one/(9*one));
2623  A(1,2) = as<Scalar>(5*one/(36*one))-sqrt15over24;
2624  A(2,0) = as<Scalar>(5*one/(36*one))+sqrt15over30;
2625  A(2,1) = as<Scalar>(2*one/(9*one))+sqrt15over15;
2626  A(2,2) = as<Scalar>(5*one/(36*one));
2627  b(0) = as<Scalar>(5*one/(18*one));
2628  b(1) = as<Scalar>(4*one/(9*one));
2629  b(2) = as<Scalar>(5*one/(18*one));
2630  c(0) = as<Scalar>(one/(2*one))-sqrt15over10;
2631  c(1) = as<Scalar>(one/(2*one));
2632  c(2) = as<Scalar>(one/(2*one))+sqrt15over10;
2633  int order = 6;
2634 
2635  this->initialize(A,b,c,order,Description.str());
2636  }
2637  virtual std::string description() const
2638  { return "RK Implicit 3 Stage 6th order Gauss"; }
2639 };
2640 
2641 
2642 // ----------------------------------------------------------------------------
2643 template<class Scalar>
2644 class Implicit1Stage1stOrderRadauA_RKBT :
2645  virtual public RKButcherTableau<Scalar>
2646 {
2647  public:
2648  Implicit1Stage1stOrderRadauA_RKBT()
2649  {
2650  std::ostringstream Description;
2651  Description << this->description() << "\n"
2652  << "A-stable\n"
2653  << "Solving Ordinary Differential Equations II:\n"
2654  << "Stiff and Differential-Algebraic Problems,\n"
2655  << "2nd Revised Edition\n"
2656  << "E. Hairer and G. Wanner\n"
2657  << "Table 5.3, pg 73\n"
2658  << "c = [ 0 ]'\n"
2659  << "A = [ 1 ]\n"
2660  << "b = [ 1 ]'" << std::endl;
2661  typedef Teuchos::ScalarTraits<Scalar> ST;
2662  int NumStages = 1;
2663  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2664  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2665  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2666  const Scalar one = ST::one();
2667  const Scalar zero = ST::zero();
2668  A(0,0) = one;
2669  b(0) = one;
2670  c(0) = zero;
2671  int order = 1;
2672 
2673  this->initialize(A,b,c,order,Description.str());
2674  }
2675  virtual std::string description() const
2676  { return "RK Implicit 1 Stage 1st order Radau left"; }
2677 };
2678 
2679 
2680 // ----------------------------------------------------------------------------
2681 template<class Scalar>
2682 class Implicit2Stage3rdOrderRadauA_RKBT :
2683  virtual public RKButcherTableau<Scalar>
2684 {
2685  public:
2686  Implicit2Stage3rdOrderRadauA_RKBT()
2687  {
2688  std::ostringstream Description;
2689  Description << this->description() << "\n"
2690  << "A-stable\n"
2691  << "Solving Ordinary Differential Equations II:\n"
2692  << "Stiff and Differential-Algebraic Problems,\n"
2693  << "2nd Revised Edition\n"
2694  << "E. Hairer and G. Wanner\n"
2695  << "Table 5.3, pg 73\n"
2696  << "c = [ 0 2/3 ]'\n"
2697  << "A = [ 1/4 -1/4 ]\n"
2698  << " [ 1/4 5/12 ]\n"
2699  << "b = [ 1/4 3/4 ]'" << std::endl;
2700  typedef Teuchos::ScalarTraits<Scalar> ST;
2701  using Teuchos::as;
2702  int NumStages = 2;
2703  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2704  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2705  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2706  const Scalar zero = ST::zero();
2707  const Scalar one = ST::one();
2708  A(0,0) = as<Scalar>(one/(4*one));
2709  A(0,1) = as<Scalar>(-one/(4*one));
2710  A(1,0) = as<Scalar>(one/(4*one));
2711  A(1,1) = as<Scalar>(5*one/(12*one));
2712  b(0) = as<Scalar>(one/(4*one));
2713  b(1) = as<Scalar>(3*one/(4*one));
2714  c(0) = zero;
2715  c(1) = as<Scalar>(2*one/(3*one));
2716  int order = 3;
2717 
2718  this->initialize(A,b,c,order,Description.str());
2719  }
2720  virtual std::string description() const
2721  { return "RK Implicit 2 Stage 3rd order Radau left"; }
2722 };
2723 
2724 
2725 // ----------------------------------------------------------------------------
2726 template<class Scalar>
2727 class Implicit3Stage5thOrderRadauA_RKBT :
2728  virtual public RKButcherTableau<Scalar>
2729 {
2730  public:
2731  Implicit3Stage5thOrderRadauA_RKBT()
2732  {
2733  std::ostringstream Description;
2734  Description << this->description() << "\n"
2735  << "A-stable\n"
2736  << "Solving Ordinary Differential Equations II:\n"
2737  << "Stiff and Differential-Algebraic Problems,\n"
2738  << "2nd Revised Edition\n"
2739  << "E. Hairer and G. Wanner\n"
2740  << "Table 5.4, pg 73\n"
2741  << "c = [ 0 (6-sqrt(6))/10 (6+sqrt(6))/10 ]'\n"
2742  << "A = [ 1/9 (-1-sqrt(6))/18 (-1+sqrt(6))/18 ]\n"
2743  << " [ 1/9 (88+7*sqrt(6))/360 (88-43*sqrt(6))/360 ]\n"
2744  << " [ 1/9 (88+43*sqrt(6))/360 (88-7*sqrt(6))/360 ]\n"
2745  << "b = [ 1/9 (16+sqrt(6))/36 (16-sqrt(6))/36 ]'"
2746  << std::endl;
2747  typedef Teuchos::ScalarTraits<Scalar> ST;
2748  using Teuchos::as;
2749  int NumStages = 3;
2750  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2751  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2752  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2753  const Scalar zero = ST::zero();
2754  const Scalar one = ST::one();
2755  A(0,0) = as<Scalar>(one/(9*one));
2756  A(0,1) = as<Scalar>( (-one-ST::squareroot(6*one))/(18*one) );
2757  A(0,2) = as<Scalar>( (-one+ST::squareroot(6*one))/(18*one) );
2758  A(1,0) = as<Scalar>(one/(9*one));
2759  A(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
2760  A(1,2) = as<Scalar>( (88*one-43*one*ST::squareroot(6*one))/(360*one) );
2761  A(2,0) = as<Scalar>(one/(9*one));
2762  A(2,1) = as<Scalar>( (88*one+43*one*ST::squareroot(6*one))/(360*one) );
2763  A(2,2) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
2764  b(0) = as<Scalar>(one/(9*one));
2765  b(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
2766  b(2) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
2767  c(0) = zero;
2768  c(1) = as<Scalar>( (6*one-ST::squareroot(6*one))/(10*one) );
2769  c(2) = as<Scalar>( (6*one+ST::squareroot(6*one))/(10*one) );
2770  int order = 5;
2771 
2772  this->initialize(A,b,c,order,Description.str());
2773  }
2774  virtual std::string description() const
2775  { return "RK Implicit 3 Stage 5th order Radau left"; }
2776 };
2777 
2778 
2779 // ----------------------------------------------------------------------------
2780 template<class Scalar>
2781 class Implicit1Stage1stOrderRadauB_RKBT :
2782  virtual public RKButcherTableau<Scalar>
2783 {
2784  public:
2785  Implicit1Stage1stOrderRadauB_RKBT()
2786  {
2787  std::ostringstream Description;
2788  Description << this->description() << "\n"
2789  << "A-stable\n"
2790  << "Solving Ordinary Differential Equations II:\n"
2791  << "Stiff and Differential-Algebraic Problems,\n"
2792  << "2nd Revised Edition\n"
2793  << "E. Hairer and G. Wanner\n"
2794  << "Table 5.5, pg 74\n"
2795  << "c = [ 1 ]'\n"
2796  << "A = [ 1 ]\n"
2797  << "b = [ 1 ]'" << std::endl;
2798  typedef Teuchos::ScalarTraits<Scalar> ST;
2799  int NumStages = 1;
2800  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2801  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2802  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2803  const Scalar one = ST::one();
2804  A(0,0) = one;
2805  b(0) = one;
2806  c(0) = one;
2807  int order = 1;
2808 
2809  this->initialize(A,b,c,order,Description.str());
2810  }
2811  virtual std::string description() const
2812  { return "RK Implicit 1 Stage 1st order Radau right"; }
2813 };
2814 
2815 
2816 // ----------------------------------------------------------------------------
2817 template<class Scalar>
2818 class Implicit2Stage3rdOrderRadauB_RKBT :
2819  virtual public RKButcherTableau<Scalar>
2820 {
2821  public:
2822  Implicit2Stage3rdOrderRadauB_RKBT()
2823  {
2824  std::ostringstream Description;
2825  Description << this->description() << "\n"
2826  << "A-stable\n"
2827  << "Solving Ordinary Differential Equations II:\n"
2828  << "Stiff and Differential-Algebraic Problems,\n"
2829  << "2nd Revised Edition\n"
2830  << "E. Hairer and G. Wanner\n"
2831  << "Table 5.5, pg 74\n"
2832  << "c = [ 1/3 1 ]'\n"
2833  << "A = [ 5/12 -1/12 ]\n"
2834  << " [ 3/4 1/4 ]\n"
2835  << "b = [ 3/4 1/4 ]'" << std::endl;
2836  typedef Teuchos::ScalarTraits<Scalar> ST;
2837  using Teuchos::as;
2838  int NumStages = 2;
2839  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2840  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2841  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2842  const Scalar one = ST::one();
2843  A(0,0) = as<Scalar>( 5*one/(12*one) );
2844  A(0,1) = as<Scalar>( -one/(12*one) );
2845  A(1,0) = as<Scalar>( 3*one/(4*one) );
2846  A(1,1) = as<Scalar>( one/(4*one) );
2847  b(0) = as<Scalar>( 3*one/(4*one) );
2848  b(1) = as<Scalar>( one/(4*one) );
2849  c(0) = as<Scalar>( one/(3*one) );
2850  c(1) = one;
2851  int order = 3;
2852 
2853  this->initialize(A,b,c,order,Description.str());
2854  }
2855  virtual std::string description() const
2856  { return "RK Implicit 2 Stage 3rd order Radau right"; }
2857 };
2858 
2859 
2860 // ----------------------------------------------------------------------------
2861 template<class Scalar>
2862 class Implicit3Stage5thOrderRadauB_RKBT :
2863  virtual public RKButcherTableau<Scalar>
2864 {
2865  public:
2866  Implicit3Stage5thOrderRadauB_RKBT()
2867  {
2868  std::ostringstream Description;
2869  Description << this->description() << "\n"
2870  << "A-stable\n"
2871  << "Solving Ordinary Differential Equations II:\n"
2872  << "Stiff and Differential-Algebraic Problems,\n"
2873  << "2nd Revised Edition\n"
2874  << "E. Hairer and G. Wanner\n"
2875  << "Table 5.6, pg 74\n"
2876  << "c = [ (4-sqrt(6))/10 (4+sqrt(6))/10 1 ]'\n"
2877  << "A = [ A1 A2 A3 ]\n"
2878  << " A1 = [ (88-7*sqrt(6))/360 ]\n"
2879  << " [ (296+169*sqrt(6))/1800 ]\n"
2880  << " [ (16-sqrt(6))/36 ]\n"
2881  << " A2 = [ (296-169*sqrt(6))/1800 ]\n"
2882  << " [ (88+7*sqrt(6))/360 ]\n"
2883  << " [ (16+sqrt(6))/36 ]\n"
2884  << " A3 = [ (-2+3*sqrt(6))/225 ]\n"
2885  << " [ (-2-3*sqrt(6))/225 ]\n"
2886  << " [ 1/9 ]\n"
2887  << "b = [ (16-sqrt(6))/36 (16+sqrt(6))/36 1/9 ]'"
2888  << std::endl;
2889  typedef Teuchos::ScalarTraits<Scalar> ST;
2890  using Teuchos::as;
2891  int NumStages = 3;
2892  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2893  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2894  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2895  const Scalar one = ST::one();
2896  A(0,0) = as<Scalar>( (88*one-7*one*ST::squareroot(6*one))/(360*one) );
2897  A(0,1) = as<Scalar>( (296*one-169*one*ST::squareroot(6*one))/(1800*one) );
2898  A(0,2) = as<Scalar>( (-2*one+3*one*ST::squareroot(6*one))/(225*one) );
2899  A(1,0) = as<Scalar>( (296*one+169*one*ST::squareroot(6*one))/(1800*one) );
2900  A(1,1) = as<Scalar>( (88*one+7*one*ST::squareroot(6*one))/(360*one) );
2901  A(1,2) = as<Scalar>( (-2*one-3*one*ST::squareroot(6*one))/(225*one) );
2902  A(2,0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
2903  A(2,1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
2904  A(2,2) = as<Scalar>( one/(9*one) );
2905  b(0) = as<Scalar>( (16*one-ST::squareroot(6*one))/(36*one) );
2906  b(1) = as<Scalar>( (16*one+ST::squareroot(6*one))/(36*one) );
2907  b(2) = as<Scalar>( one/(9*one) );
2908  c(0) = as<Scalar>( (4*one-ST::squareroot(6*one))/(10*one) );
2909  c(1) = as<Scalar>( (4*one+ST::squareroot(6*one))/(10*one) );
2910  c(2) = one;
2911  int order = 5;
2912 
2913  this->initialize(A,b,c,order,Description.str());
2914  }
2915  virtual std::string description() const
2916  { return "RK Implicit 3 Stage 5th order Radau right"; }
2917 };
2918 
2919 
2920 // ----------------------------------------------------------------------------
2921 template<class Scalar>
2922 class Implicit2Stage2ndOrderLobattoA_RKBT :
2923  virtual public RKButcherTableau<Scalar>
2924 {
2925  public:
2926  Implicit2Stage2ndOrderLobattoA_RKBT()
2927  {
2928  std::ostringstream Description;
2929  Description << this->description() << "\n"
2930  << "A-stable\n"
2931  << "Solving Ordinary Differential Equations II:\n"
2932  << "Stiff and Differential-Algebraic Problems,\n"
2933  << "2nd Revised Edition\n"
2934  << "E. Hairer and G. Wanner\n"
2935  << "Table 5.7, pg 75\n"
2936  << "c = [ 0 1 ]'\n"
2937  << "A = [ 0 0 ]\n"
2938  << " [ 1/2 1/2 ]\n"
2939  << "b = [ 1/2 1/2 ]'" << std::endl;
2940  typedef Teuchos::ScalarTraits<Scalar> ST;
2941  using Teuchos::as;
2942  int NumStages = 2;
2943  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2944  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2945  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2946  const Scalar zero = ST::zero();
2947  const Scalar one = ST::one();
2948  A(0,0) = zero;
2949  A(0,1) = zero;
2950  A(1,0) = as<Scalar>( one/(2*one) );
2951  A(1,1) = as<Scalar>( one/(2*one) );
2952  b(0) = as<Scalar>( one/(2*one) );
2953  b(1) = as<Scalar>( one/(2*one) );
2954  c(0) = zero;
2955  c(1) = one;
2956  int order = 2;
2957 
2958  this->initialize(A,b,c,order,Description.str());
2959  }
2960  virtual std::string description() const
2961  { return "RK Implicit 2 Stage 2nd order Lobatto A"; }
2962 };
2963 
2964 
2965 // ----------------------------------------------------------------------------
2966 template<class Scalar>
2967 class Implicit3Stage4thOrderLobattoA_RKBT :
2968  virtual public RKButcherTableau<Scalar>
2969 {
2970  public:
2971  Implicit3Stage4thOrderLobattoA_RKBT()
2972  {
2973  std::ostringstream Description;
2974  Description << this->description() << "\n"
2975  << "A-stable\n"
2976  << "Solving Ordinary Differential Equations II:\n"
2977  << "Stiff and Differential-Algebraic Problems,\n"
2978  << "2nd Revised Edition\n"
2979  << "E. Hairer and G. Wanner\n"
2980  << "Table 5.7, pg 75\n"
2981  << "c = [ 0 1/2 1 ]'\n"
2982  << "A = [ 0 0 0 ]\n"
2983  << " [ 5/24 1/3 -1/24 ]\n"
2984  << " [ 1/6 2/3 1/6 ]\n"
2985  << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
2986  typedef Teuchos::ScalarTraits<Scalar> ST;
2987  using Teuchos::as;
2988  int NumStages = 3;
2989  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
2990  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
2991  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
2992  const Scalar zero = ST::zero();
2993  const Scalar one = ST::one();
2994  A(0,0) = zero;
2995  A(0,1) = zero;
2996  A(0,2) = zero;
2997  A(1,0) = as<Scalar>( (5*one)/(24*one) );
2998  A(1,1) = as<Scalar>( (one)/(3*one) );
2999  A(1,2) = as<Scalar>( (-one)/(24*one) );
3000  A(2,0) = as<Scalar>( (one)/(6*one) );
3001  A(2,1) = as<Scalar>( (2*one)/(3*one) );
3002  A(2,2) = as<Scalar>( (1*one)/(6*one) );
3003  b(0) = as<Scalar>( (one)/(6*one) );
3004  b(1) = as<Scalar>( (2*one)/(3*one) );
3005  b(2) = as<Scalar>( (1*one)/(6*one) );
3006  c(0) = zero;
3007  c(1) = as<Scalar>( one/(2*one) );
3008  c(2) = one;
3009  int order = 4;
3010 
3011  this->initialize(A,b,c,order,Description.str());
3012  }
3013  virtual std::string description() const
3014  { return "RK Implicit 3 Stage 4th order Lobatto A"; }
3015 };
3016 
3017 
3018 // ----------------------------------------------------------------------------
3019 template<class Scalar>
3020 class Implicit4Stage6thOrderLobattoA_RKBT :
3021  virtual public RKButcherTableau<Scalar>
3022 {
3023  public:
3024  Implicit4Stage6thOrderLobattoA_RKBT()
3025  {
3026  using Teuchos::as;
3027  typedef Teuchos::ScalarTraits<Scalar> ST;
3028  std::ostringstream Description;
3029  Description << this->description() << "\n"
3030  << "A-stable\n"
3031  << "Solving Ordinary Differential Equations II:\n"
3032  << "Stiff and Differential-Algebraic Problems,\n"
3033  << "2nd Revised Edition\n"
3034  << "E. Hairer and G. Wanner\n"
3035  << "Table 5.8, pg 75\n"
3036  << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
3037  << "A = [ A1 A2 A3 A4 ]\n"
3038  << " A1 = [ 0 ]\n"
3039  << " [ (11+sqrt(5)/120 ]\n"
3040  << " [ (11-sqrt(5)/120 ]\n"
3041  << " [ 1/12 ]\n"
3042  << " A2 = [ 0 ]\n"
3043  << " [ (25-sqrt(5))/120 ]\n"
3044  << " [ (25+13*sqrt(5))/120 ]\n"
3045  << " [ 5/12 ]\n"
3046  << " A3 = [ 0 ]\n"
3047  << " [ (25-13*sqrt(5))/120 ]\n"
3048  << " [ (25+sqrt(5))/120 ]\n"
3049  << " [ 5/12 ]\n"
3050  << " A4 = [ 0 ]\n"
3051  << " [ (-1+sqrt(5))/120 ]\n"
3052  << " [ (-1-sqrt(5))/120 ]\n"
3053  << " [ 1/12 ]\n"
3054  << "b = [ 1/12 5/12 5/12 1/12 ]'" << std::endl;
3055  typedef Teuchos::ScalarTraits<Scalar> ST;
3056  int NumStages = 4;
3057  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3058  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3059  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3060  const Scalar zero = ST::zero();
3061  const Scalar one = ST::one();
3062  A(0,0) = zero;
3063  A(0,1) = zero;
3064  A(0,2) = zero;
3065  A(0,3) = zero;
3066  A(1,0) = as<Scalar>( (11*one+ST::squareroot(5*one))/(120*one) );
3067  A(1,1) = as<Scalar>( (25*one-ST::squareroot(5*one))/(120*one) );
3068  A(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5*one))/(120*one) );
3069  A(1,3) = as<Scalar>( (-one+ST::squareroot(5*one))/(120*one) );
3070  A(2,0) = as<Scalar>( (11*one-ST::squareroot(5*one))/(120*one) );
3071  A(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5*one))/(120*one) );
3072  A(2,2) = as<Scalar>( (25*one+ST::squareroot(5*one))/(120*one) );
3073  A(2,3) = as<Scalar>( (-one-ST::squareroot(5*one))/(120*one) );
3074  A(3,0) = as<Scalar>( one/(12*one) );
3075  A(3,1) = as<Scalar>( 5*one/(12*one) );
3076  A(3,2) = as<Scalar>( 5*one/(12*one) );
3077  A(3,3) = as<Scalar>( one/(12*one) );
3078  b(0) = as<Scalar>( one/(12*one) );
3079  b(1) = as<Scalar>( 5*one/(12*one) );
3080  b(2) = as<Scalar>( 5*one/(12*one) );
3081  b(3) = as<Scalar>( one/(12*one) );
3082  c(0) = zero;
3083  c(1) = as<Scalar>( (5*one-ST::squareroot(5))/(10*one) );
3084  c(2) = as<Scalar>( (5*one+ST::squareroot(5))/(10*one) );
3085  c(3) = one;
3086  int order = 6;
3087 
3088  this->initialize(A,b,c,order,Description.str());
3089  }
3090  virtual std::string description() const
3091  { return "RK Implicit 4 Stage 6th order Lobatto A"; }
3092 };
3093 
3094 
3095 // ----------------------------------------------------------------------------
3096 template<class Scalar>
3097 class Implicit2Stage2ndOrderLobattoB_RKBT :
3098  virtual public RKButcherTableau<Scalar>
3099 {
3100  public:
3101  Implicit2Stage2ndOrderLobattoB_RKBT()
3102  {
3103  std::ostringstream Description;
3104  Description << this->description() << "\n"
3105  << "A-stable\n"
3106  << "Solving Ordinary Differential Equations II:\n"
3107  << "Stiff and Differential-Algebraic Problems,\n"
3108  << "2nd Revised Edition\n"
3109  << "E. Hairer and G. Wanner\n"
3110  << "Table 5.9, pg 76\n"
3111  << "c = [ 0 1 ]'\n"
3112  << "A = [ 1/2 0 ]\n"
3113  << " [ 1/2 0 ]\n"
3114  << "b = [ 1/2 1/2 ]'" << std::endl;
3115  typedef Teuchos::ScalarTraits<Scalar> ST;
3116  using Teuchos::as;
3117  int NumStages = 2;
3118  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3119  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3120  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3121  const Scalar zero = ST::zero();
3122  const Scalar one = ST::one();
3123  A(0,0) = as<Scalar>( one/(2*one) );
3124  A(0,1) = zero;
3125  A(1,0) = as<Scalar>( one/(2*one) );
3126  A(1,1) = zero;
3127  b(0) = as<Scalar>( one/(2*one) );
3128  b(1) = as<Scalar>( one/(2*one) );
3129  c(0) = zero;
3130  c(1) = one;
3131  int order = 2;
3132 
3133  this->initialize(A,b,c,order,Description.str());
3134  }
3135  virtual std::string description() const
3136  { return "RK Implicit 2 Stage 2nd order Lobatto B"; }
3137 };
3138 
3139 
3140 // ----------------------------------------------------------------------------
3141 template<class Scalar>
3142 class Implicit3Stage4thOrderLobattoB_RKBT :
3143  virtual public RKButcherTableau<Scalar>
3144 {
3145  public:
3146  Implicit3Stage4thOrderLobattoB_RKBT()
3147  {
3148  std::ostringstream Description;
3149  Description << this->description() << "\n"
3150  << "A-stable\n"
3151  << "Solving Ordinary Differential Equations II:\n"
3152  << "Stiff and Differential-Algebraic Problems,\n"
3153  << "2nd Revised Edition\n"
3154  << "E. Hairer and G. Wanner\n"
3155  << "Table 5.9, pg 76\n"
3156  << "c = [ 0 1/2 1 ]'\n"
3157  << "A = [ 1/6 -1/6 0 ]\n"
3158  << " [ 1/6 1/3 0 ]\n"
3159  << " [ 1/6 5/6 0 ]\n"
3160  << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
3161  typedef Teuchos::ScalarTraits<Scalar> ST;
3162  using Teuchos::as;
3163  int NumStages = 3;
3164  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3165  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3166  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3167  const Scalar zero = ST::zero();
3168  const Scalar one = ST::one();
3169  A(0,0) = as<Scalar>( one/(6*one) );
3170  A(0,1) = as<Scalar>( -one/(6*one) );
3171  A(0,2) = zero;
3172  A(1,0) = as<Scalar>( one/(6*one) );
3173  A(1,1) = as<Scalar>( one/(3*one) );
3174  A(1,2) = zero;
3175  A(2,0) = as<Scalar>( one/(6*one) );
3176  A(2,1) = as<Scalar>( 5*one/(6*one) );
3177  A(2,2) = zero;
3178  b(0) = as<Scalar>( one/(6*one) );
3179  b(1) = as<Scalar>( 2*one/(3*one) );
3180  b(2) = as<Scalar>( one/(6*one) );
3181  c(0) = zero;
3182  c(1) = as<Scalar>( one/(2*one) );
3183  c(2) = one;
3184  int order = 4;
3185 
3186  this->initialize(A,b,c,order,Description.str());
3187  }
3188  virtual std::string description() const
3189  { return "RK Implicit 3 Stage 4th order Lobatto B"; }
3190 };
3191 
3192 
3193 // ----------------------------------------------------------------------------
3194 template<class Scalar>
3195 class Implicit4Stage6thOrderLobattoB_RKBT :
3196  virtual public RKButcherTableau<Scalar>
3197 {
3198  public:
3199  Implicit4Stage6thOrderLobattoB_RKBT()
3200  {
3201  std::ostringstream Description;
3202  Description << this->description() << "\n"
3203  << "A-stable\n"
3204  << "Solving Ordinary Differential Equations II:\n"
3205  << "Stiff and Differential-Algebraic Problems,\n"
3206  << "2nd Revised Edition\n"
3207  << "E. Hairer and G. Wanner\n"
3208  << "Table 5.10, pg 76\n"
3209  << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
3210  << "A = [ 1/12 (-1-sqrt(5))/24 (-1+sqrt(5))/24 0 ]\n"
3211  << " [ 1/12 (25+sqrt(5))/120 (25-13*sqrt(5))/120 0 ]\n"
3212  << " [ 1/12 (25+13*sqrt(5))/120 (25-sqrt(5))/120 0 ]\n"
3213  << " [ 1/12 (11-sqrt(5))/24 (11+sqrt(5))/24 0 ]\n"
3214  << "b = [ 1/12 5/12 5/12 1/12 ]'"
3215  << std::endl;
3216  typedef Teuchos::ScalarTraits<Scalar> ST;
3217  using Teuchos::as;
3218  int NumStages = 4;
3219  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3220  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3221  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3222  const Scalar zero = ST::zero();
3223  const Scalar one = ST::one();
3224  A(0,0) = as<Scalar>( one/(12*one) );
3225  A(0,1) = as<Scalar>( (-one-ST::squareroot(5))/(24*one) );
3226  A(0,2) = as<Scalar>( (-one+ST::squareroot(5))/(24*one) );
3227  A(0,3) = zero;
3228  A(1,0) = as<Scalar>( one/(12*one) );
3229  A(1,1) = as<Scalar>( (25*one+ST::squareroot(5))/(120*one) );
3230  A(1,2) = as<Scalar>( (25*one-13*one*ST::squareroot(5))/(120*one) );
3231  A(1,3) = zero;
3232  A(2,0) = as<Scalar>( one/(12*one) );
3233  A(2,1) = as<Scalar>( (25*one+13*one*ST::squareroot(5))/(120*one) );
3234  A(2,2) = as<Scalar>( (25*one-ST::squareroot(5))/(120*one) );
3235  A(2,3) = zero;
3236  A(3,0) = as<Scalar>( one/(12*one) );
3237  A(3,1) = as<Scalar>( (11*one-ST::squareroot(5*one))/(24*one) );
3238  A(3,2) = as<Scalar>( (11*one+ST::squareroot(5*one))/(24*one) );
3239  A(3,3) = zero;
3240  b(0) = as<Scalar>( one/(12*one) );
3241  b(1) = as<Scalar>( 5*one/(12*one) );
3242  b(2) = as<Scalar>( 5*one/(12*one) );
3243  b(3) = as<Scalar>( one/(12*one) );
3244  c(0) = zero;
3245  c(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
3246  c(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
3247  c(3) = one;
3248  int order = 6;
3249 
3250  this->initialize(A,b,c,order,Description.str());
3251  }
3252  virtual std::string description() const
3253  { return "RK Implicit 4 Stage 6th order Lobatto B"; }
3254 };
3255 
3256 
3257 // ----------------------------------------------------------------------------
3258 template<class Scalar>
3259 class Implicit2Stage2ndOrderLobattoC_RKBT :
3260  virtual public RKButcherTableau<Scalar>
3261 {
3262  public:
3263  Implicit2Stage2ndOrderLobattoC_RKBT()
3264  {
3265  std::ostringstream Description;
3266  Description << this->description() << "\n"
3267  << "A-stable\n"
3268  << "Solving Ordinary Differential Equations II:\n"
3269  << "Stiff and Differential-Algebraic Problems,\n"
3270  << "2nd Revised Edition\n"
3271  << "E. Hairer and G. Wanner\n"
3272  << "Table 5.11, pg 76\n"
3273  << "c = [ 0 1 ]'\n"
3274  << "A = [ 1/2 -1/2 ]\n"
3275  << " [ 1/2 1/2 ]\n"
3276  << "b = [ 1/2 1/2 ]'" << std::endl;
3277  typedef Teuchos::ScalarTraits<Scalar> ST;
3278  using Teuchos::as;
3279  int NumStages = 2;
3280  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3281  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3282  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3283  const Scalar zero = ST::zero();
3284  const Scalar one = ST::one();
3285  A(0,0) = as<Scalar>( one/(2*one) );
3286  A(0,1) = as<Scalar>( -one/(2*one) );
3287  A(1,0) = as<Scalar>( one/(2*one) );
3288  A(1,1) = as<Scalar>( one/(2*one) );
3289  b(0) = as<Scalar>( one/(2*one) );
3290  b(1) = as<Scalar>( one/(2*one) );
3291  c(0) = zero;
3292  c(1) = one;
3293  int order = 2;
3294 
3295  this->initialize(A,b,c,order,Description.str());
3296  }
3297  virtual std::string description() const
3298  { return "RK Implicit 2 Stage 2nd order Lobatto C"; }
3299 };
3300 
3301 
3302 // ----------------------------------------------------------------------------
3303 template<class Scalar>
3304 class Implicit3Stage4thOrderLobattoC_RKBT :
3305  virtual public RKButcherTableau<Scalar>
3306 {
3307  public:
3308  Implicit3Stage4thOrderLobattoC_RKBT()
3309  {
3310  std::ostringstream Description;
3311  Description << this->description() << "\n"
3312  << "A-stable\n"
3313  << "Solving Ordinary Differential Equations II:\n"
3314  << "Stiff and Differential-Algebraic Problems,\n"
3315  << "2nd Revised Edition\n"
3316  << "E. Hairer and G. Wanner\n"
3317  << "Table 5.11, pg 76\n"
3318  << "c = [ 0 1/2 1 ]'\n"
3319  << "A = [ 1/6 -1/3 1/6 ]\n"
3320  << " [ 1/6 5/12 -1/12 ]\n"
3321  << " [ 1/6 2/3 1/6 ]\n"
3322  << "b = [ 1/6 2/3 1/6 ]'" << std::endl;
3323  typedef Teuchos::ScalarTraits<Scalar> ST;
3324  using Teuchos::as;
3325  int NumStages = 3;
3326  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3327  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3328  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3329  const Scalar zero = ST::zero();
3330  const Scalar one = ST::one();
3331  A(0,0) = as<Scalar>( one/(6*one) );
3332  A(0,1) = as<Scalar>( -one/(3*one) );
3333  A(0,2) = as<Scalar>( one/(6*one) );
3334  A(1,0) = as<Scalar>( one/(6*one) );
3335  A(1,1) = as<Scalar>( 5*one/(12*one) );
3336  A(1,2) = as<Scalar>( -one/(12*one) );
3337  A(2,0) = as<Scalar>( one/(6*one) );
3338  A(2,1) = as<Scalar>( 2*one/(3*one) );
3339  A(2,2) = as<Scalar>( one/(6*one) );
3340  b(0) = as<Scalar>( one/(6*one) );
3341  b(1) = as<Scalar>( 2*one/(3*one) );
3342  b(2) = as<Scalar>( one/(6*one) );
3343  c(0) = zero;
3344  c(1) = as<Scalar>( one/(2*one) );
3345  c(2) = one;
3346  int order = 4;
3347 
3348  this->initialize(A,b,c,order,Description.str());
3349  }
3350  virtual std::string description() const
3351  { return "RK Implicit 3 Stage 4th order Lobatto C"; }
3352 };
3353 
3354 
3355 // ----------------------------------------------------------------------------
3356 template<class Scalar>
3357 class Implicit4Stage6thOrderLobattoC_RKBT :
3358  virtual public RKButcherTableau<Scalar>
3359 {
3360  public:
3361  Implicit4Stage6thOrderLobattoC_RKBT()
3362  {
3363  std::ostringstream Description;
3364  Description << this->description() << "\n"
3365  << "A-stable\n"
3366  << "Solving Ordinary Differential Equations II:\n"
3367  << "Stiff and Differential-Algebraic Problems,\n"
3368  << "2nd Revised Edition\n"
3369  << "E. Hairer and G. Wanner\n"
3370  << "Table 5.12, pg 76\n"
3371  << "c = [ 0 (5-sqrt(5))/10 (5+sqrt(5))/10 1 ]'\n"
3372  << "A = [ 1/12 -sqrt(5)/12 sqrt(5)/12 -1/12 ]\n"
3373  << " [ 1/12 1/4 (10-7*sqrt(5))/60 sqrt(5)/60 ]\n"
3374  << " [ 1/12 (10+7*sqrt(5))/60 1/4 -sqrt(5)/60 ]\n"
3375  << " [ 1/12 5/12 5/12 1/12 ]\n"
3376  << "b = [ 1/12 5/12 5/12 1/12 ]'"
3377  << std::endl;
3378  typedef Teuchos::ScalarTraits<Scalar> ST;
3379  using Teuchos::as;
3380  int NumStages = 4;
3381  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3382  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3383  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3384  const Scalar zero = ST::zero();
3385  const Scalar one = ST::one();
3386  A(0,0) = as<Scalar>( one/(12*one) );
3387  A(0,1) = as<Scalar>( -ST::squareroot(5*one)/(12*one) );
3388  A(0,2) = as<Scalar>( ST::squareroot(5*one)/(12*one) );
3389  A(0,3) = as<Scalar>( -one/(12*one) );
3390  A(1,0) = as<Scalar>( one/(12*one) );
3391  A(1,1) = as<Scalar>( one/(4*one) );
3392  A(1,2) = as<Scalar>( (10*one-7*one*ST::squareroot(5*one))/(60*one) );
3393  A(1,3) = as<Scalar>( ST::squareroot(5*one)/(60*one) );
3394  A(2,0) = as<Scalar>( one/(12*one) );
3395  A(2,1) = as<Scalar>( (10*one+7*one*ST::squareroot(5*one))/(60*one) );
3396  A(2,2) = as<Scalar>( one/(4*one) );
3397  A(2,3) = as<Scalar>( -ST::squareroot(5*one)/(60*one) );
3398  A(3,0) = as<Scalar>( one/(12*one) );
3399  A(3,1) = as<Scalar>( 5*one/(12*one) );
3400  A(3,2) = as<Scalar>( 5*one/(12*one) );
3401  A(3,3) = as<Scalar>( one/(12*one) );
3402  b(0) = as<Scalar>( one/(12*one) );
3403  b(1) = as<Scalar>( 5*one/(12*one) );
3404  b(2) = as<Scalar>( 5*one/(12*one) );
3405  b(3) = as<Scalar>( one/(12*one) );
3406  c(0) = zero;
3407  c(1) = as<Scalar>( (5*one-ST::squareroot(5*one))/(10*one) );
3408  c(2) = as<Scalar>( (5*one+ST::squareroot(5*one))/(10*one) );
3409  c(3) = one;
3410  int order = 6;
3411 
3412  this->initialize(A,b,c,order,Description.str());
3413  }
3414  virtual std::string description() const
3415  { return "RK Implicit 4 Stage 6th order Lobatto C"; }
3416 };
3417 
3418 
3419 // ----------------------------------------------------------------------------
3420 template<class Scalar>
3421 class SDIRK5Stage5thOrder_RKBT :
3422  virtual public RKButcherTableau<Scalar>
3423 {
3424  public:
3425  SDIRK5Stage5thOrder_RKBT()
3426  {
3427  std::ostringstream Description;
3428  Description << this->description() << "\n"
3429  << "A-stable\n"
3430  << "Solving Ordinary Differential Equations II:\n"
3431  << "Stiff and Differential-Algebraic Problems,\n"
3432  << "2nd Revised Edition\n"
3433  << "E. Hairer and G. Wanner\n"
3434  << "pg101 \n"
3435  << "c = [ (6-sqrt(6))/10 ]\n"
3436  << " [ (6+9*sqrt(6))/35 ]\n"
3437  << " [ 1 ]\n"
3438  << " [ (4-sqrt(6))/10 ]\n"
3439  << " [ (4+sqrt(6))/10 ]\n"
3440  << "A = [ A1 A2 A3 A4 A5 ]\n"
3441  << " A1 = [ (6-sqrt(6))/10 ]\n"
3442  << " [ (-6+5*sqrt(6))/14 ]\n"
3443  << " [ (888+607*sqrt(6))/2850 ]\n"
3444  << " [ (3153-3082*sqrt(6))/14250 ]\n"
3445  << " [ (-32583+14638*sqrt(6))/71250 ]\n"
3446  << " A2 = [ 0 ]\n"
3447  << " [ (6-sqrt(6))/10 ]\n"
3448  << " [ (126-161*sqrt(6))/1425 ]\n"
3449  << " [ (3213+1148*sqrt(6))/28500 ]\n"
3450  << " [ (-17199+364*sqrt(6))/142500 ]\n"
3451  << " A3 = [ 0 ]\n"
3452  << " [ 0 ]\n"
3453  << " [ (6-sqrt(6))/10 ]\n"
3454  << " [ (-267+88*sqrt(6))/500 ]\n"
3455  << " [ (1329-544*sqrt(6))/2500 ]\n"
3456  << " A4 = [ 0 ]\n"
3457  << " [ 0 ]\n"
3458  << " [ 0 ]\n"
3459  << " [ (6-sqrt(6))/10 ]\n"
3460  << " [ (-96+131*sqrt(6))/625 ]\n"
3461  << " A5 = [ 0 ]\n"
3462  << " [ 0 ]\n"
3463  << " [ 0 ]\n"
3464  << " [ 0 ]\n"
3465  << " [ (6-sqrt(6))/10 ]\n"
3466  << "b = [ 0 ]\n"
3467  << " [ 0 ]\n"
3468  << " [ 1/9 ]\n"
3469  << " [ (16-sqrt(6))/36 ]\n"
3470  << " [ (16+sqrt(6))/36 ]" << std::endl;
3471 
3472  this->setDescription(Description.str());
3473  this->setParameterList(Teuchos::null);
3474  }
3475 
3476  virtual std::string description() const { return "SDIRK 5 Stage 5th order"; }
3477 
3478  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
3479  {
3480  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
3481  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
3482  else pl = pList;
3483  // Can not validate because optional parameters (e.g., Solver Name).
3484  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
3485  TEUCHOS_TEST_FOR_EXCEPTION(
3486  pl->get<std::string>("Stepper Type") != this->description()
3487  ,std::runtime_error,
3488  " Stepper Type != \""+this->description()+"\"\n"
3489  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
3490 
3491  typedef Teuchos::ScalarTraits<Scalar> ST;
3492  using Teuchos::as;
3493  int NumStages = 5;
3494  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3495  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3496  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3497  const Scalar zero = ST::zero();
3498  const Scalar one = ST::one();
3499  const Scalar sqrt6 = ST::squareroot(as<Scalar>(6*one));
3500  const Scalar gamma = as<Scalar>( (6*one - sqrt6) / (10*one) ); // diagonal
3501  A(0,0) = gamma;
3502  A(0,1) = zero;
3503  A(0,2) = zero;
3504  A(0,3) = zero;
3505  A(0,4) = zero;
3506 
3507  A(1,0) = as<Scalar>( (-6*one+5*one*sqrt6)/(14*one) );
3508  A(1,1) = gamma;
3509  A(1,2) = zero;
3510  A(1,3) = zero;
3511  A(1,4) = zero;
3512 
3513  A(2,0) = as<Scalar>( (888*one+607*one*sqrt6)/(2850*one) );
3514  A(2,1) = as<Scalar>( (126*one-161*one*sqrt6)/(1425*one) );
3515  A(2,2) = gamma;
3516  A(2,3) = zero;
3517  A(2,4) = zero;
3518 
3519  A(3,0) = as<Scalar>( (3153*one-3082*one*sqrt6)/(14250*one) );
3520  A(3,1) = as<Scalar>( (3213*one+1148*one*sqrt6)/(28500*one) );
3521  A(3,2) = as<Scalar>( (-267*one+88*one*sqrt6)/(500*one) );
3522  A(3,3) = gamma;
3523  A(3,4) = zero;
3524 
3525  A(4,0) = as<Scalar>( (-32583*one+14638*one*sqrt6)/(71250*one) );
3526  A(4,1) = as<Scalar>( (-17199*one+364*one*sqrt6)/(142500*one) );
3527  A(4,2) = as<Scalar>( (1329*one-544*one*sqrt6)/(2500*one) );
3528  A(4,3) = as<Scalar>( (-96*one+131*sqrt6)/(625*one) );
3529  A(4,4) = gamma;
3530 
3531  b(0) = zero;
3532  b(1) = zero;
3533  b(2) = as<Scalar>( one/(9*one) );
3534  b(3) = as<Scalar>( (16*one-sqrt6)/(36*one) );
3535  b(4) = as<Scalar>( (16*one+sqrt6)/(36*one) );
3536 
3537  c(0) = gamma;
3538  c(1) = as<Scalar>( (6*one+9*one*sqrt6)/(35*one) );
3539  c(2) = one;
3540  c(3) = as<Scalar>( (4*one-sqrt6)/(10*one) );
3541  c(4) = as<Scalar>( (4*one+sqrt6)/(10*one) );
3542 
3543  int order = 5;
3544 
3545  this->initialize(A,b,c,order,this->getDescription());
3546  this->setMyParamList(pl);
3547  this->rkbtPL_ = pl;
3548  }
3549 
3550  Teuchos::RCP<const Teuchos::ParameterList>
3551  getValidParameters() const
3552  {
3553  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
3554  pl->setName("Default Stepper - " + this->description());
3555  pl->set<std::string>("Description", this->getDescription());
3556  pl->set<std::string>("Stepper Type", this->description());
3557  pl->set<bool>("Use Embedded", false);
3558  pl->set<std::string>("Solver Name", "",
3559  "Name of ParameterList containing the solver specifications.");
3560 
3561  return pl;
3562  }
3563 };
3564 
3565 
3566 // ----------------------------------------------------------------------------
3567 template<class Scalar>
3568 class SDIRK5Stage4thOrder_RKBT :
3569  virtual public RKButcherTableau<Scalar>
3570 {
3571  public:
3572  SDIRK5Stage4thOrder_RKBT()
3573  {
3574  std::ostringstream Description;
3575  Description << this->description() << "\n"
3576  << "L-stable\n"
3577  << "Solving Ordinary Differential Equations II:\n"
3578  << "Stiff and Differential-Algebraic Problems,\n"
3579  << "2nd Revised Edition\n"
3580  << "E. Hairer and G. Wanner\n"
3581  << "pg100 \n"
3582  << "c = [ 1/4 3/4 11/20 1/2 1 ]'\n"
3583  << "A = [ 1/4 ]\n"
3584  << " [ 1/2 1/4 ]\n"
3585  << " [ 17/50 -1/25 1/4 ]\n"
3586  << " [ 371/1360 -137/2720 15/544 1/4 ]\n"
3587  << " [ 25/24 -49/48 125/16 -85/12 1/4 ]\n"
3588  << "b = [ 25/24 -49/48 125/16 -85/12 1/4 ]'\n"
3589  << "b' = [ 59/48 -17/96 225/32 -85/12 0 ]'" << std::endl;
3590 
3591  this->setDescription(Description.str());
3592  this->setParameterList(Teuchos::null);
3593  }
3594 
3595  virtual std::string description() const { return "SDIRK 5 Stage 4th order"; }
3596 
3597  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
3598  {
3599  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
3600  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
3601  else pl = pList;
3602  // Can not validate because optional parameters (e.g., Solver Name).
3603  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
3604  TEUCHOS_TEST_FOR_EXCEPTION(
3605  pl->get<std::string>("Stepper Type") != this->description()
3606  ,std::runtime_error,
3607  " Stepper Type != \""+this->description()+"\"\n"
3608  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
3609 
3610  typedef Teuchos::ScalarTraits<Scalar> ST;
3611  using Teuchos::as;
3612  int NumStages = 5;
3613  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3614  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3615  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3616  const Scalar zero = ST::zero();
3617  const Scalar one = ST::one();
3618  const Scalar onequarter = as<Scalar>( one/(4*one) );
3619  A(0,0) = onequarter;
3620  A(0,1) = zero;
3621  A(0,2) = zero;
3622  A(0,3) = zero;
3623  A(0,4) = zero;
3624 
3625  A(1,0) = as<Scalar>( one / (2*one) );
3626  A(1,1) = onequarter;
3627  A(1,2) = zero;
3628  A(1,3) = zero;
3629  A(1,4) = zero;
3630 
3631  A(2,0) = as<Scalar>( 17*one/(50*one) );
3632  A(2,1) = as<Scalar>( -one/(25*one) );
3633  A(2,2) = onequarter;
3634  A(2,3) = zero;
3635  A(2,4) = zero;
3636 
3637  A(3,0) = as<Scalar>( 371*one/(1360*one) );
3638  A(3,1) = as<Scalar>( -137*one/(2720*one) );
3639  A(3,2) = as<Scalar>( 15*one/(544*one) );
3640  A(3,3) = onequarter;
3641  A(3,4) = zero;
3642 
3643  A(4,0) = as<Scalar>( 25*one/(24*one) );
3644  A(4,1) = as<Scalar>( -49*one/(48*one) );
3645  A(4,2) = as<Scalar>( 125*one/(16*one) );
3646  A(4,3) = as<Scalar>( -85*one/(12*one) );
3647  A(4,4) = onequarter;
3648 
3649  b(0) = as<Scalar>( 25*one/(24*one) );
3650  b(1) = as<Scalar>( -49*one/(48*one) );
3651  b(2) = as<Scalar>( 125*one/(16*one) );
3652  b(3) = as<Scalar>( -85*one/(12*one) );
3653  b(4) = onequarter;
3654 
3655  /*
3656  // Alternate version
3657  b(0) = as<Scalar>( 59*one/(48*one) );
3658  b(1) = as<Scalar>( -17*one/(96*one) );
3659  b(2) = as<Scalar>( 225*one/(32*one) );
3660  b(3) = as<Scalar>( -85*one/(12*one) );
3661  b(4) = zero;
3662  */
3663  c(0) = onequarter;
3664  c(1) = as<Scalar>( 3*one/(4*one) );
3665  c(2) = as<Scalar>( 11*one/(20*one) );
3666  c(3) = as<Scalar>( one/(2*one) );
3667  c(4) = one;
3668 
3669  int order = 4;
3670 
3671  this->initialize(A,b,c,order,this->getDescription());
3672  this->setMyParamList(pl);
3673  this->rkbtPL_ = pl;
3674  }
3675 
3676  Teuchos::RCP<const Teuchos::ParameterList>
3677  getValidParameters() const
3678  {
3679  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
3680  pl->setName("Default Stepper - " + this->description());
3681  pl->set<std::string>("Description", this->getDescription());
3682  pl->set<std::string>("Stepper Type", this->description());
3683  pl->set<bool>("Use Embedded", false);
3684  pl->set<std::string>("Solver Name", "",
3685  "Name of ParameterList containing the solver specifications.");
3686 
3687  return pl;
3688  }
3689 };
3690 
3691 
3692 // ----------------------------------------------------------------------------
3693 template<class Scalar>
3694 class SDIRK3Stage4thOrder_RKBT :
3695  virtual public RKButcherTableau<Scalar>
3696 {
3697  public:
3698  SDIRK3Stage4thOrder_RKBT()
3699  {
3700  std::ostringstream Description;
3701  Description << this->description() << "\n"
3702  << "A-stable\n"
3703  << "Solving Ordinary Differential Equations II:\n"
3704  << "Stiff and Differential-Algebraic Problems,\n"
3705  << "2nd Revised Edition\n"
3706  << "E. Hairer and G. Wanner\n"
3707  << "pg100 \n"
3708  << "gamma = (1/sqrt(3))*cos(pi/18)+1/2\n"
3709  << "delta = 1/(6*(2*gamma-1)^2)\n"
3710  << "c = [ gamma 1/2 1-gamma ]'\n"
3711  << "A = [ gamma ]\n"
3712  << " [ 1/2-gamma gamma ]\n"
3713  << " [ 2*gamma 1-4*gamma gamma ]\n"
3714  << "b = [ delta 1-2*delta delta ]'" << std::endl;
3715 
3716  this->setDescription(Description.str());
3717  this->setParameterList(Teuchos::null);
3718  }
3719 
3720  virtual std::string description() const { return "SDIRK 3 Stage 4th order"; }
3721 
3722  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
3723  {
3724  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
3725  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
3726  else pl = pList;
3727  // Can not validate because optional parameters (e.g., Solver Name).
3728  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
3729  TEUCHOS_TEST_FOR_EXCEPTION(
3730  pl->get<std::string>("Stepper Type") != this->description()
3731  ,std::runtime_error,
3732  " Stepper Type != \""+this->description()+"\"\n"
3733  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
3734 
3735  typedef Teuchos::ScalarTraits<Scalar> ST;
3736  using Teuchos::as;
3737  int NumStages = 3;
3738  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3739  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3740  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3741  const Scalar zero = ST::zero();
3742  const Scalar one = ST::one();
3743  const Scalar pi = as<Scalar>(4*one)*std::atan(one);
3744  const Scalar gamma = as<Scalar>( one/ST::squareroot(3*one)*std::cos(pi/(18*one))+one/(2*one) );
3745  const Scalar delta = as<Scalar>( one/(6*one*std::pow(2*gamma-one,2*one)) );
3746  A(0,0) = gamma;
3747  A(0,1) = zero;
3748  A(0,2) = zero;
3749 
3750  A(1,0) = as<Scalar>( one/(2*one) - gamma );
3751  A(1,1) = gamma;
3752  A(1,2) = zero;
3753 
3754  A(2,0) = as<Scalar>( 2*gamma );
3755  A(2,1) = as<Scalar>( one - 4*gamma );
3756  A(2,2) = gamma;
3757 
3758  b(0) = delta;
3759  b(1) = as<Scalar>( one-2*delta );
3760  b(2) = delta;
3761 
3762  c(0) = gamma;
3763  c(1) = as<Scalar>( one/(2*one) );
3764  c(2) = as<Scalar>( one - gamma );
3765 
3766  int order = 4;
3767 
3768  this->initialize(A,b,c,order,this->getDescription());
3769  this->setMyParamList(pl);
3770  this->rkbtPL_ = pl;
3771  }
3772 
3773  Teuchos::RCP<const Teuchos::ParameterList>
3774  getValidParameters() const
3775  {
3776  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
3777  pl->setName("Default Stepper - " + this->description());
3778  pl->set<std::string>("Description", this->getDescription());
3779  pl->set<std::string>("Stepper Type", this->description());
3780  pl->set<bool>("Use Embedded", false);
3781  pl->set<std::string>("Solver Name", "",
3782  "Name of ParameterList containing the solver specifications.");
3783 
3784  return pl;
3785  }
3786 };
3787 
3788 // ----------------------------------------------------------------------------
3789 /** \brief SDIRK 2(1) pair
3790  *
3791  * The tableau (order=2(1)) is
3792  * \f[
3793  * \begin{array}{c|c}
3794  * c & A \\ \hline
3795  * & b^T \\ \hline
3796  * & \hat{b}^T
3797  * \end{array}
3798  * \;\;\;\;\mbox{ where }\;\;\;\;
3799  * \begin{array}{c|cccc} 0 & 0 & \\
3800  * 1 & -1 & 1 \\ \hline
3801  * & 1/2 & 1/2 \\
3802  * & 1 & 0 \end{array}
3803  * \f]
3804  *
3805  */
3806 template<class Scalar>
3807 class SDIRK21_RKBT :
3808  virtual public RKButcherTableau<Scalar>
3809 {
3810  public:
3811  SDIRK21_RKBT()
3812  {
3813  std::ostringstream Description;
3814  Description << this->description() << "\n"
3815  << "c = [ 1 0 ]'\n"
3816  << "A = [ 1 ]\n"
3817  << " [ -1 1 ]\n"
3818  << "b = [ 1/2 1/2 ]\n"
3819  << "bstar = [ 1 0 ]\n" << std::endl;
3820  this->setDescription(Description.str());
3821  this->setParameterList(Teuchos::null);
3822  }
3823 
3824  virtual std::string description() const { return "SDIRK 2(1) Pair"; }
3825 
3826  void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pList)
3827  {
3828  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
3829  if (pList == Teuchos::null) *pl = *(this->getValidParameters());
3830  else pl = pList;
3831  // Can not validate because optional parameters (e.g., Solver Name).
3832  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
3833  TEUCHOS_TEST_FOR_EXCEPTION(
3834  pl->get<std::string>("Stepper Type") != this->description()
3835  ,std::runtime_error,
3836  " Stepper Type != \""+this->description()+"\"\n"
3837  " Stepper Type = " + pl->get<std::string>("Stepper Type"));
3838 
3839  typedef Teuchos::ScalarTraits<Scalar> ST;
3840  using Teuchos::as;
3841  int NumStages = 2;
3842  Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
3843  Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
3844  Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
3845  Teuchos::SerialDenseVector<int,Scalar> bstar(NumStages);
3846 
3847  const Scalar one = ST::one();
3848  const Scalar zero = ST::zero();
3849 
3850  // Fill A:
3851  A(0,0) = one;
3852  A(0,1) = zero;
3853 
3854  //A(1,0) =
3855  A(1,0) = -one;
3856  A(1,1) = one;
3857 
3858  // Fill b:
3859  b(0) = as<Scalar>(one/(2*one));
3860  b(1) = as<Scalar>(one/(2*one));
3861 
3862  // Fill c:
3863  c(0) = one;
3864  c(1) = zero;
3865 
3866  // Fill bstar
3867  bstar(0) = one;
3868  bstar(1) = zero;
3869  int order = 2;
3870 
3871  this->initialize(A,b,c,order,this->getDescription(),true,bstar);
3872  this->setMyParamList(pl);
3873  this->rkbtPL_ = pl;
3874  }
3875 
3876  Teuchos::RCP<const Teuchos::ParameterList>
3877  getValidParameters() const
3878  {
3879  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
3880  pl->setName("Default Stepper - " + this->description());
3881  pl->set<std::string>("Description", this->getDescription());
3882  pl->set<std::string>("Stepper Type", this->description());
3883  pl->set<bool>("Use Embedded", false);
3884  pl->set<std::string>("Solver Name", "",
3885  "Name of ParameterList containing the solver specifications.");
3886 
3887  return pl;
3888  }
3889 };
3890 
3891 
3892 } // namespace Tempus
3893 
3894 
3895 #endif // Tempus_RKButcherTableau_hpp
Tempus::Implicit3Stage4thOrderLobattoC_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3350
Tempus::RKButcherTableau::set_order
void set_order(const int &order)
Definition: Tempus_RKButcherTableau.hpp:196
Tempus::Implicit3Stage4thOrderLobattoB_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3188
Tempus::SDIRK2Stage2ndOrder_RKBT::gamma_
Scalar gamma_
Definition: Tempus_RKButcherTableau.hpp:1860
Tempus::SDIRK3Stage4thOrder_RKBT::SDIRK3Stage4thOrder_RKBT
SDIRK3Stage4thOrder_RKBT()
Definition: Tempus_RKButcherTableau.hpp:3698
Tempus::RKButcherTableau
Runge-Kutta methods.
Definition: Tempus_RKButcherTableau.hpp:53
Tempus::RKButcherTableau::isEmbedded
virtual bool isEmbedded() const
Return true if the RK method has embedded capabilities.
Definition: Tempus_RKButcherTableau.hpp:84
Tempus::ExplicitMerson45_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:891
Tempus::General_RKButcherTableau
Definition: Tempus_RKButcherTableau.hpp:279
Tempus::RKButcherTableau::A_
Teuchos::SerialDenseMatrix< int, Scalar > A_
Definition: Tempus_RKButcherTableau.hpp:218
Tempus::Implicit4Stage6thOrderLobattoA_RKBT::Implicit4Stage6thOrderLobattoA_RKBT
Implicit4Stage6thOrderLobattoA_RKBT()
Definition: Tempus_RKButcherTableau.hpp:3024
Tempus::RKButcherTableau::set_orderMax
void set_orderMax(const int &order)
Definition: Tempus_RKButcherTableau.hpp:198
Tempus::SDIRK1Stage1stOrder_RKBT::SDIRK1Stage1stOrder_RKBT
SDIRK1Stage1stOrder_RKBT()
Definition: Tempus_RKButcherTableau.hpp:1691
Tempus::SDIRK21_RKBT
SDIRK 2(1) pair.
Definition: Tempus_RKButcherTableau.hpp:3807
Tempus::Explicit3_8Rule_RKBT
Explicit RK 3/8th Rule Butcher Tableau.
Definition: Tempus_RKButcherTableau.hpp:916
Tempus::Implicit2Stage2ndOrderLobattoC_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3297
Tempus::ForwardEuler_RKBT
Forward Euler Runge-Kutta Butcher Tableau.
Definition: Tempus_RKButcherTableau.hpp:587
Tempus::Explicit5Stage3rdOrderKandG_RKBT
RK Explicit 5 Stage 3rd order by Kinnmark and Gray.
Definition: Tempus_RKButcherTableau.hpp:1104
Tempus::SDIRK3Stage4thOrder_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3720
Tempus::SDIRK2Stage3rdOrder_RKBT::thirdOrderAStable_default_
bool thirdOrderAStable_default_
Definition: Tempus_RKButcherTableau.hpp:2015
Tempus::RKButcherTableau::set_isImplicit
void set_isImplicit()
Definition: Tempus_RKButcherTableau.hpp:199
Tempus::SDIRK2Stage2ndOrder_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:1840
Tempus::GeneralDIRK_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:1648
Tempus::Implicit1Stage1stOrderRadauA_RKBT
Definition: Tempus_RKButcherTableau.hpp:2644
Tempus::Implicit2Stage3rdOrderRadauB_RKBT::Implicit2Stage3rdOrderRadauB_RKBT
Implicit2Stage3rdOrderRadauB_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2822
Tempus::Explicit3Stage3rdOrderTVD_RKBT::Explicit3Stage3rdOrderTVD_RKBT
Explicit3Stage3rdOrderTVD_RKBT()
Definition: Tempus_RKButcherTableau.hpp:1297
Tempus::RKButcherTableau::isImplicit_
bool isImplicit_
Definition: Tempus_RKButcherTableau.hpp:224
Tempus::Implicit2Stage2ndOrderLobattoB_RKBT::Implicit2Stage2ndOrderLobattoB_RKBT
Implicit2Stage2ndOrderLobattoB_RKBT()
Definition: Tempus_RKButcherTableau.hpp:3101
Tempus::Implicit1Stage1stOrderRadauA_RKBT::Implicit1Stage1stOrderRadauA_RKBT
Implicit1Stage1stOrderRadauA_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2648
Tempus::Implicit3Stage6thOrderKuntzmannButcher_RKBT::Implicit3Stage6thOrderKuntzmannButcher_RKBT
Implicit3Stage6thOrderKuntzmannButcher_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2128
Tempus::Implicit2Stage4thOrderHammerHollingsworth_RKBT
Definition: Tempus_RKButcherTableau.hpp:2258
Tempus::RKButcherTableau::set_orderMin
void set_orderMin(const int &order)
Definition: Tempus_RKButcherTableau.hpp:197
Tempus::RKButcherTableau::set_b
void set_b(const Teuchos::SerialDenseVector< int, Scalar > &b)
Definition: Tempus_RKButcherTableau.hpp:194
Tempus::SDIRK2Stage3rdOrder_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2012
Tempus::IRK1StageTheta_RKBT::theta_default_
Scalar theta_default_
Definition: Tempus_RKButcherTableau.hpp:2384
Tempus::Explicit2Stage2ndOrderRunge_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:1509
Tempus::BackwardEuler_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:552
Tempus::Explicit2Stage2ndOrderRunge_RKBT::Explicit2Stage2ndOrderRunge_RKBT
Explicit2Stage2ndOrderRunge_RKBT()
Definition: Tempus_RKButcherTableau.hpp:1467
Tempus_String_Utilities.hpp
Tempus::Implicit4Stage6thOrderLobattoC_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3414
Tempus::SDIRK21_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:3877
Tempus::ForwardEuler_RKBT::ForwardEuler_RKBT
ForwardEuler_RKBT()
Definition: Tempus_RKButcherTableau.hpp:591
Tempus::RKButcherTableau::getDescription
const std::string & getDescription() const
Definition: Tempus_RKButcherTableau.hpp:191
Tempus::Implicit4Stage6thOrderLobattoC_RKBT::Implicit4Stage6thOrderLobattoC_RKBT
Implicit4Stage6thOrderLobattoC_RKBT()
Definition: Tempus_RKButcherTableau.hpp:3361
Tempus::SDIRK2Stage2ndOrder_RKBT::gamma_default_
Scalar gamma_default_
Definition: Tempus_RKButcherTableau.hpp:1859
Tempus::Implicit3Stage4thOrderLobattoB_RKBT
Definition: Tempus_RKButcherTableau.hpp:3142
Tempus::Implicit3Stage6thOrderGauss_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2637
Tempus::Implicit3Stage5thOrderRadauB_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2915
Tempus::StringTokenizer
void StringTokenizer(std::vector< std::string > &tokens, const std::string &str, const std::string delimiters, bool trim)
Tokenize a string, put tokens in a vector.
Definition: Tempus_String_Utilities.cpp:31
Tempus::SDIRK5Stage4thOrder_RKBT::SDIRK5Stage4thOrder_RKBT
SDIRK5Stage4thOrder_RKBT()
Definition: Tempus_RKButcherTableau.hpp:3572
Tempus::ExplicitTrapezoidal_RKBT::ExplicitTrapezoidal_RKBT
ExplicitTrapezoidal_RKBT()
Definition: Tempus_RKButcherTableau.hpp:1534
Tempus::Implicit4Stage8thOrderKuntzmannButcher_RKBT
Definition: Tempus_RKButcherTableau.hpp:2176
Tempus::SDIRK5Stage5thOrder_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:3551
Tempus::Explicit4Stage4thOrder_RKBT
Runge-Kutta 4th order Butcher Tableau.
Definition: Tempus_RKButcherTableau.hpp:629
Tempus::RKButcherTableau::A
virtual const Teuchos::SerialDenseMatrix< int, Scalar > & A() const
Return the matrix coefficients.
Definition: Tempus_RKButcherTableau.hpp:62
Tempus::Implicit3Stage4thOrderLobattoA_RKBT::Implicit3Stage4thOrderLobattoA_RKBT
Implicit3Stage4thOrderLobattoA_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2971
Tempus::Implicit3Stage5thOrderRadauA_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2774
Tempus::SDIRK5Stage5thOrder_RKBT
Definition: Tempus_RKButcherTableau.hpp:3421
Tempus
Definition: Tempus_AdjointAuxSensitivityModelEvaluator_decl.hpp:20
Tempus::Implicit1Stage2ndOrderGauss_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2525
Tempus::ForwardEuler_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:607
Tempus::Implicit3Stage4thOrderLobattoB_RKBT::Implicit3Stage4thOrderLobattoB_RKBT
Implicit3Stage4thOrderLobattoB_RKBT()
Definition: Tempus_RKButcherTableau.hpp:3146
Tempus::BackwardEuler_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:522
Tempus::Implicit1Stage1stOrderRadauB_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2811
Tempus::Implicit4Stage6thOrderLobattoB_RKBT
Definition: Tempus_RKButcherTableau.hpp:3195
Tempus::EDIRK2Stage3rdOrder_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2066
Tempus::RKButcherTableau::orderMax
virtual int orderMax() const
Return the maximum order.
Definition: Tempus_RKButcherTableau.hpp:78
Tempus::RKButcherTableau::setParameterList
virtual void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:136
Tempus::Implicit2Stage4thOrderHammerHollingsworth_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2296
Tempus::SDIRK2Stage3rdOrder_RKBT::gamma_
Scalar gamma_
Definition: Tempus_RKButcherTableau.hpp:2020
Tempus::Explicit3Stage3rdOrderHeun_RKBT
RK Explicit 3 Stage 3rd order by Heun.
Definition: Tempus_RKButcherTableau.hpp:1381
Tempus::rKButcherTableau
Teuchos::RCP< RKButcherTableau< Scalar > > rKButcherTableau()
Definition: Tempus_RKButcherTableau.hpp:238
Tempus::Explicit4Stage4thOrder_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:697
Tempus::RKButcherTableau::b_
Teuchos::SerialDenseVector< int, Scalar > b_
Definition: Tempus_RKButcherTableau.hpp:219
Tempus::RKButcherTableau::isImplicit
virtual bool isImplicit() const
Return true if the RK method is implicit.
Definition: Tempus_RKButcherTableau.hpp:80
Tempus::RKButcherTableau::bstar_
Teuchos::SerialDenseVector< int, Scalar > bstar_
Definition: Tempus_RKButcherTableau.hpp:229
Tempus::RKButcherTableau::initialize
virtual void initialize(const Teuchos::SerialDenseMatrix< int, Scalar > &A, const Teuchos::SerialDenseVector< int, Scalar > &b, const Teuchos::SerialDenseVector< int, Scalar > &c, const int order, const std::string &longDescription, bool isEmbedded=false, const Teuchos::SerialDenseVector< int, Scalar > &bstar=Teuchos::SerialDenseVector< int, Scalar >())
Definition: Tempus_RKButcherTableau.hpp:86
Tempus::GeneralExplicit_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:459
Tempus::Implicit2Stage3rdOrderRadauB_RKBT
Definition: Tempus_RKButcherTableau.hpp:2818
Tempus::SDIRK2Stage2ndOrder_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:1798
Tempus::Implicit3Stage4thOrderLobattoA_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3013
Tempus::SDIRK1Stage1stOrder_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:1705
Tempus::Explicit3Stage3rdOrderTVD_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:1355
Tempus::RKButcherTableau::orderMax_
int orderMax_
Definition: Tempus_RKButcherTableau.hpp:223
Tempus::Implicit4Stage6thOrderLobattoB_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3252
Tempus::ExplicitBogackiShampine32_RKBT::ExplicitBogackiShampine32_RKBT
ExplicitBogackiShampine32_RKBT()
Definition: Tempus_RKButcherTableau.hpp:729
Tempus::SDIRK3Stage4thOrder_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:3774
Tempus::EDIRK2Stage3rdOrder_RKBT::EDIRK2Stage3rdOrder_RKBT
EDIRK2Stage3rdOrder_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2048
Tempus::Implicit2Stage3rdOrderRadauA_RKBT::Implicit2Stage3rdOrderRadauA_RKBT
Implicit2Stage3rdOrderRadauA_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2686
Tempus::GeneralDIRK_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:1638
Tempus::Implicit1Stage2ndOrderGauss_RKBT
Definition: Tempus_RKButcherTableau.hpp:2489
Tempus::Implicit3Stage6thOrderKuntzmannButcher_RKBT
Definition: Tempus_RKButcherTableau.hpp:2124
Tempus::Explicit5Stage3rdOrderKandG_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:1185
Tempus::Implicit3Stage4thOrderLobattoA_RKBT
Definition: Tempus_RKButcherTableau.hpp:2967
Tempus::SDIRK5Stage5thOrder_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:3478
Tempus::RKButcherTableau::description
virtual std::string description() const =0
Tempus::Implicit3Stage5thOrderRadauB_RKBT
Definition: Tempus_RKButcherTableau.hpp:2862
Tempus::RKButcherTableau::set_isDIRK
void set_isDIRK()
DIRK is defined as if a_ij = 0 for j>i and a_ii != 0 for at least one i.
Definition: Tempus_RKButcherTableau.hpp:206
Tempus::Implicit3Stage4thOrderLobattoC_RKBT
Definition: Tempus_RKButcherTableau.hpp:3304
Tempus::RKButcherTableau::bstar
virtual const Teuchos::SerialDenseVector< int, Scalar > & bstar() const
Return the vector of quadrature weights for embedded methods.
Definition: Tempus_RKButcherTableau.hpp:68
Tempus::EDIRK2Stage3rdOrder_RKBT
EDIRK 2 Stage 3rd order.
Definition: Tempus_RKButcherTableau.hpp:2044
Tempus::Implicit4Stage6thOrderLobattoA_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3090
Tempus::BackwardEuler_RKBT::BackwardEuler_RKBT
BackwardEuler_RKBT()
Definition: Tempus_RKButcherTableau.hpp:510
Tempus::SDIRK5Stage5thOrder_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3476
Tempus::SDIRK21_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:3826
Tempus::Implicit2Stage2ndOrderLobattoC_RKBT::Implicit2Stage2ndOrderLobattoC_RKBT
Implicit2Stage2ndOrderLobattoC_RKBT()
Definition: Tempus_RKButcherTableau.hpp:3263
Tempus::Implicit2Stage4thOrderGauss_RKBT::Implicit2Stage4thOrderGauss_RKBT
Implicit2Stage4thOrderGauss_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2536
Tempus::Explicit3Stage3rdOrderHeun_RKBT::Explicit3Stage3rdOrderHeun_RKBT
Explicit3Stage3rdOrderHeun_RKBT()
Definition: Tempus_RKButcherTableau.hpp:1385
Tempus::RKButcherTableau::isDIRK_
bool isDIRK_
Definition: Tempus_RKButcherTableau.hpp:225
Tempus::IRK1StageTheta_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:2363
Tempus::ExplicitMerson45_RKBT
Explicit RK Merson Butcher Tableau.
Definition: Tempus_RKButcherTableau.hpp:826
Tempus::SDIRK2Stage3rdOrder_RKBT::secondOrderLStable_
bool secondOrderLStable_
Definition: Tempus_RKButcherTableau.hpp:2018
Tempus::SDIRK2Stage2ndOrder_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:1800
Tempus::Implicit2Stage3rdOrderRadauA_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2720
Tempus::IRK1StageTheta_RKBT::theta_
Scalar theta_
Definition: Tempus_RKButcherTableau.hpp:2385
Tempus::Implicit2Stage2ndOrderLobattoA_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2960
Tempus::Implicit4Stage6thOrderLobattoA_RKBT
Definition: Tempus_RKButcherTableau.hpp:3020
Tempus::SDIRK2Stage3rdOrder_RKBT::SDIRK2Stage3rdOrder_RKBT
SDIRK2Stage3rdOrder_RKBT()
Definition: Tempus_RKButcherTableau.hpp:1896
Tempus::General_RKButcherTableau::parseGeneralPL
void parseGeneralPL(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:283
Tempus::Implicit4Stage6thOrderLobattoC_RKBT
Definition: Tempus_RKButcherTableau.hpp:3357
Tempus::ExplicitMerson45_RKBT::ExplicitMerson45_RKBT
ExplicitMerson45_RKBT()
Definition: Tempus_RKButcherTableau.hpp:830
Tempus::Explicit5Stage3rdOrderKandG_RKBT::Explicit5Stage3rdOrderKandG_RKBT
Explicit5Stage3rdOrderKandG_RKBT()
Definition: Tempus_RKButcherTableau.hpp:1108
Tempus::Implicit1Stage1stOrderRadauB_RKBT::Implicit1Stage1stOrderRadauB_RKBT
Implicit1Stage1stOrderRadauB_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2785
Tempus::IRK1StageTheta_RKBT::IRK1StageTheta_RKBT
IRK1StageTheta_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2307
Tempus::SDIRK1Stage1stOrder_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:1703
Tempus::RKButcherTableau::set_A
void set_A(const Teuchos::SerialDenseMatrix< int, Scalar > &A)
Definition: Tempus_RKButcherTableau.hpp:193
Tempus::ExplicitBogackiShampine32_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:796
Tempus::EDIRK2StageTheta_RKBT::theta_
Scalar theta_
Definition: Tempus_RKButcherTableau.hpp:2483
Tempus::EDIRK2StageTheta_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:2416
Tempus::GeneralDIRK_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:1640
Tempus::EDIRK2StageTheta_RKBT::EDIRK2StageTheta_RKBT
EDIRK2StageTheta_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2395
Tempus::RKButcherTableau::set_c
void set_c(const Teuchos::SerialDenseVector< int, Scalar > &c)
Definition: Tempus_RKButcherTableau.hpp:195
Tempus::Explicit3Stage3rdOrder_RKBT::Explicit3Stage3rdOrder_RKBT
Explicit3Stage3rdOrder_RKBT()
Definition: Tempus_RKButcherTableau.hpp:1211
Tempus::Implicit3Stage5thOrderRadauA_RKBT
Definition: Tempus_RKButcherTableau.hpp:2727
Tempus::SDIRK5Stage4thOrder_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:3597
Tempus::Implicit2Stage4thOrderHammerHollingsworth_RKBT::Implicit2Stage4thOrderHammerHollingsworth_RKBT
Implicit2Stage4thOrderHammerHollingsworth_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2262
Tempus::SDIRK3Stage4thOrder_RKBT
Definition: Tempus_RKButcherTableau.hpp:3694
Tempus::Explicit3Stage3rdOrder_RKBT
RK Explicit 3 Stage 3rd order.
Definition: Tempus_RKButcherTableau.hpp:1207
Tempus::SDIRK3Stage4thOrder_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:3722
Tempus::Explicit4Stage3rdOrderRunge_RKBT::Explicit4Stage3rdOrderRunge_RKBT
Explicit4Stage3rdOrderRunge_RKBT()
Definition: Tempus_RKButcherTableau.hpp:1015
Tempus::Implicit1Stage2ndOrderGauss_RKBT::Implicit1Stage2ndOrderGauss_RKBT
Implicit1Stage2ndOrderGauss_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2493
Tempus::SDIRK1Stage1stOrder_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:1733
Tempus::RKButcherTableau::numStages
virtual std::size_t numStages() const
Return the number of stages.
Definition: Tempus_RKButcherTableau.hpp:60
Tempus::Explicit4Stage3rdOrderRunge_RKBT
RK Explicit 4 Stage 3rd order by Runge.
Definition: Tempus_RKButcherTableau.hpp:1011
Tempus::Explicit4Stage3rdOrderRunge_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:1078
Tempus::Explicit3Stage3rdOrderHeun_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:1438
Tempus::RKButcherTableau::c_
Teuchos::SerialDenseVector< int, Scalar > c_
Definition: Tempus_RKButcherTableau.hpp:220
Tempus::Implicit2Stage2ndOrderLobattoB_RKBT
Definition: Tempus_RKButcherTableau.hpp:3097
Tempus::GeneralExplicit_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:467
Tempus::GeneralExplicit_RKBT
General Explicit Runge-Kutta Butcher Tableau.
Definition: Tempus_RKButcherTableau.hpp:428
Tempus::Implicit2Stage4thOrderGauss_RKBT
Definition: Tempus_RKButcherTableau.hpp:2532
Tempus::GeneralExplicit_RKBT::GeneralExplicit_RKBT
GeneralExplicit_RKBT()
Definition: Tempus_RKButcherTableau.hpp:432
Tempus::RKButcherTableau::order
virtual int order() const
Return the order.
Definition: Tempus_RKButcherTableau.hpp:74
Tempus::RKButcherTableau::orderMin
virtual int orderMin() const
Return the minimum order.
Definition: Tempus_RKButcherTableau.hpp:76
Tempus::Implicit3Stage5thOrderRadauB_RKBT::Implicit3Stage5thOrderRadauB_RKBT
Implicit3Stage5thOrderRadauB_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2866
Tempus::Implicit4Stage6thOrderLobattoB_RKBT::Implicit4Stage6thOrderLobattoB_RKBT
Implicit4Stage6thOrderLobattoB_RKBT()
Definition: Tempus_RKButcherTableau.hpp:3199
Tempus::Implicit3Stage6thOrderKuntzmannButcher_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2169
Tempus::SDIRK5Stage4thOrder_RKBT
Definition: Tempus_RKButcherTableau.hpp:3568
Tempus::RKButcherTableau::describe
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Definition: Tempus_RKButcherTableau.hpp:168
Tempus::IRK1StageTheta_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:2328
Tempus::IRK1StageTheta_RKBT
Definition: Tempus_RKButcherTableau.hpp:2303
Tempus::SDIRK2Stage3rdOrder_RKBT::secondOrderLStable_default_
bool secondOrderLStable_default_
Definition: Tempus_RKButcherTableau.hpp:2017
Tempus::ExplicitBogackiShampine32_RKBT
Explicit RK Bogacki-Shampine Butcher Tableau.
Definition: Tempus_RKButcherTableau.hpp:725
Tempus::SDIRK2Stage2ndOrder_RKBT::SDIRK2Stage2ndOrder_RKBT
SDIRK2Stage2ndOrder_RKBT()
Definition: Tempus_RKButcherTableau.hpp:1776
Tempus::GeneralDIRK_RKBT
General Implicit Runge-Kutta Butcher Tableau.
Definition: Tempus_RKButcherTableau.hpp:1607
Tempus::ExplicitTrapezoidal_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:1571
Tempus::SDIRK2Stage3rdOrder_RKBT::gamma_default_
Scalar gamma_default_
Definition: Tempus_RKButcherTableau.hpp:2019
Tempus::Implicit3Stage6thOrderGauss_RKBT
Definition: Tempus_RKButcherTableau.hpp:2582
Tempus::Implicit2Stage2ndOrderLobattoA_RKBT::Implicit2Stage2ndOrderLobattoA_RKBT
Implicit2Stage2ndOrderLobattoA_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2926
Tempus::Implicit4Stage8thOrderKuntzmannButcher_RKBT::Implicit4Stage8thOrderKuntzmannButcher_RKBT
Implicit4Stage8thOrderKuntzmannButcher_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2180
Tempus::GeneralExplicit_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:457
Tempus::EDIRK2Stage3rdOrder_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:2106
Tempus::RKButcherTableau::b
virtual const Teuchos::SerialDenseVector< int, Scalar > & b() const
Return the vector of quadrature weights.
Definition: Tempus_RKButcherTableau.hpp:65
Tempus::RKButcherTableau::initialize
virtual void initialize(const Teuchos::SerialDenseMatrix< int, Scalar > &A, const Teuchos::SerialDenseVector< int, Scalar > &b, const Teuchos::SerialDenseVector< int, Scalar > &c, const int order, const int orderMin, const int orderMax, const std::string &longDescription, bool isEmbedded=false, const Teuchos::SerialDenseVector< int, Scalar > &bstar=Teuchos::SerialDenseVector< int, Scalar >())
Definition: Tempus_RKButcherTableau.hpp:100
Tempus::Implicit3Stage4thOrderLobattoC_RKBT::Implicit3Stage4thOrderLobattoC_RKBT
Implicit3Stage4thOrderLobattoC_RKBT()
Definition: Tempus_RKButcherTableau.hpp:3308
Tempus::BackwardEuler_RKBT
Backward Euler Runge-Kutta Butcher Tableau.
Definition: Tempus_RKButcherTableau.hpp:506
Tempus::Implicit2Stage2ndOrderLobattoB_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3135
Tempus::RKButcherTableau::longDescription_
std::string longDescription_
Definition: Tempus_RKButcherTableau.hpp:226
Tempus::TokensToDoubles
void TokensToDoubles(std::vector< double > &values, const std::vector< std::string > &tokens)
Turn a vector of tokens into a vector of doubles.
Definition: Tempus_String_Utilities.cpp:63
Tempus::SDIRK1Stage1stOrder_RKBT
SDIRK 1 Stage 1st order.
Definition: Tempus_RKButcherTableau.hpp:1687
Tempus::Implicit2Stage2ndOrderLobattoA_RKBT
Definition: Tempus_RKButcherTableau.hpp:2922
Tempus::EDIRK2StageTheta_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:2459
Tempus::SDIRK2Stage3rdOrder_RKBT::thirdOrderAStable_
bool thirdOrderAStable_
Definition: Tempus_RKButcherTableau.hpp:2016
Tempus::SDIRK2Stage2ndOrder_RKBT
SDIRK 2 Stage 2nd order.
Definition: Tempus_RKButcherTableau.hpp:1772
Tempus::SDIRK5Stage4thOrder_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:3677
Tempus::Implicit2Stage2ndOrderLobattoC_RKBT
Definition: Tempus_RKButcherTableau.hpp:3259
Tempus::RKButcherTableau::c
virtual const Teuchos::SerialDenseVector< int, Scalar > & c() const
Return the vector of stage positions.
Definition: Tempus_RKButcherTableau.hpp:71
Tempus::Implicit2Stage3rdOrderRadauB_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2855
Tempus::RKButcherTableau::order_
int order_
Definition: Tempus_RKButcherTableau.hpp:221
Tempus::Implicit1Stage1stOrderRadauA_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2675
Tempus::Explicit3Stage3rdOrder_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:1260
Tempus::Implicit2Stage4thOrderGauss_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2575
Tempus::Explicit3_8Rule_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:985
Tempus::IRK1StageTheta_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2360
Tempus::RKButcherTableau::setDescription
virtual void setDescription(std::string longD)
Definition: Tempus_RKButcherTableau.hpp:190
Tempus::EDIRK2StageTheta_RKBT::theta_default_
Scalar theta_default_
Definition: Tempus_RKButcherTableau.hpp:2482
Tempus::Explicit4Stage4thOrder_RKBT::Explicit4Stage4thOrder_RKBT
Explicit4Stage4thOrder_RKBT()
Definition: Tempus_RKButcherTableau.hpp:633
Tempus::ExplicitTrapezoidal_RKBT
RK Explicit Trapezoidal.
Definition: Tempus_RKButcherTableau.hpp:1530
Tempus::SDIRK2Stage3rdOrder_RKBT
SDIRK 2 Stage 3rd order.
Definition: Tempus_RKButcherTableau.hpp:1892
Tempus::Implicit1Stage1stOrderRadauB_RKBT
Definition: Tempus_RKButcherTableau.hpp:2781
Tempus::SDIRK21_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3824
Tempus::EDIRK2Stage3rdOrder_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:2068
Tempus::EDIRK2StageTheta_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2479
Tempus::BackwardEuler_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:524
Tempus::SDIRK5Stage5thOrder_RKBT::SDIRK5Stage5thOrder_RKBT
SDIRK5Stage5thOrder_RKBT()
Definition: Tempus_RKButcherTableau.hpp:3425
Tempus::RKButcherTableau::rkbtPL_
Teuchos::RCP< Teuchos::ParameterList > rkbtPL_
Definition: Tempus_RKButcherTableau.hpp:232
Tempus::SDIRK2Stage3rdOrder_RKBT::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:1987
Tempus::Explicit3_8Rule_RKBT::Explicit3_8Rule_RKBT
Explicit3_8Rule_RKBT()
Definition: Tempus_RKButcherTableau.hpp:920
Tempus::Implicit3Stage5thOrderRadauA_RKBT::Implicit3Stage5thOrderRadauA_RKBT
Implicit3Stage5thOrderRadauA_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2731
Tempus::Explicit3Stage3rdOrderTVD_RKBT
RK Explicit 3 Stage 3rd order TVD.
Definition: Tempus_RKButcherTableau.hpp:1293
Tempus::SDIRK2Stage3rdOrder_RKBT::setParameterList
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &pList)
Definition: Tempus_RKButcherTableau.hpp:1925
Tempus::EDIRK2StageTheta_RKBT
Definition: Tempus_RKButcherTableau.hpp:2391
Tempus::GeneralDIRK_RKBT::GeneralDIRK_RKBT
GeneralDIRK_RKBT()
Definition: Tempus_RKButcherTableau.hpp:1611
Tempus::SDIRK5Stage4thOrder_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:3595
Tempus::Implicit3Stage6thOrderGauss_RKBT::Implicit3Stage6thOrderGauss_RKBT
Implicit3Stage6thOrderGauss_RKBT()
Definition: Tempus_RKButcherTableau.hpp:2586
Tempus::SDIRK21_RKBT::SDIRK21_RKBT
SDIRK21_RKBT()
Definition: Tempus_RKButcherTableau.hpp:3811
Tempus::RKButcherTableau::getValidParameters
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Tempus_RKButcherTableau.hpp:153
Tempus::RKButcherTableau::isEmbedded_
bool isEmbedded_
Definition: Tempus_RKButcherTableau.hpp:228
Tempus::Explicit2Stage2ndOrderRunge_RKBT
RK Explicit 2 Stage 2nd order by Runge.
Definition: Tempus_RKButcherTableau.hpp:1463
Tempus::RKButcherTableau::isDIRK
virtual bool isDIRK() const
Return true if the RK method is Diagonally Implicit.
Definition: Tempus_RKButcherTableau.hpp:82
Tempus::Implicit4Stage8thOrderKuntzmannButcher_RKBT::description
virtual std::string description() const
Definition: Tempus_RKButcherTableau.hpp:2251
Tempus::Implicit2Stage3rdOrderRadauA_RKBT
Definition: Tempus_RKButcherTableau.hpp:2682
Tempus::RKButcherTableau::orderMin_
int orderMin_
Definition: Tempus_RKButcherTableau.hpp:222