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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// 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 "Teuchos_RCP.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Thyra_VectorBase.hpp"
#include "Thyra_ModelEvaluator.hpp"
#include "Thyra_ModelEvaluatorHelpers.hpp"
#include "RTOpPack_RTOpTHelpers.hpp"

namespace Rythmos {

  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 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model);

    Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
    getModel() const;
    
    Scalar takeStep(Scalar dt, StepSizeType flag);

    const StepStatus<Scalar> getStepStatus() const;


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

    Teuchos::RCP<Teuchos::ParameterList> getParameterList();

    Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();

    std::string description() const;

    std::ostream& describe(
      std::ostream                &out
      ,const Teuchos::EVerbosityLevel      verbLevel
      ,const std::string          leadingIndent
      ,const std::string          indentSpacer
      ) 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 getPoints(
      const Array<Scalar>& time_vec
      ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec
      ,Array<Teuchos::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 computeTaylorSeriesSolution_();

    ScalarMag estimateLogRadius_();

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

    Teuchos::RCP<Teuchos::ParameterList> parameterList_;

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

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

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

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

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

    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>
    class ROpLogNormInf : public RTOpPack::ROpScalarReductionBase<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> {
      public:

        ROpLogNormInf() : RTOpPack::RTOpT<Scalar>("ROpLogInfNorm"){}

        typename Teuchos::ScalarTraits<Scalar>::magnitudeType 
          operator() (const RTOpPack::ReductTarget& reduct_obj) const { 
            return this->getRawVal(reduct_obj); 
          }


        void reduce_reduct_objs(const RTOpPack::ReductTarget& _in_reduct_obj, 
            RTOpPack::ReductTarget* _inout_reduct_obj) const
        {
          const RTOpPack::ReductTargetScalar<Scalar>& in_reduct_obj = 
            Teuchos::dyn_cast<const RTOpPack::ReductTargetScalar<Scalar> >(_in_reduct_obj);
          RTOpPack::ReductTargetScalar<Scalar>& inout_reduct_obj = 
            Teuchos::dyn_cast<RTOpPack::ReductTargetScalar<Scalar> >(*_inout_reduct_obj);
          if(in_reduct_obj.get() > inout_reduct_obj.get()) 
            inout_reduct_obj.set(in_reduct_obj.get());
        }

        void apply_op(const int num_vecs, 
            const RTOpPack::ConstSubVectorView<Scalar> sub_vecs[],
            const int num_targ_vecs, 
            const RTOpPack::SubVectorView<Scalar> targ_sub_vecs[], 
            RTOpPack::ReductTarget *_reduct_obj) const
        {
          RTOpPack::ReductTargetScalar<Scalar>& reduct_obj = 
            Teuchos::dyn_cast<RTOpPack::ReductTargetScalar<Scalar> >(*_reduct_obj);

          RTOP_APPLY_OP_1_0(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);

          typename Teuchos::ScalarTraits<Scalar>::magnitudeType norm_inf = 
            reduct_obj.get();

          // unit stride
          if(v0_s == 1) {
            for(Teuchos_Index i=0; i<subDim; i++) {
              typename Teuchos::ScalarTraits<Scalar>::magnitudeType
                mag = Teuchos::ScalarTraits<Scalar>::magnitude(*v0_val++);
              mag = std::log(Teuchos::ScalarTraits<Scalar>::one()*1e-100 + mag);
              norm_inf = mag > norm_inf ? mag : norm_inf;
            }
          } else {
            for(Teuchos_Index i=0; i<subDim; i++, v0_val+=v0_s) {
              typename Teuchos::ScalarTraits<Scalar>::magnitudeType
                mag = Teuchos::ScalarTraits<Scalar>::magnitude(*v0_val);
              mag = std::log(Teuchos::ScalarTraits<Scalar>::one()*1e-100 + mag);
              norm_inf = mag > norm_inf ? mag : norm_inf;
            }
          }
          reduct_obj.set(norm_inf);
        }

    }; // class ROpLogNormInf

  template <typename Scalar>
  typename Teuchos::ScalarTraits<Scalar>::magnitudeType
  log_norm_inf(const Thyra::VectorBase<Scalar>& x)
  {
    ROpLogNormInf<Scalar> log_norm_inf_op;
    Teuchos::RCP<RTOpPack::ReductTarget> log_norm_inf_targ = 
      log_norm_inf_op.reduct_obj_create();
    const Thyra::VectorBase<Scalar>* vecs[] = { &x };
    Thyra::applyOp<Scalar>(log_norm_inf_op,1,vecs,0,
         (Thyra::VectorBase<Scalar>**)NULL,
         log_norm_inf_targ.get());
    
    return log_norm_inf_op(*log_norm_inf_targ);
  }

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

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

  template<class Scalar>
  void ExplicitTaylorPolynomialStepper<Scalar>::setModel(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model
    )
  {
    TEST_FOR_EXCEPT(model == Teuchos::null)
    
    model_ = model;
    x_vector_ = model_->getNominalValues().get_x()->clone_v();
    x_dot_vector_ = Thyra::createMember(model_->get_x_space());
    f_vector_ = Thyra::createMember(model_->get_f_space());
  }


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

  template<class Scalar>
  Scalar 
  ExplicitTaylorPolynomialStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
  {
    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_ = Teuchos::ScalarTraits<Scalar>::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 dt = std::exp(linc_ - rho);

      // 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_;
      }

      ScalarMag local_error;

      do {

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

        // compute f( x(t_+dt), t_+dt )
        Thyra::eval_f(*model_, *x_vector_, t_+dt, f_vector_.get());

        // compute || xdot(t_+dt) - f( x(t_+dt), t_+dt ) ||
        Thyra::Vp_StV(x_dot_vector_.get(), -Teuchos::ScalarTraits<Scalar>::one(),
          *f_vector_);
        local_error = norm_inf(*x_dot_vector_);

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

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

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

      // Increment t_
      t_ += dt;

      dt_ = dt;

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

      dt_ = dt;

      return dt;
    }
  }


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

    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(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
  {
    typedef Teuchos::ScalarTraits<Scalar> ST;

    parameterList_ = paramList;

    // 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", 40);

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

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

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

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

  template<class Scalar>
  std::ostream& ExplicitTaylorPolynomialStepper<Scalar>::describe(
        std::ostream                &out
        ,const Teuchos::EVerbosityLevel      verbLevel
        ,const std::string          leadingIndent
        ,const std::string          indentSpacer
        ) const
  {
    if (verbLevel == Teuchos::VERB_EXTREME) {
      out << description() << "::describe" << std::endl;
      out << "model_ = " << std::endl;
      out << model_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
      out << "x_vector_ = " << std::endl;
      out << x_vector_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
      out << "x_dot_vector_ = " << std::endl;
      out << x_dot_vector_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
      out << "f_vector_ = " << std::endl;
      out << f_vector_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
      out << "x_poly_ = " << std::endl;
      out << x_poly_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
      out << "f_poly_ = " << std::endl;
      out << f_poly_->describe(out,verbLevel,leadingIndent,indentSpacer) << 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;
    }
    return(out);
  }


  template<class Scalar>
  void ExplicitTaylorPolynomialStepper<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_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<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec
    ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
    ,Array<ScalarMag>* accuracy_vec) const
  {
    TEST_FOR_EXCEPTION(true,std::logic_error,"Error, getPoints is not implemented for the ExplicitTaylorPolynomialStepper.\n");
  }

  template<class Scalar>
  TimeRange<Scalar> ExplicitTaylorPolynomialStepper<Scalar>::getTimeRange() const
  {
    return invalidTimeRange<Scalar>();
  }

  template<class Scalar>
  void ExplicitTaylorPolynomialStepper<Scalar>::getNodes(Array<Scalar>* time_list) const
  {
    TEST_FOR_EXCEPTION(true,std::logic_error,"Error, getNodes is not implemented for the ExplicitTaylorPolynomialStepper.\n");
  }

  template<class Scalar>
  void ExplicitTaylorPolynomialStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
  {
    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_()
  {
    Teuchos::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])
      Thyra::eval_f_poly(*model_, *x_poly_, t_, f_poly_.get());

      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.get());
      scale(Scalar(1.0)/Scalar(k), tmp.get());
    }

  }

  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
  {
    return(x_vector_->space());
  }

} // namespace Rythmos

#endif // RYTHMOS_EXPLICIT_TAYLOR_POLYNOMIAL_STEPPER_H

---