AdaptInfo.hpp

#pragma once

// std c++ headers
#include <algorithm>
#include <cmath>
#include <limits>
#include <map>
#include <string>
#include <utility>

// AMDiS includes
#include <amdis/Output.hpp>
#include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/Math.hpp>
#include <amdis/utility/TreePath.hpp>

namespace AMDiS
{

  /**
   * \ingroup Adaption
   *
   * \brief
   * Holds adapt parameters and infos about the problem.
   */
  class AdaptInfo
  {
  public:
    using Key = std::string;

  protected:
    /** \brief
     * Stores adapt infos for a scalar problem or for one component of a
     * vector valued problem.
     */
    class ScalContent
    {
    public:
      /// Constructor.
      explicit ScalContent(std::string const& prefix);

      /// Sum of all error estimates
      double est_sum = 0.0;

      /// Sum of all time error estimates
      double est_t_sum = 0.0;

      /// maximal local error estimate.
      double est_max = 0.0;

      /// Maximum of all time error estimates
      double est_t_max = 0.0;

      /// factors to combine max and integral time estimate
      double fac_max = 0.0, fac_sum = 1.0;

      /// Tolerance for the (absolute or relative) error
      double spaceTolerance = 0.0;

      /// Time tolerance.
      double timeTolerance = 0.0;

      /// Relative time tolerance
      double timeRelativeTolerance = 0.0;

      /// Lower bound for the time error.
      double timeErrLow = 0.0;

      /// true if coarsening is allowed, false otherwise.
      int coarsenAllowed = 0;

      /// true if refinement is allowed, false otherwise.
      int refinementAllowed = 1;
    };

  public:
    /// Constructor.
    explicit AdaptInfo(std::string const& name);

    /// Destructor.
    virtual ~AdaptInfo() = default;

    /// Resets all variables to zero (or something equivalent)
    void reset();

    /// Returns whether space tolerance is reached.
    virtual bool spaceToleranceReached() const
    {
      for (auto const& scalContent : scalContents)      {
        if (!(scalContent.second.est_sum < scalContent.second.spaceTolerance))
          return false;
      }

      return true;
    }

    /// Returns whether space tolerance of component associated with key is reached.
    virtual bool spaceToleranceReached(Key key) const
    {
      if (!(getScalContent(key).est_sum < getScalContent(key).spaceTolerance))
        return false;
      else
        return true;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    bool spaceToleranceReached(const TP& tp) const
    {
      return spaceToleranceReached(to_string(tp));
    }

    /// Returns whether time tolerance is reached.
    virtual bool timeToleranceReached() const
    {
      for (auto const& scalContent : scalContents)
        if (!(getTimeEstCombined(scalContent.first) < scalContent.second.timeTolerance))
          return false;

      return true;
    }

    /// Returns whether time tolerance of component associated with key is reached.
    virtual bool timeToleranceReached(Key key) const
    {
      if (!(getTimeEstCombined(key) < getScalContent(key).timeTolerance))
        return false;
      else
        return true;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    bool timeToleranceReached(const TP& tp) const
    {
      return timeToleranceReached(to_string(tp));
    }

    /// Returns whether time error is under its lower bound.
    virtual bool timeErrorLow() const
    {
      for (auto const& scalContent : scalContents)
        if (!(getTimeEstCombined(scalContent.first) < scalContent.second.timeErrLow))
          return false;

      return true;
    }

    /// Returns the time estimation as a combination
    /// of maximal and integral time error
    double getTimeEstCombined(Key key) const
    {
      return
        getScalContent(key).est_t_max * getScalContent(key).fac_max +
        getScalContent(key).est_t_sum * getScalContent(key).fac_sum;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    double getTimeEstCombined(const TP& tp) const
    {
      return getTimeEstCombined(to_string(tp));
    }

    /// Print debug information about time error and its bound.
    void printTimeErrorLowInfo() const;

    /// Returns \ref spaceIteration.
    int getSpaceIteration() const
    {
      return spaceIteration;
    }

    /// Sets \ref spaceIteration.
    void setSpaceIteration(int it)
    {
      spaceIteration = it;
    }

    /// Returns \ref maxSpaceIteration.
    int getMaxSpaceIteration() const
    {
      return maxSpaceIteration;
    }

    /// Sets \ref maxSpaceIteration.
    void setMaxSpaceIteration(int it)
    {
      maxSpaceIteration = it;
    }

    /// Increments \ref spaceIteration by 1;
    void incSpaceIteration()
    {
      spaceIteration++;
    }

    /// Sets \ref timestepIteration.
    void setTimestepIteration(int it)
    {
      timestepIteration = it;
    }

    /// Returns \ref timestepIteration.
    int getTimestepIteration() const
    {
      return timestepIteration;
    }

    /// Increments \ref timestepIteration by 1;
    void incTimestepIteration()
    {
      timestepIteration++;
    }

    /// Returns \ref maxTimestepIteration.
    int getMaxTimestepIteration() const
    {
      return maxTimestepIteration;
    }

    /// Sets \ref maxTimestepIteration.
    void setMaxTimestepIteration(int it)
    {
      maxTimestepIteration = it;
    }

    /// Sets \ref timeIteration.
    void setTimeIteration(int it)
    {
      timeIteration = it;
    }

    /// Returns \ref timeIteration.
    int getTimeIteration() const
    {
      return timeIteration;
    }

    /// Increments \ref timesIteration by 1;
    void incTimeIteration()
    {
      timeIteration++;
    }

    /// Returns \ref maxTimeIteration.
    int getMaxTimeIteration() const
    {
      return maxTimeIteration;
    }

    /// Sets \ref maxTimeIteration.
    void setMaxTimeIteration(int it)
    {
      maxTimeIteration = it;
    }

    /// Returns \ref timestepNumber.
    int getTimestepNumber() const
    {
      return timestepNumber;
    }

    /// Sets \ref timestepNumber.
    void setTimestepNumber(int num)
    {
      timestepNumber = std::min(nTimesteps, num);
    }

    /// Returns \ref nTimesteps.
    int getNumberOfTimesteps() const
    {
      return nTimesteps;
    }

    /// Sets \ref nTimesteps.
    void setNumberOfTimesteps(int num)
    {
      nTimesteps = std::max(0, num);
    }

    /// Increments \ref timestepNumber by 1;
    void incTimestepNumber()
    {
      timestepNumber++;
    }

    /// Sets \ref est_sum.
    void setEstSum(double e, Key key)
    {
      getScalContent(key).est_sum = e;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    void setEstSum(double e, const TP& tp)
    {
      setEstSum(e, to_string(tp));
    }

    /// Sets \ref est_max.
    void setEstMax(double e, Key key)
    {
      getScalContent(key).est_max = e;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    void setEstMax(double e, const TP& tp)
    {
      setEstMax(e, to_string(tp));
    }

    /// Sets \ref est_max.
    void setTimeEstMax(double e, Key key)
    {
      getScalContent(key).est_t_max = e;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    void setTimeEstMax(double e, const TP& tp)
    {
      setTimeEstMax(e, to_string(tp));
    }

    /// Sets \ref est_t_sum.
    void setTimeEstSum(double e, Key key)
    {
      getScalContent(key).est_t_sum = e;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    void setTimeEstSum(double e, const TP& tp)
    {
      setTimeEstSum(e, to_string(tp));
    }

    /// Returns \ref est_sum.
    double getEstSum(Key key) const
    {
      return getScalContent(key).est_sum;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    double getEstSum(const TP& tp)
    {
      return getEstSum(to_string(tp));
    }

    /// Returns \ref est_t_sum.
    double getEstTSum(Key key) const
    {
      return getScalContent(key).est_t_sum;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    double getEstTSum(const TP& tp)
    {
      return getEstTSum(to_string(tp));
    }

    /// Returns \ref est_max.
    double getEstMax(Key key) const
    {
      return getScalContent(key).est_max;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    double getEstMax(const TP& tp)
    {
      return getEstMax(to_string(tp));
    }

    /// Returns \ref est_max.
    double getTimeEstMax(Key key) const
    {
      return getScalContent(key).est_t_max;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    double getTimeEstmax(const TP& tp)
    {
      return getTimeEstMax(to_string(tp));
    }

    /// Returns \ref est_t_sum.
    double getTimeEstSum(Key key) const
    {
      return getScalContent(key).est_t_sum;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    double getTimeEstSum(const TP& tp)
    {
      return getTimeEstSum(to_string(tp));
    }

    /// Returns \ref est_t the estimated overall time error
    double getTimeEst() const
    {
      return est_t;
    }

    void setTimeEst(double value)
    {
      est_t = value;
    }

    /// Returns \ref spaceTolerance.
    double getSpaceTolerance(Key key) const
    {
      return getScalContent(key).spaceTolerance;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    double getSpaceTolerance(const TP& tp)
    {
      return getSpaceTolerance(to_string(tp));
    }

    /// Sets \ref spaceTolerance.
    void setSpaceTolerance(Key key, double tol)
    {
      getScalContent(key).spaceTolerance = tol;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    void setSpaceTolerance(const TP& tp, double tol)
    {
      return setSpaceTolerance(to_string(tp), tol);
    }

    /// Returns \ref timeTolerance.
    double getTimeTolerance(Key key) const
    {
      return getScalContent(key).timeTolerance;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    double getTimeTolerance(const TP& tp)
    {
      return getTimeTolerance(to_string(tp));
    }

    /// Returns \ref timeRelativeTolerance.
    double getTimeRelativeTolerance(Key key) const
    {
      return getScalContent(key).timeRelativeTolerance;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    double getTimeRelativeTolerance(const TP& tp)
    {
      return getTimeRelativeTolerance(to_string(tp));
    }

    /// Sets \ref time
    double setTime(double t)
    {
      time = t;
      if (time > endTime)
        time = endTime;
      if (time < startTime)
        time = startTime;

      return time;
    }

    /// Gets \ref time
    double getTime() const
    {
      return time;
    }

    /// Gets \ref &time
    double const* getTimePtr() const
    {
      return &time;
    }

    /// Sets \ref timestep
    double setTimestep(double t)
    {
      timestep = t;
      if (timestep > maxTimestep)
        timestep = maxTimestep;
      if (timestep < minTimestep)
        timestep = minTimestep;
      if (time + timestep > endTime)
        timestep = endTime - time;

      return timestep;
    }
    /// Gets \ref timestep
    double getTimestep() const
    {
      return timestep;
    }

    void setLastProcessedTimestep(double t)
    {
      lastProcessedTimestep = t;
    }

    double getLastProcessedTimestep() const
    {
      return lastProcessedTimestep;
    }

    /// Returns true, if the end time is reached and no more timestep
    /// computations must be done.
    bool reachedEndTime() const
    {
      if (nTimesteps > 0)
        return !(timestepNumber < nTimesteps);

      return !(std::abs(time - endTime) > std::numeric_limits<double>::epsilon());
    }


    /// Sets \ref minTimestep
    void setMinTimestep(double t)
    {
      minTimestep = t;
    }

    /// Gets \ref minTimestep
    double getMinTimestep() const
    {
      return minTimestep;
    }

    /// Sets \ref maxTimestep
    void setMaxTimestep(double t)
    {
      maxTimestep = t;
    }

    /// Gets \ref maxTimestep
    double getMaxTimestep() const
    {
      return maxTimestep;
    }

    /// Gets \ref &timestep
    double const* getTimestepPtr() const
    {
      return &timestep;
    }

    /// Sets \ref startTime = time
    void setStartTime(double time)
    {
      startTime = time;
    }

    /// Sets \ref endTime = time
    void setEndTime(double time)
    {
      endTime = time;
    }

    /// Returns \ref startTime
    double getStartTime() const
    {
      return startTime;
    }

    /// Returns \ref endTime
    double getEndTime() const
    {
      return endTime;
    }

    /// Returns \ref timeErrLow.
    double getTimeErrLow(Key key) const
    {
      return getScalContent(key).timeErrLow;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    double getTimeErrLow(const TP& tp)
    {
      return getTimeErrLow(to_string(tp));
    }

    /// Returns whether coarsening is allowed or not.
    bool isCoarseningAllowed(Key key) const
    {
      return (getScalContent(key).coarsenAllowed == 1);
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    bool isCoarseningAllowed(const TP& tp)
    {
      return isCoarseningAllowed(to_string(tp));
    }

    /// Returns whether coarsening is allowed or not.
    bool isRefinementAllowed(Key key) const
    {
      return (getScalContent(key).refinementAllowed == 1);
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    bool isRefinementAllowed(const TP& tp)
    {
      return isRefinementAllowed(to_string(tp));
    }

    ///
    void allowRefinement(bool allow, Key key)
    {
      getScalContent(key).refinementAllowed = allow;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    void allowRefinement(bool allow, const TP& tp)
    {
      return allowRefinement(allow, to_string(tp));
    }

    ///
    void allowCoarsening(bool allow, Key key)
    {
      getScalContent(key).coarsenAllowed = allow;
    }

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    void allowCoarsening(bool allow, const TP& tp)
    {
      return allowCoarsening(allow, to_string(tp));
    }

    int getSize() const
    {
      return int(scalContents.size());
    }

    // TODO: remove from AdaptInfo
    bool getRosenbrockMode() const
    {
      return rosenbrockMode;
    }

    void setSolverIterations(int it)
    {
      solverIterations = it;
    }

    int getSolverIterations() const
    {
      return solverIterations;
    }

    void setMaxSolverIterations(int it)
    {
      maxSolverIterations = it;
    }

    int getMaxSolverIterations() const
    {
      return maxSolverIterations;
    }

    void setSolverTolerance(double tol)
    {
      solverTolerance = tol;
    }

    double getSolverTolerance() const
    {
      return solverTolerance;
    }

    void setSolverResidual(double res)
    {
      solverResidual = res;
    }

    double getSolverResidual() const
    {
      return solverResidual;
    }

    void setGlobalTimeTolerance(double tol)
    {
      globalTimeTolerance = tol;
    }

    double getGlobalTimeTolerance() const
    {
      return globalTimeTolerance;
    }

    // TODO: remove from AdaptInfo
    void setRosenbrockMode(bool b)
    {
      rosenbrockMode = b;
    }

    /** \brief
     * Resets timestep, current time and time boundaries without
     * any check. Is used by the parareal algorithm.
     */
    void resetTimeValues(double newTimeStep,
                         double newStartTime,
                         double newEndTime)
    {
      time = newStartTime;
      startTime = newStartTime;
      endTime = newEndTime;
      timestep = newTimeStep;
      timestepNumber = 0;
    }

  private:
    ScalContent& getScalContent(Key key) const
    {
      auto result = scalContents.emplace(std::piecewise_construct, std::forward_as_tuple(key), std::forward_as_tuple(name_ + "[" + key + "]") );
      return result.first->second;
    }

  protected:
    /// Name.
    std::string name_;

    /// Current space iteration
    int spaceIteration = -1;

    /** \brief
     * maximal allowed number of iterations of the adaptive procedure; if
     * maxIteration <= 0, no iteration bound is used
     */
    int maxSpaceIteration = -1;

    /// Current timestep iteration
    int timestepIteration = 0;

    /// Maximal number of iterations for choosing a timestep
    int maxTimestepIteration = 30;

    /// Current time iteration
    int timeIteration = 0;

    /// Maximal number of time iterations
    int maxTimeIteration = 30;

    /// Actual time, end of time interval for current time step
    double time = 0.0;

    /// Initial time
    double startTime = 0.0;

    /// Final time
    double endTime = 1.0;

    /// Time step size to be used
    double timestep = 0.0;

    /// Last processed time step size of finished iteration
    double lastProcessedTimestep = 0.0;

    /// Minimal step size
    double minTimestep = 0.0;

    /// Maximal step size
    double maxTimestep = 1.0;

    /// Number of current time step
    int timestepNumber = 0;

    /** \brief
     * Per default this value is 0 and not used. If it is set to a non-zero value,
     * the computation of the stationary problem is done nTimesteps times with a
     * fixed timestep.
     */
    int nTimesteps = 0;

    /// number of iterations needed of linear or nonlinear solver
    int solverIterations = 0;

    /// maximal number of iterations needed of linear or nonlinear solver
    int maxSolverIterations = 0;

    ///
    double solverTolerance = 1.e-8;

    ///
    double solverResidual = 0.0;

    /// tolerance for the overall time error
    double globalTimeTolerance = 1.0;

    /// Scalar adapt infos
    mutable std::map<Key, ScalContent> scalContents;

    /// Is true, if the adaptive procedure was deserialized from a file. TODO: remove deserialization
    bool deserialized = false;

    /// Is true, if the time adaption is controlled by a Rosenbrock Method.
    bool rosenbrockMode = false;

    /// overall time error estimate
    double est_t = 0.0;
  };

} // end namespace AMDiS