gradientFlow.cpp

Go to the documentation of this file.

#include "gradientFlow.h"

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

const double GradientFlow::c4[] =
{ 0.25, -17.0 / 36.0, 8.0 / 9.0, 0.75, };

const double GradientFlow::c5[] =
{
  1.0 / 6.0,
  1.0 - 1.0 / sqrt(2.0),
  2.0 - c5[1], /* 1.0+1.0/sqrt(2.0) */
  c5[1] - 0.5,
  -4.0 * c5[2],
  6.0 - 4.0 * c5[1],
  1.0 / c5[3],
};


//- 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

//====================================================================
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, "GradientFlow: fetch error, input parameter not found.\n");
    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 GradientFlow:\n");
  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, "GradientFlow: parameter range check failed.\n");
    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 == 4) {
      //- 4th order Runge-Kutta (SA)
      gradientFlow_4th(U);
    } else if (m_order_RK == 5) {
      //- 5th order Runge-Kutta (SA)
      gradientFlow_5th(U);
    } else {
      vout.crucial(m_vl, "GradientFlow: order of Runge-Kutta is out of range.\n");
      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, t^2 E_plaq) = %.8f  %.16f\n", tt, 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_4th(Field_G& U)
{
  //- step 0
  iP = (Field_G)m_action->force();  // calculate gradient of m_action Z_0 (SA)
  update_U(c4[0] * m_Estep, iP, U); // W_1=e^{Z_0/4}*U (SA)

  //- step 1
  iP  *= c4[1];                      // -(17/36)Z_0 (SA)
  iPP  = (Field_G)m_action->force(); // Z_1 (SA)
  iPP *= c4[2];                      // (8/9)*Z_1 (SA)
  iP  += iPP;                        // calculate (8/9)*Z_1-(17/36)*Z_0 (SA)
  update_U(m_Estep, iP, U);          // W_2=e^{8*Z_1/9-17*Z_0/36}*W_1 (SA)

  //- step 2
  iPP  = (Field_G)m_action->force(); // Z_2 (SA)
  iPP *= c4[3];                      // (3/4)*Z_2 (SA)
  iPP -= iP;                         // calculate (3/4)*Z_2-(8/9)*Z_1+(17/36)*Z_0 (SA)
  update_U(m_Estep, iPP, U);         // V_out=e^{3*Z_2/4-8*Z_1/9+17*Z_0/36}*W_2 (SA)
}


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

  //- step 1
  iP1  = (Field_G)m_action->force(); // Z_1
  iPP  = iP1;
  iPP -= iP;                         // calculate Z_1-Z_0 (SA)
  update_U(c5[1] * m_Estep, iPP, U); // W_2=e^{c_1*(Z_1-Z_0)}*W_1 (SA)

  //- step 2
  iPP  = (Field_G)m_action->force(); // Z_2
  iP2  = iPP;                        // Z_2
  iPP *= c5[2];                      // c_2*Z_2
  iPP -= iP1;                        // c_2*Z_2-Z_1
  iP  *= c5[3];                      // c_3*Z_0
  iPP += iP;                         // c_2*Z_2-Z_1+c_3*Z_0
  update_U(m_Estep, iPP, U);         // W_3=e^{c_2*Z_2-Z_1+c_3*Z_0}*W_2 (SA)

  iPP  = (Field_G)m_action->force(); // Z_3
  iP2 *= c5[4];                      // c_4*Z_2
  iPP += iP2;                        // Z_3+c_4*Z_2
  iP1 *= c5[5];                      // c_5*Z_1
  iPP += iP1;                        // Z_3+c_4*Z_2+c_5*Z_1
  iP  *= c5[6];                      // Z_0
  iPP += iP;                         // Z_3+c_4*Z_2+c_5*Z_1+Z_0
  update_U(c5[0] * m_Estep, iPP, U); // V_t=e^{(Z_3+c_4*Z_2+c_5*Z_1+Z_0)/6}*W_3 (SA)
}


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