Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_StepperRKBase.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_StepperRKBase_hpp
10#define Tempus_StepperRKBase_hpp
11
12#include "Thyra_VectorBase.hpp"
13
14#include "Tempus_config.hpp"
15#include "Tempus_Stepper.hpp"
19#include "Tempus_Stepper_ErrorNorm.hpp"
20
21
22namespace Tempus {
23
30template<class Scalar>
31class StepperRKBase : virtual public Tempus::Stepper<Scalar>
32{
33
34public:
35
36 virtual Teuchos::RCP<const RKButcherTableau<Scalar> > getTableau() const
37 { return tableau_; }
38
39 virtual Scalar getOrder() const{return getTableau()->order();}
40 virtual Scalar getOrderMin() const{return getTableau()->orderMin();}
41 virtual Scalar getOrderMax() const{return getTableau()->orderMax();}
42 virtual int getNumberOfStages() const {return getTableau()->numStages();}
43
44 virtual int getStageNumber() const { return stageNumber_; }
45 virtual void setStageNumber(int s) { stageNumber_ = s; }
46
47 virtual void setUseEmbedded(bool a)
48 {
49 useEmbedded_ = a;
50 this->setEmbeddedMemory();
51 this->isInitialized_ = false;
52 }
53
54 virtual bool getUseEmbedded() const { return useEmbedded_; }
55
56 virtual void setErrorNorm(const Teuchos::RCP<Stepper_ErrorNorm<Scalar>> &errCalculator=Teuchos::null)
57 {
58 if (errCalculator != Teuchos::null) {
59 stepperErrorNormCalculator_ = errCalculator;
60 }
61 else {
62 auto er = Teuchos::rcp(new Stepper_ErrorNorm<Scalar>());
64 }
65 }
66
67
68 virtual void setAppAction(Teuchos::RCP<StepperRKAppAction<Scalar> > appAction)
69 {
70 if (appAction == Teuchos::null) {
71 // Create default appAction
73 Teuchos::rcp(new StepperRKModifierDefault<Scalar>());
74 } else {
75 stepperRKAppAction_ = appAction;
76 }
77 this->isInitialized_ = false;
78 }
79
80 virtual Teuchos::RCP<StepperRKAppAction<Scalar> > getAppAction() const
81 { return stepperRKAppAction_; }
82
84 virtual void setStepperRKValues(Teuchos::RCP<Teuchos::ParameterList> pl)
85 {
86 if (pl != Teuchos::null) {
87 pl->validateParametersAndSetDefaults(*this->getValidParameters());
88 this->setStepperValues(pl);
89 if (pl->isParameter("Use Embedded"))
90 this->setUseEmbedded(pl->get<bool>("Use Embedded"));
91 }
92 }
93
94 virtual Teuchos::RCP<RKButcherTableau<Scalar> >
95 createTableau(Teuchos::RCP<Teuchos::ParameterList> pl)
96 {
97 TEUCHOS_TEST_FOR_EXCEPTION(pl == Teuchos::null,std::runtime_error,
98 "Error parsing general tableau. ParameterList is null.\n");
99
100 Teuchos::RCP<RKButcherTableau<Scalar> > tableau;
101
102 Teuchos::RCP<Teuchos::ParameterList> tableauPL = sublist(pl,"Tableau",true);
103 std::size_t numStages = 0;
104 int order = tableauPL->get<int>("order");
105 Teuchos::SerialDenseMatrix<int,Scalar> A;
106 Teuchos::SerialDenseVector<int,Scalar> b;
107 Teuchos::SerialDenseVector<int,Scalar> c;
108 Teuchos::SerialDenseVector<int,Scalar> bstar;
109
110 // read in the A matrix
111 {
112 std::vector<std::string> A_row_tokens;
113 Tempus::StringTokenizer(A_row_tokens, tableauPL->get<std::string>("A"),
114 ";",true);
115
116 // this is the only place where numStages is set
117 numStages = A_row_tokens.size();
118
119 // allocate the matrix
120 A.shape(Teuchos::as<int>(numStages),Teuchos::as<int>(numStages));
121
122 // fill the rows
123 for(std::size_t row=0;row<numStages;row++) {
124 // parse the row (tokenize on space)
125 std::vector<std::string> tokens;
126 Tempus::StringTokenizer(tokens,A_row_tokens[row]," ",true);
127
128 std::vector<double> values;
129 Tempus::TokensToDoubles(values,tokens);
130
131 TEUCHOS_TEST_FOR_EXCEPTION(values.size()!=numStages,std::runtime_error,
132 "Error parsing A matrix, wrong number of stages in row "
133 << row << ".\n");
134
135 for(std::size_t col=0;col<numStages;col++)
136 A(row,col) = values[col];
137 }
138 }
139
140 // size b and c vectors
141 b.size(Teuchos::as<int>(numStages));
142 c.size(Teuchos::as<int>(numStages));
143
144 // read in the b vector
145 {
146 std::vector<std::string> tokens;
147 Tempus::StringTokenizer(tokens,tableauPL->get<std::string>("b")," ",true);
148 std::vector<double> values;
149 Tempus::TokensToDoubles(values,tokens);
150
151 TEUCHOS_TEST_FOR_EXCEPTION(values.size()!=numStages,std::runtime_error,
152 "Error parsing b vector, wrong number of stages.\n");
153
154 for(std::size_t i=0;i<numStages;i++)
155 b(i) = values[i];
156 }
157
158 // read in the c vector
159 {
160 std::vector<std::string> tokens;
161 Tempus::StringTokenizer(tokens,tableauPL->get<std::string>("c")," ",true);
162 std::vector<double> values;
163 Tempus::TokensToDoubles(values,tokens);
164
165 TEUCHOS_TEST_FOR_EXCEPTION(values.size()!=numStages,std::runtime_error,
166 "Error parsing c vector, wrong number of stages.\n");
167
168 for(std::size_t i=0;i<numStages;i++)
169 c(i) = values[i];
170 }
171
172 if (tableauPL->isParameter("bstar") &&
173 tableauPL->get<std::string>("bstar") != "") {
174 bstar.size(Teuchos::as<int>(numStages));
175 // read in the bstar vector
176 {
177 std::vector<std::string> tokens;
179 tokens, tableauPL->get<std::string>("bstar"), " ", true);
180 std::vector<double> values;
181 Tempus::TokensToDoubles(values,tokens);
182
183 TEUCHOS_TEST_FOR_EXCEPTION(values.size()!=numStages,std::runtime_error,
184 "Error parsing bstar vector, wrong number of stages.\n"
185 " Number of RK stages = " << numStages << "\n"
186 " Number of bstar values = " << values.size() << "\n");
187
188 for(std::size_t i=0;i<numStages;i++)
189 bstar(i) = values[i];
190 }
191 tableau = rcp(new RKButcherTableau<Scalar>(
192 "From ParameterList",A,b,c,order,order,order,bstar));
193 } else {
194 tableau = rcp(new RKButcherTableau<Scalar>(
195 "From ParameterList",A,b,c,order,order,order));
196 }
197 return tableau;
198 }
199
200
201protected:
202
203 virtual void setEmbeddedMemory() {}
204
205 Teuchos::RCP<RKButcherTableau<Scalar> > tableau_;
206
207 // For Embedded RK
209 Teuchos::RCP<Thyra::VectorBase<Scalar> > ee_;
210 Teuchos::RCP<Thyra::VectorBase<Scalar> > abs_u0;
211 Teuchos::RCP<Thyra::VectorBase<Scalar> > abs_u;
212 Teuchos::RCP<Thyra::VectorBase<Scalar> > sc;
213
214 Teuchos::RCP<Stepper_ErrorNorm<Scalar>> stepperErrorNormCalculator_;
215
218 Teuchos::RCP<StepperRKAppAction<Scalar> > stepperRKAppAction_;
219
220};
221
222} // namespace Tempus
223
224#endif // Tempus_StepperRKBase_hpp
Base class for Runge-Kutta methods, ExplicitRK, DIRK and IMEX.
Teuchos::RCP< Thyra::VectorBase< Scalar > > ee_
virtual int getStageNumber() const
virtual void setStageNumber(int s)
Teuchos::RCP< Thyra::VectorBase< Scalar > > abs_u0
virtual Scalar getOrderMax() const
Teuchos::RCP< StepperRKAppAction< Scalar > > stepperRKAppAction_
virtual void setUseEmbedded(bool a)
virtual int getNumberOfStages() const
Teuchos::RCP< Thyra::VectorBase< Scalar > > sc
virtual void setStepperRKValues(Teuchos::RCP< Teuchos::ParameterList > pl)
Set StepperRK member data from the ParameterList.
virtual Scalar getOrderMin() const
virtual Teuchos::RCP< RKButcherTableau< Scalar > > createTableau(Teuchos::RCP< Teuchos::ParameterList > pl)
virtual void setErrorNorm(const Teuchos::RCP< Stepper_ErrorNorm< Scalar > > &errCalculator=Teuchos::null)
virtual Scalar getOrder() const
virtual void setAppAction(Teuchos::RCP< StepperRKAppAction< Scalar > > appAction)
virtual Teuchos::RCP< const RKButcherTableau< Scalar > > getTableau() const
virtual bool getUseEmbedded() const
int stageNumber_
The current Runge-Kutta stage number, {0,...,s-1}. -1 indicates outside stage loop.
Teuchos::RCP< Thyra::VectorBase< Scalar > > abs_u
Teuchos::RCP< Stepper_ErrorNorm< Scalar > > stepperErrorNormCalculator_
virtual Teuchos::RCP< StepperRKAppAction< Scalar > > getAppAction() const
Teuchos::RCP< RKButcherTableau< Scalar > > tableau_
Thyra Base interface for time steppers.
bool isInitialized_
True if stepper's member data is initialized.
void setStepperValues(const Teuchos::RCP< Teuchos::ParameterList > pl)
Set Stepper member data from ParameterList.
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void TokensToDoubles(std::vector< double > &values, const std::vector< std::string > &tokens)
Turn a vector of tokens into a vector of doubles.
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.