default.hh

Go to the documentation of this file.

#pragma once
 
// standard library
#include <iomanip>
#include <sstream>
 
// external libraries
#include <spdlog/spdlog.h>
 
namespace DiFfRG
{
  template <typename NEWT> class TC_Default
  {
  public:
    TC_Default(NEWT &newton_, unsigned int alg_order_, double t_, double max_t_, double dt_, double min_dt_,
               double max_dt_, double output_dt_)
        : newton(newton_), alg_order(alg_order_), t(t_), max_t(max_t_), sug_dt(dt_), min_dt(min_dt_), max_dt(max_dt_),
          output_dt(output_dt_), cur_dt(sug_dt), last_save(t), last_t(t_), stuck(0), fin(false)
    {
    }
 
    double get_dt() const { return cur_dt; }
    double get_t() const { return t; }
    unsigned int get_stuck() const { return stuck; }
    bool finished() const { return fin; }
 
    template <typename F, typename OF> void advance(F &f, OF &of)
    {
      // first, try to do the timestep
      try {
        f(t, cur_dt);
        // on a success, possibly amend dt
        this->step_success();
      } catch (const std::exception &e) {
        // on a fail, possibly amend dt
        this->step_fail(e);
      }
 
      // check if the stepper got stuck on some time t
      if (is_close(t, last_t))
        stuck++;
      else {
        stuck = 0;
        last_t = t;
      }
 
      // check if the stepper should output data
      if (t > last_save + output_dt) {
        of(t);
        last_save = int(t / output_dt) * output_dt;
      }
 
      // should the suggested timestep be below the minimum timestep, tell the solver that it should
      // accept even timesteps which do not fulfill the convergence criterion (otherwise the evolution would just stop)
      newton.set_ignore_nonconv(false);
      if (sug_dt < min_dt) {
        cur_dt = min_dt;
        newton.set_ignore_nonconv(true);
      } else if (sug_dt > max_dt)
        cur_dt = max_dt;
      else
        cur_dt = sug_dt;
 
      // do not over-step the final time
      if (t + cur_dt > max_t) cur_dt = max_t - t;
 
      // if we are over max_t, the timestepping is finished
      if (t - 1e-4 * cur_dt >= max_t) {
        of(t);
        fin = true;
      }
      // if we got unrecoverably stuck also abort the timestepping
      if (stuck > 10) {
        of(t);
        fin = true;
        spdlog::get("log")->error("Timestepping got stuck at t = {}", t);
      }
    }
 
    void print_step_info()
    {
      std::ostringstream message;
      message << std::setw(16) << "| step from t = " << std::setw(8) << std::setprecision(5) << last_t << std::setw(8)
              << " to t = ";
    }
 
  protected:
    virtual void step_success()
    {
      t += cur_dt;
      if (t + sug_dt - 1e-4 * sug_dt > last_save + output_dt)
        cur_dt = last_save + output_dt - t;
      else
        cur_dt = sug_dt;
    }
    virtual void step_fail(const std::exception &e)
    {
      throw std::runtime_error("Timestepping failed. Error: \n" + std::string(e.what()));
    }
 
    NEWT &newton;
    unsigned int alg_order;
    double t, max_t, sug_dt, min_dt, max_dt, output_dt, cur_dt, last_save, last_t;
    unsigned int stuck;
 
    bool fin;
  };
} // namespace DiFfRG