Tempus_StepperIMEX_RK_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_StepperIMEX_RK_decl_hpp
#define Tempus_StepperIMEX_RK_decl_hpp
 
#include "Tempus_config.hpp"
#include "Tempus_StepperRKBase.hpp"
#include "Tempus_StepperImplicit.hpp"
#include "Tempus_WrapperModelEvaluatorPairIMEX_Basic.hpp"
 
namespace Tempus {
 
template <class Scalar>
class StepperIMEX_RK : virtual public Tempus::StepperImplicit<Scalar>,
                       virtual public Tempus::StepperRKBase<Scalar> {
 public:
  StepperIMEX_RK(std::string stepperType = "IMEX RK SSP2");
 
  StepperIMEX_RK(
      const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
      const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
      bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
      bool zeroInitialGuess,
      const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
      std::string stepperType,
      Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau,
      Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau,
      Scalar order);
 
 
 
  virtual Teuchos::RCP<const RKButcherTableau<Scalar> > getTableau() const
  {
    return getExplicitTableau();
  }
 
  virtual void setTableaus(
      std::string stepperType = "",
      Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau =
          Teuchos::null,
      Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau =
          Teuchos::null);
 
  virtual void setTableaus(Teuchos::RCP<Teuchos::ParameterList> stepperPL,
                           std::string stepperType);
 
  virtual Teuchos::RCP<const RKButcherTableau<Scalar> > getExplicitTableau()
      const
  {
    return explicitTableau_;
  }
 
  virtual void setExplicitTableau(
      Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau);
 
  virtual Teuchos::RCP<const RKButcherTableau<Scalar> > getImplicitTableau()
      const
  {
    return implicitTableau_;
  }
 
  virtual void setImplicitTableau(
      Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau);
 
  virtual void setModel(
      const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel);
 
  virtual Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > getModel() const
  {
    return this->wrapperModel_;
  }
 
  virtual void setModelPair(
      const Teuchos::RCP<WrapperModelEvaluatorPairIMEX_Basic<Scalar> >& mePair);
 
  virtual void setModelPair(
      const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& explicitModel,
      const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& implicitModel);
 
  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 order_; }
  virtual Scalar getOrderMax() const { return order_; }
 
  virtual bool isExplicit() const { return true; }
  virtual bool isImplicit() const { return true; }
  virtual bool isExplicitImplicit() const
  {
    return isExplicit() && isImplicit();
  }
  virtual bool isOneStepMethod() const { return true; }
  virtual bool isMultiStepMethod() const { return !isOneStepMethod(); }
  virtual OrderODE getOrderODE() const { return FIRST_ORDER_ODE; }
 
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar> > >& getStageF()
  {
    return stageF_;
  }
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar> > >& getStageG()
  {
    return stageG_;
  }
  Teuchos::RCP<Thyra::VectorBase<Scalar> >& getXTilde() { return xTilde_; }
 
  virtual Scalar getAlpha(const Scalar dt) const
  {
    const Teuchos::SerialDenseMatrix<int, Scalar>& A = implicitTableau_->A();
    return Scalar(1.0) / (dt * A(0, 0));  // Getting the first diagonal coeff!
  }
  virtual Scalar getBeta(const Scalar) const { return Scalar(1.0); }
 
  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;
 
  void evalImplicitModelExplicitly(
      const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& X, Scalar time,
      Scalar stepSize, Scalar stageNumber,
      const Teuchos::RCP<Thyra::VectorBase<Scalar> >& G) const;
 
  void evalExplicitModel(
      const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& X, Scalar time,
      Scalar stepSize, Scalar stageNumber,
      const Teuchos::RCP<Thyra::VectorBase<Scalar> >& F) const;
 
  void setOrder(Scalar order) { order_ = order; }
 
 protected:
  Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau_;
  Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau_;
 
  Scalar order_;
 
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar> > > stageF_;
  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar> > > stageG_;
 
  Teuchos::RCP<Thyra::VectorBase<Scalar> > xTilde_;
};
 
template <typename Scalar>
class StepperIMEX_RKTimeDerivative
  : virtual public Tempus::TimeDerivative<Scalar> {
 public:
  StepperIMEX_RKTimeDerivative(
      Scalar s, Teuchos::RCP<const Thyra::VectorBase<Scalar> > xTilde)
  {
    initialize(s, xTilde);
  }
 
  virtual ~StepperIMEX_RKTimeDerivative() {}
 
  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;
 
    // ith stage
    // s = 1/(dt*a_ii)
    // xOld = solution at beginning of time step
    // xTilde = xOld + dt*(Sum_{j=1}^{i-1} a_ij x_dot_j)
    // xDotTilde = - (s*x_i - s*xTilde)
    Thyra::V_StVpStV(xDot.ptr(), s_, *x, -s_, *xTilde_);
  }
 
  virtual void initialize(Scalar s,
                          Teuchos::RCP<const Thyra::VectorBase<Scalar> > xTilde)
  {
    s_      = s;
    xTilde_ = xTilde;
  }
 
 private:
  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xTilde_;
  Scalar s_;  // = 1/(dt*a_ii)
};
 
// ------------------------------------------------------------------------
template <class Scalar>
Teuchos::RCP<StepperIMEX_RK<Scalar> > createStepperIMEX_RK(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
    std::string stepperType, Teuchos::RCP<Teuchos::ParameterList> pl);
 
}  // namespace Tempus
#endif  // Tempus_StepperIMEX_RK_decl_hpp