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