Tempus_StepperDIRK_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_StepperDIRK_decl_hpp
#define Tempus_StepperDIRK_decl_hpp

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


namespace Tempus {

template<class Scalar>
class StepperDIRK : virtual public Tempus::StepperImplicit<Scalar>,
                    virtual public Tempus::StepperRKBase<Scalar>
{
public:


    virtual void initialize() override;

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

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

    virtual void setResetInitialGuess(bool reset_guess)
      { resetGuess_ = reset_guess; }
    virtual bool getResetInitialGuess() const
      { return resetGuess_; }

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

    virtual Teuchos::RCP<Tempus::StepperState<Scalar> >getDefaultStepperState() override;

    virtual bool isExplicit() const override
    {
      const int numStages = this->tableau_->numStages();
      Teuchos::SerialDenseMatrix<int,Scalar> A = this->tableau_->A();
      bool isExplicit = false;
      for (int i=0; i<numStages; ++i) if (A(i,i) == 0.0) isExplicit = true;
      return isExplicit && this->tableau_->isDIRK();
    }
    virtual bool isImplicit()         const override {return true;}
    virtual bool isExplicitImplicit() const override
      {return isExplicit() && isImplicit();}
    virtual bool isOneStepMethod()   const override {return true;}
    virtual bool isMultiStepMethod() const override {return !isOneStepMethod();}

    virtual OrderODE getOrderODE()   const override {return FIRST_ORDER_ODE;}

    virtual std::string getDescription() const = 0;

  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar> > >& getStageXDot() {return stageXDot_;}
  Teuchos::RCP<Thyra::VectorBase<Scalar> >& getXTilde() {return xTilde_;}

  virtual Scalar getAlpha(const Scalar dt) const override
  {
    const int numStages = this->tableau_->numStages();
    const Teuchos::SerialDenseMatrix<int,Scalar> & A=this->tableau_->A();
    Scalar aii = A(0,0);
    for (int i=0; i<numStages; ++i) {
      if (A(i,i) != 0.0) aii = A(i,i);
      break;
    }
    return (aii == 0.0) ? std::numeric_limits<Scalar>::infinity() : Scalar(1.0)/(dt*aii);
  }
  virtual Scalar getBeta (const Scalar   ) const override { return Scalar(1.0); }

  virtual Scalar getAlpha(const Scalar dt, int i) const
  {
    const Teuchos::SerialDenseMatrix<int,Scalar> & A=this->tableau_->A();
    return (A(i,i) == 0.0) ? std::numeric_limits<Scalar>::infinity() : Scalar(1.0)/(dt*A(i,i));
  }

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

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


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

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

  virtual void setStepperDIRKValues(Teuchos::RCP<Teuchos::ParameterList> pl)
  {
    if (pl != Teuchos::null) {
      // Can not validate because of optional Parameters, e.g., 'Solver Name'.
      //pl->validateParametersAndSetDefaults(*this->getValidParameters());
      this->setStepperValues(pl);
      if ( pl->isParameter("Use Embedded") )
        this->setUseEmbedded(pl->get<bool>("Use Embedded"));
      if ( pl->isParameter("Zero Initial Guess") )
        this->setZeroInitialGuess(pl->get<bool>("Zero Initial Guess"));
      if (pl->isParameter("Solver Name") )
        this->setStepperSolverValues(pl);
    }
  }


protected:

  virtual void setupDefault();

  virtual void setup(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& wrapperModel,
    const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
    bool useFSAL,
    std::string ICConsistency,
    bool ICConsistencyCheck,
    bool useEmbedded,
    bool zeroInitialGuess,
    const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction);

  virtual void setupTableau() = 0;

  virtual void setEmbeddedMemory() override;


  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar> > > stageXDot_;
  Teuchos::RCP<Thyra::VectorBase<Scalar> >               xTilde_;

  bool resetGuess_ = true;
};


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

  StepperDIRKTimeDerivative(
    Scalar s, Teuchos::RCP<const Thyra::VectorBase<Scalar> > xTilde)
  { initialize(s, xTilde); }

  virtual ~StepperDIRKTimeDerivative() {}

  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;
    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)
};


} // namespace Tempus

#endif // Tempus_StepperDIRK_decl_hpp