Rythmos_ExplicitTaylorPolynomialStepper.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_EXPLICIT_TAYLOR_POLYNOMIAL_STEPPER_H
#define RYTHMOS_EXPLICIT_TAYLOR_POLYNOMIAL_STEPPER_H

#include "Rythmos_StepperBase.hpp"
#include "Rythmos_StepperHelpers.hpp"
#include "Teuchos_RCP.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
#include "Thyra_VectorBase.hpp"
#include "Thyra_ModelEvaluator.hpp"
#include "Thyra_ModelEvaluatorHelpers.hpp"
#include "Thyra_PolynomialVectorTraits.hpp"
#include "RTOpPack_RTOpTHelpers.hpp"

namespace Rythmos {


RTOP_ROP_1_REDUCT_SCALAR( ROpLogNormInf,
  typename ScalarTraits<Scalar>::magnitudeType, // Reduction object type
  RTOpPack::REDUCT_TYPE_MAX // Reduction object reduction
  )
{
  using Teuchos::as;
  typedef ScalarTraits<Scalar> ST;
  typedef typename ST::magnitudeType ScalarMag;
  const ScalarMag mag = std::log(as<ScalarMag>(1e-100) + ST::magnitude(v0));
  reduct = TEUCHOS_MAX( mag, reduct );
}

template<class Scalar>
class ExplicitTaylorPolynomialStepper : virtual public StepperBase<Scalar>
{
public:

  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
    
  ExplicitTaylorPolynomialStepper();
    
  ~ExplicitTaylorPolynomialStepper();

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


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

  RCP<Teuchos::ParameterList> getNonconstParameterList();

  RCP<Teuchos::ParameterList> unsetParameterList();

  RCP<const Teuchos::ParameterList> getValidParameters() const;

  std::string description() const;

  void describe(
    Teuchos::FancyOStream &out,
    const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default
    ) 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
    );
    
  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 setRange(
    const TimeRange<Scalar>& range,
    const InterpolationBufferBase<Scalar> & IB
    );

  TimeRange<Scalar> getTimeRange() const;

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

  void removeNodes(Array<Scalar>& time_vec);

  int getOrder() const;

private:

  void defaultInitializAll_();

  void computeTaylorSeriesSolution_();

  ScalarMag estimateLogRadius_();

  RCP<const Thyra::ModelEvaluator<Scalar> > model_;

  RCP<Teuchos::ParameterList> parameterList_;

  RCP<Thyra::VectorBase<Scalar> > x_vector_;

  RCP<Thyra::VectorBase<Scalar> > x_vector_old_;

  RCP<Thyra::VectorBase<Scalar> > x_dot_vector_;

  RCP<Thyra::VectorBase<Scalar> > x_dot_vector_old_;

  RCP<Thyra::VectorBase<Scalar> > f_vector_;

  RCP<Teuchos::Polynomial<Thyra::VectorBase<Scalar> > > x_poly_;

  RCP<Teuchos::Polynomial<Thyra::VectorBase<Scalar> > > f_poly_;

  Thyra::ModelEvaluatorBase::InArgs<Scalar> basePoint_;

  bool haveInitialCondition_;

  int numSteps_;

  Scalar t_;

  Scalar dt_;

  Scalar t_initial_;

  Scalar t_final_;

  ScalarMag local_error_tolerance_;

  Scalar min_step_size_;

  Scalar max_step_size_;

  unsigned int degree_;

  Scalar linc_;
};


template <typename Scalar>
typename Teuchos::ScalarTraits<Scalar>::magnitudeType
log_norm_inf(const Thyra::VectorBase<Scalar>& x)
{
  ROpLogNormInf<Scalar> log_norm_inf_op;
  RCP<RTOpPack::ReductTarget> log_norm_inf_targ = 
    log_norm_inf_op.reduct_obj_create();
  Thyra::applyOp<Scalar>(log_norm_inf_op,
    Teuchos::tuple(Teuchos::ptrFromRef(x))(), Teuchos::null,
    log_norm_inf_targ.ptr());
  return log_norm_inf_op(*log_norm_inf_targ);
}


// Non-member constructor
template<class Scalar>
RCP<ExplicitTaylorPolynomialStepper<Scalar> > explicitTaylorPolynomialStepper()
{
  RCP<ExplicitTaylorPolynomialStepper<Scalar> > stepper = rcp(new ExplicitTaylorPolynomialStepper<Scalar>());
  return stepper;
}


template<class Scalar>
ExplicitTaylorPolynomialStepper<Scalar>::ExplicitTaylorPolynomialStepper()
{
  this->defaultInitializAll_();
  numSteps_ = 0;
}


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


template<class Scalar>
void ExplicitTaylorPolynomialStepper<Scalar>::defaultInitializAll_()
{
  typedef Teuchos::ScalarTraits<Scalar> ST;
  Scalar nan = ST::nan();
  model_ = Teuchos::null;
  parameterList_ = Teuchos::null;
  x_vector_ = Teuchos::null;
  x_vector_old_ = Teuchos::null;
  x_dot_vector_ = Teuchos::null;
  x_dot_vector_old_ = Teuchos::null;
  f_vector_ = Teuchos::null;
  x_poly_ = Teuchos::null;
  f_poly_ = Teuchos::null;
  haveInitialCondition_ = false;
  numSteps_ = -1;
  t_ = nan;
  dt_ = nan;
  t_initial_ = nan;
  t_final_ = nan;
  local_error_tolerance_ = nan;
  min_step_size_ = nan;
  max_step_size_ = nan;
  degree_ = 0;
  linc_ = nan;
}


template<class Scalar>
void ExplicitTaylorPolynomialStepper<Scalar>::setModel(
  const RCP<const Thyra::ModelEvaluator<Scalar> >& model
  )
{
  TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
  assertValidModel( *this, *model );
    
  model_ = model;
  f_vector_ = Thyra::createMember(model_->get_f_space());
}


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


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


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


template<class Scalar>
void ExplicitTaylorPolynomialStepper<Scalar>::setInitialCondition(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
  )
{
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef Thyra::ModelEvaluatorBase MEB;
  basePoint_ = initialCondition;
  if (initialCondition.supports(MEB::IN_ARG_t)) {
    t_ = initialCondition.get_t();
  } else {
    t_ = ST::zero();
  }
  dt_ = ST::zero();
  x_vector_ = initialCondition.get_x()->clone_v();
  x_dot_vector_ = x_vector_->clone_v();
  x_vector_old_ = x_vector_->clone_v();
  x_dot_vector_old_ = x_dot_vector_->clone_v();
  haveInitialCondition_ = true;
}


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


template<class Scalar>
Scalar 
ExplicitTaylorPolynomialStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
{
  typedef Teuchos::ScalarTraits<Scalar> ST;
  TEUCHOS_ASSERT( haveInitialCondition_ );
  TEUCHOS_ASSERT( !is_null(model_) );
  TEUCHOS_ASSERT( !is_null(parameterList_) ); // parameters are nan otherwise

  V_V(outArg(*x_vector_old_),*x_vector_); // x_vector_old = x_vector
  V_V(outArg(*x_dot_vector_old_),*x_dot_vector_); // x_dot_vector_old = x_dot_vector

  if (x_poly_ == Teuchos::null) {
    x_poly_ = Teuchos::rcp(new Teuchos::Polynomial<Thyra::VectorBase<Scalar> >(0,*x_vector_,degree_));
  }

  if (f_poly_ == Teuchos::null) {
    f_poly_ = Teuchos::rcp(new Teuchos::Polynomial<Thyra::VectorBase<Scalar> >(0, *f_vector_, degree_));
  }
  if (flag == STEP_TYPE_VARIABLE) {
    // If t_ > t_final_, we're done
    if (t_ > t_final_) {
      dt_ = ST::zero();
      return dt_;
    }

    // Compute a local truncated Taylor series solution to system
    computeTaylorSeriesSolution_();

    // Estimate log of radius of convergence of Taylor series
    Scalar rho = estimateLogRadius_();

    // Set step size
    Scalar shadowed_dt = std::exp(linc_ - rho);

    // If step size is too big, reduce
    if (shadowed_dt > max_step_size_) {
      shadowed_dt = max_step_size_;
    }

    // If step goes past t_final_, reduce
    if (t_+shadowed_dt > t_final_) {
      shadowed_dt = t_final_-t_;
    }

    ScalarMag local_error;

    do {

      // compute x(t_+shadowed_dt), xdot(t_+shadowed_dt)
      x_poly_->evaluate(shadowed_dt, x_vector_.get(), x_dot_vector_.get());

      // compute f( x(t_+shadowed_dt), t_+shadowed_dt )
      eval_model_explicit<Scalar>(*model_,basePoint_,*x_vector_,t_+shadowed_dt,Teuchos::outArg(*f_vector_));

      // compute || xdot(t_+shadowed_dt) - f( x(t_+shadowed_dt), t_+shadowed_dt ) ||
      Thyra::Vp_StV(x_dot_vector_.ptr(), -ST::one(),
        *f_vector_);
      local_error = norm_inf(*x_dot_vector_);

      if (local_error > local_error_tolerance_) {
        shadowed_dt *= 0.7;
      }

    } while (local_error > local_error_tolerance_ && shadowed_dt > min_step_size_);

    // Check if minimum step size was reached
    TEUCHOS_TEST_FOR_EXCEPTION(shadowed_dt < min_step_size_, 
      std::runtime_error,
      "ExplicitTaylorPolynomialStepper<Scalar>::takeStep(): " 
      << "Step size reached minimum step size " 
      << min_step_size_ << ".  Failing step." );

    // Increment t_
    t_ += shadowed_dt;

    numSteps_++;

    dt_ = shadowed_dt;

    return shadowed_dt;

  } else {

    // If t_ > t_final_, we're done
    if (t_ > t_final_) {
      dt_ = Teuchos::ScalarTraits<Scalar>::zero();
      return dt_;
    }

    // Compute a local truncated Taylor series solution to system
    computeTaylorSeriesSolution_();

    // If step size is too big, reduce
    if (dt > max_step_size_) {
      dt = max_step_size_;
    }

    // If step goes past t_final_, reduce
    if (t_+dt > t_final_) {
      dt = t_final_-t_;
    }

    // compute x(t_+dt)
    x_poly_->evaluate(dt, x_vector_.get());

    // Increment t_
    t_ += dt;

    numSteps_++;

    dt_ = dt;

    return dt;
  }
}


template<class Scalar>
const StepStatus<Scalar>
ExplicitTaylorPolynomialStepper<Scalar>::getStepStatus() const
{
  typedef Teuchos::ScalarTraits<Scalar> ST;
  StepStatus<Scalar> stepStatus;

  if (!haveInitialCondition_) {
    stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
  } 
  else if (numSteps_ == 0) {
    stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
    stepStatus.stepSize = dt_;
    stepStatus.order = this->getOrder();
    stepStatus.time = t_;
    stepStatus.solution = x_vector_;
    stepStatus.solutionDot = x_dot_vector_;
    if (!is_null(model_)) {
      stepStatus.residual = f_vector_;
    }
  } 
  else  {
    stepStatus.stepStatus = STEP_STATUS_CONVERGED;
    stepStatus.stepSize = dt_;
    stepStatus.order = this->getOrder();
    stepStatus.time = t_;
    stepStatus.solution = x_vector_;
    stepStatus.solutionDot = x_dot_vector_;
    stepStatus.residual = f_vector_;
  }
  return(stepStatus);
}


template<class Scalar>
void ExplicitTaylorPolynomialStepper<Scalar>::setParameterList(RCP<Teuchos::ParameterList> const& paramList)
{
  typedef Teuchos::ScalarTraits<Scalar> ST;

  TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
  paramList->validateParameters(*this->getValidParameters());
  parameterList_ = paramList;
  Teuchos::readVerboseObjectSublist(&*parameterList_,this);

  // Get initial time
  t_initial_ = parameterList_->get("Initial Time", ST::zero());

  // Get final time
  t_final_ = parameterList_->get("Final Time", ST::one());

  // Get local error tolerance
  local_error_tolerance_ = 
    parameterList_->get("Local Error Tolerance", ScalarMag(1.0e-10));

  // Get minimum step size
  min_step_size_ = parameterList_->get("Minimum Step Size", Scalar(1.0e-10));

  // Get maximum step size
  max_step_size_ = parameterList_->get("Maximum Step Size", Scalar(1.0));

  // Get degree_ of Taylor polynomial expansion
  degree_ = parameterList_->get("Taylor Polynomial Degree", Teuchos::as<unsigned int>(40));

  linc_ = Scalar(-16.0*std::log(10.0)/degree_);
  t_ = t_initial_;
}


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


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


template<class Scalar>
RCP<const Teuchos::ParameterList>
ExplicitTaylorPolynomialStepper<Scalar>::getValidParameters() const
{
  typedef ScalarTraits<Scalar> ST;
  static RCP<const ParameterList> validPL;
  if (is_null(validPL)) {
    RCP<ParameterList> pl = Teuchos::parameterList();

    pl->set<Scalar>("Initial Time", ST::zero());
    pl->set<Scalar>("Final Time", ST::one());
    pl->set<ScalarMag>("Local Error Tolerance", ScalarMag(1.0e-10));
    pl->set<Scalar>("Minimum Step Size", Scalar(1.0e-10));
    pl->set<Scalar>("Maximum Step Size", Scalar(1.0));
    pl->set<unsigned int>("Taylor Polynomial Degree", 40);

    Teuchos::setupVerboseObjectSublist(&*pl);
    validPL = pl;
  }
  return validPL;
}


template<class Scalar>
std::string ExplicitTaylorPolynomialStepper<Scalar>::description() const
{
  std::string name = "Rythmos::ExplicitTaylorPolynomialStepper";
  return name;
}


template<class Scalar>
void ExplicitTaylorPolynomialStepper<Scalar>::describe(
  Teuchos::FancyOStream &out,
  const Teuchos::EVerbosityLevel verbLevel
  ) const
{
  if (verbLevel == Teuchos::VERB_EXTREME) {
    out << description() << "::describe" << std::endl;
    out << "model_ = " << std::endl;
    out << Teuchos::describe(*model_, verbLevel) << std::endl;
    out << "x_vector_ = " << std::endl;
    out << Teuchos::describe(*x_vector_, verbLevel) << std::endl;
    out << "x_dot_vector_ = " << std::endl;
    out << Teuchos::describe(*x_dot_vector_, verbLevel) << std::endl;
    out << "f_vector_ = " << std::endl;
    out << Teuchos::describe(*f_vector_, verbLevel) << std::endl;
    out << "x_poly_ = " << std::endl;
    out << Teuchos::describe(*x_poly_, verbLevel) << std::endl;
    out << "f_poly_ = " << std::endl;
    out << Teuchos::describe(*f_poly_, verbLevel) << std::endl;
    out << "t_ = " << t_ << std::endl;
    out << "t_initial_ = " << t_initial_ << std::endl;
    out << "t_final_ = " << t_final_ << std::endl;
    out << "local_error_tolerance_ = " << local_error_tolerance_ << std::endl;
    out << "min_step_size_ = " << min_step_size_ << std::endl;
    out << "max_step_size_ = " << max_step_size_ << std::endl;
    out << "degree_ = " << degree_ << std::endl;
    out << "linc_ = " << linc_ << std::endl;
  }
}


template<class Scalar>
void ExplicitTaylorPolynomialStepper<Scalar>::addPoints(
  const Array<Scalar>& time_vec
  ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec
  ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
  )
{
  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for the ExplicitTaylorPolynomialStepper.\n");
}


template<class Scalar>
void ExplicitTaylorPolynomialStepper<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_-dt_,constOptInArg(*x_vector_old_),constOptInArg(*x_dot_vector_old_),
    t_,constOptInArg(*x_vector_),constOptInArg(*x_dot_vector_),
    time_vec,ptr(x_vec),ptr(xdot_vec),ptr(accuracy_vec),
    Ptr<InterpolatorBase<Scalar> >(null)
    );
}


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


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


template<class Scalar>
void ExplicitTaylorPolynomialStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
{
  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for the ExplicitTaylorPolynomialStepper.\n");
}


template<class Scalar>
int ExplicitTaylorPolynomialStepper<Scalar>::getOrder() const
{
  return degree_;
}


//
// Definitions of protected methods
//


template<class Scalar>
void
ExplicitTaylorPolynomialStepper<Scalar>::computeTaylorSeriesSolution_()
{
  RCP<Thyra::VectorBase<Scalar> > tmp;

  // Set degree_ of polynomials to 0
  x_poly_->setDegree(0);
  f_poly_->setDegree(0);

  // Set degree_ 0 coefficient
  x_poly_->setCoefficient(0, *x_vector_);

  for (unsigned int k=1; k<=degree_; k++) {

    // compute [f] = f([x])
    eval_model_explicit_poly(*model_, basePoint_, *x_poly_, t_, Teuchos::outArg(*f_poly_));

    x_poly_->setDegree(k);
    f_poly_->setDegree(k);
      
    // x[k] = f[k-1] / k
    tmp = x_poly_->getCoefficient(k);
    copy(*(f_poly_->getCoefficient(k-1)), tmp.ptr());
    scale(Scalar(1.0)/Scalar(k), tmp.ptr());
  }

}


template<class Scalar>
typename ExplicitTaylorPolynomialStepper<Scalar>::ScalarMag
ExplicitTaylorPolynomialStepper<Scalar>::estimateLogRadius_()
{
  ScalarMag rho = 0;
  ScalarMag tmp;
  for (unsigned int k=degree_/2; k<=degree_; k++) {
    tmp = log_norm_inf(*(x_poly_->getCoefficient(k))) / k;
    if (tmp > rho) {
      rho = tmp;
    }
  }
  return rho;
}


template<class Scalar>
RCP<const Thyra::VectorSpaceBase<Scalar> > ExplicitTaylorPolynomialStepper<Scalar>::get_x_space() const
{
  if (haveInitialCondition_) {
    return(x_vector_->space());
  } else {
    return Teuchos::null;
  }
}


} // namespace Rythmos

#endif // RYTHMOS_EXPLICIT_TAYLOR_POLYNOMIAL_STEPPER_H