Tempus_TimeEventComposite.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_TimeEventComposite_decl_hpp
#define Tempus_TimeEventComposite_decl_hpp

#include "Teuchos_Time.hpp"
#include "Teuchos_ParameterList.hpp"

#include "Tempus_config.hpp"
#include "Tempus_TimeEventBase.hpp"
#include "Tempus_TimeEventComposite.hpp"
#include "Tempus_TimeEventRange.hpp"
#include "Tempus_TimeEventRangeIndex.hpp"
#include "Tempus_TimeEventList.hpp"
#include "Tempus_TimeEventListIndex.hpp"

namespace Tempus {

template <class Scalar>
class TimeEventComposite : virtual public TimeEventBase<Scalar> {
 public:
  TimeEventComposite()
  {
    this->setType("Composite");
    this->setName("TimeEventComposite");
  }

  TimeEventComposite(std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > te,
                     std::string name = "TimeEventComposite")
  {
    this->setType("Composite");
    this->setName(name);
    this->setTimeEvents(te);
  }

  virtual ~TimeEventComposite() {}


  virtual std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > getTimeEvents()
      const
  {
    std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > te = timeEvents_;
    return te;
  }

  virtual void setTimeEvents(
      std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > te)
  {
    timeEvents_ = te;
  }

  virtual bool isTime(Scalar time) const
  {
    std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > timeEvents;
    return isTime(time, timeEvents);
  }

  virtual bool isTime(
      Scalar time,
      std::vector<Teuchos::RCP<TimeEventBase<Scalar> > >& timeEvents) const
  {
    timeEvents.clear();
    for (auto& e : timeEvents_) {
      if (e->isTime(time)) timeEvents.push_back(e);
    }
    return (!timeEvents.empty());
  }

  virtual Scalar timeToNextEvent(Scalar time) const
  {
    return timeOfNextEvent(time) - time;
  }

  virtual Scalar timeToNextEvent(
      Scalar time,
      std::vector<Teuchos::RCP<TimeEventBase<Scalar> > >& timeEvents) const
  {
    return timeOfNextEvent(time, timeEvents) - time;
  }

  virtual Scalar timeOfNextEvent(Scalar time) const
  {
    std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > timeEvents;
    return timeOfNextEvent(time, timeEvents);
  }

  virtual Scalar timeOfNextEvent(
      Scalar time,
      std::vector<Teuchos::RCP<TimeEventBase<Scalar> > >& timeEvents) const
  {
    timeEvents.clear();
    typedef std::pair<Scalar, Teuchos::RCP<TimeEventBase<Scalar> > > TEPAIR;
    std::vector<TEPAIR> timeEventPair;
    for (auto& e : timeEvents_)
      timeEventPair.push_back(std::make_pair(e->timeOfNextEvent(time), e));

    if (timeEventPair.empty()) return this->getDefaultTime();

    auto compare = [](TEPAIR a, TEPAIR b) { return a.first < b.first; };
    std::stable_sort(timeEventPair.begin(), timeEventPair.end(), compare);

    // The first one is the "time of next event".
    Scalar tone = timeEventPair.front().first;

    // Check if there are multiple events that match "time of next event".
    for (auto it = timeEventPair.begin(); it != timeEventPair.end(); ++it) {
      if ((*it).second->isTime(tone)) timeEvents.push_back((*it).second);
    }

    return tone;
  }

  virtual bool eventInRange(Scalar time1, Scalar time2) const
  {
    std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > timeEvents;
    return eventInRange(time1, time2, timeEvents);
  }

  virtual bool eventInRange(
      Scalar time1, Scalar time2,
      std::vector<Teuchos::RCP<TimeEventBase<Scalar> > >& timeEvents) const
  {
    typedef std::pair<Scalar, Teuchos::RCP<TimeEventBase<Scalar> > > TEPAIR;
    std::vector<TEPAIR> timeEventPair;
    for (auto& e : timeEvents_) {
      if (e->eventInRange(time1, time2))
        timeEventPair.push_back(std::make_pair(e->timeOfNextEvent(time1), e));
    }

    auto compare = [](TEPAIR a, TEPAIR b) { return a.first < b.first; };
    std::stable_sort(timeEventPair.begin(), timeEventPair.end(), compare);

    timeEvents.clear();
    for (auto& e : timeEventPair) timeEvents.push_back(e.second);

    return (!timeEvents.empty());
  }

  virtual bool isIndex(int index) const
  {
    std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > timeEvents;
    return isIndex(index, timeEvents);
  }

  virtual bool isIndex(
      int index,
      std::vector<Teuchos::RCP<TimeEventBase<Scalar> > >& timeEvents) const
  {
    timeEvents.clear();
    for (auto& e : timeEvents_) {
      if (e->isIndex(index)) timeEvents.push_back(e);
    }
    return (!timeEvents.empty());
  }

  virtual int indexToNextEvent(int index) const
  {
    return indexOfNextEvent(index) - index;
  }

  virtual int indexToNextEvent(
      int index,
      std::vector<Teuchos::RCP<TimeEventBase<Scalar> > >& timeEvents) const
  {
    return indexOfNextEvent(index, timeEvents) - index;
  }

  virtual int indexOfNextEvent(int index) const
  {
    std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > timeEvents;
    return indexOfNextEvent(index, timeEvents);
  }

  virtual int indexOfNextEvent(
      int index,
      std::vector<Teuchos::RCP<TimeEventBase<Scalar> > >& timeEvents) const
  {
    timeEvents.clear();
    typedef std::pair<int, Teuchos::RCP<TimeEventBase<Scalar> > > TEPAIR;
    std::vector<TEPAIR> timeEventPair;
    for (auto& e : timeEvents_)
      timeEventPair.push_back(std::make_pair(e->indexOfNextEvent(index), e));

    if (timeEventPair.size() == 0) return this->getDefaultIndex();

    auto compare = [](TEPAIR a, TEPAIR b) { return a.first < b.first; };
    std::stable_sort(timeEventPair.begin(), timeEventPair.end(), compare);

    // The first one is the "index of next event".
    int ione = timeEventPair.front().first;

    // Check if there are multiple events that match "index of next event".
    for (auto it = timeEventPair.begin(); it != timeEventPair.end(); ++it) {
      if ((*it).second->isIndex(ione)) timeEvents.push_back((*it).second);
    }

    return ione;
  }

  virtual bool eventInRangeIndex(int index1, int index2) const
  {
    std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > timeEvents;
    return eventInRangeIndex(index1, index2, timeEvents);
  }

  virtual bool eventInRangeIndex(
      int index1, int index2,
      std::vector<Teuchos::RCP<TimeEventBase<Scalar> > >& timeEvents) const
  {
    typedef std::pair<int, Teuchos::RCP<TimeEventBase<Scalar> > > TEPAIR;
    std::vector<TEPAIR> timeEventPair;
    for (auto& e : timeEvents_) {
      if (e->eventInRangeIndex(index1, index2))
        timeEventPair.push_back(std::make_pair(e->indexOfNextEvent(index1), e));
    }

    auto compare = [](TEPAIR a, TEPAIR b) { return a.first < b.first; };
    std::stable_sort(timeEventPair.begin(), timeEventPair.end(), compare);

    timeEvents.clear();
    for (auto& e : timeEventPair) timeEvents.push_back(e.second);

    return (!timeEvents.empty());
  }

  virtual Scalar getAbsTol() const
  {
    std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > timeEvents;
    return getAbsTol(timeEvents);
  }

  virtual Scalar getAbsTol(
      std::vector<Teuchos::RCP<TimeEventBase<Scalar> > >& timeEvents) const
  {
    timeEvents.clear();
    Scalar largestAbsTol = timeEvents_.front()->getAbsTol();
    timeEvents.push_back(timeEvents_.front());
    for (auto& e : timeEvents_)
      if (e->getAbsTol() > largestAbsTol) largestAbsTol = e->getAbsTol();

    for (auto& e : timeEvents_)
      if (e->getAbsTol() - largestAbsTol < largestAbsTol * 1.0e-14)
        timeEvents.push_back(e);

    return largestAbsTol;
  }

  virtual bool getLandOnExactly() const
  {
    std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > timeEvents;
    return getLandOnExactly(timeEvents);
  }

  virtual bool getLandOnExactly(
      std::vector<Teuchos::RCP<TimeEventBase<Scalar> > >& timeEvents) const
  {
    timeEvents.clear();
    for (auto& e : timeEvents_)
      if (e->getLandOnExactly()) timeEvents.push_back(e);

    return (!timeEvents.empty());
  }

  void add(Teuchos::RCP<TimeEventBase<Scalar> > timeEvent)
  {
    std::string name                        = timeEvent->getName();
    Teuchos::RCP<TimeEventBase<Scalar> > te = find(name);
    if (!te.is_null()) {
      // auto l_out = Teuchos::fancyOStream( Teuchos::rcpFromRef(std::cout) );
      // Teuchos::OSTab ostab(*l_out, 2, "TimeEventComposite::add");
      // l_out->setOutputToRootOnly(0);

      //*l_out << "TimeEventComposite::add: Replacing Time Event, "
      //       << name << "." << "\n";
      remove(name);
    }
    timeEvents_.push_back(timeEvent);
  }

  void remove(std::string name)
  {
    for (std::size_t i = 0; i < timeEvents_.size(); ++i) {
      if (timeEvents_[i]->getName() == name) {
        timeEvents_.erase(timeEvents_.begin() + i);
        break;
      }
    }
    // Did not find 'name', so did nothing.
  }

  Teuchos::RCP<TimeEventBase<Scalar> > find(std::string name)
  {
    for (std::size_t i = 0; i < timeEvents_.size(); ++i)
      if (timeEvents_[i]->getName() == name) return timeEvents_[i];

    return Teuchos::null;
  }

  void clear() { timeEvents_.clear(); }

  std::size_t getSize() const { return timeEvents_.size(); }

  std::string getTimeEventNames() const
  {
    std::stringstream tecList;
    for (std::size_t i = 0; i < timeEvents_.size(); ++i) {
      tecList << timeEvents_[i]->getName();
      if (i < timeEvents_.size() - 1) tecList << ", ";
    }
    return tecList.str();
  }

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

    *l_out << "TimeEventComposite:"
           << "\n"
           << "  name                 = " << this->getName() << "\n"
           << "  Type                 = " << this->getType() << "\n"
           << "  Number of TimeEvents = " << this->getSize() << "\n"
           << "  Time Events          = " << this->getTimeEventNames()
           << std::endl;
    *l_out << "--------------------------------------------" << std::endl;
    for (auto& e : timeEvents_) {
      (*e).describe(*l_out, verbLevel);
      *l_out << "--------------------------------------------" << std::endl;
    }
  }

  virtual Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const
  {
    Teuchos::RCP<Teuchos::ParameterList> pl =
        Teuchos::parameterList("Time Event Composite");

    pl->setName(this->getName());
    pl->set("Name", this->getName());
    pl->set("Type", this->getType());
    pl->set<std::string>("Time Events", this->getTimeEventNames());

    for (auto& s : timeEvents_) pl->set(s->getName(), *s->getValidParameters());

    return pl;
  }

 protected:
  std::vector<Teuchos::RCP<TimeEventBase<Scalar> > > timeEvents_;
};

// Nonmember constructor - ParameterList
// ------------------------------------------------------------------------

template <class Scalar>
Teuchos::RCP<TimeEventComposite<Scalar> > createTimeEventComposite(
    Teuchos::RCP<Teuchos::ParameterList> const& pList)
{
  using Teuchos::ParameterList;
  using Teuchos::RCP;

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

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

  tec->setName(pList->get("Name", "From createTimeEventComposite"));

  // string tokenizer
  std::vector<std::string> teList;
  teList.clear();
  std::string str = pList->get<std::string>("Time Events");
  std::string delimiters(",");
  const char* WhiteSpace = " \t\v\r\n";
  // 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);

    std::size_t start = token.find_first_not_of(WhiteSpace);
    std::size_t end   = token.find_last_not_of(WhiteSpace);
    token =
        (start == end ? std::string() : token.substr(start, end - start + 1));

    teList.push_back(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
  }

  // For each sublist name tokenized, add the TimeEvent.
  for (auto teName : teList) {
    RCP<ParameterList> pl =
        Teuchos::rcp(new ParameterList(pList->sublist(teName)));

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

  if (tec->getSize() == 0) {
    RCP<Teuchos::FancyOStream> out =
        Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
    out->setOutputToRootOnly(0);
    Teuchos::OSTab ostab(out, 1, "createTimeEventComposite()");
    *out << "Warning -- createTimeEventComposite() - Did not\n"
         << "           find/recognize any TimeEvents to create!\n"
         << "           If there is a app-specific TimeEvent,\n"
         << "           explicitly add it to this TimeEventComposite.\n"
         << std::endl;
  }

  return tec;
}

}  // namespace Tempus

#endif  // Tempus_TimeEventComposite_decl_hpp