Rythmos_ExplicitRKStepper_def.hpp

//@HEADER
// ***********************************************************************
//
//                           Rythmos Package
//                 Copyright (2006) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
// USA
// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
//
// ***********************************************************************
//@HEADER

#ifndef Rythmos_ExplicitRK_STEPPER_DEF_H
#define Rythmos_ExplicitRK_STEPPER_DEF_H

#include "Rythmos_ExplicitRKStepper_decl.hpp"

#include "Rythmos_RKButcherTableau.hpp"
#include "Rythmos_RKButcherTableauHelpers.hpp"
#include "Rythmos_RKButcherTableauBuilder.hpp"
#include "Rythmos_StepperHelpers.hpp"
#include "Rythmos_LinearInterpolator.hpp"
#include "Rythmos_InterpolatorBaseHelpers.hpp"

#include "Teuchos_VerboseObjectParameterListHelpers.hpp"

#include "Thyra_ModelEvaluatorHelpers.hpp"
#include "Thyra_MultiVectorStdOps.hpp"
#include "Thyra_VectorStdOps.hpp"


namespace Rythmos {

// Non-member constructors
template<class Scalar>
RCP<ExplicitRKStepper<Scalar> > explicitRKStepper()
{
  RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
  return stepper;
}

template<class Scalar>
RCP<ExplicitRKStepper<Scalar> > explicitRKStepper(
    const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model
    )
{
  RCP<RKButcherTableauBase<Scalar> > rkbt = createRKBT<Scalar>("Explicit 4 Stage");
  //RCP<RKButcherTableauBase<Scalar> > rkbt = rcp(new Explicit4Stage4thOrder_RKBT<Scalar>());
  RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
  stepper->setModel(model);
  stepper->setRKButcherTableau(rkbt);
  return stepper;
}

template<class Scalar>
RCP<ExplicitRKStepper<Scalar> > explicitRKStepper(
    const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
    const RCP<const RKButcherTableauBase<Scalar> >& rkbt
    )
{
  RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
  stepper->setModel(model);
  stepper->setRKButcherTableau(rkbt);
  return stepper;
}

template<class Scalar>
ExplicitRKStepper<Scalar>::ExplicitRKStepper()
{
  this->defaultInitializeAll_();
  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
  out->precision(15);
  erkButcherTableau_ = rKButcherTableau<Scalar>();
  numSteps_ = 0;
}

template<class Scalar>
void ExplicitRKStepper<Scalar>::defaultInitializeAll_()
{
  model_ = Teuchos::null;
  solution_vector_ = Teuchos::null;
  solution_vector_old_ = Teuchos::null;
  //k_vector_;
  ktemp_vector_ = Teuchos::null;
  //basePoint_;
  erkButcherTableau_ = Teuchos::null;
  t_ = ST::nan();
  t_old_ = ST::nan();
  dt_ = ST::nan();
  numSteps_ = -1;
  parameterList_ = Teuchos::null;
  isInitialized_ = false;
  haveInitialCondition_ = false;
}

template<class Scalar>
void ExplicitRKStepper<Scalar>::setRKButcherTableau(const RCP<const RKButcherTableauBase<Scalar> > &rkbt)
{
  TEUCHOS_ASSERT( !is_null(rkbt) );
  validateERKButcherTableau(*rkbt);
  int numStages_old = erkButcherTableau_->numStages();
  int numStages_new = rkbt->numStages();
  TEUCHOS_TEST_FOR_EXCEPTION( numStages_new == 0, std::logic_error,
      "Error!  The Runge-Kutta Butcher tableau has no stages!"
      );
  if (!is_null(model_)) {
    int numNewStages = numStages_new - numStages_old;
    if ( numNewStages > 0 ) {
      k_vector_.reserve(numStages_new);
      for (int i=0 ; i<numNewStages ; ++i) {
        k_vector_.push_back(Thyra::createMember(model_->get_f_space()));
      }
    }
  }
  erkButcherTableau_ = rkbt;
}

template<class Scalar>
RCP<const RKButcherTableauBase<Scalar> > ExplicitRKStepper<Scalar>::getRKButcherTableau() const
{
  return erkButcherTableau_;
}

template<class Scalar>
void ExplicitRKStepper<Scalar>::initialize_()
{
  if (!isInitialized_) {
    TEUCHOS_ASSERT( !is_null(model_) );
    TEUCHOS_ASSERT( !is_null(erkButcherTableau_) );
    TEUCHOS_ASSERT( haveInitialCondition_ );
    TEUCHOS_TEST_FOR_EXCEPTION( erkButcherTableau_->numStages() == 0, std::logic_error,
        "Error!  The Runge-Kutta Butcher tableau has no stages!"
        );
    ktemp_vector_ = Thyra::createMember(model_->get_f_space());
    // Initialize the stage vectors
    int numStages = erkButcherTableau_->numStages();
    k_vector_.reserve(numStages);
    for (int i=0 ; i<numStages ; ++i) {
      k_vector_.push_back(Thyra::createMember(model_->get_f_space()));
    }
  }
#ifdef HAVE_RYTHMOS_DEBUG
  THYRA_ASSERT_VEC_SPACES(
    "Rythmos::ExplicitRKStepper::initialize_(...)",
    *solution_vector_->space(), *model_->get_x_space() );
#endif // HAVE_RYTHMOS_DEBUG
  isInitialized_ = true;
}


template<class Scalar>
ExplicitRKStepper<Scalar>::~ExplicitRKStepper()
{
}

template<class Scalar>
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > ExplicitRKStepper<Scalar>::get_x_space() const
{
  TEUCHOS_ASSERT( !is_null(model_) );
  return(model_->get_x_space());
}

template<class Scalar>
Scalar ExplicitRKStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
{

  RCP<FancyOStream> out = this->getOStream();
  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
  Teuchos::OSTab ostab(out,1,"takeStep");
  Scalar stepSizeTaken;


  // not needed for this
  int desiredOrder;

  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
    *out
      << "\nEntering "
      << Teuchos::TypeNameTraits<ExplicitRKStepper<Scalar> >::name()
      << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n";
  }


  if (stepSizeType == STEP_TYPE_FIXED) {
    stepSizeTaken = takeFixedStep_(dt , stepSizeType);
    return stepSizeTaken;
  }  else {
      isVariableStep_ = true;

      stepControl_->setOStream(out);
      stepControl_->setVerbLevel(verbLevel);

      rkNewtonConvergenceStatus_ = -1;

      while (rkNewtonConvergenceStatus_ < 0){

         stepControl_->setRequestedStepSize(*this, dt, stepSizeType);
         stepControl_->nextStepSize(*this, &dt, &stepSizeType, &desiredOrder);

         stepSizeTaken = takeVariableStep_(dt, stepSizeType);

      }
     return stepSizeTaken;
  }


}


template<class Scalar>
Scalar ExplicitRKStepper<Scalar>::takeVariableStep_(Scalar dt, StepSizeType /* stepSizeType */)
{
  typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag;
  this->initialize_();

  // Store old solution & old time
  V_V(solution_vector_old_.ptr(), *solution_vector_);
  t_old_ = t_;
  Scalar current_dt = dt;
  Scalar t = timeRange_.upper();
  Scalar dt_to_return;
  AttemptedStepStatusFlag status;
  Scalar dt_old = dt;

  dt_ = dt;

  int stages = erkButcherTableau_->numStages();
  Teuchos::SerialDenseMatrix<int,Scalar> A = erkButcherTableau_->A();
  Teuchos::SerialDenseVector<int,Scalar> b = erkButcherTableau_->b();
  Teuchos::SerialDenseVector<int,Scalar> c = erkButcherTableau_->c();

 
  // Compute stage solutions
  for (int s=0 ; s < stages ; ++s) {
    Thyra::assign(ktemp_vector_.ptr(), *solution_vector_); // ktemp = solution_vector
    for (int j=0 ; j < s ; ++j) { // assuming Butcher matix is strictly lower triangular
      if (A(s,j) != ST::zero()) {
        Thyra::Vp_StV(ktemp_vector_.ptr(), A(s,j), *k_vector_[j]); // ktemp = ktemp + a_{s+1,j+1}*k_{j+1}
      }
    }
    TScalarMag ts = t_ + c(s)*dt;
    TScalarMag scaled_dt = (s == 0)? Scalar(dt/stages) : c(s)*dt;

    // need to check here the status of the solver (the linear solve)
    //eval_model_explicit<Scalar>(*model_,basePoint_,*ktemp_vector_,ts,Teuchos::outArg(*k_vector_[s]));
    eval_model_explicit<Scalar>(*model_,basePoint_,*ktemp_vector_,ts,Teuchos::outArg(*k_vector_[s]), scaled_dt, c(s));
    Thyra::Vt_S(k_vector_[s].ptr(),dt); // k_s = k_s*dt
  }
  // Sum for solution:
  for (int s=0 ; s < stages ; ++s) {
    if (b(s) != ST::zero()) {
      Thyra::Vp_StV(solution_vector_.ptr(), b(s), *k_vector_[s]); // solution_vector += b_{s+1}*k_{s+1}
    }
  }

   // cheat and say that the solver converged ( although no solver is needed for explicit method )
   rkNewtonConvergenceStatus_ = 0; 
   stepControl_->setCorrection(*this, solution_vector_, Teuchos::null , rkNewtonConvergenceStatus_);
   bool stepPass = stepControl_->acceptStep(*this, &LETvalue_);

   if (erkButcherTableau_->isEmbeddedMethod() ){ 
       
       Teuchos::SerialDenseVector<int,Scalar> bhat = erkButcherTableau_->bhat();

        // Sum for Embedded solution:
        for (int s=0 ; s < stages ; ++s) {
          if (bhat(s) != ST::zero()) {
             
              // solution_hat_vector += bhat_{s+1}*k_{s+1}
            Thyra::Vp_StV(solution_hat_vector_.ptr(), bhat(s), *k_vector_[s]); 
          }
        }

        // ee_ = (solution_vector__ - solution_hat_vector_)
        Thyra::V_VmV(ee_.ptr(), *solution_vector_, *solution_hat_vector_);
        
        stepControl_->setCorrection(*this, solution_vector_, ee_ , rkNewtonConvergenceStatus_);
   } 
   else {
        stepControl_->setCorrection(*this, solution_vector_, Teuchos::null, rkNewtonConvergenceStatus_);
   }

    stepPass = stepControl_->acceptStep(*this, &LETvalue_);

    if (!stepPass) { // stepPass = false
       stepLETStatus_ = STEP_LET_STATUS_FAILED;
       rkNewtonConvergenceStatus_ = -1; // just making sure here
    } else { // stepPass = true
       stepLETStatus_ = STEP_LET_STATUS_PASSED;
       rkNewtonConvergenceStatus_ = 0; // just making sure here
    }

   if (rkNewtonConvergenceStatus_ == 0) {

     // Update time range
     timeRange_ = timeRange(t,t+current_dt);
     numSteps_++;

     // completeStep only if the none of the stage solution's failed to converged
     stepControl_->completeStep(*this);

     dt_to_return = current_dt;

     } else {
        
        rkNewtonConvergenceStatus_ = -1;
        status = stepControl_-> rejectStep(*this); // reject the stage value
        (void) status; // avoid "set but not used" build warning
        dt_to_return = dt_old;
     }

  // update current time:
  t_ = t_ + dt;
  numSteps_++;
  return( dt_to_return );
}


template<class Scalar>
Scalar ExplicitRKStepper<Scalar>::takeFixedStep_(Scalar dt, StepSizeType flag)
{
  typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag;
  this->initialize_();
#ifdef HAVE_RYTHMOS_DEBUG
    TEUCHOS_TEST_FOR_EXCEPTION( flag == STEP_TYPE_VARIABLE, std::logic_error,
        "Error!  ExplicitRKStepper does not support variable time steps at this time."
        );
#endif // HAVE_RYTHMOS_DEBUG
  if ((flag == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
    return(Scalar(-ST::one()));
  }
  // Store old solution & old time
  V_V(solution_vector_old_.ptr(), *solution_vector_);
  t_old_ = t_;

  dt_ = dt;

  int stages = erkButcherTableau_->numStages();
  Teuchos::SerialDenseMatrix<int,Scalar> A = erkButcherTableau_->A();
  Teuchos::SerialDenseVector<int,Scalar> b = erkButcherTableau_->b();
  Teuchos::SerialDenseVector<int,Scalar> c = erkButcherTableau_->c();
  // Compute stage solutions
  for (int s=0 ; s < stages ; ++s) {
    Thyra::assign(ktemp_vector_.ptr(), *solution_vector_); // ktemp = solution_vector
    for (int j=0 ; j < s ; ++j) { // assuming Butcher matix is strictly lower triangular
      if (A(s,j) != ST::zero()) {
        Thyra::Vp_StV(ktemp_vector_.ptr(), A(s,j), *k_vector_[j]); // ktemp = ktemp + a_{s+1,j+1}*k_{j+1}
      }
    }
    TScalarMag ts = t_ + c(s)*dt;
    TScalarMag scaled_dt = (s == 0)? Scalar(dt/stages) : c(s)*dt;

    eval_model_explicit<Scalar>(*model_,basePoint_,*ktemp_vector_,ts,Teuchos::outArg(*k_vector_[s]), scaled_dt, c(s));
    Thyra::Vt_S(k_vector_[s].ptr(),dt); // k_s = k_s*dt
  }
  // Sum for solution:
  for (int s=0 ; s < stages ; ++s) {
    if (b(s) != ST::zero()) {
      Thyra::Vp_StV(solution_vector_.ptr(), b(s), *k_vector_[s]); // solution_vector += b_{s+1}*k_{s+1}
    }
  }

  // update current time:
  t_ = t_ + dt;

  numSteps_++;

  return(dt);
}

template<class Scalar>
const StepStatus<Scalar> ExplicitRKStepper<Scalar>::getStepStatus() const
{
  StepStatus<Scalar> stepStatus;

  if (!haveInitialCondition_) {
    stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
  } else if (numSteps_ == 0) {
    stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
    stepStatus.order = erkButcherTableau_->order();
    stepStatus.time = t_;
    stepStatus.solution = solution_vector_;
  } else {
    stepStatus.stepStatus = STEP_STATUS_CONVERGED;
    stepStatus.stepSize = dt_;
    stepStatus.order = erkButcherTableau_->order();
    stepStatus.time = t_;
    stepStatus.solution = solution_vector_;
  }

  return(stepStatus);
}

template<class Scalar>
void ExplicitRKStepper<Scalar>::describe(
    Teuchos::FancyOStream &out,
    const Teuchos::EVerbosityLevel verbLevel
    ) const
{
  if ( (static_cast<int>(verbLevel) == static_cast<int>(Teuchos::VERB_DEFAULT) ) ||
       (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)     )
     ) {
    out << this->description() << "::describe" << std::endl;
    out << "model = " << model_->description() << std::endl;
    out << erkButcherTableau_->numStages() << " stage Explicit RK method" << std::endl;
  } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)) {
    out << "solution_vector = " << std::endl;
    out << Teuchos::describe(*solution_vector_,verbLevel);
  } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM)) {
  } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_HIGH)) {
    out << "model = " << std::endl;
    out << Teuchos::describe(*model_,verbLevel);
    int stages = erkButcherTableau_->numStages();
    for (int i=0 ; i<stages ; ++i) {
      out << "k_vector[" << i << "] = " << std::endl;
      out << Teuchos::describe(*k_vector_[i],verbLevel);
    }
    out << "ktemp_vector = " << std::endl;
    out << Teuchos::describe(*ktemp_vector_,verbLevel);
    out << "ERK Butcher Tableau A matrix: " << printMat(erkButcherTableau_->A()) << std::endl;
    out << "ERK Butcher Tableau b vector: " << printMat(erkButcherTableau_->b()) << std::endl;
    out << "ERK Butcher Tableau c vector: " << printMat(erkButcherTableau_->c()) << std::endl;
    out << "t = " << t_ << std::endl;
  }
}

template<class Scalar>
void ExplicitRKStepper<Scalar>::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 */
    )
{
  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for ExplicitRKStepper at this time.\n");
}

template<class Scalar>
TimeRange<Scalar> ExplicitRKStepper<Scalar>::getTimeRange() const
{
  if (!haveInitialCondition_) {
    return(invalidTimeRange<Scalar>());
  } else {
    return(TimeRange<Scalar>(t_old_,t_));
  }
}

template<class Scalar>
void ExplicitRKStepper<Scalar>::getPoints(
  const Array<Scalar>& time_vec,
  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
  Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
  Array<ScalarMag>* accuracy_vec
  ) const
{
  TEUCHOS_ASSERT( haveInitialCondition_ );
  using Teuchos::constOptInArg;
  using Teuchos::null;
  defaultGetPoints<Scalar>(
      t_old_, constOptInArg(*solution_vector_old_),
      Ptr<const VectorBase<Scalar> >(null),
      t_, constOptInArg(*solution_vector_),
      Ptr<const VectorBase<Scalar> >(null),
      time_vec,ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
      Ptr<InterpolatorBase<Scalar> >(null)
      );
}

template<class Scalar>
void ExplicitRKStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
{
  TEUCHOS_ASSERT( time_vec != NULL );
  time_vec->clear();
  if (!haveInitialCondition_) {
    return;
  }
  time_vec->push_back(t_old_);
  if (t_ != t_old_) {
    time_vec->push_back(t_);
  }
}

template<class Scalar>
void ExplicitRKStepper<Scalar>::removeNodes(Array<Scalar>& /* time_vec */)
{
  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ExplicitRKStepper at this time.\n");
}

template<class Scalar>
int ExplicitRKStepper<Scalar>::getOrder() const
{
  return(erkButcherTableau_->order());
}

template <class Scalar>
void ExplicitRKStepper<Scalar>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
{
  TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
  parameterList_ = paramList;
  Teuchos::readVerboseObjectSublist(&*parameterList_,this);
}

template <class Scalar>
Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::getNonconstParameterList()
{
  return(parameterList_);
}

template <class Scalar>
Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::unsetParameterList()
{
  Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
  parameterList_ = Teuchos::null;
  return(temp_param_list);
}

template<class Scalar>
RCP<const Teuchos::ParameterList>
ExplicitRKStepper<Scalar>::getValidParameters() const
{
  using Teuchos::ParameterList;
  static RCP<const ParameterList> validPL;
  if (is_null(validPL)) {
    RCP<ParameterList> pl = Teuchos::parameterList();
    Teuchos::setupVerboseObjectSublist(&*pl);
    validPL = pl;
  }
  return validPL;
}

template<class Scalar>
void ExplicitRKStepper<Scalar>::setModel(const RCP<const Thyra::ModelEvaluator<Scalar> >& model)
{
  TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
  TEUCHOS_TEST_FOR_EXCEPT( !is_null(model_) ); // For now you can only call this once.
  assertValidModel( *this, *model );
  model_ = model;
}


template<class Scalar>
void ExplicitRKStepper<Scalar>::setNonconstModel(const RCP<Thyra::ModelEvaluator<Scalar> >& model)
{
  this->setModel(model); // TODO 09/09/09 tscoffe:  use ConstNonconstObjectContainer!
}


template<class Scalar>
Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
ExplicitRKStepper<Scalar>::getModel() const
{
  return model_;
}


template<class Scalar>
Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >
ExplicitRKStepper<Scalar>::getNonconstModel()
{
  return Teuchos::null;
}


template<class Scalar>
void ExplicitRKStepper<Scalar>::setInitialCondition(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
  )
{

  // typedef Thyra::ModelEvaluatorBase MEB; // unused

  basePoint_ = initialCondition;

  // x

  RCP<const Thyra::VectorBase<Scalar> >
    x_init = initialCondition.get_x();

#ifdef HAVE_RYTHMOS_DEBUG
  TEUCHOS_TEST_FOR_EXCEPTION(
    is_null(x_init), std::logic_error,
    "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
    "then x can not be null!" );
#endif

  solution_vector_ = x_init->clone_v();
  solution_vector_old_ = x_init->clone_v();

 // for the embedded RK method
  solution_hat_vector_ = x_init->clone_v();
  ee_ = x_init->clone_v();

  // t

  t_ = initialCondition.get_t();
  t_old_ = t_;

  haveInitialCondition_ = true;

}


template<class Scalar>
Thyra::ModelEvaluatorBase::InArgs<Scalar>
ExplicitRKStepper<Scalar>::getInitialCondition() const
{
  return basePoint_;
}

template<class Scalar>
bool ExplicitRKStepper<Scalar>::supportsCloning() const
{
  return true;
}

template<class Scalar>
RCP<StepperBase<Scalar> > ExplicitRKStepper<Scalar>::cloneStepperAlgorithm() const
{
  // Just use the interface to clone the algorithm in a basically
  // uninitialized state
  RCP<ExplicitRKStepper<Scalar> >
    stepper = Teuchos::rcp(new ExplicitRKStepper<Scalar>());

  if (!is_null(model_)) {
    stepper->setModel(model_); // Shallow copy is okay!
  }

  if (!is_null(erkButcherTableau_)) {
    // 06/16/09 tscoffe:  should we clone the RKBT here?
    stepper->setRKButcherTableau(erkButcherTableau_);
  }

  if (!is_null(parameterList_)) {
    stepper->setParameterList(Teuchos::parameterList(*parameterList_));
  }

  return stepper;

}

template<class Scalar>
RCP<StepControlStrategyBase<Scalar> > ExplicitRKStepper<Scalar>::getNonconstStepControlStrategy()
{
  return(stepControl_);
}

template<class Scalar>
RCP<const StepControlStrategyBase<Scalar> > ExplicitRKStepper<Scalar>::getStepControlStrategy() const
{
  return(stepControl_);
}

template<class Scalar>
void ExplicitRKStepper<Scalar>::setStepControlStrategy(const RCP<StepControlStrategyBase<Scalar> >& stepControl)
{
  TEUCHOS_TEST_FOR_EXCEPTION(stepControl == Teuchos::null,std::logic_error,"Error, stepControl == Teuchos::null!\n");
  stepControl_ = stepControl;
}
//
// Explicit Instantiation macro
//
// Must be expanded from within the Rythmos namespace!
//

#define RYTHMOS_EXPLICIT_RK_STEPPER_INSTANT(SCALAR) \
  \
  template class ExplicitRKStepper< SCALAR >; \
  \
  template RCP< ExplicitRKStepper< SCALAR > > \
  explicitRKStepper();  \
  \
  template RCP< ExplicitRKStepper< SCALAR > > \
  explicitRKStepper( \
    const RCP<Thyra::ModelEvaluator< SCALAR > >& model \
      ); \
  \
  template RCP< ExplicitRKStepper< SCALAR > > \
  explicitRKStepper( \
    const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
    const RCP<const RKButcherTableauBase< SCALAR > >& rkbt \
      ); \

} // namespace Rythmos

#endif //Rythmos_ExplicitRK_STEPPER_DEF_H