Rythmos - Transient Integration for Differential Equations  Version of the Day
Rythmos_ExplicitRKStepper_def.hpp
1 //@HEADER
2 // ***********************************************************************
3 //
4 // Rythmos Package
5 // Copyright (2006) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25 //
26 // ***********************************************************************
27 //@HEADER
28 
29 #ifndef Rythmos_ExplicitRK_STEPPER_DEF_H
30 #define Rythmos_ExplicitRK_STEPPER_DEF_H
31 
32 #include "Rythmos_ExplicitRKStepper_decl.hpp"
33 
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"
40 
41 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
42 
43 #include "Thyra_ModelEvaluatorHelpers.hpp"
44 #include "Thyra_MultiVectorStdOps.hpp"
45 #include "Thyra_VectorStdOps.hpp"
46 
47 
48 namespace Rythmos {
49 
50 // Non-member constructors
51 template<class Scalar>
52 RCP<ExplicitRKStepper<Scalar> > explicitRKStepper()
53 {
54  RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
55  return stepper;
56 }
57 
58 template<class Scalar>
59 RCP<ExplicitRKStepper<Scalar> > explicitRKStepper(
60  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model
61  )
62 {
63  RCP<RKButcherTableauBase<Scalar> > rkbt = createRKBT<Scalar>("Explicit 4 Stage");
64  //RCP<RKButcherTableauBase<Scalar> > rkbt = rcp(new Explicit4Stage4thOrder_RKBT<Scalar>());
65  RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
66  stepper->setModel(model);
67  stepper->setRKButcherTableau(rkbt);
68  return stepper;
69 }
70 
71 template<class Scalar>
72 RCP<ExplicitRKStepper<Scalar> > explicitRKStepper(
73  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
74  const RCP<const RKButcherTableauBase<Scalar> >& rkbt
75  )
76 {
77  RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
78  stepper->setModel(model);
79  stepper->setRKButcherTableau(rkbt);
80  return stepper;
81 }
82 
83 template<class Scalar>
85 {
86  this->defaultInitializeAll_();
87  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
88  out->precision(15);
89  erkButcherTableau_ = rKButcherTableau<Scalar>();
90  numSteps_ = 0;
91 }
92 
93 template<class Scalar>
95 {
96  model_ = Teuchos::null;
97  solution_vector_ = Teuchos::null;
98  solution_vector_old_ = Teuchos::null;
99  //k_vector_;
100  ktemp_vector_ = Teuchos::null;
101  //basePoint_;
102  erkButcherTableau_ = Teuchos::null;
103  t_ = ST::nan();
104  t_old_ = ST::nan();
105  dt_ = ST::nan();
106  numSteps_ = -1;
107  parameterList_ = Teuchos::null;
108  isInitialized_ = false;
109  haveInitialCondition_ = false;
110 }
111 
112 template<class Scalar>
113 void ExplicitRKStepper<Scalar>::setRKButcherTableau(const RCP<const RKButcherTableauBase<Scalar> > &rkbt)
114 {
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!"
121  );
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()));
128  }
129  }
130  }
131  erkButcherTableau_ = rkbt;
132 }
133 
134 template<class Scalar>
135 RCP<const RKButcherTableauBase<Scalar> > ExplicitRKStepper<Scalar>::getRKButcherTableau() const
136 {
137  return erkButcherTableau_;
138 }
139 
140 template<class Scalar>
142 {
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!"
149  );
150  ktemp_vector_ = Thyra::createMember(model_->get_f_space());
151  // Initialize the stage vectors
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()));
156  }
157  }
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;
164 }
165 
166 
167 template<class Scalar>
169 {
170 }
171 
172 template<class Scalar>
173 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > ExplicitRKStepper<Scalar>::get_x_space() const
174 {
175  TEUCHOS_ASSERT( !is_null(model_) );
176  return(model_->get_x_space());
177 }
178 
179 template<class Scalar>
180 Scalar ExplicitRKStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
181 {
182  typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag;
183  this->initialize_();
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."
187  );
188 #endif // HAVE_RYTHMOS_DEBUG
189  if ((flag == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
190  return(Scalar(-ST::one()));
191  }
192  // Store old solution & old time
193  V_V(solution_vector_old_.ptr(), *solution_vector_);
194  t_old_ = t_;
195 
196  dt_ = dt;
197 
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();
202  // Compute stage solutions
203  for (int s=0 ; s < stages ; ++s) {
204  Thyra::assign(ktemp_vector_.ptr(), *solution_vector_); // ktemp = solution_vector
205  for (int j=0 ; j < s ; ++j) { // assuming Butcher matix is strictly lower triangular
206  if (A(s,j) != ST::zero()) {
207  Thyra::Vp_StV(ktemp_vector_.ptr(), A(s,j), *k_vector_[j]); // ktemp = ktemp + a_{s+1,j+1}*k_{j+1}
208  }
209  }
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); // k_s = k_s*dt
213  }
214  // Sum for solution:
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]); // solution_vector += b_{s+1}*k_{s+1}
218  }
219  }
220 
221  // update current time:
222  t_ = t_ + dt;
223 
224  numSteps_++;
225 
226  return(dt);
227 }
228 
229 template<class Scalar>
231 {
232  StepStatus<Scalar> stepStatus;
233 
234  if (!haveInitialCondition_) {
235  stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
236  } else if (numSteps_ == 0) {
237  stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
238  stepStatus.order = erkButcherTableau_->order();
239  stepStatus.time = t_;
240  stepStatus.solution = solution_vector_;
241  } else {
242  stepStatus.stepStatus = STEP_STATUS_CONVERGED;
243  stepStatus.stepSize = dt_;
244  stepStatus.order = erkButcherTableau_->order();
245  stepStatus.time = t_;
246  stepStatus.solution = solution_vector_;
247  }
248 
249  return(stepStatus);
250 }
251 
252 template<class Scalar>
254  Teuchos::FancyOStream &out,
255  const Teuchos::EVerbosityLevel verbLevel
256  ) const
257 {
258  if ( (static_cast<int>(verbLevel) == static_cast<int>(Teuchos::VERB_DEFAULT) ) ||
259  (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW) )
260  ) {
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);
275  }
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;
282  }
283 }
284 
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
290  )
291 {
292  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for ExplicitRKStepper at this time.\n");
293 }
294 
295 template<class Scalar>
297 {
298  if (!haveInitialCondition_) {
299  return(invalidTimeRange<Scalar>());
300  } else {
301  return(TimeRange<Scalar>(t_old_,t_));
302  }
303 }
304 
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
311  ) const
312 {
313  TEUCHOS_ASSERT( haveInitialCondition_ );
314  using Teuchos::constOptInArg;
315  using Teuchos::null;
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)
323  );
324 }
325 
326 template<class Scalar>
327 void ExplicitRKStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
328 {
329  TEUCHOS_ASSERT( time_vec != NULL );
330  time_vec->clear();
331  if (!haveInitialCondition_) {
332  return;
333  }
334  time_vec->push_back(t_old_);
335  if (t_ != t_old_) {
336  time_vec->push_back(t_);
337  }
338 }
339 
340 template<class Scalar>
341 void ExplicitRKStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
342 {
343  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ExplicitRKStepper at this time.\n");
344 }
345 
346 template<class Scalar>
348 {
349  return(erkButcherTableau_->order());
350 }
351 
352 template <class Scalar>
353 void ExplicitRKStepper<Scalar>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
354 {
355  TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
356  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
357  parameterList_ = paramList;
358  Teuchos::readVerboseObjectSublist(&*parameterList_,this);
359 }
360 
361 template <class Scalar>
362 Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::getNonconstParameterList()
363 {
364  return(parameterList_);
365 }
366 
367 template <class Scalar>
368 Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::unsetParameterList()
369 {
370  Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
371  parameterList_ = Teuchos::null;
372  return(temp_param_list);
373 }
374 
375 template<class Scalar>
376 RCP<const Teuchos::ParameterList>
378 {
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);
384  validPL = pl;
385  }
386  return validPL;
387 }
388 
389 template<class Scalar>
390 void ExplicitRKStepper<Scalar>::setModel(const RCP<const Thyra::ModelEvaluator<Scalar> >& model)
391 {
392  TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
393  TEUCHOS_TEST_FOR_EXCEPT( !is_null(model_) ); // For now you can only call this once.
394  assertValidModel( *this, *model );
395  model_ = model;
396 }
397 
398 
399 template<class Scalar>
400 void ExplicitRKStepper<Scalar>::setNonconstModel(const RCP<Thyra::ModelEvaluator<Scalar> >& model)
401 {
402  this->setModel(model); // TODO 09/09/09 tscoffe: use ConstNonconstObjectContainer!
403 }
404 
405 
406 template<class Scalar>
407 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
409 {
410  return model_;
411 }
412 
413 
414 template<class Scalar>
415 Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >
417 {
418  return Teuchos::null;
419 }
420 
421 
422 template<class Scalar>
424  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
425  )
426 {
427 
428  // typedef Thyra::ModelEvaluatorBase MEB; // unused
429 
430  basePoint_ = initialCondition;
431 
432  // x
433 
434  RCP<const Thyra::VectorBase<Scalar> >
435  x_init = initialCondition.get_x();
436 
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!" );
442 #endif
443 
444  solution_vector_ = x_init->clone_v();
445  solution_vector_old_ = x_init->clone_v();
446 
447  // t
448 
449  t_ = initialCondition.get_t();
450  t_old_ = t_;
451 
452  haveInitialCondition_ = true;
453 
454 }
455 
456 
457 template<class Scalar>
458 Thyra::ModelEvaluatorBase::InArgs<Scalar>
460 {
461  return basePoint_;
462 }
463 
464 template<class Scalar>
466 {
467  return true;
468 }
469 
470 template<class Scalar>
471 RCP<StepperBase<Scalar> > ExplicitRKStepper<Scalar>::cloneStepperAlgorithm() const
472 {
473  // Just use the interface to clone the algorithm in a basically
474  // uninitialized state
475  RCP<ExplicitRKStepper<Scalar> >
476  stepper = Teuchos::rcp(new ExplicitRKStepper<Scalar>());
477 
478  if (!is_null(model_)) {
479  stepper->setModel(model_); // Shallow copy is okay!
480  }
481 
482  if (!is_null(erkButcherTableau_)) {
483  // 06/16/09 tscoffe: should we clone the RKBT here?
484  stepper->setRKButcherTableau(erkButcherTableau_);
485  }
486 
487  if (!is_null(parameterList_)) {
488  stepper->setParameterList(Teuchos::parameterList(*parameterList_));
489  }
490 
491  return stepper;
492 
493 }
494 
495 
496 
497 //
498 // Explicit Instantiation macro
499 //
500 // Must be expanded from within the Rythmos namespace!
501 //
502 
503 #define RYTHMOS_EXPLICIT_RK_STEPPER_INSTANT(SCALAR) \
504  \
505  template class ExplicitRKStepper< SCALAR >; \
506  \
507  template RCP< ExplicitRKStepper< SCALAR > > \
508  explicitRKStepper(); \
509  \
510  template RCP< ExplicitRKStepper< SCALAR > > \
511  explicitRKStepper( \
512  const RCP<Thyra::ModelEvaluator< SCALAR > >& model \
513  ); \
514  \
515  template RCP< ExplicitRKStepper< SCALAR > > \
516  explicitRKStepper( \
517  const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
518  const RCP<const RKButcherTableauBase< SCALAR > >& rkbt \
519  ); \
520 
521 } // namespace Rythmos
522 
523 #endif //Rythmos_ExplicitRK_STEPPER_DEF_H
524 
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 &paramList)
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()
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)