Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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 "Teuchos_SerialDenseMatrix.hpp"
18#include "Teuchos_SerialDenseVector.hpp"
19
20#include "Tempus_config.hpp"
22
23
24namespace Tempus {
25
26
44template<class Scalar>
46 virtual public Teuchos::Describable,
47 virtual public Teuchos::VerboseObject<RKButcherTableau<Scalar> >
48{
49 public:
50
52 std::string stepperType,
53 const Teuchos::SerialDenseMatrix<int,Scalar>& A,
54 const Teuchos::SerialDenseVector<int,Scalar>& b,
55 const Teuchos::SerialDenseVector<int,Scalar>& c,
56 const int order,
57 const int orderMin,
58 const int orderMax,
59 const Teuchos::SerialDenseVector<int,Scalar>&
60 bstar = Teuchos::SerialDenseVector<int,Scalar>(),
61 bool checkC = true)
62 : description_(stepperType)
63 {
64 const int numStages = A.numRows();
65 TEUCHOS_ASSERT_EQUALITY( A.numCols(), numStages );
66 TEUCHOS_ASSERT_EQUALITY( b.length(), numStages );
67 TEUCHOS_ASSERT_EQUALITY( c.length(), numStages );
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
79 typedef Teuchos::ScalarTraits<Scalar> ST;
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 ) {
110 TEUCHOS_ASSERT_EQUALITY( bstar.length(), numStages );
111 isEmbedded_ = true;
112 } else {
113 isEmbedded_ = false;
114 }
115 bstar_ = bstar;
116 }
117
119 virtual std::size_t numStages() const { return A_.numRows(); }
121 virtual const Teuchos::SerialDenseMatrix<int,Scalar>& A() const
122 { return A_; }
124 virtual const Teuchos::SerialDenseVector<int,Scalar>& b() const
125 { return b_; }
127 virtual const Teuchos::SerialDenseVector<int,Scalar>& bstar() const
128 { return bstar_ ; }
130 virtual const Teuchos::SerialDenseVector<int,Scalar>& c() const
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
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
242 Teuchos::SerialDenseMatrix<int,Scalar> A_;
243 Teuchos::SerialDenseVector<int,Scalar> b_;
244 Teuchos::SerialDenseVector<int,Scalar> c_;
251 bool isTVD_ = false;
252 Teuchos::SerialDenseVector<int,Scalar> bstar_;
253 Scalar tvdCoeff_ = 0;
254};
255
256
257} // namespace Tempus
258
259
260#endif // Tempus_RKButcherTableau_hpp
bool operator!=(const RKButcherTableau &t) const
virtual bool isImplicit() const
Return true if the RK method is implicit.
Teuchos::SerialDenseVector< int, Scalar > b_
virtual const Teuchos::SerialDenseMatrix< int, Scalar > & A() const
Return the matrix coefficients.
virtual int orderMin() const
Return the minimum order.
virtual int orderMax() const
Return the maximum order.
virtual void setTVDCoeff(const Scalar a)
set TVD coefficient of RK method
virtual bool isTVD() const
Return true if the RK method is TVD.
bool operator==(const RKButcherTableau &t) const
virtual Scalar getTVDCoeff() const
Return TVD coefficient of RK method.
Teuchos::SerialDenseMatrix< int, Scalar > A_
virtual const Teuchos::SerialDenseVector< int, Scalar > & c() const
Return the vector of stage positions.
virtual int order() const
Return the order.
virtual const Teuchos::SerialDenseVector< int, Scalar > & b() const
Return the vector of quadrature weights.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual std::size_t numStages() const
Return the number of stages.
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 bool isEmbedded() const
Return true if the RK method has embedded capabilities.
void set_isDIRK()
DIRK is defined as if a_ij = 0 for j>i and a_ii != 0 for at least one i.
virtual const Teuchos::SerialDenseVector< int, Scalar > & bstar() const
Return the vector of quadrature weights for embedded methods.
Teuchos::SerialDenseVector< int, Scalar > c_
virtual bool isDIRK() const
Return true if the RK method is Diagonally Implicit.
virtual std::string description() const
Teuchos::SerialDenseVector< int, Scalar > bstar_