Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
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 stepSizeType)
181 {
182 
183  RCP<FancyOStream> out = this->getOStream();
184  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
185  Teuchos::OSTab ostab(out,1,"takeStep");
186  Scalar stepSizeTaken;
187 
188 
189  // not needed for this
190  int desiredOrder;
191 
192  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
193  *out
194  << "\nEntering "
195  << Teuchos::TypeNameTraits<ExplicitRKStepper<Scalar> >::name()
196  << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n";
197  }
198 
199 
200  if (stepSizeType == STEP_TYPE_FIXED) {
201  stepSizeTaken = takeFixedStep_(dt , stepSizeType);
202  return stepSizeTaken;
203  } else {
204  isVariableStep_ = true;
205 
206  stepControl_->setOStream(out);
207  stepControl_->setVerbLevel(verbLevel);
208 
209  rkNewtonConvergenceStatus_ = -1;
210 
211  while (rkNewtonConvergenceStatus_ < 0){
212 
213  stepControl_->setRequestedStepSize(*this, dt, stepSizeType);
214  stepControl_->nextStepSize(*this, &dt, &stepSizeType, &desiredOrder);
215 
216  stepSizeTaken = takeVariableStep_(dt, stepSizeType);
217 
218  }
219  return stepSizeTaken;
220  }
221 
222 
223 }
224 
225 
226 template<class Scalar>
227 Scalar ExplicitRKStepper<Scalar>::takeVariableStep_(Scalar dt, StepSizeType /* stepSizeType */)
228 {
229  typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag;
230  this->initialize_();
231 
232  // Store old solution & old time
233  V_V(solution_vector_old_.ptr(), *solution_vector_);
234  t_old_ = t_;
235  Scalar current_dt = dt;
236  Scalar t = timeRange_.upper();
237  Scalar dt_to_return;
238  AttemptedStepStatusFlag status;
239  Scalar dt_old = dt;
240 
241  dt_ = dt;
242 
243  int stages = erkButcherTableau_->numStages();
244  Teuchos::SerialDenseMatrix<int,Scalar> A = erkButcherTableau_->A();
245  Teuchos::SerialDenseVector<int,Scalar> b = erkButcherTableau_->b();
246  Teuchos::SerialDenseVector<int,Scalar> c = erkButcherTableau_->c();
247 
248 
249  // Compute stage solutions
250  for (int s=0 ; s < stages ; ++s) {
251  Thyra::assign(ktemp_vector_.ptr(), *solution_vector_); // ktemp = solution_vector
252  for (int j=0 ; j < s ; ++j) { // assuming Butcher matix is strictly lower triangular
253  if (A(s,j) != ST::zero()) {
254  Thyra::Vp_StV(ktemp_vector_.ptr(), A(s,j), *k_vector_[j]); // ktemp = ktemp + a_{s+1,j+1}*k_{j+1}
255  }
256  }
257  TScalarMag ts = t_ + c(s)*dt;
258  TScalarMag scaled_dt = (s == 0)? Scalar(dt/stages) : c(s)*dt;
259 
260  // need to check here the status of the solver (the linear solve)
261  //eval_model_explicit<Scalar>(*model_,basePoint_,*ktemp_vector_,ts,Teuchos::outArg(*k_vector_[s]));
262  eval_model_explicit<Scalar>(*model_,basePoint_,*ktemp_vector_,ts,Teuchos::outArg(*k_vector_[s]), scaled_dt, c(s));
263  Thyra::Vt_S(k_vector_[s].ptr(),dt); // k_s = k_s*dt
264  }
265  // Sum for solution:
266  for (int s=0 ; s < stages ; ++s) {
267  if (b(s) != ST::zero()) {
268  Thyra::Vp_StV(solution_vector_.ptr(), b(s), *k_vector_[s]); // solution_vector += b_{s+1}*k_{s+1}
269  }
270  }
271 
272  // cheat and say that the solver converged ( although no solver is needed for explicit method )
273  rkNewtonConvergenceStatus_ = 0;
274  stepControl_->setCorrection(*this, solution_vector_, Teuchos::null , rkNewtonConvergenceStatus_);
275  bool stepPass = stepControl_->acceptStep(*this, &LETvalue_);
276 
277  if (erkButcherTableau_->isEmbeddedMethod() ){
278 
279  Teuchos::SerialDenseVector<int,Scalar> bhat = erkButcherTableau_->bhat();
280 
281  // Sum for Embedded solution:
282  for (int s=0 ; s < stages ; ++s) {
283  if (bhat(s) != ST::zero()) {
284 
285  // solution_hat_vector += bhat_{s+1}*k_{s+1}
286  Thyra::Vp_StV(solution_hat_vector_.ptr(), bhat(s), *k_vector_[s]);
287  }
288  }
289 
290  // ee_ = (solution_vector__ - solution_hat_vector_)
291  Thyra::V_VmV(ee_.ptr(), *solution_vector_, *solution_hat_vector_);
292 
293  stepControl_->setCorrection(*this, solution_vector_, ee_ , rkNewtonConvergenceStatus_);
294  }
295  else {
296  stepControl_->setCorrection(*this, solution_vector_, Teuchos::null, rkNewtonConvergenceStatus_);
297  }
298 
299  stepPass = stepControl_->acceptStep(*this, &LETvalue_);
300 
301  if (!stepPass) { // stepPass = false
302  stepLETStatus_ = STEP_LET_STATUS_FAILED;
303  rkNewtonConvergenceStatus_ = -1; // just making sure here
304  } else { // stepPass = true
305  stepLETStatus_ = STEP_LET_STATUS_PASSED;
306  rkNewtonConvergenceStatus_ = 0; // just making sure here
307  }
308 
309  if (rkNewtonConvergenceStatus_ == 0) {
310 
311  // Update time range
312  timeRange_ = timeRange(t,t+current_dt);
313  numSteps_++;
314 
315  // completeStep only if the none of the stage solution's failed to converged
316  stepControl_->completeStep(*this);
317 
318  dt_to_return = current_dt;
319 
320  } else {
321 
322  rkNewtonConvergenceStatus_ = -1;
323  status = stepControl_-> rejectStep(*this); // reject the stage value
324  (void) status; // avoid "set but not used" build warning
325  dt_to_return = dt_old;
326  }
327 
328  // update current time:
329  t_ = t_ + dt;
330  numSteps_++;
331  return( dt_to_return );
332 }
333 
334 
335 template<class Scalar>
336 Scalar ExplicitRKStepper<Scalar>::takeFixedStep_(Scalar dt, StepSizeType flag)
337 {
338  typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag;
339  this->initialize_();
340 #ifdef HAVE_RYTHMOS_DEBUG
341  TEUCHOS_TEST_FOR_EXCEPTION( flag == STEP_TYPE_VARIABLE, std::logic_error,
342  "Error! ExplicitRKStepper does not support variable time steps at this time."
343  );
344 #endif // HAVE_RYTHMOS_DEBUG
345  if ((flag == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
346  return(Scalar(-ST::one()));
347  }
348  // Store old solution & old time
349  V_V(solution_vector_old_.ptr(), *solution_vector_);
350  t_old_ = t_;
351 
352  dt_ = dt;
353 
354  int stages = erkButcherTableau_->numStages();
355  Teuchos::SerialDenseMatrix<int,Scalar> A = erkButcherTableau_->A();
356  Teuchos::SerialDenseVector<int,Scalar> b = erkButcherTableau_->b();
357  Teuchos::SerialDenseVector<int,Scalar> c = erkButcherTableau_->c();
358  // Compute stage solutions
359  for (int s=0 ; s < stages ; ++s) {
360  Thyra::assign(ktemp_vector_.ptr(), *solution_vector_); // ktemp = solution_vector
361  for (int j=0 ; j < s ; ++j) { // assuming Butcher matix is strictly lower triangular
362  if (A(s,j) != ST::zero()) {
363  Thyra::Vp_StV(ktemp_vector_.ptr(), A(s,j), *k_vector_[j]); // ktemp = ktemp + a_{s+1,j+1}*k_{j+1}
364  }
365  }
366  TScalarMag ts = t_ + c(s)*dt;
367  TScalarMag scaled_dt = (s == 0)? Scalar(dt/stages) : c(s)*dt;
368 
369  eval_model_explicit<Scalar>(*model_,basePoint_,*ktemp_vector_,ts,Teuchos::outArg(*k_vector_[s]), scaled_dt, c(s));
370  Thyra::Vt_S(k_vector_[s].ptr(),dt); // k_s = k_s*dt
371  }
372  // Sum for solution:
373  for (int s=0 ; s < stages ; ++s) {
374  if (b(s) != ST::zero()) {
375  Thyra::Vp_StV(solution_vector_.ptr(), b(s), *k_vector_[s]); // solution_vector += b_{s+1}*k_{s+1}
376  }
377  }
378 
379  // update current time:
380  t_ = t_ + dt;
381 
382  numSteps_++;
383 
384  return(dt);
385 }
386 
387 template<class Scalar>
389 {
390  StepStatus<Scalar> stepStatus;
391 
392  if (!haveInitialCondition_) {
393  stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
394  } else if (numSteps_ == 0) {
395  stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
396  stepStatus.order = erkButcherTableau_->order();
397  stepStatus.time = t_;
398  stepStatus.solution = solution_vector_;
399  } else {
400  stepStatus.stepStatus = STEP_STATUS_CONVERGED;
401  stepStatus.stepSize = dt_;
402  stepStatus.order = erkButcherTableau_->order();
403  stepStatus.time = t_;
404  stepStatus.solution = solution_vector_;
405  }
406 
407  return(stepStatus);
408 }
409 
410 template<class Scalar>
412  Teuchos::FancyOStream &out,
413  const Teuchos::EVerbosityLevel verbLevel
414  ) const
415 {
416  if ( (static_cast<int>(verbLevel) == static_cast<int>(Teuchos::VERB_DEFAULT) ) ||
417  (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW) )
418  ) {
419  out << this->description() << "::describe" << std::endl;
420  out << "model = " << model_->description() << std::endl;
421  out << erkButcherTableau_->numStages() << " stage Explicit RK method" << std::endl;
422  } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)) {
423  out << "solution_vector = " << std::endl;
424  out << Teuchos::describe(*solution_vector_,verbLevel);
425  } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM)) {
426  } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_HIGH)) {
427  out << "model = " << std::endl;
428  out << Teuchos::describe(*model_,verbLevel);
429  int stages = erkButcherTableau_->numStages();
430  for (int i=0 ; i<stages ; ++i) {
431  out << "k_vector[" << i << "] = " << std::endl;
432  out << Teuchos::describe(*k_vector_[i],verbLevel);
433  }
434  out << "ktemp_vector = " << std::endl;
435  out << Teuchos::describe(*ktemp_vector_,verbLevel);
436  out << "ERK Butcher Tableau A matrix: " << erkButcherTableau_->A() << std::endl;
437  out << "ERK Butcher Tableau b vector: " << erkButcherTableau_->b() << std::endl;
438  out << "ERK Butcher Tableau c vector: " << erkButcherTableau_->c() << std::endl;
439  out << "t = " << t_ << std::endl;
440  }
441 }
442 
443 template<class Scalar>
445  const Array<Scalar>& /* time_vec */
446  ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& /* x_vec */
447  ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& /* xdot_vec */
448  )
449 {
450  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for ExplicitRKStepper at this time.\n");
451 }
452 
453 template<class Scalar>
455 {
456  if (!haveInitialCondition_) {
457  return(invalidTimeRange<Scalar>());
458  } else {
459  return(TimeRange<Scalar>(t_old_,t_));
460  }
461 }
462 
463 template<class Scalar>
465  const Array<Scalar>& time_vec,
466  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
467  Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
468  Array<ScalarMag>* accuracy_vec
469  ) const
470 {
471  TEUCHOS_ASSERT( haveInitialCondition_ );
472  using Teuchos::constOptInArg;
473  using Teuchos::null;
474  defaultGetPoints<Scalar>(
475  t_old_, constOptInArg(*solution_vector_old_),
476  Ptr<const VectorBase<Scalar> >(null),
477  t_, constOptInArg(*solution_vector_),
478  Ptr<const VectorBase<Scalar> >(null),
479  time_vec,ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
480  Ptr<InterpolatorBase<Scalar> >(null)
481  );
482 }
483 
484 template<class Scalar>
485 void ExplicitRKStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
486 {
487  TEUCHOS_ASSERT( time_vec != NULL );
488  time_vec->clear();
489  if (!haveInitialCondition_) {
490  return;
491  }
492  time_vec->push_back(t_old_);
493  if (t_ != t_old_) {
494  time_vec->push_back(t_);
495  }
496 }
497 
498 template<class Scalar>
499 void ExplicitRKStepper<Scalar>::removeNodes(Array<Scalar>& /* time_vec */)
500 {
501  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ExplicitRKStepper at this time.\n");
502 }
503 
504 template<class Scalar>
506 {
507  return(erkButcherTableau_->order());
508 }
509 
510 template <class Scalar>
511 void ExplicitRKStepper<Scalar>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
512 {
513  TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
514  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
515  parameterList_ = paramList;
516  Teuchos::readVerboseObjectSublist(&*parameterList_,this);
517 }
518 
519 template <class Scalar>
520 Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::getNonconstParameterList()
521 {
522  return(parameterList_);
523 }
524 
525 template <class Scalar>
526 Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::unsetParameterList()
527 {
528  Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
529  parameterList_ = Teuchos::null;
530  return(temp_param_list);
531 }
532 
533 template<class Scalar>
534 RCP<const Teuchos::ParameterList>
536 {
537  using Teuchos::ParameterList;
538  static RCP<const ParameterList> validPL;
539  if (is_null(validPL)) {
540  RCP<ParameterList> pl = Teuchos::parameterList();
541  Teuchos::setupVerboseObjectSublist(&*pl);
542  validPL = pl;
543  }
544  return validPL;
545 }
546 
547 template<class Scalar>
548 void ExplicitRKStepper<Scalar>::setModel(const RCP<const Thyra::ModelEvaluator<Scalar> >& model)
549 {
550  TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
551  TEUCHOS_TEST_FOR_EXCEPT( !is_null(model_) ); // For now you can only call this once.
552  assertValidModel( *this, *model );
553  model_ = model;
554 }
555 
556 
557 template<class Scalar>
558 void ExplicitRKStepper<Scalar>::setNonconstModel(const RCP<Thyra::ModelEvaluator<Scalar> >& model)
559 {
560  this->setModel(model); // TODO 09/09/09 tscoffe: use ConstNonconstObjectContainer!
561 }
562 
563 
564 template<class Scalar>
565 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
567 {
568  return model_;
569 }
570 
571 
572 template<class Scalar>
573 Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >
575 {
576  return Teuchos::null;
577 }
578 
579 
580 template<class Scalar>
582  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
583  )
584 {
585 
586  // typedef Thyra::ModelEvaluatorBase MEB; // unused
587 
588  basePoint_ = initialCondition;
589 
590  // x
591 
592  RCP<const Thyra::VectorBase<Scalar> >
593  x_init = initialCondition.get_x();
594 
595 #ifdef HAVE_RYTHMOS_DEBUG
596  TEUCHOS_TEST_FOR_EXCEPTION(
597  is_null(x_init), std::logic_error,
598  "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
599  "then x can not be null!" );
600 #endif
601 
602  solution_vector_ = x_init->clone_v();
603  solution_vector_old_ = x_init->clone_v();
604 
605  // for the embedded RK method
606  solution_hat_vector_ = x_init->clone_v();
607  ee_ = x_init->clone_v();
608 
609  // t
610 
611  t_ = initialCondition.get_t();
612  t_old_ = t_;
613 
614  haveInitialCondition_ = true;
615 
616 }
617 
618 
619 template<class Scalar>
620 Thyra::ModelEvaluatorBase::InArgs<Scalar>
622 {
623  return basePoint_;
624 }
625 
626 template<class Scalar>
628 {
629  return true;
630 }
631 
632 template<class Scalar>
633 RCP<StepperBase<Scalar> > ExplicitRKStepper<Scalar>::cloneStepperAlgorithm() const
634 {
635  // Just use the interface to clone the algorithm in a basically
636  // uninitialized state
637  RCP<ExplicitRKStepper<Scalar> >
638  stepper = Teuchos::rcp(new ExplicitRKStepper<Scalar>());
639 
640  if (!is_null(model_)) {
641  stepper->setModel(model_); // Shallow copy is okay!
642  }
643 
644  if (!is_null(erkButcherTableau_)) {
645  // 06/16/09 tscoffe: should we clone the RKBT here?
646  stepper->setRKButcherTableau(erkButcherTableau_);
647  }
648 
649  if (!is_null(parameterList_)) {
650  stepper->setParameterList(Teuchos::parameterList(*parameterList_));
651  }
652 
653  return stepper;
654 
655 }
656 
657 template<class Scalar>
658 RCP<StepControlStrategyBase<Scalar> > ExplicitRKStepper<Scalar>::getNonconstStepControlStrategy()
659 {
660  return(stepControl_);
661 }
662 
663 template<class Scalar>
664 RCP<const StepControlStrategyBase<Scalar> > ExplicitRKStepper<Scalar>::getStepControlStrategy() const
665 {
666  return(stepControl_);
667 }
668 
669 template<class Scalar>
671 {
672  TEUCHOS_TEST_FOR_EXCEPTION(stepControl == Teuchos::null,std::logic_error,"Error, stepControl == Teuchos::null!\n");
673  stepControl_ = stepControl;
674 }
675 //
676 // Explicit Instantiation macro
677 //
678 // Must be expanded from within the Rythmos namespace!
679 //
680 
681 #define RYTHMOS_EXPLICIT_RK_STEPPER_INSTANT(SCALAR) \
682  \
683  template class ExplicitRKStepper< SCALAR >; \
684  \
685  template RCP< ExplicitRKStepper< SCALAR > > \
686  explicitRKStepper(); \
687  \
688  template RCP< ExplicitRKStepper< SCALAR > > \
689  explicitRKStepper( \
690  const RCP<Thyra::ModelEvaluator< SCALAR > >& model \
691  ); \
692  \
693  template RCP< ExplicitRKStepper< SCALAR > > \
694  explicitRKStepper( \
695  const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
696  const RCP<const RKButcherTableauBase< SCALAR > >& rkbt \
697  ); \
698 
699 } // namespace Rythmos
700 
701 #endif //Rythmos_ExplicitRK_STEPPER_DEF_H
702 
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Scalar takeStep(Scalar dt, StepSizeType flag)
void setStepControlStrategy(const RCP< StepControlStrategyBase< Scalar > > &stepControlStrategy)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
RCP< const Teuchos::ParameterList > getValidParameters() const
The member functions in the StepControlStrategyBase move you between these states in the following fa...
void getNodes(Array< Scalar > *time_vec) const
Get interpolation nodes.
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Redefined from Teuchos::ParameterListAcceptor.
RCP< const StepControlStrategyBase< Scalar > > getStepControlStrategy() const
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< StepControlStrategyBase< Scalar > > getNonconstStepControlStrategy()
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)