Rythmos_ImplicitBDFStepper.hpp

Go to the documentation of this file.

//@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_IMPLICITBDF_STEPPER_H
#define Rythmos_IMPLICITBDF_STEPPER_H

#include "Rythmos_Stepper.hpp"
#include "Teuchos_RefCountPtr.hpp"
#include "Thyra_VectorBase.hpp"
#include "Thyra_ModelEvaluator.hpp"
#include "Thyra_ModelEvaluatorHelpers.hpp"
#include "Thyra_NonlinearSolverBase.hpp"
#include "Thyra_SingleResidSSDAEModelEvaluator.hpp"
#include "Thyra_SolveSupportTypes.hpp"

namespace Rythmos {

enum BDFactionFlag { ACTION_UNSET, ACTION_LOWER, ACTION_MAINTAIN, ACTION_RAISE };
enum BDFstatusFlag { PREDICT_AGAIN, CONTINUE_ANYWAY, REP_ERR_FAIL, REP_CONV_FAIL };

template<class Scalar>
class ImplicitBDFStepper : virtual public Stepper<Scalar>
{
  public:

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

    ImplicitBDFStepper(
      const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> >  &model
      ,const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> >  &solver
      );

    ImplicitBDFStepper(
      const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> >  &model
      ,const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> >  &solver
      ,Teuchos::ParameterList &parameterList
      );

    void InitializeStepper(Teuchos::ParameterList &implicitBDFParameters);

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

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

    Scalar TakeStep(Scalar dt);
   
    Scalar TakeStep();

    Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > get_solution() const;

    Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > get_residual() const;

    std::string description() const;

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

    bool ErrWtVecSet(
      Thyra::VectorBase<Scalar> *w, 
      const Thyra::VectorBase<Scalar> &y
      );

    Scalar WRMSNorm(
      const Thyra::VectorBase<Scalar> &w
      , const Thyra::VectorBase<Scalar> &y
      ) const;
    
    bool SetPoints(
      const std::vector<Scalar>& time_list
      ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& x_list
      ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& xdot_list);
    
    bool GetPoints(
      const std::vector<Scalar>& time_list
      ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* x_list
      ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* xdot_list
      ,std::vector<ScalarMag>* accuracy_list) const;

    bool SetRange(
      const Scalar& time_lower
      ,const Scalar& time_upper
      ,const InterpolationBuffer<Scalar> & IB);

    bool GetNodes(std::vector<Scalar>* time_list) const;

    bool RemoveNodes(std::vector<Scalar>& time_list) const;

    int GetOrder() const;


  private:


    void obtainPredictor();
    void obtainResidual();
    void obtainJacobian();
    // bool interpolateSolution(Scalar time, Thyra::VectorBase<Scalar> &tmpSolVector);
    void updateHistory();
    void restoreHistory();
    void updateCoeffs();
    void initialize();
    Scalar checkReduceOrder();
    BDFstatusFlag rejectStep();
    void completeStep();

    void setDefaultMagicNumbers(Teuchos::ParameterList &magicNumberList);

    // 05/05/06 tscoffe:  I hate the underscores for private variables!
    Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > model;
    Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > solver;
    Thyra::SingleResidSSDAEModelEvaluator<Scalar>   neModel;

    Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > xn0;
    Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > xpn0;
    Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > x_dot_base;
    std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > > xHistory;
    Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > ee;
    Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > delta;
    Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > residual;
    Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > errWtVec;

    Scalar time;


    //typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag ScalarMag;
    ScalarMag relErrTol; // relative error tolerance
    ScalarMag absErrTol; // absolute error tolerance
    Scalar hh;        // Current step-size
    int currentOrder; // Current order of integration
    int oldOrder;     // previous order of integration
    int maxOrder;     // maximum order = min(5,user option maxord) - see below.
    int usedOrder;    // order used in current step (used after currentOrder is updated)
    Scalar alpha_s;    // $\alpha_s$ fixed-leading coefficient of this BDF method
    vector<Scalar> alpha;    // $\alpha_j(n)=h_n/\psi_j(n)$ coefficient used in local error test
                      // note:   $h_n$ = current step size, n = current time step
    Scalar alpha_0;     // $-\sum_{j=1}^k \alpha_j(n)$ coefficient used in local error test
    Scalar cj ;        // $-\alpha_s/h_n$ coefficient used in local error test
    Scalar ck ;        // local error coefficient gamma[0] = 0; // $\gamma_j(n)=\sum_{l=1}^{j-1}1/\psi_l(n)$ coefficient used to
    vector<Scalar> gamma;    // calculate time derivative of history array for predictor 
    vector<Scalar> beta;     // coefficients used to evaluate predictor from history array
    vector<Scalar> psi;      // $\psi_j(n) = t_n-t_{n-j}$ intermediary variable used to 
                      // compute $\beta_j(n)$
    vector<Scalar> sigma;    // $\sigma_j(n) = \frac{h_n^j(j-1)!}{\psi_1(n)*\cdots *\psi_j(n)}$
    int numberOfSteps;// number of total time integration steps taken
    int  nef;
    Scalar usedStep;
    int nscsco;
    Scalar Ek;
    Scalar Ekm1;
    Scalar Ekm2;
    Scalar Ekp1;
    Scalar Est;
    Scalar Tk;
    Scalar Tkm1;
    Scalar Tkm2;
    Scalar Tkp1;
    int newOrder;
    bool initialPhase;
    Scalar stopTime;
    bool constantStepSize;

    // Magic Numbers:
    Scalar h0_safety;
    Scalar h0_max_factor;
    Scalar h_phase0_incr;
    Scalar h_max_inv;
    Scalar Tkm1_Tk_safety;
    Scalar Tkp1_Tk_safety;
    Scalar r_factor;
    Scalar r_safety;
    Scalar r_fudge;
    Scalar r_min;
    Scalar r_max;
    Scalar r_hincr_test;
    Scalar r_hincr;
    int    max_LET_fail;
    Scalar minTimeStep;
    Scalar maxTimeStep;

    int newtonConvergenceStatus;

#ifdef Rythmos_DEBUG
    int debugLevel;
    Teuchos::RefCountPtr<Teuchos::FancyOStream> debug_out;
#endif // Rythmos_DEBUG


};

// ////////////////////////////
// Defintions

template<class Scalar>
void ImplicitBDFStepper<Scalar>::InitializeStepper(
    Teuchos::ParameterList &implicitBDFParameters
    )
{
  typedef Teuchos::ScalarTraits<Scalar> ST;
  // Initialize algorithm coefficients
  maxOrder = implicitBDFParameters.get("maxOrder",int(5)); // maximum order
  maxOrder = max(1,min(maxOrder,5)); // 1 <= maxOrder <= 5
  currentOrder=1; // Current order of integration
  oldOrder=1; // previous order of integration
  usedOrder=1;  // order used in current step (used after currentOrder is updated)
  alpha_s=Scalar(-ST::one());  // $\alpha_s$ fixed-leading coefficient of this BDF method
  alpha.reserve(maxOrder+1);  // $\alpha_j(n)=h_n/\psi_j(n)$ coefficient used in local error test
                  // note:   $h_n$ = current step size, n = current time step
  gamma.reserve(maxOrder+1);  // calculate time derivative of history array for predictor 
  beta.reserve(maxOrder+1);   // coefficients used to evaluate predictor from history array
  psi.reserve(maxOrder+1);    // $\psi_j(n) = t_n-t_{n-j}$ intermediary variable used to 
                  // compute $\beta_j(n;$
  sigma.reserve(maxOrder+1);  // $\sigma_j(n) = \frac{h_n^j(j-1)!}{\psi_1(n)*\cdots *\psi_j(n)}$
  for (int i=0 ; i<maxOrder ; ++i)
  {
    alpha.push_back(ST::zero());
    beta.push_back(ST::zero());
    gamma.push_back(ST::zero());
    psi.push_back(ST::zero());
    sigma.push_back(ST::zero());
  }
  alpha_0=ST::zero();   // $-\sum_{j=1}^k \alpha_j(n)$ coefficient used in local error test
  cj=ST::zero();      // $-\alpha_s/h_n$ coefficient used in local error test
  ck=ST::zero();      // local error coefficient gamma_[0] = 0; // $\gamma_j(n)=\sum_{l=1}^{j-1}1/\psi_l(n)$ coefficient used to
  hh=ST::zero();
  numberOfSteps=0;   // number of total time integration steps taken
  nef=0;
  usedStep=ST::zero();
  nscsco=0;
  Ek=ST::zero();
  Ekm1=ST::zero();
  Ekm2=ST::zero();
  Ekp1=ST::zero();
  Est=ST::zero();
  Tk=ST::zero();
  Tkm1=ST::zero();
  Tkm2=ST::zero();
  Tkp1=ST::zero();
  newOrder=1;
  initialPhase=true;

  relErrTol = implicitBDFParameters.get( "relErrTol", Scalar(1.0e-4) );
  absErrTol = implicitBDFParameters.get( "absErrTol", Scalar(1.0e-6) );
  constantStepSize = implicitBDFParameters.get( "constantStepSize", false );
  stopTime = implicitBDFParameters.get( "stopTime", Scalar(10.0) );


#ifdef Rythmos_DEBUG
  debugLevel = implicitBDFParameters.get( "debugLevel", int(1) );
  debug_out = Teuchos::VerboseObjectBase::getDefaultOStream();
  debug_out->precision(15);
  debug_out->setMaxLenLinePrefix(28);
  debug_out->pushLinePrefix("Rythmos::ImplicitBDFStepper");
  debug_out->setShowLinePrefix(true);
  debug_out->setTabIndentStr("    ");
#endif // Rythmos_DEBUG

  setDefaultMagicNumbers(implicitBDFParameters.sublist("magicNumbers"));

#ifdef Rythmos_DEBUG
  Teuchos::OSTab ostab(debug_out,1,"InitializeStepper");
  if (debugLevel > 1)
  {
    *debug_out << "maxOrder = " << maxOrder << endl;
    *debug_out << "currentOrder = " << currentOrder << endl;
    *debug_out << "oldOrder = " << oldOrder << endl;
    *debug_out << "usedOrder = " << usedOrder << endl;
    *debug_out << "alpha_s = " << alpha_s << endl;
    for (int i=0 ; i<maxOrder ; ++i)
    {
      *debug_out << "alpha[" << i << "] = " << alpha[i] << endl;
      *debug_out << "beta[" << i << "] = " << beta[i] << endl;
      *debug_out << "gamma[" << i << "] = " << gamma[i] << endl;
      *debug_out << "psi[" << i << "] = " << psi[i] << endl;
      *debug_out << "sigma[" << i << "] = " << sigma[i] << endl;
    }
    *debug_out << "alpha_0 = " << alpha_0 << endl;
    *debug_out << "cj = " << cj << endl;
    *debug_out << "ck = " << ck << endl;
    *debug_out << "numberOfSteps = " << numberOfSteps << endl;
    *debug_out << "nef = " << nef << endl;
    *debug_out << "usedStep = " << usedStep << endl;
    *debug_out << "nscsco = " << nscsco << endl;
    *debug_out << "Ek = " << Ek << endl;
    *debug_out << "Ekm1 = " << Ekm1 << endl;
    *debug_out << "Ekm2 = " << Ekm2 << endl;
    *debug_out << "Ekp1 = " << Ekp1 << endl;
    *debug_out << "Est = " << Est << endl;
    *debug_out << "Tk = " << Tk << endl;
    *debug_out << "Tkm1 = " << Tkm1 << endl;
    *debug_out << "Tkm2 = " << Tkm2 << endl;
    *debug_out << "Tkp1 = " << Tkp1 << endl;
    *debug_out << "newOrder = " << newOrder << endl;
    *debug_out << "initialPhase = " << initialPhase << endl;
    *debug_out << "relErrTol  = " << relErrTol  << endl;
    *debug_out << "absErrTol  = " << absErrTol  << endl;
    *debug_out << "constantStepSize  = " << constantStepSize  << endl;
    *debug_out << "stopTime  = " << stopTime  << endl;
  }
#endif // Rythmos_DEBUG
}

template<class Scalar>
ImplicitBDFStepper<Scalar>::ImplicitBDFStepper(
  const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model
  ,const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > &solver
  ,Teuchos::ParameterList &parameterList
  )
{
  InitializeStepper(parameterList);

  // Now we instantiate the model and the solver
  setModel(model);
  setSolver(solver);
}

template<class Scalar>
ImplicitBDFStepper<Scalar>::ImplicitBDFStepper(
  const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model
  ,const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > &solver
  )
{
  Teuchos::ParameterList emptyList;
  InitializeStepper(emptyList);

  // Now we instantiate the model and the solver
  setModel(model);
  setSolver(solver);
}

template<class Scalar>
void ImplicitBDFStepper<Scalar>::setDefaultMagicNumbers(
    Teuchos::ParameterList &magicNumberList)
{
  // Magic Number Defaults:
  h0_safety      = magicNumberList.get( "h0_safety",      Scalar(2.0)     );
  h0_max_factor  = magicNumberList.get( "h0_max_factor",  Scalar(0.001)   );
  h_phase0_incr  = magicNumberList.get( "h_phase0_incr",  Scalar(2.0)     );
  h_max_inv      = magicNumberList.get( "h_max_inv",      Scalar(0.0)     );
  Tkm1_Tk_safety = magicNumberList.get( "Tkm1_Tk_safety", Scalar(2.0)     );
  Tkp1_Tk_safety = magicNumberList.get( "Tkp1_Tk_safety", Scalar(0.5)     );
  r_factor       = magicNumberList.get( "r_factor",       Scalar(0.9)     );
  r_safety       = magicNumberList.get( "r_safety",       Scalar(2.0)     );
  r_fudge        = magicNumberList.get( "r_fudge",        Scalar(0.0001)  );
  r_min          = magicNumberList.get( "r_min",          Scalar(0.125)   );
  r_max          = magicNumberList.get( "r_max",          Scalar(0.9)     );
  r_hincr_test   = magicNumberList.get( "r_hincr_test",   Scalar(2.0)     );
  r_hincr        = magicNumberList.get( "r_hincr",        Scalar(2.0)     );
  max_LET_fail   = magicNumberList.get( "max_LET_fail",   int(15)         );
  minTimeStep    = magicNumberList.get( "minTimeStep",    Scalar(0.0)     );
  maxTimeStep    = magicNumberList.get( "maxTimeStep",    Scalar(10.0)    ); 

#ifdef Rythmos_DEBUG
  Teuchos::OSTab ostab(debug_out,1,"setDefaultMagicNumbers");
  if (debugLevel > 1)
  {
    *debug_out << "h0_safety = " << h0_safety << endl;
    *debug_out << "h0_max_factor = " << h0_max_factor << endl;
    *debug_out << "h_phase0_incr = " << h_phase0_incr << endl;
    *debug_out << "h_max_inv = " << h_max_inv << endl;
    *debug_out << "Tkm1_Tk_safety = " << Tkm1_Tk_safety  << endl;
    *debug_out << "Tkp1_Tk_safety = " << Tkp1_Tk_safety << endl;
    *debug_out << "r_factor = " << r_factor << endl;
    *debug_out << "r_safety = " << r_safety << endl;
    *debug_out << "r_fudge = " << r_fudge << endl;
    *debug_out << "r_min = " << r_min << endl;
    *debug_out << "r_max = " << r_max << endl;
    *debug_out << "r_hincr_test = " << r_hincr_test << endl;
    *debug_out << "r_hincr = " << r_hincr << endl;
    *debug_out << "max_LET_fail = " << max_LET_fail << endl;
    *debug_out << "minTimeStep = " << minTimeStep << endl;
    *debug_out << "maxTimeStep = " << maxTimeStep << endl;
  }
#endif // Rythmos_DEBUG
}

template<class Scalar>
void ImplicitBDFStepper<Scalar>::setModel(const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model_)
{
  typedef Teuchos::ScalarTraits<Scalar> ST;
  model = model_;
  time = ST::zero();
  xn0 = model->getNominalValues().get_x()->clone_v();
  xpn0 = model->getNominalValues().get_x_dot()->clone_v(); 
  x_dot_base = model->getNominalValues().get_x_dot()->clone_v();
  ee = model->getNominalValues().get_x()->clone_v();
  delta = model->getNominalValues().get_x()->clone_v();
  residual = Thyra::createMember(model->get_f_space());
  errWtVec = xn0->clone_v();
  ErrWtVecSet(&*errWtVec,*xn0);
  xHistory.push_back(xn0->clone_v());
  xHistory.push_back(xpn0->clone_v());
  for (int i=2 ; i<=maxOrder ; ++i) // Store maxOrder+1 vectors
  {
    xHistory.push_back(xn0->clone_v()); 
    V_S(&*xHistory[i],ST::zero());
  }
  initialize(); // Now that we have the model, we can do our initialization 
}

template<class Scalar>
void ImplicitBDFStepper<Scalar>::setSolver(const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > &solver_)
{
  solver = solver_;
}

template<class Scalar>
Scalar ImplicitBDFStepper<Scalar>::TakeStep()
{
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag ScalarMag;
#ifdef Rythmos_DEBUG
  Teuchos::OSTab ostab(debug_out,1,"TakeStep");
#endif // Rythmos_DEBUG
  BDFstatusFlag status;
  while (1)
  {
    // Set up problem coefficients (and handle first step)
    updateCoeffs();
    // Set Error weight vector
    ErrWtVecSet(&*errWtVec,*xn0);
    // compute predictor
    obtainPredictor();
    // solve nonlinear problem (as follows)
    
    //
    // Setup the nonlinear equations:
    //
    //   f_bar( x_dot_coeff * x_bar + x_dot_base, x_coeff * x_bar + x_base, t_base ) = 0
    //   x_dot_coeff = -alpha_s/dt
    //   x_dot_base = x_prime_pred + (alpha_s/dt) * x_pred
    //   x_coeff = 1
    //   x_base = 0
    //   t_base = tn+dt
    //
    Scalar coeff_x_dot = Scalar(-ST::one())*alpha_s/hh;
    V_StVpStV( &*x_dot_base, ST::one(), *xpn0, alpha_s/hh, *xn0 );
    neModel.initialize(model,coeff_x_dot,x_dot_base,ST::one(),Teuchos::null,time+hh,xn0);
    //
    // Solve the implicit nonlinear system to a tolerance of ???
    // 
    // 05/08/06 tscoffe:  I really need to get the update, not the solution from
    // the nonlinear solver.
    if(solver->getModel().get()!=&neModel)
      solver->setModel( Teuchos::rcp(&neModel,false) );
    /* // Thyra::TimeStepNewtonNonlinearSolver uses a built in solveCriteria, so you can't pass one in.
       // I believe this is the correct solveCriteria for IDA though.
    Thyra::SolveMeasureType nonlinear_solve_measure_type(Thyra::SOLVE_MEASURE_NORM_RESIDUAL,Thyra::SOLVE_MEASURE_ONE); 
    ScalarMag tolerance = relErrTol/ScalarMag(10.0); // This should be changed to match the condition in IDA
    Thyra::SolveCriteria<Scalar> nonlinearSolveCriteria(nonlinear_solve_measure_type, tolerance);
    Thyra::SolveStatus<Scalar> nonlinearSolveStatus = solver->solve( &*xn0, &nonlinearSolveCriteria, &*delta ); 
    */
    //Thyra::assign(&*xn0,ST::zero()); // 08/10/06 tscoffe:  what is this doing here?  It hoses the solve.
    Thyra::SolveStatus<Scalar> nonlinearSolveStatus = solver->solve( &*xn0, NULL, &*ee ); 
    if (nonlinearSolveStatus.solveStatus == Thyra::SOLVE_STATUS_CONVERGED) 
      newtonConvergenceStatus = 0;
    else 
      newtonConvergenceStatus = -1;

    // check error and evaluate LTE
    Scalar enorm = checkReduceOrder();
    
#ifdef Rythmos_DEBUG
    if (debugLevel > 1)
    {
      *debug_out << "xn0 = " << std::endl;
      xn0->describe(*debug_out,Teuchos::VERB_EXTREME);
      *debug_out << "ee = " << std::endl;
      ee->describe(*debug_out,Teuchos::VERB_EXTREME);
      *debug_out << "delta = " << std::endl;
      delta->describe(*debug_out,Teuchos::VERB_EXTREME);
      for (int i=0; i<max(2,maxOrder); ++i)
      {
        *debug_out << "xHistory[" << i << "] = " << std::endl;
        xHistory[i]->describe(*debug_out,Teuchos::VERB_EXTREME);
      }
      *debug_out << "ck = " << ck << endl;
      *debug_out << "enorm = " << enorm << endl;
      *debug_out << "Local Truncation Error Check: (ck*enorm) < 1:  (" << ck*enorm << ") <?= 1" << endl;
    }
#endif // Rythmos_DEBUG
    // Check LTE here:
    if ((ck*enorm) > ST::one())
      status = rejectStep(); 
    else 
      break;
    if (status == CONTINUE_ANYWAY)
      break;
    if (status == REP_ERR_FAIL)
      return(Scalar(-ST::one()));
  }

  completeStep();  
  return(usedStep);
}

template<class Scalar>
Scalar ImplicitBDFStepper<Scalar>::TakeStep(Scalar dt)
{
  constantStepSize = true;
  hh = dt;
  dt = TakeStep();
  return(dt);
}

template<class Scalar>
Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > ImplicitBDFStepper<Scalar>::get_solution() const
{
  return(xHistory[0]);
}

template<class Scalar>
Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > ImplicitBDFStepper<Scalar>::get_residual() const
{
  return(residual);
}

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

template<class Scalar>
void ImplicitBDFStepper<Scalar>::describe(
      Teuchos::FancyOStream                &out
      ,const Teuchos::EVerbosityLevel      verbLevel
      ) const
{
  if (verbLevel == Teuchos::VERB_EXTREME)
  {
    out << description() << "::describe" << std::endl;
    out << "model = " << std::endl;
    model->describe(out,verbLevel);
    out << "solver = " << std::endl;
    solver->describe(out,verbLevel);
    out << "xn0 = " << std::endl;
    xn0->describe(out,verbLevel);
    out << "xpn0 = " << std::endl;
    xpn0->describe(out,verbLevel);
    out << "x_dot_base = " << std::endl;
    x_dot_base->describe(out,verbLevel);
    out << "xHistory = " << std::endl;
    for (int i=0 ; i < max(2,maxOrder) ; ++i)
    {
      out << "xHistory[" << i << "] = " << std::endl;
      xHistory[i]->describe(out,verbLevel);
    }
    out << "ee = " << std::endl;
    ee->describe(out,verbLevel);
    out << "delta = " << std::endl;
    delta->describe(out,verbLevel);
    out << "residual = " << std::endl;
    residual->describe(out,verbLevel);
    out << "errWtVec = " << std::endl;
    errWtVec->describe(out,verbLevel);
    out << "time = " << time << std::endl;
    out << "hh = " << hh << std::endl;
    out << "currentOrder = " << currentOrder << std::endl;
    out << "neModel = " << std::endl;
    neModel.describe(out,verbLevel);
  }
}

template<class Scalar>
void ImplicitBDFStepper<Scalar>::obtainPredictor()
{
  typedef Teuchos::ScalarTraits<Scalar> ST;

#ifdef Rythmos_DEBUG
  Teuchos::OSTab ostab(debug_out,1,"obtainPredictor");
  if (debugLevel > 1)
  {
    *debug_out << "currentOrder = " << currentOrder << std::endl;
  }
#endif // Rythmos_DEBUG
  
  // prepare history array for prediction
  for (int i=nscsco;i<=currentOrder;++i)
  {
    Vt_S(&*xHistory[i],beta[i]);
  }
  
  // evaluate predictor
  V_V(&*xn0,*xHistory[0]);
  V_S(&*xpn0,ST::zero());
  for (int i=1;i<=currentOrder;++i)
  {
    Vp_V(&*xn0,*xHistory[i]);
    Vp_StV(&*xpn0,gamma[i],*xHistory[i]);
  }
#ifdef Rythmos_DEBUG
  if (debugLevel > 1)
  {
    *debug_out << "xn0 = " << std::endl;
    xn0->describe(*debug_out,Teuchos::VERB_EXTREME);
    *debug_out << "xpn0 = " << std::endl;
    xpn0->describe(*debug_out,Teuchos::VERB_EXTREME);
  }
#endif // Rythmos_DEBUG
}

/* // This function is not ready yet.
template<class Scalar>
void ImplicitBDFStepper<Scalar>::interpolateSolution(Scalar timepoint, 
      Thyra::VectorBase<Scalar>  &tmpSolVector)
{
  typedef Teuchos::ScalarTraits<Scalar> ST;

// 03/23/04 tscoffe:  Currently this code is nearly identical to the IDA code
// for interpolating to an output time.  Either we acknowledge the copyright,
// the list of conditions in the license and the disclaimer or we rewrite this
// function.  The IDA license is included after this routine.
  Scalar tfuzz;   // fuzz factor to check for valid output time
  Scalar tp;      // approximately t{n-1}
  Scalar delt;    // distance between timepoint and time
  Scalar c = Scalar(1.0); // coefficient for interpolation
  Scalar gam;     // coefficient for interpolation
  int kord;       // order of interpolation
  Scalar tn = time;
  Scalar hh = hh;
  Scalar hused = usedStep;
  int kused = usedOrder;
  Scalar uround = ST::zero();  // unit round-off (set to zero for now)

  tfuzz = 100 * uround * (tn + hh);
  tp = tn - hused - tfuzz;
  if ( (timepoint - tp)*hh < ST::zero() ) 
    return false;

  assign(&*tmpSolVector,*xHistory[0]),
  kord = kused;
  if ( (kused == 0) || (timepoint == tn) ) 
    kord = 1;

  delt = timepoint - tn;
  gam = delt/psi[0];
  for (int j=1 ; j <= kord ; ++j)
  {
    c = c*gam;
    gam = (delt + psi[j-1])/psi[j];
    Vp_StV(&*tmpSolVector,c,*xHistory[j]);
  }
  return true;
}
*/

template<class Scalar>
void ImplicitBDFStepper<Scalar>::updateHistory()
{

  // Save Newton correction for potential order increase on next step.
  if (usedOrder < maxOrder)  
  {
    assign( &*xHistory[usedOrder+1], *ee );
  }
  // Update history arrays
  Vp_V( &*xHistory[usedOrder], *ee );
  for (int j=usedOrder-1;j>=0;j--) 
  {
    Vp_V( &*xHistory[j], *xHistory[j+1] );
  }
#ifdef Rythmos_DEBUG
  Teuchos::OSTab ostab(debug_out,1,"updateHistory");
  if (debugLevel > 1)
  {
    for (int i=0;i<max(2,maxOrder);++i)
    {
      *debug_out << "xHistory[" << i << "] = " << endl;
      xHistory[i]->describe(*debug_out,Teuchos::VERB_EXTREME);
    }
  }
#endif // Rythmos_DEBUG

}

template<class Scalar>
void ImplicitBDFStepper<Scalar>::restoreHistory()
{
  typedef Teuchos::ScalarTraits<Scalar> ST;

  // undo preparation of history array for prediction
  for (int i=nscsco;i<=currentOrder;++i)
  {
    Vt_S( &*xHistory[i], ST::one()/beta[i] );
  }
  for (int i=1;i<=currentOrder;++i)
  {
    psi[i-1] = psi[i] - hh;
  }
#ifdef Rythmos_DEBUG
  Teuchos::OSTab ostab(debug_out,1,"restoreHistory");
  if (debugLevel > 1)
  {
    for (int i=0;i<maxOrder;++i)
      *debug_out << "psi[" << i << "] = " << psi[i] << endl;
    for (int i=0;i<maxOrder;++i)
    {
      *debug_out << "xHistory[" << i << "] = " << endl;
      xHistory[i]->describe(*debug_out,Teuchos::VERB_EXTREME);
    }
  }
#endif // Rythmos_DEBUG
} 

template<class Scalar>
void ImplicitBDFStepper<Scalar>::updateCoeffs()
{
  typedef Teuchos::ScalarTraits<Scalar> ST;
  // If the number of steps taken with constant order and constant stepsize is
  // more than the current order + 1 then we don't bother to update the
  // coefficients because we've reached a constant step-size formula.  When
  // this is is not true, then we update the coefficients for the variable
  // step-sizes. 
  if ((hh != usedStep) || (currentOrder != usedOrder))
    nscsco = 0;
  nscsco = min(nscsco+1,usedOrder+2);
  if (currentOrder+1 >= nscsco)
  {
    beta[0] = ST::one();
    alpha[0] = ST::one();
    Scalar temp1 = hh;
    sigma[0] = ST::one();
    gamma[0] = ST::zero();
    for (int i=1;i<=currentOrder;++i)
    {
      Scalar temp2 = psi[i-1];
      psi[i-1] = temp1;
      beta[i] = beta[i-1]*psi[i-1]/temp2;
      temp1 = temp2 + hh;
      alpha[i] = hh/temp1;
      sigma[i] = Scalar(i+1)*sigma[i-1]*alpha[i];
      gamma[i] = gamma[i-1]+alpha[i-1]/hh;
    }
    psi[currentOrder] = temp1;
    alpha_s = ST::zero();
    alpha_0 = ST::zero();
    for (int i=0;i<currentOrder;++i)
    {
      alpha_s = alpha_s - Scalar(ST::one()/(i+ST::one()));
      alpha_0 = alpha_0 - alpha[i];
    }
    cj = -alpha_s/hh;
    ck = abs(alpha[currentOrder]+alpha_s-alpha_0);
    ck = max(ck,alpha[currentOrder]);
  }
#ifdef Rythmos_DEBUG
  Teuchos::OSTab ostab(debug_out,1,"updateCoeffs");
  if (debugLevel > 1)
  {
    for (int i=0;i<=maxOrder;++i)
    {
      *debug_out << "alpha[" << i << "] = " << alpha[i] << endl;
      *debug_out << "beta[" << i << "] = " << beta[i] << endl;
      *debug_out << "sigma[" << i << "] = " << sigma[i] << endl;
      *debug_out << "gamma[" << i << "] = " << gamma[i] << endl;
      *debug_out << "psi[" << i << "] = " << psi[i] << endl;
      *debug_out << "alpha_s = " << alpha_s << endl;
      *debug_out << "alpha_0 = " << alpha_0 << endl;
      *debug_out << "cj = " << cj << endl;
      *debug_out << "ck = " << ck << endl;
    }
  }
#endif // Rythmos_DEBUG

}

template<class Scalar>
void ImplicitBDFStepper<Scalar>::initialize()
{
  typedef Teuchos::ScalarTraits<Scalar> ST;

  // Choose initial step-size
  Scalar time_to_stop = stopTime - time;
  Scalar currentTimeStep;
  if (constantStepSize)
  {
    currentTimeStep = hh;
    //currentTimeStep = 0.1 * time_to_stop;
    //currentTimeStep = min(hh, currentTimeStep);
  }
  else
  {
    // compute an initial step-size based on rate of change in the solution initially
    Scalar ypnorm = WRMSNorm(*errWtVec,*xHistory[1]);
    if (ypnorm > ST::zero())  // time-dependent DAE
    {
      currentTimeStep = min(h0_max_factor*abs(time_to_stop),sqrt(2.0)/(h0_safety*ypnorm));
    } 
    else  // non-time-dependent DAE
    {
      currentTimeStep = h0_max_factor*abs(time_to_stop);
    }
    // choose min of user specified value and our value:
    if (hh > ST::zero())
      currentTimeStep = min(hh, currentTimeStep);
    // check for maximum step-size:
    Scalar rh = abs(currentTimeStep)*h_max_inv; 
    if (rh>1.0) currentTimeStep = currentTimeStep/rh;
  }
  hh = currentTimeStep;
#ifdef Rythmos_DEBUG
  Teuchos::OSTab ostab(debug_out,1,"initialize");
  if (debugLevel > 1)
  {
    *debug_out << "hh = " << hh << endl;
  }
#endif // Rythmos_DEBUG

  // x history
  assign(&*xHistory[0],*xn0);
  V_S(&*xHistory[1],ST::zero());

  // Coefficient initialization 
  numberOfSteps = 0;    // number of total time integration steps taken
  currentOrder = 1;
  usedOrder = 1;
  psi[0] = hh;
  cj = 1/psi[0];
  nscsco = 0;
}

template<class Scalar>
Scalar ImplicitBDFStepper<Scalar>::checkReduceOrder()
{
// This routine puts its output in newOrder_

// This routine changes the following variables:
//    Ek, Tk, Est, newOrder, dsDae.delta_x, dsDae.delta_q,
//    Ekm1, Tkm1, Ekm2, Tkm2 

// This routine reads but does not change the following variables:
//    currentOrder, sigma, dsDae.newtonCorrectionPtr, dsDae.qNewtonCorrectionPtr,
//    dsDae.errWtVecPtr, dsDae.qErrWtVecPtr, dsDae.xHistory, dsDae.qHistory

  Scalar enorm = WRMSNorm(*errWtVec,*ee);
  Ek = sigma[currentOrder]*enorm;
  Tk = Scalar(currentOrder+1)*Ek;
  Est = Ek;
  newOrder = currentOrder;
#ifdef Rythmos_DEBUG
  Teuchos::OSTab ostab(debug_out,1,"checkReduceOrder");
  if (debugLevel > 1)
  {
    *debug_out << "currentOrder = " << currentOrder << std::endl;
    *debug_out << "Ek = " << Ek << std::endl;
    *debug_out << "Tk = " << Tk << std::endl;
    *debug_out << "enorm = " << enorm << std::endl;
  }
#endif // Rythmos_DEBUG
  if (currentOrder>1)
  {
    V_VpV(&*delta,*xHistory[currentOrder],*ee);
    Ekm1 = sigma[currentOrder-1]*WRMSNorm(*errWtVec,*delta);