Tempus_StepperRKBase.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_StepperRKBase_hpp
#define Tempus_StepperRKBase_hpp

#include "Thyra_VectorBase.hpp"

#include "Tempus_config.hpp"
#include "Tempus_Stepper.hpp"
#include "Tempus_RKButcherTableau.hpp"
#include "Tempus_StepperRKAppAction.hpp"
#include "Tempus_StepperRKModifierDefault.hpp"
#include "Tempus_Stepper_ErrorNorm.hpp"

namespace Tempus {

template <class Scalar>
class StepperRKBase : virtual public Tempus::Stepper<Scalar> {
 public:
  virtual Teuchos::RCP<const RKButcherTableau<Scalar>> getTableau() const
  {
    return tableau_;
  }

  virtual Scalar getOrder() const { return getTableau()->order(); }
  virtual Scalar getOrderMin() const { return getTableau()->orderMin(); }
  virtual Scalar getOrderMax() const { return getTableau()->orderMax(); }
  virtual int getNumberOfStages() const { return getTableau()->numStages(); }

  virtual int getStageNumber() const { return stageNumber_; }
  virtual void setStageNumber(int s) { stageNumber_ = s; }

  virtual void setUseEmbedded(bool a)
  {
    useEmbedded_ = a;
    this->setEmbeddedMemory();
    this->isInitialized_ = false;
  }

  virtual bool getUseEmbedded() const { return useEmbedded_; }

  virtual void setErrorNorm(const Teuchos::RCP<Stepper_ErrorNorm<Scalar>>
                                &errCalculator = Teuchos::null)
  {
    if (errCalculator != Teuchos::null) {
      stepperErrorNormCalculator_ = errCalculator;
    }
    else {
      auto er                     = Teuchos::rcp(new Stepper_ErrorNorm<Scalar>());
      stepperErrorNormCalculator_ = er;
    }
  }

  virtual void setAppAction(Teuchos::RCP<StepperRKAppAction<Scalar>> appAction)
  {
    if (appAction == Teuchos::null) {
      // Create default appAction
      stepperRKAppAction_ =
          Teuchos::rcp(new StepperRKModifierDefault<Scalar>());
    }
    else {
      stepperRKAppAction_ = appAction;
    }
    this->isInitialized_ = false;
  }

  virtual Teuchos::RCP<StepperRKAppAction<Scalar>> getAppAction() const
  {
    return stepperRKAppAction_;
  }

  virtual void setStepperRKValues(Teuchos::RCP<Teuchos::ParameterList> pl)
  {
    if (pl != Teuchos::null) {
      pl->validateParametersAndSetDefaults(*this->getValidParameters());
      this->setStepperValues(pl);
      if (pl->isParameter("Use Embedded"))
        this->setUseEmbedded(pl->get<bool>("Use Embedded"));
    }
  }

  virtual Teuchos::RCP<RKButcherTableau<Scalar>> createTableau(
      Teuchos::RCP<Teuchos::ParameterList> pl)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(
        pl == Teuchos::null, std::runtime_error,
        "Error parsing general tableau.  ParameterList is null.\n");

    Teuchos::RCP<RKButcherTableau<Scalar>> tableau;

    Teuchos::RCP<Teuchos::ParameterList> tableauPL =
        sublist(pl, "Tableau", true);
    std::size_t numStages = 0;
    int order             = tableauPL->get<int>("order");
    Teuchos::SerialDenseMatrix<int, Scalar> A;
    Teuchos::SerialDenseVector<int, Scalar> b;
    Teuchos::SerialDenseVector<int, Scalar> c;
    Teuchos::SerialDenseVector<int, Scalar> bstar;

    // read in the A matrix
    {
      std::vector<std::string> A_row_tokens;
      Tempus::StringTokenizer(A_row_tokens, tableauPL->get<std::string>("A"),
                              ";", true);

      // this is the only place where numStages is set
      numStages = A_row_tokens.size();

      // allocate the matrix
      A.shape(Teuchos::as<int>(numStages), Teuchos::as<int>(numStages));

      // fill the rows
      for (std::size_t row = 0; row < numStages; row++) {
        // parse the row (tokenize on space)
        std::vector<std::string> tokens;
        Tempus::StringTokenizer(tokens, A_row_tokens[row], " ", true);

        std::vector<double> values;
        Tempus::TokensToDoubles(values, tokens);

        TEUCHOS_TEST_FOR_EXCEPTION(
            values.size() != numStages, std::runtime_error,
            "Error parsing A matrix, wrong number of stages in row "
                << row << ".\n");

        for (std::size_t col = 0; col < numStages; col++)
          A(row, col) = values[col];
      }
    }

    // size b and c vectors
    b.size(Teuchos::as<int>(numStages));
    c.size(Teuchos::as<int>(numStages));

    // read in the b vector
    {
      std::vector<std::string> tokens;
      Tempus::StringTokenizer(tokens, tableauPL->get<std::string>("b"), " ",
                              true);
      std::vector<double> values;
      Tempus::TokensToDoubles(values, tokens);

      TEUCHOS_TEST_FOR_EXCEPTION(
          values.size() != numStages, std::runtime_error,
          "Error parsing b vector, wrong number of stages.\n");

      for (std::size_t i = 0; i < numStages; i++) b(i) = values[i];
    }

    // read in the c vector
    {
      std::vector<std::string> tokens;
      Tempus::StringTokenizer(tokens, tableauPL->get<std::string>("c"), " ",
                              true);
      std::vector<double> values;
      Tempus::TokensToDoubles(values, tokens);

      TEUCHOS_TEST_FOR_EXCEPTION(
          values.size() != numStages, std::runtime_error,
          "Error parsing c vector, wrong number of stages.\n");

      for (std::size_t i = 0; i < numStages; i++) c(i) = values[i];
    }

    if (tableauPL->isParameter("bstar") &&
        tableauPL->get<std::string>("bstar") != "") {
      bstar.size(Teuchos::as<int>(numStages));
      // read in the bstar vector
      {
        std::vector<std::string> tokens;
        Tempus::StringTokenizer(tokens, tableauPL->get<std::string>("bstar"),
                                " ", true);
        std::vector<double> values;
        Tempus::TokensToDoubles(values, tokens);

        TEUCHOS_TEST_FOR_EXCEPTION(
            values.size() != numStages, std::runtime_error,
            "Error parsing bstar vector, wrong number of stages.\n"
                << "      Number of RK stages    = " << numStages
                << "\n      Number of bstar values = "
                << values.size() << "\n");

        for (std::size_t i = 0; i < numStages; i++) bstar(i) = values[i];
      }
      tableau = rcp(new RKButcherTableau<Scalar>("From ParameterList", A, b, c,
                                                 order, order, order, bstar));
    }
    else {
      tableau = rcp(new RKButcherTableau<Scalar>("From ParameterList", A, b, c,
                                                 order, order, order));
    }
    return tableau;
  }

 protected:
  virtual void setEmbeddedMemory() {}

  Teuchos::RCP<RKButcherTableau<Scalar>> tableau_;

  // For Embedded RK
  bool useEmbedded_;
  Teuchos::RCP<Thyra::VectorBase<Scalar>> ee_;
  Teuchos::RCP<Thyra::VectorBase<Scalar>> abs_u0;
  Teuchos::RCP<Thyra::VectorBase<Scalar>> abs_u;
  Teuchos::RCP<Thyra::VectorBase<Scalar>> sc;

  Teuchos::RCP<Stepper_ErrorNorm<Scalar>> stepperErrorNormCalculator_;

  int stageNumber_;
  Teuchos::RCP<StepperRKAppAction<Scalar>> stepperRKAppAction_;
};

}  // namespace Tempus

#endif  // Tempus_StepperRKBase_hpp