Tempus_TimeStepControl_impl.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_TimeStepControl_impl_hpp
#define Tempus_TimeStepControl_impl_hpp

#include "Teuchos_TimeMonitor.hpp"

#include "Tempus_TimeStepControlStrategyConstant.hpp"
#include "Tempus_TimeStepControlStrategyComposite.hpp"
#include "Tempus_TimeStepControlStrategyBasicVS.hpp"
#include "Tempus_TimeStepControlStrategyIntegralController.hpp"
#include "Tempus_TimeEventRange.hpp"
#include "Tempus_TimeEventRangeIndex.hpp"
#include "Tempus_TimeEventList.hpp"
#include "Tempus_TimeEventListIndex.hpp"


namespace Tempus {

template<class Scalar>
TimeStepControl<Scalar>::TimeStepControl()
  : isInitialized_      (false),
    initTime_           (0.0),
    finalTime_          (1.0e+99),
    minTimeStep_        (0.0),
    initTimeStep_       (1.0e+99),
    maxTimeStep_        (1.0e+99),
    initIndex_          (0),
    finalIndex_         (1000000),
    maxAbsError_        (1.0e-08),
    maxRelError_        (1.0e-08),
    maxFailures_        (10),
    maxConsecFailures_  (5),
    numTimeSteps_       (-1),
    printDtChanges_     (true),
    teAdjustedDt_       (false),
    dtAfterTimeEvent_   (0.0)
{
  setTimeStepControlStrategy();
  setTimeEvents();
  this->initialize();
}


template<class Scalar>
TimeStepControl<Scalar>::TimeStepControl(
  Scalar              initTime,
  Scalar              finalTime,
  Scalar              minTimeStep,
  Scalar              initTimeStep,
  Scalar              maxTimeStep,
  int                 initIndex,
  int                 finalIndex,
  Scalar              maxAbsError,
  Scalar              maxRelError,
  int                 maxFailures,
  int                 maxConsecFailures,
  int                 numTimeSteps,
  bool                printDtChanges,
  bool                outputExactly,
  std::vector<int>    outputIndices,
  std::vector<Scalar> outputTimes,
  int                 outputIndexInterval,
  Scalar              outputTimeInterval,
  Teuchos::RCP<TimeEventComposite<Scalar> > timeEvent,
  Teuchos::RCP<TimeStepControlStrategy<Scalar>> stepControlStrategy)
  : isInitialized_      (false              ),
    initTime_           (initTime           ),
    finalTime_          (finalTime          ),
    minTimeStep_        (minTimeStep        ),
    initTimeStep_       (initTimeStep       ),
    maxTimeStep_        (maxTimeStep        ),
    initIndex_          (initIndex          ),
    finalIndex_         (finalIndex         ),
    maxAbsError_        (maxAbsError        ),
    maxRelError_        (maxRelError        ),
    maxFailures_        (maxFailures        ),
    maxConsecFailures_  (maxConsecFailures  ),
    printDtChanges_     (printDtChanges     ),
    teAdjustedDt_       (false              ),
    dtAfterTimeEvent_   (0.0                )
{
  using Teuchos::rcp;

  setTimeStepControlStrategy(stepControlStrategy);
  setNumTimeSteps(numTimeSteps);

  auto tec = rcp(new TimeEventComposite<Scalar>());
  tec->setName("Time Step Control Events");
  auto tes = timeEvent->getTimeEvents();
  for(auto& e : tes) tec->add(e);

  // Add a range of times.
  auto teRange = rcp(new TimeEventRange<Scalar>(
    initTime, finalTime, outputTimeInterval, "Output Time Interval", outputExactly));
  tec->add(teRange);

  // Add a list of times.
  if (!outputTimes.empty())
    tec->add(rcp(new TimeEventList<Scalar>(
      outputTimes, "Output Time List", outputExactly)));

  // Add a range of indices.
  tec->add(rcp(new TimeEventRangeIndex<Scalar>(
    initIndex, finalIndex, outputIndexInterval, "Output Index Interval")));

  // Add a list of indices.
  if (!outputTimes.empty())
    tec->add(rcp(new TimeEventListIndex<Scalar>(
      outputIndices, "Output Index List")));

  setTimeEvents(tec);

  this->initialize();
}


template<class Scalar>
void TimeStepControl<Scalar>::initialize() const
{
  TEUCHOS_TEST_FOR_EXCEPTION(
    (getInitTime() > getFinalTime() ), std::logic_error,
    "Error - Inconsistent time range.\n"
    "    (timeMin = "<<getInitTime()<<") > (timeMax = "<<getFinalTime()<<")\n");

  TEUCHOS_TEST_FOR_EXCEPTION(
    (getMinTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
    std::logic_error,
    "Error - Negative minimum time step.  dtMin = "<<getMinTimeStep()<<")\n");
  TEUCHOS_TEST_FOR_EXCEPTION(
    (getMaxTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
    std::logic_error,
    "Error - Negative maximum time step.  dtMax = "<<getMaxTimeStep()<<")\n");
  TEUCHOS_TEST_FOR_EXCEPTION(
    (getMinTimeStep() > getMaxTimeStep() ), std::logic_error,
    "Error - Inconsistent time step range.\n"
    "  (dtMin = "<<getMinTimeStep()<<") > (dtMax = "<<getMaxTimeStep()<<")\n");
  TEUCHOS_TEST_FOR_EXCEPTION(
    (getInitTimeStep() < Teuchos::ScalarTraits<Scalar>::zero() ),
    std::logic_error,
    "Error - Negative initial time step.  dtInit = "<<getInitTimeStep()<<")\n");
  TEUCHOS_TEST_FOR_EXCEPTION(
    (getInitTimeStep() < getMinTimeStep() ||
     getInitTimeStep() > getMaxTimeStep() ),
    std::out_of_range,
    "Error - Initial time step is out of range.\n"
    << "    [dtMin, dtMax] = [" << getMinTimeStep() << ", "
                                << getMaxTimeStep() << "]\n"
    << "    dtInit = " << getInitTimeStep() << "\n");

  TEUCHOS_TEST_FOR_EXCEPTION(
    (getInitIndex() > getFinalIndex() ), std::logic_error,
    "Error - Inconsistent time index range.\n"
    "  (iStepMin = "<<getInitIndex()<<") > (iStepMax = "
    <<getFinalIndex()<<")\n");

  TEUCHOS_TEST_FOR_EXCEPTION(
    (getMaxAbsError() < Teuchos::ScalarTraits<Scalar>::zero() ),
    std::logic_error,
    "Error - Negative maximum time step.  errorMaxAbs = "
    <<getMaxAbsError()<<")\n");
  TEUCHOS_TEST_FOR_EXCEPTION(
    (getMaxRelError() < Teuchos::ScalarTraits<Scalar>::zero() ),
    std::logic_error,
    "Error - Negative maximum time step.  errorMaxRel = "
    <<getMaxRelError()<<")\n");

  TEUCHOS_TEST_FOR_EXCEPTION(
    (stepControlStrategy_ == Teuchos::null), std::logic_error,
    "Error - Strategy is unset!\n");

  stepControlStrategy_->initialize();

  TEUCHOS_TEST_FOR_EXCEPTION(
    (getStepType() != "Constant" && getStepType() != "Variable"),
    std::out_of_range,
      "Error - 'Step Type' does not equal one of these:\n"
    << "  'Constant' - Integrator will take constant time step sizes.\n"
    << "  'Variable' - Integrator will allow changes to the time step size.\n"
    << "  stepType = " << getStepType() << "\n");

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


template<class Scalar>
void TimeStepControl<Scalar>::
printDtChanges(int istep, Scalar dt_old, Scalar dt_new, std::string reason) const
{
  if (!getPrintDtChanges()) return;

  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
  out->setOutputToRootOnly(0);
  Teuchos::OSTab ostab(out,0,"printDtChanges");

  std::stringstream message;
  message << std::scientific
                   <<std::setw(6)<<std::setprecision(3)<<istep
    << " *  (dt = "<<std::setw(9)<<std::setprecision(3)<<dt_old
    <<   ", new = "<<std::setw(9)<<std::setprecision(3)<<dt_new
    << ")  " << reason << std::endl;
  *out << message.str();
}


template<class Scalar>
void TimeStepControl<Scalar>::checkInitialized()
{
  if ( !isInitialized_ ) {
    this->describe( *(this->getOStream()), Teuchos::VERB_MEDIUM);
    TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
      "Error - " << this->description() << " is not initialized!");
  }
}


template<class Scalar>
void TimeStepControl<Scalar>::setNextTimeStep(
  const Teuchos::RCP<SolutionHistory<Scalar> > & solutionHistory,
  Status & integratorStatus)
{
  using Teuchos::RCP;

  checkInitialized();

  TEMPUS_FUNC_TIME_MONITOR("Tempus::TimeStepControl::setNextTimeStep()");
  {
    RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
    const Scalar lastTime = solutionHistory->getCurrentState()->getTime();
    const int iStep = workingState->getIndex();
    Scalar dt   = workingState->getTimeStep();
    Scalar time = workingState->getTime();

    RCP<StepperState<Scalar> > stepperState = workingState->getStepperState();

    // If last time step was adjusted for time event, reinstate previous dt.
    if (getStepType() == "Variable") {
      if (teAdjustedDt_ == true) {
        printDtChanges(iStep, dt, dtAfterTimeEvent_, "Reset dt after time event.");
        dt = dtAfterTimeEvent_;
        time = lastTime + dt;
        teAdjustedDt_ = false;
        dtAfterTimeEvent_ = 0.0;
      }

      if (dt <= 0.0) {
        printDtChanges(iStep, dt, getInitTimeStep(), "Reset dt to initial dt.");
        dt = getInitTimeStep();
        time = lastTime + dt;
      }

      if (dt < getMinTimeStep()) {
        printDtChanges(iStep, dt, getMinTimeStep(), "Reset dt to minimum dt.");
        dt = getMinTimeStep();
        time = lastTime + dt;
      }
    }

    // Update dt for the step control strategy to be informed
    workingState->setTimeStep(dt);
    workingState->setTime(time);

    // Call the step control strategy (to update dt if needed)
    stepControlStrategy_->setNextTimeStep(*this, solutionHistory,
                                          integratorStatus);

    // Get the dt (It was probably changed by stepControlStrategy_.)
    dt = workingState->getTimeStep();
    time = workingState->getTime();

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


    // Check if this index is a TimeEvent and whether it is an output event.
    bool outputDueToIndex = false;
    std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > constrainingTEs;
    if (timeEvent_->isIndex(iStep, constrainingTEs)) {
      for(auto& e : constrainingTEs) {
        if (e->getName() == "Output Index Interval" ||
            e->getName() == "Output Index List"        )
        { outputDueToIndex = true; }
      }
    }

    // Check if during this time step is there a TimeEvent.
    Scalar endTime = lastTime+dt;
    // For "Variable", need to check t+dt+t_min in case we need to split
    // the time step into two steps (e.g., time step lands within dt_min).
    if (getStepType() == "Variable") endTime = lastTime+dt+getMinTimeStep();

    bool teThisStep = timeEvent_->eventInRange(lastTime, endTime, constrainingTEs);

    bool outputDueToTime = false;
    Scalar tone = endTime;
    bool landOnExactly = false;
    if (teThisStep) {
      for (auto& e : constrainingTEs) {
        if (e->getName() == "Output Time Interval" ||
            e->getName() == "Output Time List"        )
        { outputDueToTime = true; }

        if (e->getLandOnExactly() == true) {
          landOnExactly = true;
          tone = e->timeOfNextEvent(lastTime);
          break;
        }
      }
    }

    Scalar reltol = 1.0e-6;
    if (teThisStep && getStepType() == "Variable") {
      if (landOnExactly == true) {
        // Adjust time step to hit TimeEvent.
        if ( time > tone ) {
          // Next TimeEvent is not near next time.  It is more than
          // getMinTimeStep() away from it.
          printDtChanges(iStep, dt, tone - lastTime,
            "Adjusting dt to hit the next TimeEvent.");
          teAdjustedDt_ = true;
          dtAfterTimeEvent_ = dt;
          dt = tone - lastTime;
          time = lastTime + dt;
        } else if (std::fabs(time-tone) < reltol*std::fabs(time)) {
          // Next TimeEvent IS VERY near next time.  It is less than
          // reltol away from it, e.g., adjust for numerical roundoff.
          printDtChanges(iStep, dt, tone - lastTime,
            "Adjusting dt for numerical roundoff to hit the next TimeEvent.");
          teAdjustedDt_ = true;
          dtAfterTimeEvent_ = dt;
          dt = tone - lastTime;
          time = lastTime + dt;
        } else  if (std::fabs(time + getMinTimeStep() - tone)
                    < reltol*std::fabs(tone)) {
          // Next TimeEvent IS near next time.  It is less than
          // getMinTimeStep() away from it.  Take two time steps
          // to get to next TimeEvent.
          printDtChanges(iStep, dt, (tone - lastTime)/2.0,
            "The next TimeEvent is within the minimum dt of the next time. "
            "Adjusting dt to take two steps.");
          dt = (tone - lastTime)/2.0;
          time = lastTime + dt;
          // This time step is no longer an output step due the time constraint.
          outputDueToTime = false;
        }
      }
    }

    // Adjust time step to hit final time or correct for small
    // numerical differences.
    if (getStepType() == "Variable") {
      if ( (lastTime+dt > getFinalTime()) ||
           (std::fabs((lastTime+dt-getFinalTime()))
            < reltol*std::fabs(lastTime+dt)) ) {
        printDtChanges(iStep, dt, getFinalTime() - lastTime,
          "Adjusting dt to hit final time.");
        dt = getFinalTime() - lastTime;
        time = lastTime + dt;
      }
    }

    // Check for negative time step.
    TEUCHOS_TEST_FOR_EXCEPTION( dt <= Scalar(0.0), std::out_of_range,
      "Error - Time step is not positive.  dt = " << dt <<"\n");

    // Time step always needs to keep time within range.
    TEUCHOS_TEST_FOR_EXCEPTION(
      (lastTime + dt < getInitTime()), std::out_of_range,
      "Error - Time step does not move time INTO time range.\n"
      "    [timeMin, timeMax] = [" << getInitTime() << ", "
      << getFinalTime() << "]\n"
      "    T + dt = " << lastTime <<" + "<< dt <<" = " << lastTime + dt <<"\n");

    if (getStepType() == "Variable") {
      TEUCHOS_TEST_FOR_EXCEPTION(
        (lastTime + dt > getFinalTime()), std::out_of_range,
        "Error - Time step move time OUT OF time range.\n"
        "    [timeMin, timeMax] = [" << getInitTime() << ", "
        << getFinalTime() << "]\n"
        "    T + dt = " << lastTime <<" + "<< dt <<" = " << lastTime + dt <<"\n");
    }

    workingState->setTimeStep(dt);
    workingState->setTime(time);
    workingState->setOutput(outputDueToIndex || outputDueToTime);
  }
  return;
}


template<class Scalar>
bool TimeStepControl<Scalar>::timeInRange(const Scalar time) const
{
  // Get absolute tolerance 1.0e-(i+14), i.e., 14 digits of accuracy.
  const int relTol = 14;
  const int i =
    (getInitTime() == 0) ? 0 : 1 + (int)std::floor(std::log10(std::fabs(getInitTime()) ) );
  const Scalar absTolInit = std::pow(10, i-relTol);
  const int j =
    (getFinalTime() == 0) ? 0 : 1 + (int)std::floor(std::log10(std::fabs(getFinalTime()) ) );
  const Scalar absTolFinal = std::pow(10, j-relTol);

  const bool test1 = getInitTime() - absTolInit <= time;
  const bool test2 = time < getFinalTime() - absTolFinal;

  return (test1 && test2);
}


template<class Scalar>
bool TimeStepControl<Scalar>::indexInRange(const int iStep) const
{
  bool iir = (getInitIndex() <= iStep && iStep < getFinalIndex());
  return iir;
}


template<class Scalar>
std::string TimeStepControl<Scalar>::getStepType() const
{
  TEUCHOS_TEST_FOR_EXCEPTION(
    (stepControlStrategy_ == Teuchos::null), std::logic_error,
    "Error - getStepType() - Strategy is unset!\n");
  return stepControlStrategy_->getStepType();
}


template<class Scalar>
bool TimeStepControl<Scalar>::getOutputExactly() const
{
  auto event = timeEvent_->find("Output Time Interval");
  if (event != Teuchos::null) return event->getLandOnExactly();

  event = timeEvent_->find("Output Time List");
  if (event != Teuchos::null) return event->getLandOnExactly();

  return true;
}


template<class Scalar>
std::vector<int> TimeStepControl<Scalar>::getOutputIndices() const
{
  std::vector<int> indices;
  if ( timeEvent_->find("Output Index List") == Teuchos::null)
    return indices;
  auto teli = Teuchos::rcp_dynamic_cast<TimeEventListIndex<Scalar> >(
    timeEvent_->find("Output Index List"));
  return teli->getIndexList();
}


template<class Scalar>
std::vector<Scalar> TimeStepControl<Scalar>::getOutputTimes() const
{
  std::vector<Scalar> times;
  if ( timeEvent_->find("Output Time List") == Teuchos::null)
    return times;
  auto tel = Teuchos::rcp_dynamic_cast<TimeEventList<Scalar> >(
    timeEvent_->find("Output Time List"));
  return tel->getTimeList();
}


template<class Scalar>
int TimeStepControl<Scalar>::getOutputIndexInterval() const
{
  if ( timeEvent_->find("Output Index Interval") == Teuchos::null)
    return 1000000;
  auto teri = Teuchos::rcp_dynamic_cast<TimeEventRangeIndex<Scalar> >(
    timeEvent_->find("Output Index Interval"));
  return teri->getIndexStride();
}


template<class Scalar>
Scalar TimeStepControl<Scalar>::getOutputTimeInterval() const
{
  if ( timeEvent_->find("Output Time Interval") == Teuchos::null)
    return 1.0e+99;
  auto ter = Teuchos::rcp_dynamic_cast<TimeEventRange<Scalar> >(
    timeEvent_->find("Output Time Interval"));
  return ter->getTimeStride();
}


template<class Scalar>
void TimeStepControl<Scalar>::setNumTimeSteps(int numTimeSteps)
{
  TEUCHOS_TEST_FOR_EXCEPTION( (getStepType() != "Constant"), std::out_of_range,
    "Error - Can only use setNumTimeSteps() when 'Step Type' == 'Constant'.\n");

  numTimeSteps_ = numTimeSteps;
  if (numTimeSteps_ >= 0) {
    setFinalIndex(getInitIndex() + numTimeSteps_);
    Scalar initTimeStep;
    if (numTimeSteps_ == 0)
      initTimeStep = Scalar(0.0);
    else
      initTimeStep = (getFinalTime() - getInitTime())/numTimeSteps_;
    setInitTimeStep(initTimeStep);
    setMinTimeStep (initTimeStep);
    setMaxTimeStep (initTimeStep);

    Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
    out->setOutputToRootOnly(0);
    Teuchos::OSTab ostab(out,1,"setNumTimeSteps");
    *out << "Warning - setNumTimeSteps() Setting 'Number of Time Steps' = " << getNumTimeSteps()
         << "  Set the following parameters: \n"
         << "  'Final Time Index'     = " << getFinalIndex() << "\n"
         << "  'Initial Time Step'    = " << getInitTimeStep() << "\n"
         << "  'Step Type'            = " << getStepType() << std::endl;

    isInitialized_ = false;
  }
}


template<class Scalar>
void TimeStepControl<Scalar>::setOutputExactly(bool b)
{
  auto event = timeEvent_->find("Output Time Interval");
  if (event != Teuchos::null) {
    auto ter = Teuchos::rcp_dynamic_cast<TimeEventRange<Scalar> >(event);
    ter->setLandOnExactly(b);
  }
  event = timeEvent_->find("Output Time List");
  if (event != Teuchos::null) {
    auto tel = Teuchos::rcp_dynamic_cast<TimeEventList<Scalar> >(event);
    tel->setLandOnExactly(b);
  }
}


template<class Scalar>
void TimeStepControl<Scalar>::setOutputIndices(std::vector<int> outputIndices)
{
  timeEvent_->add(Teuchos::rcp(new Tempus::TimeEventListIndex<Scalar>(
    outputIndices, "Output Index List")));
  isInitialized_ = false;
}


template<class Scalar>
void TimeStepControl<Scalar>::setOutputTimes(std::vector<Scalar> outputTimes)
{
  timeEvent_->add(Teuchos::rcp(new Tempus::TimeEventList<Scalar>(
    outputTimes, "Output Time List", getOutputExactly())));
  isInitialized_ = false;
}


template<class Scalar>
void TimeStepControl<Scalar>::setOutputIndexInterval(int i)
{
  timeEvent_->add(Teuchos::rcp(new TimeEventRangeIndex<Scalar>(
    getInitIndex(), getFinalIndex(), i, "Output Index Interval")));
  isInitialized_ = false;
}


template<class Scalar>
void TimeStepControl<Scalar>::setOutputTimeInterval(Scalar t)
{
  timeEvent_->add(Teuchos::rcp(new TimeEventRange<Scalar>(
    getInitTime(), getFinalTime(), t, "Output Time Interval", getOutputExactly())));
  isInitialized_ = false;
}


template<class Scalar>
void TimeStepControl<Scalar>::setTimeEvents(
  Teuchos::RCP<TimeEventComposite<Scalar> > timeEvent)
{
  using Teuchos::rcp;

  if ( timeEvent != Teuchos::null ) {
    timeEvent_ = timeEvent;
  } else {
    auto tec = rcp(new TimeEventComposite<Scalar>());
    tec->setName("Time Step Control Events");
    tec->add(rcp(new TimeEventRangeIndex<Scalar>(
      getInitIndex(), getFinalIndex(), 1000000, "Output Index Interval")));
    tec->add(rcp(new TimeEventRange<Scalar>(
      getInitTime(), getFinalTime(), 1.0e+99, "Output Time Interval", true)));
    timeEvent_ = tec;
  }
  isInitialized_ = false;
}


template<class Scalar>
void TimeStepControl<Scalar>::setTimeStepControlStrategy(
  Teuchos::RCP<TimeStepControlStrategy<Scalar> > tscs)
{
  using Teuchos::rcp;

  if ( tscs != Teuchos::null ) {
    stepControlStrategy_ = tscs;
  } else {
    stepControlStrategy_ =
      rcp(new TimeStepControlStrategyConstant<Scalar>(getInitTimeStep()));
  }
  isInitialized_ = false;
}


template<class Scalar>
std::string TimeStepControl<Scalar>::description() const
{
  std::string name = "Tempus::TimeStepControl";
  return(name);
}


template<class Scalar>
void TimeStepControl<Scalar>::describe(
   Teuchos::FancyOStream               &out,
   const Teuchos::EVerbosityLevel      verbLevel) const
{
  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)) {
    std::vector<int> idx = getOutputIndices();
    std::ostringstream listIdx;
    if (!idx.empty()) {
      for(std::size_t i = 0; i < idx.size()-1; ++i) listIdx << idx[i] << ", ";
      listIdx << idx[idx.size()-1];
    }

    std::vector<Scalar> times = getOutputTimes();
    std::ostringstream listTimes;
    if (!times.empty()) {
      for(std::size_t i = 0; i < times.size()-1; ++i)
        listTimes << times[i] << ", ";
      listTimes << times[times.size()-1];
    }

    *l_out << "  stepType           = " << getStepType()            << std::endl
           << "  initTime           = " << getInitTime()            << std::endl
           << "  finalTime          = " << getFinalTime()           << std::endl
           << "  minTimeStep        = " << getMinTimeStep()         << std::endl
           << "  initTimeStep       = " << getInitTimeStep()        << std::endl
           << "  maxTimeStep        = " << getMaxTimeStep()         << std::endl
           << "  initIndex          = " << getInitIndex()           << std::endl
           << "  finalIndex         = " << getFinalIndex()          << std::endl
           << "  maxAbsError        = " << getMaxAbsError()         << std::endl
           << "  maxRelError        = " << getMaxRelError()         << std::endl
           << "  maxFailures        = " << getMaxFailures()         << std::endl
           << "  maxConsecFailures  = " << getMaxConsecFailures()   << std::endl
           << "  numTimeSteps       = " << getNumTimeSteps()        << std::endl
           << "  printDtChanges     = " << getPrintDtChanges()      << std::endl
           << "  outputExactly      = " << getOutputExactly()       << std::endl
           << "  outputIndices      = " << listIdx.str()            << std::endl
           << "  outputTimes        = " << listTimes.str()          << std::endl
           << "  outputIndexInterval= " << getOutputIndexInterval() << std::endl
           << "  outputTimeInterval = " << getOutputTimeInterval()  << std::endl
           << "  outputAdjustedDt   = " << teAdjustedDt_            << std::endl
           << "  dtAfterOutput      = " << dtAfterTimeEvent_        << std::endl;
           stepControlStrategy_->describe(*l_out, verbLevel);
           timeEvent_->describe(*l_out, verbLevel);
  }
  *l_out << std::string(this->description().length()+8, '-') <<std::endl;

}


template<class Scalar>
Teuchos::RCP<const Teuchos::ParameterList>
TimeStepControl<Scalar>::getValidParameters() const
{
  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList("Time Step Control");

  pl->set<double>("Initial Time"          , getInitTime()    , "Initial time");
  pl->set<double>("Final Time"            , getFinalTime()   , "Final time");
  pl->set<double>("Minimum Time Step"     , getMinTimeStep() , "Minimum time step size");
  pl->set<double>("Initial Time Step"     , getInitTimeStep(), "Initial time step size");
  pl->set<double>("Maximum Time Step"     , getMaxTimeStep() , "Maximum time step size");
  pl->set<int>   ("Initial Time Index"    , getInitIndex()   , "Initial time index");
  pl->set<int>   ("Final Time Index"      , getFinalIndex()  , "Final time index");
  pl->set<int>   ("Number of Time Steps", getNumTimeSteps(),
    "The number of constant time steps.  The actual step size gets computed\n"
    "on the fly given the size of the time domain.  Overides and resets\n"
    "  'Final Time Index'     = 'Initial Time Index' + 'Number of Time Steps'\n"
    "  'Initial Time Step'    = "
    "('Final Time' - 'Initial Time')/'Number of Time Steps'\n");
  pl->set<double>("Maximum Absolute Error", getMaxAbsError() , "Maximum absolute error");
  pl->set<double>("Maximum Relative Error", getMaxRelError() , "Maximum relative error");

  pl->set<bool>  ("Print Time Step Changes", getPrintDtChanges(),
    "Print timestep size when it changes");

  auto event = timeEvent_->find("Output Index Interval");
  if (event != Teuchos::null) {
    auto teri = Teuchos::rcp_dynamic_cast<TimeEventRangeIndex<Scalar> >(event);
    pl->set<int>   ("Output Index Interval", teri->getIndexStride(), "Output Index Interval");
  }

  event = timeEvent_->find("Output Time Interval");
  if (event != Teuchos::null) {
    auto ter = Teuchos::rcp_dynamic_cast<TimeEventRange<Scalar> >(event);
    pl->set<double>("Output Time Interval", ter->getTimeStride(), "Output time interval");
  }

  pl->set<bool>("Output Exactly On Output Times", getOutputExactly(),
    "This determines if the timestep size will be adjusted to exactly land\n"
    "on the output times for 'Variable' timestepping (default=true).\n"
    "When set to 'false' or for 'Constant' time stepping, the timestep\n"
    "following the output time will be flagged for output.\n");

  {
    event = timeEvent_->find("Output Index List");
    std::ostringstream list;
    if (event != Teuchos::null) {
      auto teli = Teuchos::rcp_dynamic_cast<TimeEventListIndex<Scalar> >(event);
      std::vector<int> idx = teli->getIndexList();
      if (!idx.empty()) {
        for(std::size_t i = 0; i < idx.size()-1; ++i) list << idx[i] << ", ";
        list << idx[idx.size()-1];
      }
    }
    pl->set<std::string>("Output Index List", list.str(),
      "Comma deliminated list of output indices");
  }

  {
    event = timeEvent_->find("Output Time List");
    std::ostringstream list;
    if (event != Teuchos::null) {
      auto teli = Teuchos::rcp_dynamic_cast<TimeEventList<Scalar> >(event);
      std::vector<Scalar> times = teli->getTimeList();
      if (!times.empty()) {
        for(std::size_t i = 0; i < times.size()-1; ++i) list << times[i] << ", ";
        list << times[times.size()-1];
      }
    }
    pl->set<std::string>("Output Time List", list.str(),
      "Comma deliminated list of output times");
  }

  { // Do not duplicate the above "Output" events in "Time Step Control Events".
    auto tecTmp = Teuchos::rcp(new TimeEventComposite<Scalar>());
    tecTmp->setTimeEvents(timeEvent_->getTimeEvents());
    tecTmp->remove("Output Index Interval");
    tecTmp->remove("Output Time Interval");
    tecTmp->remove("Output Index List");
    tecTmp->remove("Output Time List");
    if (tecTmp->getSize() > 0)
      pl->set("Time Step Control Events", *tecTmp->getValidParameters());
  }

  pl->set<int>   ("Maximum Number of Stepper Failures", getMaxFailures(),
    "Maximum number of Stepper failures");
  pl->set<int>   ("Maximum Number of Consecutive Stepper Failures", getMaxConsecFailures(),
    "Maximum number of consecutive Stepper failures");

  pl->set("Time Step Control Strategy", *stepControlStrategy_->getValidParameters());

  return pl;
}


// Nonmember constructor - ParameterList
// ------------------------------------------------------------------------
template <class Scalar>
Teuchos::RCP<TimeStepControl<Scalar> > createTimeStepControl(
  Teuchos::RCP<Teuchos::ParameterList> const& pList, bool runInitialize)
{
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::ParameterList;

  auto tsc = Teuchos::rcp(new TimeStepControl<Scalar>());
  if (pList == Teuchos::null || pList->numParams() == 0) return tsc;

  auto tscValidPL = Teuchos::rcp_const_cast<ParameterList>(tsc->getValidParameters());
  // Handle optional "Time Step Control Events" sublist.
  if (pList->isSublist("Time Step Control Events")) {
    RCP<ParameterList> tsctePL =
      Teuchos::sublist(pList, "Time Step Control Events", true);
    tscValidPL->set("Time Step Control Events", *tsctePL);
  }

  pList->validateParametersAndSetDefaults(*tscValidPL, 0);

  tsc->setInitTime(         pList->get<double>("Initial Time"));
  tsc->setFinalTime(        pList->get<double>("Final Time"));
  tsc->setMinTimeStep(      pList->get<double>("Minimum Time Step"));
  tsc->setInitTimeStep(     pList->get<double>("Initial Time Step"));
  tsc->setMaxTimeStep(      pList->get<double>("Maximum Time Step"));
  tsc->setInitIndex(        pList->get<int>   ("Initial Time Index"));
  tsc->setFinalIndex(       pList->get<int>   ("Final Time Index"));
  tsc->setMaxAbsError(      pList->get<double>("Maximum Absolute Error"));
  tsc->setMaxRelError(      pList->get<double>("Maximum Relative Error"));
  tsc->setMaxFailures(      pList->get<int>   ("Maximum Number of Stepper Failures"));
  tsc->setMaxConsecFailures(pList->get<int>   ("Maximum Number of Consecutive Stepper Failures"));
  tsc->setPrintDtChanges(   pList->get<bool>  ("Print Time Step Changes"));
  tsc->setNumTimeSteps(     pList->get<int>   ("Number of Time Steps"));


  auto tec = rcp(new TimeEventComposite<Scalar>());
  tec->setName("Time Step Control Events");

  { // Parse output indices, "Output Index List".
    std::vector<int> outputIndices;
    outputIndices.clear();
    std::string str = pList->get<std::string>("Output Index List");
    std::string delimiters(",");
    std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
    std::string::size_type pos     = str.find_first_of(delimiters, lastPos);
    while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
      std::string token = str.substr(lastPos,pos-lastPos);
      outputIndices.push_back(int(std::stoi(token)));
      if(pos==std::string::npos) break;

      lastPos = str.find_first_not_of(delimiters, pos);
      pos = str.find_first_of(delimiters, lastPos);
    }

    if (!outputIndices.empty()) {
      // order output indices
      std::sort(outputIndices.begin(),outputIndices.end());
      outputIndices.erase(std::unique(outputIndices.begin(),
                                      outputIndices.end()   ),
                                      outputIndices.end()     );

      tec->add(rcp(new Tempus::TimeEventListIndex<Scalar>(
                   outputIndices, "Output Index List")));
    }
  }

  // Parse output index internal
  tec->add(rcp(new TimeEventRangeIndex<Scalar>(
    tsc->getInitIndex(), tsc->getFinalIndex(),
    pList->get<int>("Output Index Interval"),
    "Output Index Interval")));

  auto outputExactly = pList->get<bool>("Output Exactly On Output Times",true);

  { // Parse output times, "Output Time List".
    std::vector<Scalar> outputTimes;
    outputTimes.clear();
    std::string str = pList->get<std::string>("Output Time List");
    std::string delimiters(",");
    // Skip delimiters at the beginning
    std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
    // Find the first delimiter
    std::string::size_type pos     = str.find_first_of(delimiters, lastPos);
    while ((pos != std::string::npos) || (lastPos != std::string::npos)) {
      // Found a token, add it to the vector
      std::string token = str.substr(lastPos,pos-lastPos);
      outputTimes.push_back(Scalar(std::stod(token)));
      if(pos==std::string::npos) break;

      lastPos = str.find_first_not_of(delimiters, pos); // Skip delimiters
      pos = str.find_first_of(delimiters, lastPos);     // Find next delimiter
    }

    if (!outputTimes.empty()) {
      // order output times
      std::sort(outputTimes.begin(),outputTimes.end());
      outputTimes.erase(std::unique(outputTimes.begin(),
                                    outputTimes.end()   ),
                                    outputTimes.end()     );

      tec->add(rcp(new Tempus::TimeEventList<Scalar>(
                   outputTimes, "Output Time List", outputExactly)));
    }
  }

  // Parse output time interval
  tec->add(rcp(new TimeEventRange<Scalar>(
    tsc->getInitTime(), tsc->getFinalTime(),
    pList->get<double>("Output Time Interval"),
    "Output Time Interval", outputExactly)));

  if (pList->isSublist("Time Step Control Events")) {
    RCP<ParameterList> tsctePL =
      Teuchos::sublist(pList, "Time Step Control Events", true);

    auto timeEventType = tsctePL->get<std::string>("Type", "Unknown");
    if (timeEventType == "Range") {
      tec->add(createTimeEventRange<Scalar>(tsctePL));
    } else if (timeEventType == "Range Index") {
      tec->add(createTimeEventRangeIndex<Scalar>(tsctePL));
    } else if (timeEventType == "List") {
      tec->add(createTimeEventList<Scalar>(tsctePL));
    } else if (timeEventType == "List Index") {
      tec->add(createTimeEventListIndex<Scalar>(tsctePL));
    } else if (timeEventType == "Composite") {
      auto tecTmp = createTimeEventComposite<Scalar>(tsctePL);
      auto timeEvents = tecTmp->getTimeEvents();
      for(auto& e : timeEvents) tec->add(e);
    } else {
      RCP<Teuchos::FancyOStream> out =
        Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
      out->setOutputToRootOnly(0);
      Teuchos::OSTab ostab(out,1, "createTimeStepControl()");
      *out << "Warning -- Unknown Time Event Type!\n"
           << "'Type' = '" << timeEventType << "'\n"
           << "Should call add() with this "
           << "(app-specific?) Time Event.\n" << std::endl;
    }
  }

  tsc->setTimeEvents(tec);

  if ( !pList->isParameter("Time Step Control Strategy") ) {

    tsc->setTimeStepControlStrategy();  // i.e., set default Constant timestep strategy.

  } else {

    RCP<ParameterList> tscsPL =
      Teuchos::sublist(pList, "Time Step Control Strategy", true);

    auto strategyType = tscsPL->get<std::string>("Strategy Type");
    if (strategyType == "Constant") {
      tsc->setTimeStepControlStrategy(
        createTimeStepControlStrategyConstant<Scalar>(tscsPL));
    } else if (strategyType == "Basic VS") {
      tsc->setTimeStepControlStrategy(
        createTimeStepControlStrategyBasicVS<Scalar>(tscsPL));
    } else if (strategyType == "Integral Controller") {
      tsc->setTimeStepControlStrategy(
        createTimeStepControlStrategyIntegralController<Scalar>(tscsPL));
    } else if (strategyType == "Composite") {
      tsc->setTimeStepControlStrategy(
        createTimeStepControlStrategyComposite<Scalar>(tscsPL));
    } else {
      RCP<Teuchos::FancyOStream> out =
        Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
      out->setOutputToRootOnly(0);
      Teuchos::OSTab ostab(out,1, "createTimeStepControl()");
      *out << "Warning -- Did not find a Tempus strategy to create!\n"
           << "'Strategy Type' = '" << strategyType << "'\n"
           << "Should call setTimeStepControlStrategy() with this\n"
           << "(app-specific?) strategy, and initialize().\n" << std::endl;
    }
  }

  if (runInitialize) tsc->initialize();
  return tsc;
}


} // namespace Tempus
#endif // Tempus_TimeStepControl_impl_hpp