29 #ifndef Rythmos_ExplicitRK_STEPPER_DEF_H 30 #define Rythmos_ExplicitRK_STEPPER_DEF_H 32 #include "Rythmos_ExplicitRKStepper_decl.hpp" 34 #include "Rythmos_RKButcherTableau.hpp" 35 #include "Rythmos_RKButcherTableauHelpers.hpp" 36 #include "Rythmos_RKButcherTableauBuilder.hpp" 37 #include "Rythmos_StepperHelpers.hpp" 38 #include "Rythmos_LinearInterpolator.hpp" 39 #include "Rythmos_InterpolatorBaseHelpers.hpp" 41 #include "Teuchos_VerboseObjectParameterListHelpers.hpp" 43 #include "Thyra_ModelEvaluatorHelpers.hpp" 44 #include "Thyra_MultiVectorStdOps.hpp" 45 #include "Thyra_VectorStdOps.hpp" 51 template<
class Scalar>
52 RCP<ExplicitRKStepper<Scalar> > explicitRKStepper()
54 RCP<ExplicitRKStepper<Scalar> > stepper = rcp(
new ExplicitRKStepper<Scalar>());
58 template<
class Scalar>
59 RCP<ExplicitRKStepper<Scalar> > explicitRKStepper(
60 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model
63 RCP<RKButcherTableauBase<Scalar> > rkbt = createRKBT<Scalar>(
"Explicit 4 Stage");
65 RCP<ExplicitRKStepper<Scalar> > stepper = rcp(
new ExplicitRKStepper<Scalar>());
66 stepper->setModel(model);
67 stepper->setRKButcherTableau(rkbt);
71 template<
class Scalar>
72 RCP<ExplicitRKStepper<Scalar> > explicitRKStepper(
73 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
74 const RCP<
const RKButcherTableauBase<Scalar> >& rkbt
77 RCP<ExplicitRKStepper<Scalar> > stepper = rcp(
new ExplicitRKStepper<Scalar>());
78 stepper->setModel(model);
79 stepper->setRKButcherTableau(rkbt);
83 template<
class Scalar>
86 this->defaultInitializeAll_();
87 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
89 erkButcherTableau_ = rKButcherTableau<Scalar>();
93 template<
class Scalar>
96 model_ = Teuchos::null;
97 solution_vector_ = Teuchos::null;
98 solution_vector_old_ = Teuchos::null;
100 ktemp_vector_ = Teuchos::null;
102 erkButcherTableau_ = Teuchos::null;
107 parameterList_ = Teuchos::null;
108 isInitialized_ =
false;
109 haveInitialCondition_ =
false;
112 template<
class Scalar>
115 TEUCHOS_ASSERT( !is_null(rkbt) );
116 validateERKButcherTableau(*rkbt);
117 int numStages_old = erkButcherTableau_->numStages();
118 int numStages_new = rkbt->numStages();
119 TEUCHOS_TEST_FOR_EXCEPTION( numStages_new == 0, std::logic_error,
120 "Error! The Runge-Kutta Butcher tableau has no stages!" 122 if (!is_null(model_)) {
123 int numNewStages = numStages_new - numStages_old;
124 if ( numNewStages > 0 ) {
125 k_vector_.reserve(numStages_new);
126 for (
int i=0 ; i<numNewStages ; ++i) {
127 k_vector_.push_back(Thyra::createMember(model_->get_f_space()));
131 erkButcherTableau_ = rkbt;
134 template<
class Scalar>
137 return erkButcherTableau_;
140 template<
class Scalar>
143 if (!isInitialized_) {
144 TEUCHOS_ASSERT( !is_null(model_) );
145 TEUCHOS_ASSERT( !is_null(erkButcherTableau_) );
146 TEUCHOS_ASSERT( haveInitialCondition_ );
147 TEUCHOS_TEST_FOR_EXCEPTION( erkButcherTableau_->numStages() == 0, std::logic_error,
148 "Error! The Runge-Kutta Butcher tableau has no stages!" 150 ktemp_vector_ = Thyra::createMember(model_->get_f_space());
152 int numStages = erkButcherTableau_->numStages();
153 k_vector_.reserve(numStages);
154 for (
int i=0 ; i<numStages ; ++i) {
155 k_vector_.push_back(Thyra::createMember(model_->get_f_space()));
158 #ifdef HAVE_RYTHMOS_DEBUG 159 THYRA_ASSERT_VEC_SPACES(
160 "Rythmos::ExplicitRKStepper::initialize_(...)",
161 *solution_vector_->space(), *model_->get_x_space() );
162 #endif // HAVE_RYTHMOS_DEBUG 163 isInitialized_ =
true;
167 template<
class Scalar>
172 template<
class Scalar>
175 TEUCHOS_ASSERT( !is_null(model_) );
176 return(model_->get_x_space());
179 template<
class Scalar>
182 typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag;
184 #ifdef HAVE_RYTHMOS_DEBUG 185 TEUCHOS_TEST_FOR_EXCEPTION( flag == STEP_TYPE_VARIABLE, std::logic_error,
186 "Error! ExplicitRKStepper does not support variable time steps at this time." 188 #endif // HAVE_RYTHMOS_DEBUG 189 if ((flag == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
190 return(Scalar(-ST::one()));
193 V_V(solution_vector_old_.ptr(), *solution_vector_);
198 int stages = erkButcherTableau_->numStages();
199 Teuchos::SerialDenseMatrix<int,Scalar> A = erkButcherTableau_->A();
200 Teuchos::SerialDenseVector<int,Scalar> b = erkButcherTableau_->b();
201 Teuchos::SerialDenseVector<int,Scalar> c = erkButcherTableau_->c();
203 for (
int s=0 ; s < stages ; ++s) {
204 Thyra::assign(ktemp_vector_.ptr(), *solution_vector_);
205 for (
int j=0 ; j < s ; ++j) {
206 if (A(s,j) != ST::zero()) {
207 Thyra::Vp_StV(ktemp_vector_.ptr(), A(s,j), *k_vector_[j]);
210 TScalarMag ts = t_ + c(s)*dt;
211 eval_model_explicit<Scalar>(*model_,basePoint_,*ktemp_vector_,ts,Teuchos::outArg(*k_vector_[s]));
212 Thyra::Vt_S(k_vector_[s].ptr(),dt);
215 for (
int s=0 ; s < stages ; ++s) {
216 if (b(s) != ST::zero()) {
217 Thyra::Vp_StV(solution_vector_.ptr(), b(s), *k_vector_[s]);
229 template<
class Scalar>
234 if (!haveInitialCondition_) {
235 stepStatus.
stepStatus = STEP_STATUS_UNINITIALIZED;
236 }
else if (numSteps_ == 0) {
238 stepStatus.
order = erkButcherTableau_->order();
239 stepStatus.
time = t_;
240 stepStatus.
solution = solution_vector_;
242 stepStatus.
stepStatus = STEP_STATUS_CONVERGED;
244 stepStatus.
order = erkButcherTableau_->order();
245 stepStatus.
time = t_;
246 stepStatus.
solution = solution_vector_;
252 template<
class Scalar>
254 Teuchos::FancyOStream &out,
255 const Teuchos::EVerbosityLevel verbLevel
258 if ( (static_cast<int>(verbLevel) == static_cast<int>(Teuchos::VERB_DEFAULT) ) ||
259 (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW) )
261 out << this->description() <<
"::describe" << std::endl;
262 out <<
"model = " << model_->description() << std::endl;
263 out << erkButcherTableau_->numStages() <<
" stage Explicit RK method" << std::endl;
264 }
else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)) {
265 out <<
"solution_vector = " << std::endl;
266 out << Teuchos::describe(*solution_vector_,verbLevel);
267 }
else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM)) {
268 }
else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_HIGH)) {
269 out <<
"model = " << std::endl;
270 out << Teuchos::describe(*model_,verbLevel);
271 int stages = erkButcherTableau_->numStages();
272 for (
int i=0 ; i<stages ; ++i) {
273 out <<
"k_vector[" << i <<
"] = " << std::endl;
274 out << Teuchos::describe(*k_vector_[i],verbLevel);
276 out <<
"ktemp_vector = " << std::endl;
277 out << Teuchos::describe(*ktemp_vector_,verbLevel);
278 out <<
"ERK Butcher Tableau A matrix: " << erkButcherTableau_->A() << std::endl;
279 out <<
"ERK Butcher Tableau b vector: " << erkButcherTableau_->b() << std::endl;
280 out <<
"ERK Butcher Tableau c vector: " << erkButcherTableau_->c() << std::endl;
281 out <<
"t = " << t_ << std::endl;
285 template<
class Scalar>
287 const Array<Scalar>& time_vec
288 ,
const Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >& x_vec
289 ,
const Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >& xdot_vec
292 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Error, addPoints is not implemented for ExplicitRKStepper at this time.\n");
295 template<
class Scalar>
298 if (!haveInitialCondition_) {
299 return(invalidTimeRange<Scalar>());
305 template<
class Scalar>
307 const Array<Scalar>& time_vec,
308 Array<RCP<
const Thyra::VectorBase<Scalar> > >* x_vec,
309 Array<RCP<
const Thyra::VectorBase<Scalar> > >* xdot_vec,
310 Array<ScalarMag>* accuracy_vec
313 TEUCHOS_ASSERT( haveInitialCondition_ );
314 using Teuchos::constOptInArg;
316 defaultGetPoints<Scalar>(
317 t_old_, constOptInArg(*solution_vector_old_),
318 Ptr<const VectorBase<Scalar> >(null),
319 t_, constOptInArg(*solution_vector_),
320 Ptr<const VectorBase<Scalar> >(null),
321 time_vec,ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
322 Ptr<InterpolatorBase<Scalar> >(null)
326 template<
class Scalar>
329 TEUCHOS_ASSERT( time_vec != NULL );
331 if (!haveInitialCondition_) {
334 time_vec->push_back(t_old_);
336 time_vec->push_back(t_);
340 template<
class Scalar>
343 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Error, removeNodes is not implemented for ExplicitRKStepper at this time.\n");
346 template<
class Scalar>
349 return(erkButcherTableau_->order());
352 template <
class Scalar>
355 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
356 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
357 parameterList_ = paramList;
358 Teuchos::readVerboseObjectSublist(&*parameterList_,
this);
361 template <
class Scalar>
364 return(parameterList_);
367 template <
class Scalar>
370 Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
371 parameterList_ = Teuchos::null;
372 return(temp_param_list);
375 template<
class Scalar>
376 RCP<const Teuchos::ParameterList>
379 using Teuchos::ParameterList;
380 static RCP<const ParameterList> validPL;
381 if (is_null(validPL)) {
382 RCP<ParameterList> pl = Teuchos::parameterList();
383 Teuchos::setupVerboseObjectSublist(&*pl);
389 template<
class Scalar>
392 TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
393 TEUCHOS_TEST_FOR_EXCEPT( !is_null(model_) );
394 assertValidModel( *
this, *model );
399 template<
class Scalar>
402 this->setModel(model);
406 template<
class Scalar>
407 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
414 template<
class Scalar>
415 Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >
418 return Teuchos::null;
422 template<
class Scalar>
424 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
430 basePoint_ = initialCondition;
434 RCP<const Thyra::VectorBase<Scalar> >
435 x_init = initialCondition.get_x();
437 #ifdef HAVE_RYTHMOS_DEBUG 438 TEUCHOS_TEST_FOR_EXCEPTION(
439 is_null(x_init), std::logic_error,
440 "Error, if the client passes in an intial condition to setInitialCondition(...),\n" 441 "then x can not be null!" );
444 solution_vector_ = x_init->clone_v();
445 solution_vector_old_ = x_init->clone_v();
449 t_ = initialCondition.get_t();
452 haveInitialCondition_ =
true;
457 template<
class Scalar>
458 Thyra::ModelEvaluatorBase::InArgs<Scalar>
464 template<
class Scalar>
470 template<
class Scalar>
475 RCP<ExplicitRKStepper<Scalar> >
478 if (!is_null(model_)) {
479 stepper->setModel(model_);
482 if (!is_null(erkButcherTableau_)) {
484 stepper->setRKButcherTableau(erkButcherTableau_);
487 if (!is_null(parameterList_)) {
488 stepper->setParameterList(Teuchos::parameterList(*parameterList_));
503 #define RYTHMOS_EXPLICIT_RK_STEPPER_INSTANT(SCALAR) \ 505 template class ExplicitRKStepper< SCALAR >; \ 507 template RCP< ExplicitRKStepper< SCALAR > > \ 508 explicitRKStepper(); \ 510 template RCP< ExplicitRKStepper< SCALAR > > \ 512 const RCP<Thyra::ModelEvaluator< SCALAR > >& model \ 515 template RCP< ExplicitRKStepper< SCALAR > > \ 517 const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \ 518 const RCP<const RKButcherTableauBase< SCALAR > >& rkbt \ 523 #endif //Rythmos_ExplicitRK_STEPPER_DEF_H Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Scalar takeStep(Scalar dt, StepSizeType flag)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
RCP< const Teuchos::ParameterList > getValidParameters() const
void getNodes(Array< Scalar > *time_vec) const
Get interpolation nodes.
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const ¶mList)
Redefined from Teuchos::ParameterListAcceptor.
const StepStatus< Scalar > getStepStatus() const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
int getOrder() const
Get order of interpolation.
void setRKButcherTableau(const RCP< const RKButcherTableauBase< Scalar > > &rkbt)
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
void getPoints(const Array< Scalar > &time_vec, Array< RCP< const VectorBase< Scalar > > > *x_vec, Array< RCP< const VectorBase< Scalar > > > *xdot_vec, Array< ScalarMag > *accuracy_vec) const
Get values from buffer.
void removeNodes(Array< Scalar > &time_vec)
Remove interpolation nodes.
Thyra::ModelEvaluatorBase::InArgs< Scalar > getInitialCondition() const
RCP< StepperBase< Scalar > > cloneStepperAlgorithm() const
RCP< const RKButcherTableauBase< Scalar > > getRKButcherTableau() const
RCP< const Thyra::VectorBase< Scalar > > solution
RCP< Thyra::ModelEvaluator< Scalar > > getNonconstModel()
bool supportsCloning() const
void setNonconstModel(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model)
TimeRange< Scalar > getTimeRange() const
void addPoints(const Array< Scalar > &time_vec, const Array< Teuchos::RCP< const Thyra::VectorBase< Scalar > > > &x_vec, const Array< Teuchos::RCP< const Thyra::VectorBase< Scalar > > > &xdot_vec)
void setInitialCondition(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &initialCondition)