Tempus_StepperBDF2_decl.hpp

Go to the documentation of this file.

// @HEADER
// ****************************************************************************
//                Tempus: Copyright (2017) Sandia Corporation
//
// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
// ****************************************************************************
// @HEADER

#ifndef Tempus_StepperBDF2_decl_hpp
#define Tempus_StepperBDF2_decl_hpp

#include "Tempus_config.hpp"
#include "Tempus_StepperImplicit.hpp"
#include "Tempus_WrapperModelEvaluator.hpp"
#include "Tempus_StepperBDF2AppAction.hpp"

namespace Tempus {

template<class Scalar>
class StepperBDF2 : virtual public Tempus::StepperImplicit<Scalar>
{
public:

  StepperBDF2();

  StepperBDF2(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    const Teuchos::RCP<Stepper<Scalar> >& startUpStepper,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperBDF2AppAction<Scalar> >& stepperBDF2AppAction);


  virtual void setModel(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel);

  virtual void setAppAction(
    Teuchos::RCP<StepperBDF2AppAction<Scalar> > appAction);

  virtual Teuchos::RCP<StepperBDF2AppAction<Scalar> > getAppAction() const
  { return stepperBDF2AppAction_; }

  void setStartUpStepper(std::string startupStepperType);
  void setStartUpStepper(Teuchos::RCP<Stepper<Scalar> > startupStepper);

  virtual void initialize();

  virtual void setInitialConditions (
    const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory);

  virtual void takeStep(
    const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory);

  virtual Teuchos::RCP<Tempus::StepperState<Scalar> > getDefaultStepperState();
  virtual Scalar getOrder() const {return order_;}
  virtual Scalar getOrderMin() const {return 1.0;}
  virtual Scalar getOrderMax() const {return 2.0;}

  virtual bool isExplicit()         const {return false;}
  virtual bool isImplicit()         const {return true;}
  virtual bool isExplicitImplicit() const
  {return isExplicit() && isImplicit();}
  virtual bool isOneStepMethod()   const {return false;}
  virtual bool isMultiStepMethod() const {return !isOneStepMethod();}
  virtual OrderODE getOrderODE()   const {return FIRST_ORDER_ODE;}

  virtual Scalar getAlpha(const Scalar dt) const {return getAlpha(dt,dt);}
  virtual Scalar getAlpha(const Scalar dt, const Scalar dtOld) const
  { return (Scalar(2.0)*dt + dtOld)/(dt*(dt + dtOld)); }
  virtual Scalar getBeta (const Scalar   ) const { return Scalar(1.0); }

  virtual void computeStartUp(
    const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory);

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


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

  virtual bool isValidSetup(Teuchos::FancyOStream & out) const;

private:

  Teuchos::RCP<Stepper<Scalar> >              startUpStepper_;
  Teuchos::RCP<StepperBDF2AppAction<Scalar> > stepperBDF2AppAction_;
  Scalar                                      order_ = Scalar(2.0);
};


template <typename Scalar>
class StepperBDF2TimeDerivative
  : virtual public Tempus::TimeDerivative<Scalar>
{
public:

  StepperBDF2TimeDerivative(
    Scalar dt, Scalar dtOld, Teuchos::RCP<const Thyra::VectorBase<Scalar> > xOld,
    Teuchos::RCP<const Thyra::VectorBase<Scalar> > xOldOld)
  { initialize(dt, dtOld, xOld, xOldOld); }

  virtual ~StepperBDF2TimeDerivative() {}

  virtual void compute(
    Teuchos::RCP<const Thyra::VectorBase<Scalar> > x,
    Teuchos::RCP<      Thyra::VectorBase<Scalar> > xDot,
    Teuchos::RCP<      Thyra::VectorBase<Scalar> > xDotDot = Teuchos::null)
  {
    xDotDot = Teuchos::null;
    // Calculate the BDF2 x dot vector
    const Scalar a = ((Scalar(2.0)*dt_ + dtOld_)/(dt_ + dtOld_))/dt_;
    const Scalar b = (                       dt_/(dt_ + dtOld_))/dtOld_;
    //xDot = a*(x_n - x_{n-1}) - b*(x_{n-1} - x_{n-2})
    Thyra::V_StVpStV(xDot.ptr(), a, *x, -(a+b), *xOld_);
    Thyra::Vp_StV(xDot.ptr(), b, *xOldOld_);
  }

  virtual void initialize(Scalar dt, Scalar dtOld,
    Teuchos::RCP<const Thyra::VectorBase<Scalar> > xOld,
    Teuchos::RCP<const Thyra::VectorBase<Scalar> > xOldOld)
    { dt_ = dt; dtOld_ = dtOld; xOld_ = xOld; xOldOld_ = xOldOld;}

private:

  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xOld_;
  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xOldOld_;
  Scalar                                         dt_;    // = t_n - t_{n-1}
  Scalar                                         dtOld_; // = t_{n-1} - t_{n-2}
};


// ------------------------------------------------------------------------
template<class Scalar>
Teuchos::RCP<StepperBDF2<Scalar> >
createStepperBDF2(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  Teuchos::RCP<Teuchos::ParameterList> pl);


} // namespace Tempus

#endif // Tempus_StepperBDF2_decl_hpp