Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 
19 
20 #include "Tempus_config.hpp"
22 
23 
24 namespace Tempus {
25 
26 
44 template<class Scalar>
46  virtual public Teuchos::Describable,
47  virtual public Teuchos::VerboseObject<RKButcherTableau<Scalar> >
48 {
49  public:
50 
52  std::string stepperType,
56  const int order,
57  const int orderMin,
58  const int orderMax,
61  bool checkC = true)
62  : description_(stepperType)
63  {
64  const int numStages = A.numRows();
68  TEUCHOS_ASSERT( order > 0 );
69  A_ = A;
70  b_ = b;
71  c_ = c;
72  order_ = order;
75  this->set_isImplicit();
76  this->set_isDIRK();
77 
78  // Consistency check on b
80  Scalar sumb = ST::zero();
81  for (size_t i = 0; i < this->numStages(); i++) sumb += b_(i);
82  TEUCHOS_TEST_FOR_EXCEPTION( std::abs(ST::one()-sumb) > 1.0e-08,
83  std::runtime_error,
84  "Error - Butcher Tableau b fails to satisfy Sum(b_i) = 1.\n"
85  << " Sum(b_i) = " << sumb << "\n");
86 
87  // Consistency check on c
88  if (checkC) {
89  for (size_t i = 0; i < this->numStages(); i++) {
90  Scalar sumai = ST::zero();
91  for (size_t j = 0; j < this->numStages(); j++) sumai += A_(i,j);
92  bool failed = false;
93  if (std::abs(sumai) > 1.0e-08)
94  failed = (std::abs((sumai-c_(i))/sumai) > 1.0e-08);
95  else
96  failed = (std::abs(c_(i)) > 1.0e-08);
97 
98  TEUCHOS_TEST_FOR_EXCEPTION( failed, std::runtime_error,
99  "Error - Butcher Tableau fails to satisfy c_i = Sum_j(a_ij).\n"
100  << " Stage i = " << i << "\n"
101  << " c_i = " << c_(i) << "\n"
102  << " Sum_j(a_ij) = " << sumai << "\n"
103  << " This may be OK as some valid tableaus do not satisfy\n"
104  << " this condition. If so, construct this RKButcherTableau\n"
105  << " with checkC = false.\n");
106  }
107  }
108 
109  if ( bstar.length() > 0 ) {
111  isEmbedded_ = true;
112  } else {
113  isEmbedded_ = false;
114  }
115  bstar_ = bstar;
116  }
117 
119  virtual std::size_t numStages() const { return A_.numRows(); }
122  { return A_; }
125  { return b_; }
128  { return bstar_ ; }
131  { return c_; }
133  virtual int order() const { return order_; }
135  virtual int orderMin() const { return orderMin_; }
137  virtual int orderMax() const { return orderMax_; }
139  virtual bool isImplicit() const { return isImplicit_; }
141  virtual bool isDIRK() const { return isDIRK_; }
143  virtual bool isEmbedded() const { return isEmbedded_; }
145  virtual bool isTVD() const { return isTVD_; }
147  virtual Scalar getTVDCoeff() const { return tvdCoeff_; }
149  virtual void setTVDCoeff(const Scalar a) { tvdCoeff_ = a; }
150  virtual void setTVD(bool a) { isTVD_ = a; }
151  virtual void setEmbedded(bool a) { isEmbedded_ = a; }
152 
153 
154  /* \brief Redefined from Teuchos::Describable */
156  virtual std::string description() const { return description_; }
157 
158  virtual void describe( Teuchos::FancyOStream &out,
159  const Teuchos::EVerbosityLevel verbLevel) const
160  {
161  out.setOutputToRootOnly(0);
162 
163  if (verbLevel != Teuchos::VERB_NONE) {
164  out << this->description() << std::endl;
165  out << "number of Stages = " << this->numStages() << std::endl;
166  out << "A = " << printMat(this->A()) << std::endl;
167  out << "b = " << printMat(this->b()) << std::endl;
168  out << "c = " << printMat(this->c()) << std::endl;
169  out << "bstar = " << printMat(this->bstar()) << std::endl;
170  out << "order = " << this->order() << std::endl;
171  out << "orderMin = " << this->orderMin() << std::endl;
172  out << "orderMax = " << this->orderMax() << std::endl;
173  out << "isImplicit = " << this->isImplicit() << std::endl;
174  out << "isDIRK = " << this->isDIRK() << std::endl;
175  out << "isEmbedded = " << this->isEmbedded() << std::endl;
176  if (this->isTVD())
177  out << "TVD Coeff = " << this->getTVDCoeff() << std::endl;
178  }
179  }
181 
182  bool operator == (const RKButcherTableau & t) const {
183  const Scalar relTol = 1.0e-15;
184  if ( A_->numRows() != t.A_->numRows() ||
185  A_->numCols() != t.A_->numCols() ) {
186  return false;
187  } else {
188  int i, j;
189  for(i = 0; i < A_->numRows(); i++) {
190  for(j = 0; j < A_->numCols(); j++) {
191  if(std::abs((t.A_(i,j) - A_(i,j))/A_(i,j)) > relTol) return false;
192  }
193  }
194  }
195 
196  if ( b_->length() != t.b_->length() ||
197  b_->length() != t.c_->length() ||
198  b_->length() != t.bstar_->length() ) {
199  return false;
200  } else {
201  int i;
202  for(i = 0; i < A_->numRows(); i++) {
203  if(std::abs((t.b_(i) - b_(i))/b_(i)) > relTol) return false;
204  if(std::abs((t.c_(i) - c_(i))/c_(i)) > relTol) return false;
205  if(std::abs((t.bstar_(i) - bstar_(i))/bstar_(i)) > relTol) return false;
206  }
207  }
208  return true;
209  }
210 
211  bool operator != (const RKButcherTableau & t) const {
212  return !((*this) == t);
213  }
214 
215  private:
216 
218 
219  protected:
220 
221  void set_isImplicit() {
222  isImplicit_ = false;
223  for (size_t i = 0; i < this->numStages(); i++)
224  for (size_t j = i; j < this->numStages(); j++)
225  if (A_(i,j) != 0.0) isImplicit_ = true;
226  }
227 
229  void set_isDIRK() {
230  isDIRK_ = true;
231  bool nonZero = false;
232  for (size_t i = 0; i < this->numStages(); i++) {
233  if (A_(i,i) != 0.0) nonZero = true;
234  for (size_t j = i+1; j < this->numStages(); j++)
235  if (A_(i,j) != 0.0) isDIRK_ = false;
236  }
237  if (nonZero == false) isDIRK_ = false;
238  }
239 
240  std::string description_;
241 
245  int order_;
249  bool isDIRK_;
251  bool isTVD_ = false;
253  Scalar tvdCoeff_ = 0;
254 };
255 
256 
257 } // namespace Tempus
258 
259 
260 #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)
bool operator!=(const RKButcherTableau &t) const
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
virtual const Teuchos::SerialDenseVector< int, Scalar > & bstar() const
Return the vector of quadrature weights for embedded methods.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void setTVDCoeff(const Scalar a)
set TVD coefficient of RK method
virtual bool isTVD() const
Return true if the RK method is TVD.
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.
virtual Scalar getTVDCoeff() const
Return TVD coefficient of RK method.
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
OrdinalType length() const
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
OrdinalType numCols() const
virtual std::size_t numStages() const
Return the number of stages.
Teuchos::SerialDenseVector< int, Scalar > c_
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
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.
bool operator==(const RKButcherTableau &t) const
OrdinalType numRows() const