gradientFlow.cpp

Go to the documentation of this file.

#include "gradientFlow.h"

#ifdef USE_PARAMETERS_FACTORY
#include "parameters_factory.h"
#endif


//- 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("number_of_steps", 0);
    param.Register_int("order_of_approx_for_exp_iP", 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";

//====================================================================
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    Nstep, Nprec;

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

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


  set_parameters(order_RK, Estep, Nstep, Nprec);
}

//====================================================================
void GradientFlow::set_parameters(const int order_RK,
                                  const double Estep, const int Nstep, const int Nprec)
{
  //- 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, "  Nstep    = %d\n", Nstep);
  vout.general(m_vl, "  Nprec    = %d\n", Nprec);

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

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

  //- store values
  m_order_RK = order_RK; // order of Runge-Kutta
  m_Estep    = Estep;    // step size (SA)
  m_Nstep    = Nstep;    // a number of steps (SA)
  m_Nprec    = Nprec;    // order of approximation for e^{iP} (SA)
}

//====================================================================
double GradientFlow::evolve(Field_G& U)
{
  double Eplaq = 0.0;  // superficial initialization
  double tt;

  m_action->set_config(&U); // give gauge conf pointer &U to d_action (SA)

  Staples wl;
  double  plaq = wl.plaquette(U);                   // calculate plaqutte (SA)

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

  //- time evolution (SA)
  for (int i = 0; i < m_Nstep; ++i) {
    if (m_order_RK == 1) {
      gradientFlow_1st(U);
    } else if (m_order_RK == 2) {
      gradientFlow_2nd(U);
    } else if (m_order_RK == 3) {
      gradientFlow_3rd(U);
    } else if (m_order_RK == 4) {
      gradientFlow_4th(U);
    } else {
      vout.crucial(m_vl, "%s: order of Runge-Kutta is out of range.\n", class_name.c_str());
      abort();
    }

    plaq  = wl.plaquette(U);     // calculate plaqutte (SA)
    Eplaq = 36.0 * (1.0 - plaq); // E=36*(1-Plaq)
    tt    = (i + 1) * m_Estep;
    Eplaq = tt * tt * Eplaq;

    vout.general(m_vl, "  (t, plaq, t^2 E_plaq) = %.8f %.16f %.16f\n", tt, plaq, Eplaq); // write-out
  }

  double result = Eplaq;

  return result;
}

//====================================================================
void GradientFlow::update_U(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), u1(Nc), u2(Nc);
  Mat_SU_N h1(Nc);

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

      for (int iprec = 0; iprec < m_Nprec; ++iprec) {
        double exf = estep / (m_Nprec - iprec); // step_size/(N-k) (SA)

        u2  = h1 * u1;                          // iP*u1 (SA)
        u2 *= exf;                              // step_size*iP*u1/(N-k) (SA)
        u1  = u2;
        u1 += u0;                               // U + step_size*iP*u1/(N-k) (SA)
      }

      // u1 = sum_{k=0}^{N-1} (step_size * iP)^k/ k!*U,  N=m_Nprec (SA)
      u1.reunit();             // reunitarize u1 (SA)
      U.set_mat(site, ex, u1); // U=u1 (SA)
    }
  }

  m_action->notify_linkv(); // notify action about update of U (SA)
}

//====================================================================
void GradientFlow::gradientFlow_1st(Field_G& U)
{
  //- step 0
  m_iP = (Field_G)m_action->force();  // calculate gradient of m_action Z_0 (SA)
  update_U(m_Estep, m_iP, U);         // W_1=e^{Z_0}*U
}

//====================================================================
void GradientFlow::gradientFlow_2nd(Field_G& U)
{
  const double c2 = 0.5;

  copy(m_U_0, U);

  //- step 0
  m_iP = (Field_G)m_action->force();  // calculate gradient of m_action Z_0 (SA)
  update_U(c2 * m_Estep, m_iP, U);    // W_1=e^{c2*Z_0}*U

  //- step 1
  m_iPP = (Field_G)m_action->force(); // Z_1
  copy(U, m_U_0);
  update_U(m_Estep, m_iPP, U);        // V_out=e^{Z_1}*W_0
}

//====================================================================
void GradientFlow::gradientFlow_3rd(Field_G& U)
{
  const double c3[] =
  { 0.25, -17.0 / 36.0, 8.0 / 9.0, 0.75, };

  //- step 0
  m_iP = (Field_G)m_action->force();  // calculate gradient of m_action Z_0 (SA)
  update_U(c3[0] * m_Estep, m_iP, U); // W_1=e^{Z_0/4}*U (SA)

  //- step 1
  scal(m_iP, c3[1]);                  // -(17/36)Z_0 (SA)
  m_iPP = (Field_G)m_action->force(); // Z_1 (SA)
  scal(m_iPP, c3[2]);                 // (8/9)*Z_1 (SA)
  axpy(m_iP, 1.0, m_iPP);             // (8/9)*Z_1-(17/36)*Z_0 (SA)
  update_U(m_Estep, m_iP, U);         // W_2=e^{8*Z_1/9-17*Z_0/36}*W_1 (SA)

  //- step 2
  m_iPP = (Field_G)m_action->force(); // Z_2 (SA)
  scal(m_iPP, c3[3]);                 // (3/4)*Z_2 (SA)
  axpy(m_iPP, -1.0, m_iP);            // (3/4)*Z_2-(8/9)*Z_1+(17/36)*Z_0 (SA)
  update_U(m_Estep, m_iPP, U);        // V_out=e^{3*Z_2/4-8*Z_1/9+17*Z_0/36}*W_2 (SA)
}

//====================================================================
void GradientFlow::gradientFlow_4th(Field_G& U)
{
  const double c4[] =
  {
    0.5,
    -0.5,
    3.0,
    2.0,
    -1.0,
    1.0 / 12.0,
  };

  copy(m_U_0, U);  // m_U_0 = U

  //- step 0
  m_iP = (Field_G)m_action->force();   // calculate gradient of m_action Z_0 (SA)
  update_U(c4[0] * m_Estep, m_iP, U);  // W_1=e^{Z_0/2}*U (SA)

  copy(m_W_1, U);                      // W_1 = U

  //- step 1
  m_iP1 = (Field_G)m_action->force();  // Z_1
  copy(U, m_U_0);                      // U = m_U_0
  update_U(c4[0] * m_Estep, m_iP1, U); // W_2=e^{(1/2)*Z_1}*W_0

  //- step 2
  m_iP2 = (Field_G)m_action->force();  // Z_2
  copy(m_iPP, m_iP2);                  // m_iPP = m_iP2
  axpy(m_iPP, c4[1], m_iP);            // Z_2 - (1/2)*Z_0
  copy(U, m_W_1);                      // U = W_1  NB. not m_U_0
  update_U(m_Estep, m_iPP, U);         // W_3=e^{Z_2 - (1/2)*Z_0}*W_1  NB. not W_0

  //- step 3
  m_iP3 = (Field_G)m_action->force();  // Z_3

  copy(m_iPP, m_iP);                   // m_iPP = Z_0
  scal(m_iPP, c4[2]);                  // 3*Z_0
  axpy(m_iPP, c4[3], m_iP1);           // 3*Z_0+2*Z_1
  axpy(m_iPP, c4[3], m_iP2);           // 3*Z_0+2*Z_1+2*Z_2
  axpy(m_iPP, c4[4], m_iP3);           // 3*Z_0+2*Z_1+2*Z_2-Z_3
  copy(U, m_U_0);                      // U = m_U_0
  update_U(c4[5] * m_Estep, m_iPP, U); // W_4=e^{(1/12)*(3*Z_0+2*Z_1+2*Z_2-Z_3)}*W_0

  //- step 4 (extra step)
  copy(m_iPP, m_iP);                   // m_iPP = Z_0
  scal(m_iPP, c4[4]);                  // -Z_0
  axpy(m_iPP, c4[3], m_iP1);           // -Z_0+2*Z_1
  axpy(m_iPP, c4[3], m_iP2);           // -Z_0+2*Z_1+2*Z_2
  axpy(m_iPP, c4[2], m_iP3);           // -Z_0+2*Z_1+2*Z_2+3*Z_3
  update_U(c4[5] * m_Estep, m_iPP, U); // V_t=e^{(1/12)*(-Z_0+2*Z_1+2*Z_2+3*Z_3)}*W_4  NB. not W_0
}

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