Tempus_StepperBDF2_decl.hpp

Go to the documentation of this file.

//@HEADER
// *****************************************************************************
//          Tempus: Time Integration and Sensitivity Analysis Package
//
// Copyright 2017 NTESS and the Tempus contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
//@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