Tempus_TimeStepControlStrategyBasicVS.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_BasicVS_hpp
#define Tempus_TimeStepControlStrategy_BasicVS_hpp

#include "Thyra_VectorStdOps.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 TimeStepControlStrategyBasicVS
  : virtual public TimeStepControlStrategy<Scalar>
{
public:

  TimeStepControlStrategyBasicVS()
    : rho_(1.75), sigma_(0.5), minEta_(0.0), maxEta_(1.0e+16)
  {
    this->setStrategyType("Basic VS");
    this->setStepType("Variable");
    this->setName("Basic VS");
    this->initialize();
  }

  TimeStepControlStrategyBasicVS(Scalar rho, Scalar sigma, Scalar minEta,
    Scalar maxEta, std::string name = "Basic VS")
    : rho_(rho), sigma_(sigma), minEta_(minEta), maxEta_(maxEta)
  {
    this->setStrategyType("Basic VS");
    this->setStepType("Variable");
    this->setName(name);
    this->initialize();
  }

  virtual ~TimeStepControlStrategyBasicVS(){}

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

    this->checkInitialized();

    RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
    const Scalar errorAbs = workingState->getErrorAbs();
    const Scalar errorRel = workingState->getErrorRel();
    const int iStep = workingState->getIndex();
    Scalar dt       = workingState->getTimeStep();

    Scalar rho   = getAmplFactor();
    Scalar sigma = getReductFactor();
    Scalar eta   = solutionHistory->getCurrentState()->getDxNormL2Rel();
    if (iStep == 1) eta = getMinEta();  // For first step use initial dt.

    // General rule: only increase/decrease dt once for any given reason.
    if (workingState->getSolutionStatus() == Status::FAILED) {
      tsc.printDtChanges(iStep, dt, dt*sigma,
                          "Stepper failure - Decreasing dt.");
      dt *= sigma;
    }
    else { //Stepper passed
      if (eta < getMinEta()) { // increase dt
        tsc.printDtChanges(iStep, dt, dt*rho,
          "Change too small ("
          + std::to_string(eta) + " < " + std::to_string(getMinEta())
          + ").  Increasing dt.");
        dt *= rho;
      }
      else if (eta > getMaxEta()) { // reduce dt
        tsc.printDtChanges(iStep, dt, dt*sigma,
          "Change too large ("
          + std::to_string(eta) + " > " + std::to_string(getMaxEta())
          + ").  Decreasing dt.");
        dt *= sigma;
      }
      else if (errorAbs > tsc.getMaxAbsError()) { // reduce dt
        tsc.printDtChanges(iStep, dt, dt*sigma,
          "Absolute error is too large ("
          + std::to_string(errorAbs)+" > "+std::to_string(tsc.getMaxAbsError())
          + ").  Decreasing dt.");
        dt *= sigma;
      }
      else if (errorRel > tsc.getMaxRelError()) { // reduce dt
        tsc.printDtChanges(iStep, dt, dt*sigma,
          "Relative error is too large ("
          + std::to_string(errorRel)+" > "+std::to_string(tsc.getMaxRelError())
          + ").  Decreasing dt.");
        dt *= sigma;
      }
    }

    if (dt < tsc.getMinTimeStep()) { // decreased below minimum dt
      tsc.printDtChanges(iStep, dt, tsc.getMinTimeStep(),
        "dt is too small.  Resetting to minimum dt.");
      dt = tsc.getMinTimeStep();
    }
    if (dt > tsc.getMaxTimeStep()) { // increased above maximum dt
      tsc.printDtChanges(iStep, dt, tsc.getMaxTimeStep(),
        "dt is too large.  Resetting to maximum dt.");
      dt = tsc.getMaxTimeStep();
    }

    workingState->setTimeStep(dt);
    workingState->setTime(solutionHistory->getCurrentState()->getTime() + dt);
    workingState->setComputeNorms(true);
  }


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

    void describe(Teuchos::FancyOStream          &out,
                  const Teuchos::EVerbosityLevel verbLevel) const override
    {
      auto l_out = Teuchos::fancyOStream( out.getOStream() );
      Teuchos::OSTab ostab(*l_out, 2, this->description());
      l_out->setOutputToRootOnly(0);

      *l_out << "\n--- " << this->description() << " ---" << std::endl;

      if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_MEDIUM)) {
        *l_out << "  StrategyType                      = " << this->getStrategyType()<< std::endl
               << "  Step Type                         = " << this->getStepType() << std::endl
               << "  Amplification Factor              = " << getAmplFactor()   << std::endl
               << "  Reduction Factor                  = " << getReductFactor() << std::endl
               << "  Minimum Value Monitoring Function = " << getMinEta()       << std::endl
               << "  Maximum Value Monitoring Function = " << getMaxEta()       << std::endl;
        *l_out << std::string(this->description().length()+8, '-') <<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(), "Basic VS");
     pl->set<double>("Amplification Factor", getAmplFactor(), "Amplification factor");
     pl->set<double>("Reduction Factor"    , getReductFactor() , "Reduction factor");
     pl->set<double>("Minimum Value Monitoring Function", getMinEta(), "Min value eta");
     pl->set<double>("Maximum Value Monitoring Function", getMaxEta(), "Max value eta");
     return pl;
  }


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

    TEUCHOS_TEST_FOR_EXCEPTION(getReductFactor() >= 1.0, std::out_of_range,
      "Error - Invalid value of Reduction Factor = " << getReductFactor()
      << "!  \n" << "Reduction Factor must be < 1.0.\n");

    TEUCHOS_TEST_FOR_EXCEPTION(getMinEta() > getMaxEta(), std::out_of_range,
      "Error - Invalid values of 'Minimum Value Monitoring Function' = "
      << getMinEta() << "\n and 'Maximum Value Monitoring Function' = "
      << getMaxEta() <<"! \n Mininum Value cannot be > Maximum Value! \n");

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


  virtual Scalar getAmplFactor() const { return rho_; }
  virtual Scalar getReductFactor() const { return sigma_;}
  virtual Scalar getMinEta() const { return minEta_; }
  virtual Scalar getMaxEta() const { return maxEta_; }

  virtual void setAmplFactor(Scalar rho) { rho_ = rho; this->isInitialized_ = false; }
  virtual void setReductFactor(Scalar sigma) { sigma_ = sigma; this->isInitialized_ = false; }
  virtual void setMinEta(Scalar minEta) { minEta_ = minEta; this->isInitialized_ = false; }
  virtual void setMaxEta(Scalar maxEta) { maxEta_ = maxEta; this->isInitialized_ = false; }


private:

  Scalar rho_;      
  Scalar sigma_;    
  Scalar minEta_;   
  Scalar maxEta_;   

};


template <class Scalar>
Teuchos::RCP<TimeStepControlStrategyBasicVS<Scalar> >
createTimeStepControlStrategyBasicVS(
  const Teuchos::RCP<Teuchos::ParameterList> & pList,
  std::string name = "Basic VS")
{
  auto tscs = Teuchos::rcp(new TimeStepControlStrategyBasicVS<Scalar>());
  if (pList == Teuchos::null || pList->numParams() == 0) return tscs;

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

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

  tscs->setAmplFactor(  pList->get<double>("Amplification Factor"));
  tscs->setReductFactor(pList->get<double>("Reduction Factor"));
  tscs->setMinEta(pList->get<double>("Minimum Value Monitoring Function"));
  tscs->setMaxEta(pList->get<double>("Maximum Value Monitoring Function"));

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

  return tscs;
}


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


} // namespace Tempus
#endif // Tempus_TimeStepControlStrategy_BasicVS_hpp