Tempus_TimeStepControlStrategyIntegralController.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_TimeStepControlStrategy_IntegralController_hpp
#define Tempus_TimeStepControlStrategy_IntegralController_hpp

#include "Tempus_config.hpp"
#include "Tempus_TimeStepControlStrategy.hpp"
#include "Tempus_SolutionState.hpp"
#include "Tempus_SolutionHistory.hpp"
#include "Tempus_StepperState.hpp"


namespace Tempus {

template<class Scalar>
class TimeStepControlStrategyIntegralController
  : virtual public TimeStepControlStrategy<Scalar>
{
public:

  TimeStepControlStrategyIntegralController()
    : controller_("PID"), KI_(0.58), KP_(0.21), KD_(0.10),
      safetyFactor_(0.90), safetyFactorAfterReject_(0.9),
      facMax_(5.0), facMin_(0.5)
  {
    errN_ = Scalar(1.0);
    errNm1_ = Scalar(1.0);
    errNm2_ = Scalar(1.0);
    facMaxINPUT_ = facMax_;

    this->setStrategyType("Integral Controller");
    this->setStepType("Variable");
    this->setName("Integral Controller");
    this->initialize();
  }

  TimeStepControlStrategyIntegralController( std::string controller,
    Scalar KI, Scalar KP, Scalar KD,
    Scalar safetyFactor, Scalar safetyFactorAfterReject,
    Scalar facMax, Scalar facMin, std::string name = "Integral Controller")
    : controller_(controller), KI_(KI), KP_(KP), KD_(KD),
      safetyFactor_(safetyFactor),
      safetyFactorAfterReject_(safetyFactorAfterReject),
      facMax_(facMax), facMin_(facMin)
  {
    errN_ = Scalar(1.0);
    errNm1_ = Scalar(1.0);
    errNm2_ = Scalar(1.0);
    facMaxINPUT_ = facMax_;

    this->setStrategyType("Integral Controller");
    this->setStepType("Variable");
    this->setName(name);
    this->initialize();
  }


  virtual ~TimeStepControlStrategyIntegralController(){}

  virtual void setNextTimeStep(const TimeStepControl<Scalar> & tsc,
    Teuchos::RCP<SolutionHistory<Scalar> > solutionHistory,
    Status & /* integratorStatus */) override
  {
    using Teuchos::RCP;

    this->checkInitialized();

    // Take first step with initial time step provided
    if (!firstSuccessfulStep_){
      firstSuccessfulStep_ = true;
      return;
    }

    RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
    const Scalar errorRel = workingState->getErrorRel();
    Scalar beta = 1.0;

    // assumes the embedded solution is the low order solution
    int order = workingState->getOrder() - 1;
    Scalar dt = workingState->getTimeStep();

    // update errors
    errNm2_ = errNm1_;
    errNm1_ = errN_;
    errN_ = errorRel;

    Scalar k1 = Teuchos::as<Scalar>(-KI_ / order);
    Scalar k2 = Teuchos::as<Scalar>( KP_ / order);
    Scalar k3 = Teuchos::as<Scalar>(-KD_ / order);

    k1 = std::pow(errN_, k1);
    k2 = std::pow(errNm1_, k2);
    k3 = std::pow(errNm2_, k3);

    if (controller_ == "I")
       beta = safetyFactor_*k1;
    else if (controller_ == "PI")
       beta = safetyFactor_ *k1*k2;
    else // (controller_ == "PID")
       beta = safetyFactor_*k1*k2*k3;

    beta = std::max(facMin_, beta);
    beta = std::min(facMax_, beta);

    // new (optimal) suggested time step
    dt = beta * dt;

    if (workingState->getSolutionStatus() == Status::PASSED ||
        workingState->getSolutionStatus() == Status::WORKING) {
      if(lastStepRejected_){
         dt = std::min(dt, workingState->getTimeStep());
      } else {
         facMax_ = facMaxINPUT_;
      }
      lastStepRejected_ = false;
    } else {
      facMax_ = safetyFactorAfterReject_;
      lastStepRejected_ = true;
    }

    // update dt
    workingState->setTimeStep(dt);
    workingState->setTime(solutionHistory->getCurrentState()->getTime() + dt);
  }



    std::string description() const override
    { return "Tempus::TimeStepControlStrategyIntegralController"; }

    void describe(Teuchos::FancyOStream          &out,
                  const Teuchos::EVerbosityLevel verbLevel) const override
    {
      Teuchos::OSTab ostab(out,2,"describe");
      out << description() << "::describe:" << std::endl
          << "Strategy Type                      = " << this->getStrategyType() << std::endl
          << "Step Type                          = " << this->getStepType()     << std::endl
          << "Controller Type                    = " << getController()         << std::endl
          << "KI                                 = " << getKI()                 << std::endl
          << "KP                                 = " << getKP()                 << std::endl
          << "KD                                 = " << getKD()                 << std::endl
          << "errN_                              = " << errN_                   << std::endl
          << "errNm1_                            = " << errNm1_                 << std::endl
          << "errNm2_                            = " << errNm2_                 << std::endl
          << "Safety Factor                      = " << getSafetyFactor()       << std::endl
          << "Safety Factor After Step Rejection = " << getSafetyFactorAfterReject() << std::endl
          << "Maximum Safety Factor (INPUT)      = " << facMaxINPUT_            << std::endl
          << "Maximum Safety Factor              = " << getFacMax()             << std::endl
          << "Minimum Safety Factor              = " << getFacMin()             << std::endl;
    }

  virtual Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override
  {
    Teuchos::RCP<Teuchos::ParameterList> pl =
      Teuchos::parameterList("Time Step Control Strategy");

    pl->set<std::string>("Strategy Type", this->getStrategyType(), "Integral Controller");
    pl->set<std::string>("Controller Type", getController(),
                         "Proportional-Integral-Derivative");
    pl->set<Scalar>("KI" , getKI(), "Integral gain");
    pl->set<Scalar>("KP" , getKP(), "Proportional gain");
    pl->set<Scalar>("KD" , getKD(), "Derivative gain");
    pl->set<Scalar>("Safety Factor" , getSafetyFactor(), "Safety Factor");
    pl->set<Scalar>("Safety Factor After Step Rejection",
                    getSafetyFactorAfterReject(),
                    "Safety Factor Following Step Rejection");
    pl->set<Scalar>("Maximum Safety Factor" , getFacMax(), "Maximum Safety Factor");
    pl->set<Scalar>("Minimum Safety Factor" , getFacMin(), "Minimum Safety Factor");
    return pl;
  }


  virtual void initialize() const override
  {
    TEUCHOS_TEST_FOR_EXCEPTION(safetyFactor_ <= 0.0, std::out_of_range,
    "Error - Invalid value of Safety Factory= " << safetyFactor_ << "!  \n"
    << "Safety Factor must be > 0.0.\n");

    TEUCHOS_TEST_FOR_EXCEPTION(facMax_ <= 0.0, std::out_of_range,
    "Error - Invalid value of Maximum Safety Factory= " << facMax_ << "!  \n"
    << "Maximum Safety Factor must be > 0.0.\n");

    TEUCHOS_TEST_FOR_EXCEPTION(facMax_<= 0.0, std::out_of_range,
    "Error - Invalid value of Minimum Safety Factory= " << facMin_ << "!  \n"
    << "Minimum Safety Factor must be > 0.0.\n");

    TEUCHOS_TEST_FOR_EXCEPTION(((controller_ != "I") &&
                                (controller_ != "PI") &&
                                (controller_ != "PID")), std::invalid_argument,
    "Error - Invalid choice of Controller Type = " << controller_ << "!  \n"
    << "Valid Choice are ['I', 'PI', 'PID'].\n");

    this->isInitialized_ = true;   // Only place where this is set to true!
  }


  virtual std::string getController() const { return controller_; }
  virtual Scalar getKI() const { return KI_; }
  virtual Scalar getKP() const { return KP_; }
  virtual Scalar getKD() const { return KD_; }
  virtual Scalar getSafetyFactor() const { return safetyFactor_; }
  virtual Scalar getSafetyFactorAfterReject() const { return safetyFactorAfterReject_; }
  virtual Scalar getFacMax() const { return facMax_; }
  virtual Scalar getFacMin() const { return facMin_; }

  virtual void setController(std::string c) { controller_ = c; this->isInitialized_ = false; }
  virtual void setKI(Scalar k) { KI_ = k; this->isInitialized_ = false; }
  virtual void setKP(Scalar k) { KP_ = k; this->isInitialized_ = false; }
  virtual void setKD(Scalar k) { KD_ = k; this->isInitialized_ = false; }
  virtual void setSafetyFactor(Scalar f) { safetyFactor_ = f; this->isInitialized_ = false; }
  virtual void setSafetyFactorAfterReject(Scalar f) { safetyFactorAfterReject_ = f; this->isInitialized_ = false; }
  virtual void setFacMax(Scalar f) { facMax_ = f; facMaxINPUT_ = f; this->isInitialized_ = false; }
  virtual void setFacMin(Scalar f) { facMin_ = f; this->isInitialized_ = false; }

private:

  std::string controller_;          
  Scalar KI_;                       
  Scalar KP_;                       
  Scalar KD_;                       
  Scalar errN_;
  Scalar errNm1_;
  Scalar errNm2_;
  Scalar safetyFactor_;             
  Scalar safetyFactorAfterReject_;  
  Scalar facMaxINPUT_;              
  Scalar facMax_;                   
  Scalar facMin_;                   
  bool firstSuccessfulStep_ = false;
  bool lastStepRejected_ = false;

};


// Nonmember constructor.
template <class Scalar>
Teuchos::RCP<TimeStepControlStrategyIntegralController<Scalar> >
createTimeStepControlStrategyIntegralController(
  const Teuchos::RCP<Teuchos::ParameterList> pList,
  std::string name = "Integral Controller")
{
  using Teuchos::rcp;
  auto tscs = rcp(new TimeStepControlStrategyIntegralController<Scalar>());
  if (pList == Teuchos::null) return tscs;

  TEUCHOS_TEST_FOR_EXCEPTION(
    pList->get<std::string>("Strategy Type", "Integral Controller") !=
      "Integral Controller", std::logic_error,
    "Error - Strategy Type != 'Integral Controller'.  (='"
    +pList->get<std::string>("Strategy Type")+"')\n");

  pList->validateParametersAndSetDefaults(*tscs->getValidParameters());

  tscs->setController  (pList->get<std::string>("Controller Type"));
  tscs->setKI          (pList->get<Scalar>("KI"));
  tscs->setKP          (pList->get<Scalar>("KP"));
  tscs->setKD          (pList->get<Scalar>("KD"));
  tscs->setSafetyFactor(pList->get<Scalar>("Safety Factor"));
  tscs->setSafetyFactorAfterReject(pList->get<Scalar>("Safety Factor After Step Rejection"));
  tscs->setFacMax      (pList->get<Scalar>("Maximum Safety Factor"));
  tscs->setFacMin      (pList->get<Scalar>("Minimum Safety Factor"));

  tscs->setName(name);
  tscs->initialize();

  return tscs;
}


template<class Scalar>
Teuchos::RCP<Teuchos::ParameterList> getTimeStepControlStrategyIntegralControllerPL()
{
  auto t = rcp(new Tempus::TimeStepControlStrategyIntegralController<Scalar>());
  return Teuchos::rcp_const_cast<Teuchos::ParameterList> (t->getValidParameters());
}


} // namespace Tempus
#endif // Tempus_TimeStepControlStrategy_IntegralController_hpp