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 
18 
19 #include "Teuchos_Assert.hpp"
20 #include "Teuchos_as.hpp"
21 #include "Teuchos_Describable.hpp"
22 #include "Teuchos_VerboseObject.hpp"
23 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
24 #include "Teuchos_SerialDenseMatrix.hpp"
25 #include "Teuchos_SerialDenseVector.hpp"
26 #include "Thyra_MultiVectorStdOps.hpp"
27 
28 
29 namespace Tempus {
30 
31 
32 /** \brief Runge-Kutta methods.
33  *
34  * This base class specifies the Butcher tableau which defines the
35  * Runge-Kutta (RK) method. Both explicit and implicit RK methods
36  * can be specified here, and of arbitrary number of stages and orders.
37  * Embedded methods are also supported.
38  *
39  * Since this is a generic RK class, no low-storage methods are
40  * incorporated here, however any RK method with a Butcher tableau
41  * can be created with the base class.
42  *
43  * There are over 40 derived RK methods that have been implemented,
44  * ranging from first order and eight order, and from single stage
45  * to 5 stages.
46  *
47  * This class was taken and modified from Rythmos' RKButcherTableau class.
48  */
49 template<class Scalar>
51  virtual public Teuchos::Describable,
52  virtual public Teuchos::VerboseObject<RKButcherTableau<Scalar> >
53 {
54  public:
55 
57  std::string stepperType,
58  const Teuchos::SerialDenseMatrix<int,Scalar>& A,
59  const Teuchos::SerialDenseVector<int,Scalar>& b,
60  const Teuchos::SerialDenseVector<int,Scalar>& c,
61  const int order,
62  const int orderMin,
63  const int orderMax,
64  const Teuchos::SerialDenseVector<int,Scalar>&
65  bstar = Teuchos::SerialDenseVector<int,Scalar>(),
66  bool checkC = true)
67  : description_(stepperType)
68  {
69  const int numStages = A.numRows();
70  TEUCHOS_ASSERT_EQUALITY( A.numCols(), numStages );
71  TEUCHOS_ASSERT_EQUALITY( b.length(), numStages );
72  TEUCHOS_ASSERT_EQUALITY( c.length(), numStages );
73  TEUCHOS_ASSERT( order > 0 );
74  A_ = A;
75  b_ = b;
76  c_ = c;
77  order_ = order;
80  this->set_isImplicit();
81  this->set_isDIRK();
82 
83  // Consistency check on b
84  typedef Teuchos::ScalarTraits<Scalar> ST;
85  Scalar sumb = ST::zero();
86  for (size_t i = 0; i < this->numStages(); i++) sumb += b_(i);
87  TEUCHOS_TEST_FOR_EXCEPTION( std::abs(ST::one()-sumb) > 1.0e-08,
88  std::runtime_error,
89  "Error - Butcher Tableau b fails to satisfy Sum(b_i) = 1.\n"
90  << " Sum(b_i) = " << sumb << "\n");
91 
92  // Consistency check on c
93  if (checkC) {
94  for (size_t i = 0; i < this->numStages(); i++) {
95  Scalar sumai = ST::zero();
96  for (size_t j = 0; j < this->numStages(); j++) sumai += A_(i,j);
97  bool failed = false;
98  if (std::abs(sumai) > 1.0e-08)
99  failed = (std::abs((sumai-c_(i))/sumai) > 1.0e-08);
100  else
101  failed = (std::abs(c_(i)) > 1.0e-08);
102 
103  TEUCHOS_TEST_FOR_EXCEPTION( failed, std::runtime_error,
104  "Error - Butcher Tableau fails to satisfy c_i = Sum_j(a_ij).\n"
105  << " Stage i = " << i << "\n"
106  << " c_i = " << c_(i) << "\n"
107  << " Sum_j(a_ij) = " << sumai << "\n"
108  << " This may be OK as some valid tableaus do not satisfy\n"
109  << " this condition. If so, construct this RKButcherTableau\n"
110  << " with checkC = false.\n");
111  }
112  }
113 
114  if ( bstar.length() > 0 ) {
115  TEUCHOS_ASSERT_EQUALITY( bstar.length(), numStages );
116  isEmbedded_ = true;
117  } else {
118  isEmbedded_ = false;
119  }
120  bstar_ = bstar;
121  }
122 
123  /** \brief Return the number of stages */
124  virtual std::size_t numStages() const { return A_.numRows(); }
125  /** \brief Return the matrix coefficients */
126  virtual const Teuchos::SerialDenseMatrix<int,Scalar>& A() const
127  { return A_; }
128  /** \brief Return the vector of quadrature weights */
129  virtual const Teuchos::SerialDenseVector<int,Scalar>& b() const
130  { return b_; }
131  /** \brief Return the vector of quadrature weights for embedded methods */
132  virtual const Teuchos::SerialDenseVector<int,Scalar>& bstar() const
133  { return bstar_ ; }
134  /** \brief Return the vector of stage positions */
135  virtual const Teuchos::SerialDenseVector<int,Scalar>& c() const
136  { return c_; }
137  /** \brief Return the order */
138  virtual int order() const { return order_; }
139  /** \brief Return the minimum order */
140  virtual int orderMin() const { return orderMin_; }
141  /** \brief Return the maximum order */
142  virtual int orderMax() const { return orderMax_; }
143  /** \brief Return true if the RK method is implicit */
144  virtual bool isImplicit() const { return isImplicit_; }
145  /** \brief Return true if the RK method is Diagonally Implicit */
146  virtual bool isDIRK() const { return isDIRK_; }
147  /** \brief Return true if the RK method has embedded capabilities */
148  virtual bool isEmbedded() const { return isEmbedded_; }
149 
150  /* \brief Redefined from Teuchos::Describable */
151  //@{
152  virtual std::string description() const { return description_; }
153 
154  virtual void describe( Teuchos::FancyOStream &out,
155  const Teuchos::EVerbosityLevel verbLevel) const
156  {
157  if (verbLevel != Teuchos::VERB_NONE) {
158  out << this->description() << std::endl;
159  out << "number of Stages = " << this->numStages() << std::endl;
160  out << "A = " << printMat(this->A()) << std::endl;
161  out << "b = " << printMat(this->b()) << std::endl;
162  out << "c = " << printMat(this->c()) << std::endl;
163  out << "bstar = " << printMat(this->bstar()) << std::endl;
164  out << "order = " << this->order() << std::endl;
165  out << "orderMin = " << this->orderMin() << std::endl;
166  out << "orderMax = " << this->orderMax() << std::endl;
167  out << "isImplicit = " << this->isImplicit() << std::endl;
168  out << "isDIRK = " << this->isDIRK() << std::endl;
169  out << "isEmbedded = " << this->isEmbedded() << std::endl;
170  }
171  }
172  //@}
173 
174  private:
175 
177 
178  protected:
179 
180  void set_isImplicit() {
181  isImplicit_ = false;
182  for (size_t i = 0; i < this->numStages(); i++)
183  for (size_t j = i; j < this->numStages(); j++)
184  if (A_(i,j) != 0.0) isImplicit_ = true;
185  }
186 
187  /// DIRK is defined as if a_ij = 0 for j>i and a_ii != 0 for at least one i.
188  void set_isDIRK() {
189  isDIRK_ = true;
190  bool nonZero = false;
191  for (size_t i = 0; i < this->numStages(); i++) {
192  if (A_(i,i) != 0.0) nonZero = true;
193  for (size_t j = i+1; j < this->numStages(); j++)
194  if (A_(i,j) != 0.0) isDIRK_ = false;
195  }
196  if (nonZero == false) isDIRK_ = false;
197  }
198 
199  std::string description_;
200 
201  Teuchos::SerialDenseMatrix<int,Scalar> A_;
202  Teuchos::SerialDenseVector<int,Scalar> b_;
203  Teuchos::SerialDenseVector<int,Scalar> c_;
204  int order_;
208  bool isDIRK_;
210  Teuchos::SerialDenseVector<int,Scalar> bstar_;
211 };
212 
213 
214 } // namespace Tempus
215 
216 
217 #endif // Tempus_RKButcherTableau_hpp
Teuchos::SerialDenseMatrix< int, Scalar > A_
virtual const Teuchos::SerialDenseVector< int, Scalar > & b() const
Return the vector of quadrature weights.
virtual const Teuchos::SerialDenseVector< int, Scalar > & c() const
Return the vector of stage positions.
RKButcherTableau(std::string stepperType, 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 Teuchos::SerialDenseVector< int, Scalar > &bstar=Teuchos::SerialDenseVector< int, Scalar >(), bool checkC=true)
virtual const Teuchos::SerialDenseVector< int, Scalar > & bstar() const
Return the vector of quadrature weights for embedded methods.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Teuchos::SerialDenseVector< int, Scalar > bstar_
virtual bool isImplicit() const
Return true if the RK method is implicit.
virtual bool isEmbedded() const
Return true if the RK method has embedded capabilities.
virtual bool isDIRK() const
Return true if the RK method is Diagonally Implicit.
virtual int orderMax() const
Return the maximum order.
void set_isDIRK()
DIRK is defined as if a_ij = 0 for j&gt;i and a_ii != 0 for at least one i.
Teuchos::SerialDenseVector< int, Scalar > b_
virtual std::string description() const
virtual std::size_t numStages() const
Return the number of stages.
Teuchos::SerialDenseVector< int, Scalar > c_
virtual const Teuchos::SerialDenseMatrix< int, Scalar > & A() const
Return the matrix coefficients.
virtual int order() const
Return the order.
virtual int orderMin() const
Return the minimum order.