Tempus_StepperDIRK_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_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 i         = this->getStageNumber();
    const int numStages = this->tableau_->numStages();
    TEUCHOS_TEST_FOR_EXCEPTION(
        (i < 0 || i > numStages - 1), std::logic_error,
        "Invalid stage number, "
            << i << ".  Valid range is 0 to " << (numStages - 1) << "\n"
            << "(-1 indicates not currently in the stage loop).  Try calling\n"
            << "getAlpha(dt, i) directly with a valid stage number.\n");

    return getAlpha(dt, i);
  }

  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