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/typetree/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 (!(scalContent(key).est_sum < scalContent(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 (!(timeEstCombined(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 (!(timeEstCombined(key) < scalContent(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 (!(timeEstCombined(scalContent.first) < scalContent.second.timeErrLow))
          return false;

      return true;
    }

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

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

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

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

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

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

    /// Sets \ref maxSpaceIteration_.
    void maxSpaceIteration(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 timestepIteration() const
    {
      return timestepIteration_;
    }

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

    /// Returns \ref maxTimestepIteration_.
    int maxTimestepIteration() 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 timeIteration() const
    {
      return timeIteration_;
    }

    /// Increments \ref timesIteration_ by 1;
    void incTimeIteration()
    {
      timeIteration_++;
    }

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

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

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

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

    /// Returns \ref nTimesteps_.
    int numberOfTimesteps() 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)
    {
      scalContent(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)
    {
      scalContent(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)
    {
      scalContent(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)
    {
      scalContent(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 estSum(Key key) const
    {
      return scalContent(key).est_sum;
    }

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

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

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

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

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

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

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

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

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

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

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

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

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

    /// Sets \ref spaceTolerance.
    void setSpaceTolerance(Key key, double tol)
    {
      scalContent(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 timeTolerance(Key key) const
    {
      return scalContent(key).timeTolerance;
    }

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

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

    template <class TP, REQUIRES( Concepts::PreTreePath<TP> )>
    double timeRelativeTolerance(const TP& tp)
    {
      return timeRelativeTolerance(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 const& time() 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 const& timestep() const
    {
      return timestep_;
    }

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

    double lastProcessedTimestep() 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 minTimestep() const
    {
      return minTimestep_;
    }

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

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

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

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

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

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

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

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

    /// Returns whether coarsening is allowed or not.
    bool isCoarseningAllowed(Key key) const
    {
      return (scalContent(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 (scalContent(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)
    {
      scalContent(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)
    {
      scalContent(key).coarsenAllowed = allow;
    }

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

    int size() const
    {
      return int(scalContents_.size());
    }

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

    int solverIterations() const
    {
      return solverIterations_;
    }

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

    int maxSolverIterations() const
    {
      return maxSolverIterations_;
    }

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

    double solverTolerance() const
    {
      return solverTolerance_;
    }

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

    double solverResidual() const
    {
      return solverResidual_;
    }

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

    double globalTimeTolerance() const
    {
      return globalTimeTolerance_;
    }


    /** \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& scalContent(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_;

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

} // end namespace AMDiS