gradientFlow.cpp

Go to the documentation of this file.

#include "gradientFlow.h"

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_int("order_of_RungeKutta", 0);
    param.Register_double("step_size", 0.0);
    param.Register_int("order_of_approx_for_exp_iP", 0);
    param.Register_double("tolerance", 0.0);
    param.Register_double("safety_factor", 0.0);
    param.Register_int("adaptive", 0);

    param.Register_string("verbose_level", "NULL");
  }


#ifdef USE_PARAMETERS_FACTORY
  bool init_param = ParametersFactory::Register("GradientFlow", append_entry);
#endif
}
//- end

//- parameters class
Parameters_GradientFlow::Parameters_GradientFlow() { append_entry(*this); }
//- end

const std::string GradientFlow::class_name = "GradientFlow";

//====================================================================
GradientFlow::GradientFlow(Action *action)
  : m_vl(CommonParameters::Vlevel()),
    m_action(action),
    m_impl(0)
{
  initialize();
}


//====================================================================
GradientFlow::GradientFlow(unique_ptr<Action>& action)
  : m_vl(CommonParameters::Vlevel()),
    m_action(action.get()),
    m_impl(0)
{
  initialize();
}


//====================================================================
void GradientFlow::initialize()
{
  m_Nprec = 0;
  m_Estep = 0.0;
}


//====================================================================
void GradientFlow::set_parameters(const Parameters& params)
{
  const string str_vlevel = params.get_string("verbose_level");

  m_vl = vout.set_verbose_level(str_vlevel);

  //- fetch and check input parameters
  int    order_RK;
  double Estep;
  int    Nprec;
  double tol_RK, safety;
  int    is_adaptive;

  int err = 0;
  err += params.fetch_int("order_of_RungeKutta", order_RK);
  err += params.fetch_double("step_size", Estep);
  err += params.fetch_int("order_of_approx_for_exp_iP", Nprec);

  err += params.fetch_int("adaptive", is_adaptive);

  if (is_adaptive != 0) {  // mandatory only if adaptive is on.
    err += params.fetch_double("tolerance", tol_RK);
    err += params.fetch_double("safety_factor", safety);
  } else {
    params.fetch_double("tolerance", tol_RK);
    params.fetch_double("safety_factor", safety);
  }

  if (err) {
    vout.crucial(m_vl, "%s: fetch error, input parameter not found.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }

  set_parameters(order_RK, Estep, Nprec, tol_RK, safety, is_adaptive);
}


//====================================================================
void GradientFlow::set_parameters(
  const int order_RK,
  const double Estep, const int Nprec,
  const double tol_RK, const double safety,
  const int is_adaptive)
{
  //- print input parameters
  vout.general(m_vl, "Parameters of %s:\n", class_name.c_str());
  vout.general(m_vl, "  order_RK   = %d\n", order_RK);
  vout.general(m_vl, "  Estep      = %10.6f\n", Estep);
  vout.general(m_vl, "  Nprec      = %d\n", Nprec);
  vout.general(m_vl, "  tol_RK     = %10.6e\n", tol_RK);
  vout.general(m_vl, "  safety     = %10.6f\n", safety);
  vout.general(m_vl, "  adaptive   = %d\n", is_adaptive);

  //- range check
  int err = 0;
  err += ParameterCheck::non_negative(order_RK);
  err += ParameterCheck::square_non_zero(Estep);
  err += ParameterCheck::non_negative(Nprec);
  err += ParameterCheck::non_negative(tol_RK);
  err += ParameterCheck::non_negative(safety);

  if (err) {
    vout.crucial(m_vl, "%s: parameter range check failed.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }

  //- store values
  m_order_RK = order_RK; // order of Runge-Kutta
  m_Estep    = Estep;    // step size (SA)
  m_Nprec    = Nprec;    // order of approximation for e^{iP} (SA)
  m_tol_RK   = tol_RK;
  m_safety   = safety;
  m_adaptive = (is_adaptive == 0) ? false : true;

  set_parameter_order(order_RK, m_adaptive);
}


//====================================================================
class GradientFlow::RungeKutta
{
 public:
  RungeKutta(Action *action, int nprec, Bridge::VerboseLevel vl)
    : m_action(action),
      m_Nprec(nprec),
      m_Ndim(CommonParameters::Ndim()),
      m_Nvol(CommonParameters::Nvol()),
      m_vl(vl)
  {}

  virtual ~RungeKutta() {}

  virtual void flow(double& t, double& Estep, Field_G& U) = 0;

  void update_u(Field_G& W, const double Estep, Field_G& iP, Field_G& U);

  virtual int order() const = 0;

 protected:
  Action *m_action;
  int    m_Nprec;

  int m_Ndim;
  int m_Nvol;

  Bridge::VerboseLevel m_vl;
};

//====================================================================
namespace {
  // w = exp(alpha iv) u
  void exponential(Mat_SU_N& w, const double alpha, const Mat_SU_N& iv, const Mat_SU_N& u, const int nprec)
  {
    Mat_SU_N p(CommonParameters::Nc());

    w = u;

    for (int iprec = 0; iprec < nprec; ++iprec) {
      double exf = alpha / (nprec - iprec);

      p  = iv * w;
      p *= exf;
      w  = p;
      w += u;  // w' = u + 1/(n-k) alpha iv ( w )
    }

    w.reunit();
  }
}

//====================================================================
void GradientFlow::RungeKutta::update_u(Field_G& W, const double Estep, Field_G& iP, Field_G& U)
{
  int Nvol = U.nvol();
  int Nex  = U.nex();
  int Nc   = CommonParameters::Nc();

  Mat_SU_N u0(Nc), v0(Nc), w0(Nc);

  for (int ex = 0; ex < Nex; ++ex) {
    for (int site = 0; site < Nvol; ++site) {
      u0 = U.mat(site, ex);
      v0 = iP.mat(site, ex);

      exponential(w0, Estep, v0, u0, m_Nprec);

      W.set_mat(site, ex, w0);
    }
  }
}


//====================================================================
class GradientFlow::RungeKutta_1st : public GradientFlow::RungeKutta
{
 public:
  RungeKutta_1st(Action *action, int nprec, Bridge::VerboseLevel vl)
    : RungeKutta(action, nprec, vl)
  {
    z0.reset(m_Nvol, m_Ndim);
  }

  void flow(double& t, double& Estep, Field_G& U);

  int order() const { return 1; }

 private:

  Field_G z0;
};

//====================================================================
void GradientFlow::RungeKutta_1st::flow(double& t, double& Estep, Field_G& U)
{
  Field_G& w0 = U;
  Field_G& w1 = U;

  // Field_G z0 = w0.clone();

  //- step 0
  // calculate gradient of m_action Z_0 (SA)
  m_action->force(z0, &w0);

  // W_1=e^{Z_0}*U
  update_u(w1, Estep, z0, w0);  // assume it is ok for w = u case.

  t += Estep;
}


//====================================================================
class GradientFlow::RungeKutta_2nd : public GradientFlow::RungeKutta
{
 public:
  RungeKutta_2nd(Action *action, int nprec, Bridge::VerboseLevel vl)
    : RungeKutta(action, nprec, vl)
  {
    w1.reset(m_Nvol, m_Ndim);
    z0.reset(m_Nvol, m_Ndim);
    z1.reset(m_Nvol, m_Ndim);
  }

  void flow(double& t, double& Estep, Field_G& U);

  int order() const { return 2; }

 private:
  Field_G w1;
  Field_G z0, z1;
};

//====================================================================
void GradientFlow::RungeKutta_2nd::flow(double& t, double& Estep, Field_G& U)
{
  Field_G& w0 = U;
  // Field_G w1 = w0.clone();
  Field_G& w2 = U;

  // Field_G z0 = w0.clone();
  // Field_G z1 = w0.clone();
  // Field_G zt = w0.clone();

  //- step 0
  // calculate gradient of m_action Z_0 (SA)
  m_action->force(z0, &w0);

  // W_1=e^{c2*Z_0}*U
  update_u(w1, 0.5 * Estep, z0, w0);

  //- step 1
  // Z_1
  m_action->force(z1, &w1);

  // V_out=e^{Z_1}*W_0
  update_u(w2, Estep, z1, w0);

  t += Estep;
}


//====================================================================
class GradientFlow::RungeKutta_3rd : public GradientFlow::RungeKutta
{
 public:
  RungeKutta_3rd(Action *action, int nprec, Bridge::VerboseLevel vl)
    : RungeKutta(action, nprec, vl)
  {
    w1.reset(m_Nvol, m_Ndim);
    w2.reset(m_Nvol, m_Ndim);

    z0.reset(m_Nvol, m_Ndim);
    z1.reset(m_Nvol, m_Ndim);
    z2.reset(m_Nvol, m_Ndim);

    zt.reset(m_Nvol, m_Ndim);
  }

  void flow(double& t, double& Estep, Field_G& U);

  int order() const { return 3; }

 private:
  Field_G w1, w2;
  Field_G z0, z1, z2;
  Field_G zt;
};

//====================================================================
void GradientFlow::RungeKutta_3rd::flow(double& t, double& Estep, Field_G& U)
{
  const double a0 = 0.25;

  const double b0 = -17.0 / 36.0;
  const double b1 = 8.0 / 9.0;

  const double c0 = 17.0 / 36.0;
  const double c1 = -8.0 / 9.0;
  const double c2 = 0.75;

  Field_G& w0 = U;
  // Field_G w1 = w0.clone();
  // Field_G w2 = w0.clone();
  Field_G& w3 = U;

  // Field_G z0 = w0.clone();
  // Field_G z1 = w0.clone();
  // Field_G z2 = w0.clone();

  // Field_G zt = w0.clone();

  //- step 0
  // calculate gradient of m_action Z_0 (SA)
  m_action->force(z0, &w0);

  // zt = a0 * z0 = 1/4 z0
  copy(zt, z0);
  scal(zt, a0);

  // W_1=e^{Z_0/4}*U (SA)
  update_u(w1, Estep, zt, w0);

  //- step 1
  // Z_1 (SA)
  m_action->force(z1, &w1);

  // zt = b1 * z1 + b0 * z0 = 8/9 z1 - 17/36 z0
  copy(zt, z1);
  scal(zt, b1);
  axpy(zt, b0, z0);

  // W_2=e^{8*Z_1/9-17*Z_0/36}*W_1 (SA)
  update_u(w2, Estep, zt, w1);

  //- step 2
  // Z_2 (SA)
  m_action->force(z2, &w2);

  // zt = c2 * z2 + c1 * z1 * c0 * z0
  //    = 3/4 z2 - 8/9 z1 + 17/36 z0
  //    = 3/4 z2 - zt_prev
  copy(zt, z2);
  scal(zt, c2);
  axpy(zt, c1, z1);
  axpy(zt, c0, z0);

  // V_out=e^{3*Z_2/4-8*Z_1/9+17*Z_0/36}*W_2 (SA)
  update_u(w3, Estep, zt, w2);

  t += Estep;
}


//====================================================================
class GradientFlow::RungeKutta_4th : public GradientFlow::RungeKutta
{
 public:
  RungeKutta_4th(Action *action, int nprec, Bridge::VerboseLevel vl)
    : RungeKutta(action, nprec, vl)
  {
    w1.reset(m_Nvol, m_Ndim);
    w2.reset(m_Nvol, m_Ndim);
    w3.reset(m_Nvol, m_Ndim);
    w4.reset(m_Nvol, m_Ndim);

    z0.reset(m_Nvol, m_Ndim);
    z1.reset(m_Nvol, m_Ndim);
    z2.reset(m_Nvol, m_Ndim);
    z3.reset(m_Nvol, m_Ndim);
    z4.reset(m_Nvol, m_Ndim);

    zt.reset(m_Nvol, m_Ndim);
  }

  void flow(double& t, double& Estep, Field_G& U);

  int order() const { return 4; }

 private:
  Field_G w1, w2, w3, w4;
  Field_G z0, z1, z2, z3, z4;
  Field_G zt;
};

//====================================================================
void GradientFlow::RungeKutta_4th::flow(double& t, double& Estep, Field_G& U)
{
  Field_G& w0 = U;
  // Field_G w1 = w0.clone();
  // Field_G w2 = w0.clone();
  // Field_G w3 = w0.clone();
  // Field_G w4 = w0.clone();
  Field_G& w5 = U;

  // Field_G z0 = w0.clone();
  // Field_G z1 = w0.clone();
  // Field_G z2 = w0.clone();
  // Field_G z3 = w0.clone();
  // Field_G z4 = w0.clone();

  // Field_G zt = w0.clone();

  //- step 0
  // calculate gradient of m_action Z_0 (SA)
  m_action->force(z0, &w0);

  copy(zt, z0);
  scal(zt, 0.5);

  // W_1=e^{Z_0/2}*U (SA)
  update_u(w1, Estep, zt, w0);

  //- step 1
  // Z_1
  m_action->force(z1, &w1);

  copy(zt, z1);
  scal(zt, 0.5);

  // W_2=e^{(1/2)*Z_1}*W_0
  update_u(w2, Estep, zt, w0);

  //- step 2
  // Z_2
  m_action->force(z2, &w2);

  // Z_2 - (1/2)*Z_0
  copy(zt, z2);
  axpy(zt, -0.5, z0);

  // W_3=e^{Z_2 - (1/2)*Z_0}*W_1  NB. not W_0
  update_u(w3, Estep, zt, w1);

  //- step 3
  // Z_3
  m_action->force(z3, &w3);

  // 3*Z_0+2*Z_1+2*Z_2-Z_3
  copy(zt, z0);
  scal(zt, 0.25);
  axpy(zt, 1.0 / 6.0, z1);
  axpy(zt, 1.0 / 6.0, z2);
  axpy(zt, -1.0 / 12.0, z3);

  // W_4=e^{(1/12)*(3*Z_0+2*Z_1+2*Z_2-Z_3)}*W_0
  update_u(w4, Estep, zt, w0);

  //- step 4 (extra step)

  // 1/12*(-Z_0+2*Z_1+2*Z_2+3*Z_3)
  copy(zt, z0);
  scal(zt, -1.0 / 12.0);
  axpy(zt, 1.0 / 6.0, z1);
  axpy(zt, 1.0 / 6.0, z2);
  axpy(zt, 0.25, z3);

  // V_t=e^{(1/12)*(-Z_0+2*Z_1+2*Z_2+3*Z_3)}*W_4  NB. not W_0
  update_u(w5, Estep, zt, w4);

  t += Estep;
}


//====================================================================
class GradientFlow::AdaptiveRungeKutta : public RungeKutta
{
 public:
  AdaptiveRungeKutta(RungeKutta *rk, double tol, double safety, Bridge::VerboseLevel vl)
    : RungeKutta(0, 0, vl),
      m_rk(rk),
      m_tol_RK(tol),
      m_safety(safety)
  {
    u_rough.reset(m_Nvol, m_Ndim);
    u_fine.reset(m_Nvol, m_Ndim);
  }

  ~AdaptiveRungeKutta() { if (m_rk) delete m_rk; }

  void flow(double& t, double& Estep, Field_G& U);

  int order() const { return m_rk ? m_rk->order() : 0; }

  double max_diff(const Field_G& u1, const Field_G& u2) const;

 protected:
  RungeKutta *m_rk;

  double m_tol_RK;
  double m_safety;

  Field_G u_rough, u_fine;
};

//====================================================================
void GradientFlow::AdaptiveRungeKutta::flow(double& t, double& Estep, Field_G& U)
{
  copy(u_rough, U);
  copy(u_fine, U);

  while (true)
  {
    //- rough step
    double t2     = t;
    double estep2 = 2.0 * Estep;

    m_rk->flow(t2, estep2, u_rough);

    //- two original steps
    double t1     = t;
    double estep1 = Estep;

    m_rk->flow(t1, estep1, u_fine);
    m_rk->flow(t1, estep1, u_fine);

    double diff = max_diff(u_rough, u_fine);

    vout.general(m_vl, "  Estep,diff,m_tol_RK = %e %e %e\n", Estep, diff, m_tol_RK);

    if (diff < m_tol_RK) {
      copy(U, u_fine);
      t += estep2;

      // extend stepsize
      Estep *= m_safety * pow(m_tol_RK / diff, 1.0 / (order() + 1));

      break;  // exit loop
    } else {
      // shrink stepsize
      Estep *= m_safety * pow(m_tol_RK / diff, 1.0 / (order() + 1));

      if (Estep < m_tol_RK) {
        vout.crucial(m_vl, "%s: too small Estep = %e\n", class_name.c_str(), Estep);
        exit(EXIT_FAILURE);
      }

      // proceed to next trial
    }
  }
}


//====================================================================
double GradientFlow::AdaptiveRungeKutta::max_diff(const Field_G& f1, const Field_G& f2) const
{
  int Nvol = f1.nvol();
  int Nex  = f2.nex();
  int Nc   = CommonParameters::Nc();

  Mat_SU_N u0(Nc), u1(Nc), u2(Nc);

  double max_norm = 0.0;

  for (int ex = 0; ex < Nex; ++ex) {
    for (int site = 0; site < Nvol; ++site) {
      f1.mat(u0, site, ex);
      f2.mat(u1, site, ex);

      u2 = u0 - u1;

      double norm = sqrt(u2.norm2()) / Nc;

      if (norm > max_norm) max_norm = norm;
    }
  }

  max_norm = Communicator::reduce_max(max_norm);

  return max_norm;
}


//====================================================================
GradientFlow::~GradientFlow()
{
  if (m_impl) delete m_impl;
}


//====================================================================
void GradientFlow::set_parameter_order(const int order, bool is_adaptive)
{
  m_order_RK = order;

  switch (order)
  {
  case 1:
    m_impl = new RungeKutta_1st(m_action, m_Nprec, m_vl);
    break;

  case 2:
    m_impl = new RungeKutta_2nd(m_action, m_Nprec, m_vl);
    break;

  case 3:
    m_impl = new RungeKutta_3rd(m_action, m_Nprec, m_vl);
    break;

  case 4:
    m_impl = new RungeKutta_4th(m_action, m_Nprec, m_vl);
    break;

  default:
    vout.crucial(m_vl, "%s: order of Runge-Kutta is out of range.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }

  if (is_adaptive) {
    m_impl = new AdaptiveRungeKutta(m_impl, m_tol_RK, m_safety, m_vl);
  }
}


//====================================================================
double GradientFlow::evolve(double& t_flow, Field_G& U)
{
  double plaq = m_staple.plaquette(U);

  vout.general(m_vl, "  plaq_org = %.16f\n", plaq); // write-out

  //- time evolution (SA)
  m_impl->flow(t_flow, m_Estep, U);

//  t_flow += m_Estep;

  plaq = m_staple.plaquette(U);
  double Eplaq = t_flow * t_flow * 36.0 * (1.0 - plaq); // t^2*E=t^2*36*(1-Plaq)

  vout.general(m_vl, "  (t, plaq, t^2 E_plaq) = %.8f %.16f %.16f\n", t_flow, plaq, Eplaq);

  double result = Eplaq;

  return result;
}


//====================================================================
//============================================================END=====