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