Rythmos_ImplicitRKStepper_decl.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_IMPLICIT_RK_STEPPER_DECL_H
#define Rythmos_IMPLICIT_RK_STEPPER_DECL_H

#include "Rythmos_Types.hpp"
#include "Rythmos_StepperBase.hpp"
#include "Rythmos_DataStore.hpp"
#include "Rythmos_SolverAcceptingStepperBase.hpp"
#include "Rythmos_RKButcherTableauAcceptingStepperBase.hpp"
#include "Rythmos_RKButcherTableauBase.hpp"

#include "Thyra_ModelEvaluator.hpp"
#include "Thyra_ProductVectorBase.hpp"
#include "Thyra_NonlinearSolverBase.hpp"

namespace Rythmos {


template<class Scalar>
class ImplicitRKStepper : 
  virtual public SolverAcceptingStepperBase<Scalar>,
  virtual public RKButcherTableauAcceptingStepperBase<Scalar>
{
public:
  
  typedef typename ScalarTraits<Scalar>::magnitudeType ScalarMag;
  

  ImplicitRKStepper();
  
  void set_W_factory( const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &irk_W_factory );

  RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
  
  
  void setRKButcherTableau( const RCP<const RKButcherTableauBase<Scalar> > &rkButcherTableau );

  RCP<const RKButcherTableauBase<Scalar> > getRKButcherTableau() const;


  // This function is mostly for testing purposes to explicitely over-ride the
  // internal RKBT detection to allow testing of 1-stage RKBTs as both fully
  // implicit RK and as DIRK methods.
  void setDirk(bool isDirk);

  

  void setSolver(
    const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
    );

  RCP<Thyra::NonlinearSolverBase<Scalar> >
  getNonconstSolver();

  RCP<const Thyra::NonlinearSolverBase<Scalar> >
  getSolver() const;



  bool isImplicit() const;
 
  bool supportsCloning() const;

  RCP<StepperBase<Scalar> > cloneStepperAlgorithm() 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();

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

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

  Scalar takeStep(Scalar dt, StepSizeType flag);
  
  const StepStatus<Scalar> getStepStatus() const;
  


  RCP<const Thyra::VectorSpaceBase<Scalar> >
  get_x_space() const;

  void addPoints(
    const Array<Scalar>& time_vec,
    const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
    const Array<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 setParameterList(RCP<ParameterList> const& paramList);
  
  RCP<ParameterList> getNonconstParameterList();
  
  RCP<ParameterList> unsetParameterList();
  
  RCP<const ParameterList> getValidParameters() const;
 

  
  void describe(
    FancyOStream  &out,
    const Teuchos::EVerbosityLevel verbLevel
    ) const;


private:

  // ///////////////////////
  // Private date members

  bool isInitialized_;
  RCP<const Thyra::ModelEvaluator<Scalar> > model_;
  RCP<Thyra::NonlinearSolverBase<Scalar> > solver_;
  RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > irk_W_factory_;
  RCP<ParameterList> paramList_;

  Thyra::ModelEvaluatorBase::InArgs<Scalar> basePoint_;
  RCP<Thyra::VectorBase<Scalar> > x_;
  RCP<Thyra::VectorBase<Scalar> > x_old_;
  RCP<Thyra::VectorBase<Scalar> > x_dot_;

  TimeRange<Scalar> timeRange_;

  RCP<Thyra::ModelEvaluator<Scalar> > irkModel_;
  RCP<const RKButcherTableauBase<Scalar> > irkButcherTableau_;

  bool isDirk_; // Used for Diagonal Implicit RK 

  int numSteps_;

  bool haveInitialCondition_;

  // Cache
  RCP<Thyra::ProductVectorBase<Scalar> > x_stage_bar_;

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

  void defaultInitializeAll_();
  void initialize_();

};


template<class Scalar>
RCP<ImplicitRKStepper<Scalar> >
implicitRKStepper();

template<class Scalar>
RCP<ImplicitRKStepper<Scalar> >
implicitRKStepper(
  const RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
  const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& irk_W_factory,
  const RCP<const RKButcherTableauBase<Scalar> >& irkbt
  );


} // namespace Rythmos

#endif //Rythmos_IMPLICIT_RK_STEPPER_DECL_H