Rythmos_ForwardSensitivityStepper.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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
//
// ***********************************************************************
//@HEADER

#ifndef RYTHMOS_FORWARD_SENSITIVITY_STEPPER_HPP
#define RYTHMOS_FORWARD_SENSITIVITY_STEPPER_HPP


#include "Rythmos_StepperBase.hpp"
#include "Rythmos_StepperHelpers.hpp"
#include "Rythmos_ForwardSensitivityModelEvaluatorBase.hpp"
#include "Rythmos_ForwardSensitivityImplicitModelEvaluator.hpp"
#include "Rythmos_ForwardSensitivityExplicitModelEvaluator.hpp"
#include "Rythmos_StateAndForwardSensitivityModelEvaluator.hpp"
#include "Rythmos_SolverAcceptingStepperBase.hpp"
#include "Rythmos_IntegratorBase.hpp"
#include "Rythmos_SingleResidualModelEvaluatorBase.hpp"
#include "Thyra_ModelEvaluatorHelpers.hpp"
#include "Thyra_LinearNonlinearSolver.hpp"
#include "Thyra_ProductVectorBase.hpp"
#include "Thyra_AssertOp.hpp"
#include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
#include "Teuchos_ConstNonconstObjectContainer.hpp"
#include "Teuchos_Assert.hpp"
#include "Teuchos_as.hpp"


namespace Rythmos {


template<class Scalar> 
class ForwardSensitivityStepper
  : virtual public StepperBase<Scalar>,
    virtual public Teuchos::ParameterListAcceptorDefaultBase
{
public:

  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;


  ForwardSensitivityStepper();

  void initializeSyncedSteppers(
    const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
    const int p_index,
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
    const RCP<StepperBase<Scalar> > &stateStepper,
    const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
    const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
    const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
    );

  void initializeSyncedSteppersInitCondOnly(
    const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
    const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space,
    const Thyra::ModelEvaluatorBase::InArgs<Scalar>& stateBasePoint,
    const RCP<StepperBase<Scalar> >& stateStepper,
    const RCP<Thyra::NonlinearSolverBase<Scalar> >& stateTimeStepSolver,
    const RCP<StepperBase<Scalar> >& sensStepper = Teuchos::null,
    const RCP<Thyra::NonlinearSolverBase<Scalar> >& sensTimeStepSolver = Teuchos::null
    );

  void initializeDecoupledSteppers(
    const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
    const int p_index,
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
    const RCP<StepperBase<Scalar> > &stateStepper,
    const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
    const RCP<IntegratorBase<Scalar> > &stateIntegrator,
    const Scalar &finalTime,
    const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
    const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
    );

  bool stateModelIsConst() const;

  RCP<const Thyra::ModelEvaluator<Scalar> >
  getStateModel() const;

  RCP<StepperBase<Scalar> >
  getNonconstStateStepper();

  RCP<const ForwardSensitivityModelEvaluatorBase<Scalar> >
  getFwdSensModel() const;

  RCP<const StateAndForwardSensitivityModelEvaluator<Scalar> >
  getStateAndFwdSensModel() const;



  void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
  RCP<const Teuchos::ParameterList> getValidParameters() const;



  bool acceptsModel() const;

  void setModel(
    const RCP<const Thyra::ModelEvaluator<Scalar> >& model
    );

  void setNonconstModel(
    const RCP<Thyra::ModelEvaluator<Scalar> >& model
    );

  RCP<const Thyra::ModelEvaluator<Scalar> > getModel() const;

  RCP<Thyra::ModelEvaluator<Scalar> > getNonconstModel();

  // RAB: ToDo: 2007/05/15: I need to talk with Todd about potentially
  // removing the setModel() and getModel() functions from the StepperBase
  // interface.  In the case of this forward sensitivity solver, I am not sure
  // that it makes a lot of sense to define a model.  This surely will not be
  // the model that a generic client would expect.  The assumption I am sure
  // would be that this model has the same space for x as the interpolation
  // buffer but that is not true in this case.

  void setInitialCondition(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
    );

  Thyra::ModelEvaluatorBase::InArgs<Scalar> getInitialCondition() const;

  Scalar takeStep( Scalar dt, StepSizeType stepType );

  const StepStatus<Scalar> getStepStatus() const;



  RCP<const Thyra::VectorSpaceBase<Scalar> >
  get_x_space() 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
    );

  TimeRange<Scalar> getTimeRange() const;

  void 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;

  void getNodes(Array<Scalar>* time_vec) const;

  void removeNodes(Array<Scalar>& time_vec);

  int getOrder() const;



  void initialize(
    const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
    const int p_index,
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
    const RCP<StepperBase<Scalar> > &stateStepper,
    const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
    const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
    const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
    )
    {
      initializeSyncedSteppers(
        stateModel, p_index, stateBasePoint, stateStepper, stateTimeStepSolver,
        sensStepper, sensTimeStepSolver
        );
    }


private:
  // ///////////////////
  // Private types

  typedef Teuchos::ConstNonconstObjectContainer<Thyra::ModelEvaluator<Scalar> > CNCME;

  // /////////////////////////
  // Private data members

  bool forceUpToDateW_;
  CNCME stateModel_;
  Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
  RCP<StepperBase<Scalar> > stateStepper_;
  RCP<Thyra::NonlinearSolverBase<Scalar> > stateTimeStepSolver_;
  RCP<IntegratorBase<Scalar> > stateIntegrator_;
  Scalar finalTime_;
  Thyra::ModelEvaluatorBase::InArgs<Scalar> stateAndSensBasePoint_;
  RCP<StepperBase<Scalar> > sensStepper_;
  RCP<Thyra::NonlinearSolverBase<Scalar> > sensTimeStepSolver_;

  bool isSingleResidualStepper_;
  RCP<ForwardSensitivityModelEvaluatorBase<Scalar> > sensModel_;
  RCP<StateAndForwardSensitivityModelEvaluator<Scalar> > stateAndSensModel_;
  Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_t_;

  static const std::string forceUpToDateW_name_;
  static const bool forceUpToDateW_default_;

  // /////////////////////////
  // Private member functions
  

  // Common initialization helper
  //
  // Preconditions:
  // (*) p_index >=0 or nonnull(p_space) == true
  //
  void initializeCommon(
    const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
    const int p_index,
    const RCP<const Thyra::VectorSpaceBase<Scalar> > &p_space,
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
    const RCP<StepperBase<Scalar> > &stateStepper,
    const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
    const RCP<StepperBase<Scalar> > &sensStepper,
    const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
    );

  Scalar takeSyncedStep( Scalar dt, StepSizeType stepType );

  Scalar takeDecoupledStep( Scalar dt, StepSizeType stepType );

};


// 2009/09/05: rabartl: ToDo: To fix the const and non-const handling of the
// stateModel in this class is going to be a lot of work but here is what
// needs to be done:
//
// (*) Duplicate each function that sets the stateModel, one for const and one
// for non-const.
//
// (*) Create a single a private version for each of these functions that
// accepts a Teuchos::ConstNonconstObjectContainer<> object and will implement
// the guts of the set up same as the existing functions.
//
// (*) Get all of the concrete StepperBase subclasses to implement the
// setModel(const RCP<const ME>&) and modelIsConst() functions and get them to
// use the Teuchos::ConstNonconstObjectContainer<> class as described above.
// This should be pretty easy as the only function that needs to be addressed
// in most cases is just the setModel(...) function.
//


template<class Scalar> 
inline
RCP<ForwardSensitivityStepper<Scalar> >
forwardSensitivityStepper()
{
  return Teuchos::rcp(new ForwardSensitivityStepper<Scalar>());
}


template<class Scalar> 
inline
RCP<ForwardSensitivityStepper<Scalar> >
forwardSensitivityStepper(
  const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
  const int p_index,
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
  const RCP<StepperBase<Scalar> > &stateStepper,
  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
  const RCP<StepperBase<Scalar> > &sensStepper = Teuchos::null,
  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver = Teuchos::null
  )
{
  RCP<ForwardSensitivityStepper<Scalar> >
    fwdSensStepper = Teuchos::rcp(new ForwardSensitivityStepper<Scalar>());
  fwdSensStepper->initializeSyncedSteppers(
    stateModel, p_index, stateBasePoint, stateStepper, stateTimeStepSolver );
  return fwdSensStepper;
}


template<class Scalar>
int getParameterIndex(
  const ForwardSensitivityStepper<Scalar> &fwdSensStepper
  )
{
  return fwdSensStepper.getFwdSensModel()->get_p_index();
}


template<class Scalar> 
inline
Thyra::ModelEvaluatorBase::InArgs<Scalar>
createStateAndSensInitialCondition(
  const ForwardSensitivityStepper<Scalar> &fwdSensStepper,
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_ic,
  const RCP<const Thyra::MultiVectorBase<Scalar> > S_init = Teuchos::null,
  const RCP<const Thyra::MultiVectorBase<Scalar> > S_dot_init = Teuchos::null
  )
{

  using Teuchos::outArg;
  using Thyra::assign;
  typedef Thyra::ModelEvaluatorBase MEB;
  
  RCP<const Thyra::VectorBase<Scalar> > s_bar_init;
  if (nonnull(S_init)) {
    s_bar_init = create_s_bar_given_S(*fwdSensStepper.getFwdSensModel(), S_init);
  }
  else {
    RCP<Thyra::VectorBase<Scalar> > s_bar_init_loc =
      createMember(fwdSensStepper.getFwdSensModel()->get_x_space());
    assign( outArg(*s_bar_init_loc), 0.0 );
    s_bar_init = s_bar_init_loc;
  }

  RCP<const Thyra::VectorBase<Scalar> > s_bar_dot_init;
  if (nonnull(S_dot_init)) {
    s_bar_dot_init = create_s_bar_given_S(*fwdSensStepper.getFwdSensModel(), S_dot_init);
  }
  else {
    RCP<Thyra::VectorBase<Scalar> > s_bar_dot_init_loc =
      createMember(fwdSensStepper.getFwdSensModel()->get_x_space());
    assign( outArg(*s_bar_dot_init_loc), 0.0 );
    s_bar_dot_init = s_bar_dot_init_loc;
  }
  
  RCP<const Rythmos::StateAndForwardSensitivityModelEvaluator<Scalar> >
    stateAndSensModel = fwdSensStepper.getStateAndFwdSensModel();
  
  MEB::InArgs<Scalar>
    state_and_sens_ic = fwdSensStepper.getModel()->createInArgs();
  
  // Copy time, parameters etc.
  state_and_sens_ic.setArgs(state_ic);
  // Set initial condition for x_bar = [ x; s_bar ]
  state_and_sens_ic.set_x(
    stateAndSensModel->create_x_bar_vec(state_ic.get_x(), s_bar_init)
    );
  // Set initial condition for x_bar_dot = [ x_dot; s_bar_dot ]
  state_and_sens_ic.set_x_dot(
    stateAndSensModel->create_x_bar_vec(state_ic.get_x_dot(), s_bar_dot_init)
    );

  return state_and_sens_ic;

}


template<class Scalar> 
inline
Thyra::ModelEvaluatorBase::InArgs<Scalar>
extractStateInitialCondition(
  const ForwardSensitivityStepper<Scalar> &fwdSensStepper,
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
  )
{

  using Thyra::productVectorBase;
  typedef Thyra::ModelEvaluatorBase MEB;

  MEB::InArgs<Scalar>
    state_ic = fwdSensStepper.getStateModel()->createInArgs();
  
  // Copy time, parameters etc.
  state_ic.setArgs(state_and_sens_ic);
  state_ic.set_x(
    productVectorBase(state_and_sens_ic.get_x())->getVectorBlock(0));
  state_ic.set_x_dot(
    productVectorBase(state_and_sens_ic.get_x_dot())->getVectorBlock(0));
  
  return state_ic;
  
}


//
// Implementation
//


// Static members


template<class Scalar> 
const std::string ForwardSensitivityStepper<Scalar>::forceUpToDateW_name_
= "Force Up-To-Date Jacobian";

template<class Scalar> 
const bool ForwardSensitivityStepper<Scalar>::forceUpToDateW_default_
= true;


// Constructors, Intializers, Misc.


template<class Scalar> 
ForwardSensitivityStepper<Scalar>::ForwardSensitivityStepper()
  :forceUpToDateW_(false),
   isSingleResidualStepper_(false)
{}


template<class Scalar> 
void ForwardSensitivityStepper<Scalar>::initializeSyncedSteppers(
  const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
  const int p_index,
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
  const RCP<StepperBase<Scalar> > &stateStepper,
  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
  const RCP<StepperBase<Scalar> > &sensStepper,
  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
  )
  
{
  initializeCommon( stateModel, p_index, Teuchos::null, stateBasePoint, stateStepper,
    stateTimeStepSolver, sensStepper, sensTimeStepSolver );
}


template<class Scalar> 
void ForwardSensitivityStepper<Scalar>::initializeSyncedSteppersInitCondOnly(
  const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
  const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space,
  const Thyra::ModelEvaluatorBase::InArgs<Scalar>& stateBasePoint,
  const RCP<StepperBase<Scalar> >& stateStepper,
  const RCP<Thyra::NonlinearSolverBase<Scalar> >& stateTimeStepSolver,
  const RCP<StepperBase<Scalar> >& sensStepper,
  const RCP<Thyra::NonlinearSolverBase<Scalar> >& sensTimeStepSolver
  )
{
  initializeCommon(stateModel, -1, p_space, stateBasePoint, stateStepper,
    stateTimeStepSolver, sensStepper, sensTimeStepSolver );
}


template<class Scalar> 
void ForwardSensitivityStepper<Scalar>::initializeDecoupledSteppers(
  const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
  const int p_index,
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
  const RCP<StepperBase<Scalar> > &stateStepper,
  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
  const RCP<IntegratorBase<Scalar> > &stateIntegrator,
  const Scalar &finalTime,
  const RCP<StepperBase<Scalar> > &sensStepper,
  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
  )
{
  TEUCHOS_ASSERT(nonnull(stateIntegrator));
  initializeCommon( stateModel, p_index, Teuchos::null, stateBasePoint, stateStepper,
    stateTimeStepSolver, sensStepper, sensTimeStepSolver );
  stateIntegrator_ = stateIntegrator;
  finalTime_ = finalTime;
}


template<class Scalar> 
bool ForwardSensitivityStepper<Scalar>::stateModelIsConst() const
{
  return stateModel_.isConst();
}


template<class Scalar> 
RCP<const Thyra::ModelEvaluator<Scalar> >
ForwardSensitivityStepper<Scalar>::getStateModel() const
{
  return stateModel_.getConstObj();
}

  
template<class Scalar> 
RCP<StepperBase<Scalar> >
ForwardSensitivityStepper<Scalar>::getNonconstStateStepper()
{
  return stateStepper_;
}

  
template<class Scalar> 
RCP<const ForwardSensitivityModelEvaluatorBase<Scalar> >
ForwardSensitivityStepper<Scalar>::getFwdSensModel() const
{
  return sensModel_;
}

  
template<class Scalar> 
RCP<const StateAndForwardSensitivityModelEvaluator<Scalar> >
ForwardSensitivityStepper<Scalar>::getStateAndFwdSensModel() const
{
  return stateAndSensModel_;
}


// Overridden from Teuchos::ParameterListAcceptor


template<class Scalar> 
void ForwardSensitivityStepper<Scalar>::setParameterList(
  RCP<Teuchos::ParameterList> const& paramList
  )
{
  TEST_FOR_EXCEPT(is_null(paramList));
  paramList->validateParameters(*getValidParameters());
  this->setMyParamList(paramList);
  forceUpToDateW_ = paramList->get(forceUpToDateW_name_,forceUpToDateW_default_);
  Teuchos::readVerboseObjectSublist(&*paramList,this);
}


template<class Scalar> 
RCP<const Teuchos::ParameterList>
ForwardSensitivityStepper<Scalar>::getValidParameters() const
{
  static RCP<const ParameterList> validPL;
  if (is_null(validPL) ) {
    RCP<ParameterList> pl = Teuchos::parameterList();
    pl->set( forceUpToDateW_name_, forceUpToDateW_default_,
      "If set to true, then the Jacobian matrix W used in the\n"
      "state timestep equation will be forced to be up to date\n"
      "with the final value for x for the nonlinear solve.  If\n"
      "you are willing to live with slightly less accurate sensitivities\n"
      "then set this to false."
      );
    Teuchos::setupVerboseObjectSublist(&*pl);
    validPL = pl;
  }
  return validPL;
}


// Overridden from StepperBase

template<class Scalar> 
bool ForwardSensitivityStepper<Scalar>::acceptsModel() const
{
  return false;
}

template<class Scalar> 
void ForwardSensitivityStepper<Scalar>::setModel(
  const RCP<const Thyra::ModelEvaluator<Scalar> >& model
  )
{
  TEST_FOR_EXCEPT_MSG( true,
    "Error, this stepper subclass does not accept a model"
    " as defined by the StepperBase interface!");
}


template<class Scalar> 
void ForwardSensitivityStepper<Scalar>::setNonconstModel(
  const RCP<Thyra::ModelEvaluator<Scalar> >& model
  )
{
  TEST_FOR_EXCEPT_MSG( true,
    "Error, this stepper subclass does not accept a model"
    " as defined by the StepperBase interface!");
}


template<class Scalar> 
RCP<const Thyra::ModelEvaluator<Scalar> >
ForwardSensitivityStepper<Scalar>::getModel() const
{
  return stateAndSensModel_;
}


template<class Scalar> 
RCP<Thyra::ModelEvaluator<Scalar> >
ForwardSensitivityStepper<Scalar>::getNonconstModel() 
{
  return stateAndSensModel_;
}


template<class Scalar> 
void ForwardSensitivityStepper<Scalar>::setInitialCondition(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &state_and_sens_ic
  )
{
  
  typedef Thyra::ModelEvaluatorBase MEB;

  stateAndSensBasePoint_ = state_and_sens_ic;

  // Get the product vectors for x_bar = [ x; s_bar ] and x_bar_dot

  TEST_FOR_EXCEPTION(
    is_null(state_and_sens_ic.get_x()), std::logic_error,
    "Error, the initial condition for x_bar = [ x; s_bar ] can not be null!" );

  const RCP<const Thyra::ProductVectorBase<Scalar> >
    x_bar_init = Thyra::productVectorBase<Scalar>(
      state_and_sens_ic.get_x()
      );

  RCP<const Thyra::ProductVectorBase<Scalar> > x_bar_dot_init;
  if (state_and_sens_ic.supports(MEB::IN_ARG_x_dot)) {
      x_bar_dot_init = Thyra::productVectorBase<Scalar>(
        state_and_sens_ic.get_x_dot()
        );
  }

  // Remove x and x_dot from state_and_sens_ic_in to avoid cloning x and x dot!
  
  Thyra::ModelEvaluatorBase::InArgs<Scalar>
    state_and_sens_ic_no_x = state_and_sens_ic;
  state_and_sens_ic_no_x.set_x(Teuchos::null);
  if (state_and_sens_ic_no_x.supports(MEB::IN_ARG_x_dot)) {
    state_and_sens_ic_no_x.set_x_dot(Teuchos::null);
  }

  // Set initial condition for the state

  MEB::InArgs<Scalar> state_ic = stateModel_->createInArgs();
  state_ic.setArgs(state_and_sens_ic_no_x,true,true); // Set time, parameters etc.
  state_ic.set_x(x_bar_init->getVectorBlock(0)->clone_v());
  if (state_ic.supports(MEB::IN_ARG_x_dot)) {
    state_ic.set_x_dot(
        !is_null(x_bar_dot_init)
        ? x_bar_dot_init->getVectorBlock(0)->clone_v()
        : Teuchos::null
        );
  }
  stateStepper_->setInitialCondition(state_ic);

  // Set up the integrator if needed
  //if (!is_null(stateIntegrator_)) {
  //  stateIntegrator_->setStepper( stateStepper_, finalTime_ );
  //  sensModel_->setStateIntegrator( stateIntegrator_, state_ic );
  //}

  // Set initial condition for the sensitivities
  
  MEB::InArgs<Scalar> sens_ic = sensModel_->createInArgs();
  sens_ic.setArgs(state_and_sens_ic_no_x,true,true); // Set time etc.
  sens_ic.set_x(x_bar_init->getVectorBlock(1)->clone_v());
  if (sens_ic.supports(MEB::IN_ARG_x_dot)) {
    sens_ic.set_x_dot(
        !is_null(x_bar_dot_init)
        ? x_bar_dot_init->getVectorBlock(1)->clone_v()
        : Teuchos::null
        );
  }
  sensStepper_->setInitialCondition(sens_ic);

}


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

template<class Scalar> 
Scalar
ForwardSensitivityStepper<Scalar>::takeStep(
  Scalar dt, StepSizeType stepType
  )
{

#ifdef ENABLE_RYTHMOS_TIMERS
  TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep");
#endif

  if (!is_null(stateIntegrator_)) {
    return takeDecoupledStep(dt,stepType);
  }

  return takeSyncedStep(dt,stepType);

}


template<class Scalar> 
const StepStatus<Scalar>
ForwardSensitivityStepper<Scalar>::getStepStatus() const
{

  const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
  StepStatus<Scalar> stepStatus;
  
  stepStatus.message = sensStepStatus.message;
  stepStatus.stepStatus = sensStepStatus.stepStatus;
  stepStatus.stepLETStatus = sensStepStatus.stepLETStatus;
  stepStatus.stepSize = sensStepStatus.stepSize;
  stepStatus.order = sensStepStatus.order;
  stepStatus.time = sensStepStatus.time;
  stepStatus.stepLETValue = sensStepStatus.stepLETValue;
  stepStatus.extraParameters = sensStepStatus.extraParameters;

  if (is_null(stateIntegrator_)) {
    const StepStatus<Scalar> 
      stateStepStatus = stateStepper_->getStepStatus();
    if (!is_null(stateStepStatus.solution) && !is_null(sensStepStatus.solution))
      stepStatus.solution = stateAndSensModel_->create_x_bar_vec(
        stateStepStatus.solution, sensStepStatus.solution
        );
    if (!is_null(stateStepStatus.solutionDot) && !is_null(sensStepStatus.solutionDot))
      stepStatus.solutionDot = stateAndSensModel_->create_x_bar_vec(
        stateStepStatus.solutionDot, sensStepStatus.solutionDot
        );
  }

  return stepStatus;

}


// Overridden from InterpolationBufferBase


template<class Scalar> 
RCP<const Thyra::VectorSpaceBase<Scalar> >
ForwardSensitivityStepper<Scalar>::get_x_space() const
{
  return stateAndSensModel_->get_x_space();
}


template<class Scalar> 
void ForwardSensitivityStepper<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
  )
{
  TEST_FOR_EXCEPT("Not implemented addPoints(...) yet but we could if we wanted!");
}


template<class Scalar> 
TimeRange<Scalar>
ForwardSensitivityStepper<Scalar>::getTimeRange() const
{
  return sensStepper_->getTimeRange();
}


template<class Scalar> 
void ForwardSensitivityStepper<Scalar>::getPoints(
  const Array<Scalar>& time_vec,
  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_vec,
  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_dot_vec,
  Array<ScalarMag>* accuracy_vec
  ) const
{

  using Teuchos::as;

#ifdef RYTHMOS_DEBUG
  TEST_FOR_EXCEPT( as<int>(time_vec.size()) == 0 );
#endif

  const int numTimePoints = time_vec.size();

  if (x_bar_vec)
    x_bar_vec->clear();

  if (x_bar_dot_vec)
    x_bar_dot_vec->clear();
  
  Array<RCP<const Thyra::VectorBase<Scalar> > >
    x_vec, x_dot_vec;

  if (!is_null(stateIntegrator_)) {
    stateIntegrator_->getPoints(
      time_vec,
      x_bar_vec ? &x_vec: 0,
      x_bar_dot_vec ? &x_dot_vec: 0,
      0 // Ignoring accuracy from state for now!
      );
  }
  else {
    stateStepper_->getPoints(
      time_vec,
      x_bar_vec ? &x_vec: 0,
      x_bar_dot_vec ? &x_dot_vec: 0,
      0 // Ignoring accuracy from state for now!
      );
  }
  
  Array<RCP<const Thyra::VectorBase<Scalar> > >
    s_bar_vec, s_bar_dot_vec;

  sensStepper_->getPoints(
    time_vec,
    x_bar_vec ? &s_bar_vec: 0,
    x_bar_dot_vec ? &s_bar_dot_vec: 0,
    accuracy_vec
    );
  
  if ( x_bar_vec ) {
    for ( int i = 0; i < numTimePoints; ++i ) {
      x_bar_vec->push_back(
        stateAndSensModel_->create_x_bar_vec(x_vec[i],s_bar_vec[i])
        );
    }
  }
  
  if ( x_bar_dot_vec ) {
    for ( int i = 0; i < numTimePoints; ++i ) {
      x_bar_dot_vec->push_back(
        stateAndSensModel_->create_x_bar_vec(x_dot_vec[i],s_bar_dot_vec[i])
        );
    }
  }

}


template<class Scalar>
void ForwardSensitivityStepper<Scalar>::getNodes(
  Array<Scalar>* time_vec
  ) const
{
  TEUCHOS_ASSERT( time_vec != NULL );
  time_vec->clear();
  if (is_null(stateIntegrator_) && is_null(stateStepper_)) {
    return;
  }
  if (!is_null(stateIntegrator_)) {
    stateIntegrator_->getNodes(time_vec);
  }
  else {
    stateStepper_->getNodes(time_vec);
  }
}


template<class Scalar> 
void ForwardSensitivityStepper<Scalar>::removeNodes(
  Array<Scalar>& time_vec
  )
{
  TEST_FOR_EXCEPT("Not implemented yet but we can!");
}


template<class Scalar> 
int ForwardSensitivityStepper<Scalar>::getOrder() const
{
  return sensStepper_->getOrder();
  // Note: This assumes that stateStepper will have the same order!
}


// private


template<class Scalar> 
void ForwardSensitivityStepper<Scalar>::initializeCommon(
  const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
  const int p_index,
  const RCP<const Thyra::VectorSpaceBase<Scalar> > &p_space,
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
  const RCP<StepperBase<Scalar> > &stateStepper,
  const RCP<Thyra::NonlinearSolverBase<Scalar> > &stateTimeStepSolver,
  const RCP<StepperBase<Scalar> > &sensStepper,
  const RCP<Thyra::NonlinearSolverBase<Scalar> > &sensTimeStepSolver
  )
{

  using Teuchos::rcp_implicit_cast;
  using Teuchos::rcp_dynamic_cast;

  typedef Thyra::ModelEvaluatorBase MEB;

  //
  // Validate input
  //

  TEUCHOS_ASSERT( p_index >= 0 || nonnull(p_space) );
  if (nonnull(p_space)) {
    TEUCHOS_ASSERT_EQUALITY(p_index, -1);
  }
  if (p_index >= 0) {
    TEUCHOS_ASSERT(is_null(p_space));
  }
  TEST_FOR_EXCEPT( is_null(stateModel) );
  TEST_FOR_EXCEPT( is_null(stateStepper) );
  if (stateStepper->isImplicit()) {
    TEST_FOR_EXCEPT( is_null(stateTimeStepSolver) ); // allow to be null for explicit methods
  }

  //
  // Create the sensModel which will do some more validation
  //
  
  RCP<ForwardSensitivityModelEvaluatorBase<Scalar> > sensModel;
  MEB::InArgs<Scalar> stateModelInArgs = stateModel->createInArgs();
  if (stateModelInArgs.supports(MEB::IN_ARG_x_dot)) {
    // Implicit DE formulation
    sensModel = Teuchos::rcp(new ForwardSensitivityImplicitModelEvaluator<Scalar>);
  }
  else {
    // Explicit DE formulation
    sensModel = Teuchos::rcp(new ForwardSensitivityExplicitModelEvaluator<Scalar>);
  }

  if (p_index >= 0) {
    sensModel->initializeStructure(stateModel, p_index);
  }
  else {
    sensModel->initializeStructureInitCondOnly(stateModel, p_space);
  }
  
  //
  // Get the input objects
  //

  stateModel_.initialize(stateModel);

  stateBasePoint_ = stateBasePoint;

  stateStepper_ = stateStepper;
  
  stateTimeStepSolver_ = stateTimeStepSolver;

  sensModel_ = sensModel;

  stateAndSensModel_ = Teuchos::rcp(new StateAndForwardSensitivityModelEvaluator<Scalar>);
  stateAndSensModel_->initializeStructure(sensModel_);

  if (!is_null(sensStepper)) {
    sensStepper_ = sensStepper;
  }
  else {
    sensStepper_ = stateStepper_->cloneStepperAlgorithm();
    TEST_FOR_EXCEPTION(
      is_null(sensStepper_), std::logic_error,
      "Error, if the client does not pass in a stepper for the senitivity\n"
      "equations then the stateStepper object must support cloning to create\n"
      "the sensitivity stepper!"
      );
  }

  if (!is_null(sensTimeStepSolver)) {
    sensTimeStepSolver_ = sensTimeStepSolver;
  }
  else {
    RCP<Thyra::LinearNonlinearSolver<Scalar> >
      linearNonlinearSolver(new Thyra::LinearNonlinearSolver<Scalar>);
    // ToDo: Set tolerance on the nonlinear solver???
    sensTimeStepSolver_ = linearNonlinearSolver;
  }

  //
  // Setup the steppers
  //

  isSingleResidualStepper_ = true; // ToDo: Add dynamic cast on
                                   // stateTimeStepSolver to check this!

  setStepperModel(Teuchos::inOutArg(*stateStepper_),stateModel_);
  if (stateStepper_->isImplicit()) {
    rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
        stateStepper_,true)->setSolver(stateTimeStepSolver_);
  }
  sensStepper_->setModel(sensModel_); 
  if (sensStepper_->isImplicit()) {
    rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
        sensStepper_,true)->setSolver(sensTimeStepSolver_);
  }

  stateBasePoint_t_ = stateModel_->createInArgs();

  // 2007/05/18: rabartl: ToDo: Move the above initialization code to give
  // setInitializeCondition(...) a chance to set the initial condition.

}


template<class Scalar> 
Scalar ForwardSensitivityStepper<Scalar>::takeSyncedStep(
  Scalar dt, StepSizeType stepType
  )
{

#ifdef ENABLE_RYTHMOS_TIMERS
  TEUCHOS_FUNC_TIME_MONITOR_DIFF("Rythmos:ForwardSensitivityStepper::takeStep: synced",
    TopLevel);
#endif

  using Teuchos::as;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSIBB;
  typedef Thyra::ModelEvaluatorBase MEB;

  RCP<Teuchos::FancyOStream> out = this->getOStream();
  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
  const bool lowTrace =
    ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) );
  const bool mediumTrace =
    ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) );
  Teuchos::OSTab tab(out);

  if (lowTrace) {
    *out
      << "\nEntering " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
      << "::takeSyncedStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
  }

  //
  // A) Compute the state timestep
  //

  if (lowTrace) {
    *out
      << "\nTaking state step using stepper : "
      << stateStepper_->description() << "\n";
  }

  Scalar state_dt = -1.0;
  {
#ifdef ENABLE_RYTHMOS_TIMERS
    TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: stateStep");
#endif
    VOTSIBB stateStepper_outputTempState(stateStepper_,out,verbLevel);
    state_dt = stateStepper_->takeStep(dt,stepType);
  }

  if (state_dt < Scalar(-ST::one())) {
    if (lowTrace)
      *out << "\nThe state stepper has failed so return a failed timestep!\n";
    return state_dt;
  }

  {
#ifdef ENABLE_RYTHMOS_TIMERS
    TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: updateSensModel");
#endif
    // Set up the sensitivity model for this timestep
    sensModel_->initializePointState(Teuchos::inOutArg(*stateStepper_),forceUpToDateW_);
  } 

  //
  // C) Compute the sensitivity timestep for the exact same timestep as was
  // used for the state solve.
  //

  if (lowTrace) {
    *out
      << "\nTaking sensitivity step using stepper : "
      << sensStepper_->description() << "\n";
  }

  Scalar sens_dt = -1.0;
  {
#ifdef ENABLE_RYTHMOS_TIMERS
    TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: sensStep");
#endif
    // Copy the step control data to make sure that the sensStepper takes the
    // same type of step that the statStepper took.  This is needed to ensure
    // that the W matrix is the same for one.
    sensStepper_->setStepControlData(*stateStepper_);
    VOTSIBB sensStepper_outputTempState(sensStepper_,out,verbLevel);
    sens_dt = sensStepper_->takeStep(state_dt,STEP_TYPE_FIXED);
  }

  if (mediumTrace) {
    const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
    *out << "\nSensitivity step status:\n" << sensStepStatus;
  }
  
  TEST_FOR_EXCEPTION(
    sens_dt != state_dt, std::logic_error,
    "Error, the sensitivity step failed for some reason.  We should\n"
    "just return a negative step size and reject the step but currently\n"
    "there is no way to roll back the state timestep it for back to\n"
    "the status before this function was called!"
    );

  // 2007/05/18: rabartl: ToDo: If stepType == STEP_TYPE_VARIABLE and the state
  // timestep sucessed but the sensitivity timestep failed, then we need to
  // really throw an excpetion because there is nothing that we can really do
  // here!

  // 2007/05/18: rabartl: ToDo: Replace the above std::logic_error type with
  // a Rythmos::CatastrophicFailure or just use Thyra::CatastrophicFailure!

  if (lowTrace) {
    *out
      << "\nLeaving " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
      << "::takeSyncedStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
  }

  return state_dt;

}


template<class Scalar> 
Scalar ForwardSensitivityStepper<Scalar>::takeDecoupledStep(
  Scalar dt, StepSizeType stepType
  )
{

#ifdef ENABLE_RYTHMOS_TIMERS
  TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityStepper::takeStep: decoupled");
#endif

  using Teuchos::as;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSIBB;
  typedef Thyra::ModelEvaluatorBase MEB;

  RCP<Teuchos::FancyOStream> out = this->getOStream();
  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
  const bool lowTrace =
    ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) );
  const bool mediumTrace =
    ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) );
  Teuchos::OSTab tab(out);

  if (lowTrace) {
    *out
      << "\nEntering " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
      << "::takeDecoupledStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
  }
  
  //
  // A) Take the sens timestep
  //

  if (lowTrace) {
    *out
      << "\nTaking sensitivity step using stepper : "
      << sensStepper_->description() << "\n";
  }

  Scalar sens_dt = -1.0;
  VOTSIBB sensStepper_outputTempState(sensStepper_,out,verbLevel);
  sens_dt = sensStepper_->takeStep(dt,stepType);

  if (mediumTrace) {
    const StepStatus<Scalar> sensStepStatus = sensStepper_->getStepStatus();
    *out << "\nSensitivity step status:\n" << sensStepStatus;
  }

  //
  // B) Wipe out all state interp buffer info before this sens timestep
  //
  
  //TEST_FOR_EXCEPT(true);

  if (lowTrace) {
    *out
      << "\nLeaving " << TypeNameTraits<ForwardSensitivityStepper<Scalar> >::name()
      << "::takeDecoupledStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
  }

  return sens_dt;
  
}


} // namespace Rythmos


#endif //RYTHMOS_FORWARD_SENSITIVITY_STEPPER_HPP