9#ifndef Tempus_RKButcherTableau_hpp
10#define Tempus_RKButcherTableau_hpp
14#pragma clang system_header
17#include "Teuchos_SerialDenseMatrix.hpp"
18#include "Teuchos_SerialDenseVector.hpp"
20#include "Tempus_config.hpp"
46 virtual public Teuchos::Describable,
47 virtual public Teuchos::VerboseObject<RKButcherTableau<Scalar> >
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,
59 const Teuchos::SerialDenseVector<int,Scalar>&
60 bstar = Teuchos::SerialDenseVector<int,Scalar>(),
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 );
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,
84 "Error - Butcher Tableau b fails to satisfy Sum(b_i) = 1.\n"
85 <<
" Sum(b_i) = " << sumb <<
"\n");
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);
93 if (std::abs(sumai) > 1.0e-08)
94 failed = (std::abs((sumai-
c_(i))/sumai) > 1.0e-08);
96 failed = (std::abs(
c_(i)) > 1.0e-08);
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");
109 if (
bstar.length() > 0 ) {
121 virtual const Teuchos::SerialDenseMatrix<int,Scalar>&
A()
const
124 virtual const Teuchos::SerialDenseVector<int,Scalar>&
b()
const
127 virtual const Teuchos::SerialDenseVector<int,Scalar>&
bstar()
const
130 virtual const Teuchos::SerialDenseVector<int,Scalar>&
c()
const
159 const Teuchos::EVerbosityLevel verbLevel)
const
161 out.setOutputToRootOnly(0);
163 if (verbLevel != Teuchos::VERB_NONE) {
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;
177 out <<
"TVD Coeff = " << this->
getTVDCoeff() << std::endl;
183 const Scalar relTol = 1.0e-15;
184 if (
A_->numRows() != t.
A_->numRows() ||
185 A_->numCols() != t.
A_->numCols() ) {
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;
196 if (
b_->length() != t.
b_->length() ||
197 b_->length() != t.
c_->length() ||
198 b_->length() != t.
bstar_->length() ) {
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;
212 return !((*this) == t);
223 for (
size_t i = 0; i < this->
numStages(); i++)
224 for (
size_t j = i; j < this->
numStages(); j++)
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++)
237 if (nonZero ==
false)
isDIRK_ =
false;
242 Teuchos::SerialDenseMatrix<int,Scalar>
A_;
243 Teuchos::SerialDenseVector<int,Scalar>
b_;
244 Teuchos::SerialDenseVector<int,Scalar>
c_;
252 Teuchos::SerialDenseVector<int,Scalar>
bstar_;
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.
virtual void setTVD(bool a)
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
virtual void setEmbedded(bool a)
Teuchos::SerialDenseVector< int, Scalar > bstar_