staple_SF.cpp

Go to the documentation of this file.

 
#include "staple_SF.h"
 
const std::string Staple_SF::class_name = "Staple_SF";
 
//====================================================================
void Staple_SF::set_parameters(const Parameters& params)
{
  std::string vlevel;
  if (!params.fetch_string("verbose_level", vlevel)) {
    m_vl = vout.set_verbose_level(vlevel);
  }
 
  //- fetch and check input parameters
  std::vector<double> phi, phipr, p_omega;
 
  int err = 0;
  err += params.fetch_double_vector("phi", phi);
  err += params.fetch_double_vector("phipr", phipr);
  err += params.fetch_double_vector("p_omega", p_omega);
 
  if (err) {
    vout.crucial(m_vl, "Error at %s: input parameter not found.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
 
 
  set_parameters(phi, phipr, p_omega);  // call std::vector version
}
 
 
//====================================================================
void Staple_SF::get_parameters(Parameters& params) const
{
  params.set_double_vector("phi", m_phi);
  params.set_double_vector("phipr", m_phipr);
  params.set_double_vector("p_omega", m_p_omega);
 
  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
}
 
 
//====================================================================
 
void Staple_SF::set_parameters(const std::vector<double>& phi,
                               const std::vector<double>& phipr,
                               const std::vector<double>& p_omega)
{
  //- print input parameters
  vout.general(m_vl, "%s:\n", class_name.c_str());
  vout.general(m_vl, "  phi1     = %12.6f\n", phi[0]);
  vout.general(m_vl, "  phi2     = %12.6f\n", phi[1]);
  vout.general(m_vl, "  phi3     = %12.6f\n", phi[2]);
  vout.general(m_vl, "  phipr1   = %12.6f\n", phipr[0]);
  vout.general(m_vl, "  phipr2   = %12.6f\n", phipr[1]);
  vout.general(m_vl, "  phipr3   = %12.6f\n", phipr[2]);
  vout.general(m_vl, "  p_omega1 = %12.6f\n", p_omega[0]);
  vout.general(m_vl, "  p_omega2 = %12.6f\n", p_omega[1]);
  vout.general(m_vl, "  p_omega3 = %12.6f\n", p_omega[2]);
 
  //- range check
  // NB. phi,phipr,p_omega == 0 is allowed.
  assert(phi.size() == 3);
  assert(phipr.size() == 3);
  assert(p_omega.size() == 3);
 
  m_i_omega0.zero();
  m_i_omega0.set(0, 0, 0.0, p_omega[0]);
  m_i_omega0.set(1, 1, 0.0, p_omega[1]);
  m_i_omega0.set(2, 2, 0.0, p_omega[2]);
 
  m_initialized = 1;
 
  // store raw parameter values
  m_phi.resize(3);
  m_phi[0] = phi[0];
  m_phi[1] = phi[1];
  m_phi[2] = phi[2];
 
  m_phipr.resize(3);
  m_phipr[0] = phipr[0];
  m_phipr[1] = phipr[1];
  m_phipr[2] = phipr[2];
 
  m_p_omega.resize(3);
  m_p_omega[0] = p_omega[0];
  m_p_omega[1] = p_omega[1];
  m_p_omega[2] = p_omega[2];
 
  Field_SF::set_boundary_matrix(m_wk, m_phi);
  Field_SF::set_boundary_matrix(m_wkpr, m_phipr);
}
 
 
//====================================================================
void Staple_SF::set_parameters(const std::vector<double>& phi,
                               const std::vector<double>& phipr)
{
  m_i_omega0.zero();
  m_i_omega0.set(0, 0, 0.0, 1.0);
  m_i_omega0.set(1, 1, 0.0, -0.5);
  m_i_omega0.set(2, 2, 0.0, -0.5);
 
  m_initialized = 1;
 
  // store raw parameter values
  m_phi.resize(3);
  m_phi[0] = phi[0];
  m_phi[1] = phi[1];
  m_phi[2] = phi[2];
 
  m_phipr.resize(3);
  m_phipr[0] = phipr[0];
  m_phipr[1] = phipr[1];
  m_phipr[2] = phipr[2];
 
  Field_SF::set_boundary_matrix(m_wk, m_phi);
  Field_SF::set_boundary_matrix(m_wkpr, m_phipr);
}
 
 
//====================================================================
 
double Staple_SF::sf_coupling_plaq(const Field_G& U, const double ct)
{
  if (!m_initialized) {
    vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
 
  const int Nc   = CommonParameters::Nc();
  const int Ndim = CommonParameters::Ndim();
 
  const int Nx   = CommonParameters::Nx();
  const int Ny   = CommonParameters::Ny();
  const int Nz   = CommonParameters::Nz();
  const int Nt   = CommonParameters::Nt();
  const int NPEt = CommonParameters::NPEt();
 
  const int Lx = CommonParameters::Lx();
 
  double plaq   = 0.0;
  double plaqt0 = 0.0;
  double plaqtT = 0.0;
 
  for (int nu = 0; nu < Ndim - 1; nu++) {
    Field_G staple;
    upper(staple, U, 3, nu);
 
    for (int z = 0; z < Nz; z++) {
      for (int y = 0; y < Ny; y++) {
        for (int x = 0; x < Nx; x++) {
          // boundary
          if (Communicator::ipe(3) == 0) {
            int t    = 0;
            int site = m_index.site(x, y, z, t);
 
            Mat_SU_N up(Nc);
            up = staple.mat(site) * U.mat_dag(site, 3);
            double scr = ReTr(m_i_omega0 * up);
 
            /*
            up.unit();
            up *= staple.mat(site);
            up *= U->mat_dag(site,3);
            up *= m_i_omega0;
            scr = ReTr( up );
            */
            plaq   -= scr;
            plaqt0 += scr;
          }
          // boundary
          if (Communicator::ipe(3) == NPEt - 1) {
            int t    = Nt - 1;
            int site = m_index.site(x, y, z, t);
 
            Mat_SU_N up(Nc);
            up = staple.mat_dag(site) * U.mat(site, 3);
            double scr = ReTr(m_i_omega0 * up);
 
            /*
            up.unit();
            up *= staple.mat_dag(site);
            up *= U->mat(site,3);
            up *= m_i_omega0;
            scr = ReTr( up );
            */
            plaq   += scr;
            plaqtT += scr;
          }
        }
      }
    }
  }
  plaq   = Communicator::reduce_sum(plaq);
  plaqt0 = Communicator::reduce_sum(plaqt0);
  plaqtT = Communicator::reduce_sum(plaqtT);
 
  plaq   *= ct / Lx;
  plaqt0 *= ct / Lx;
  plaqtT *= ct / Lx;
 
  vout.general(m_vl, "SF_delSg_plaq, from 0, from T = %.8f %.8f %.8f\n", plaq, plaqt0, plaqtT);
 
  return plaq;
}
 
 
//====================================================================
 
double Staple_SF::sf_coupling_rect(const Field_G& U, const double ctr)
{
  if (!m_initialized) {
    vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
 
  const int Nc   = CommonParameters::Nc();
  const int Ndim = CommonParameters::Ndim();
  const int Nvol = CommonParameters::Nvol();
 
  const int Nx   = CommonParameters::Nx();
  const int Ny   = CommonParameters::Ny();
  const int Nz   = CommonParameters::Nz();
  const int Nt   = CommonParameters::Nt();
  const int NPEt = CommonParameters::NPEt();
 
  const int Lx = CommonParameters::Lx();
 
  const int nu = 3;
 
  double rect01 = 0.0;
  double rect02 = 0.0;
  double rect03 = 0.0;
  double rectt1 = 0.0;
  double rectt2 = 0.0;
  double rectt3 = 0.0;
 
  for (int mu = 0; mu < Ndim - 1; mu++) {
    // rect01
    //      <---<---+
    //      |       |
    //  t=0 x--->---+
    //  omega0
 
    // rect02
    //      <---<---+
    //      |       |
    //  t=0 +---x---+
    //       omega0
 
    // rectt1
    // omega0
    // t=Nt x--->---+
    //      |       |
    //      +---<---+
 
    // rectt2
    //       omega0
    // t=Nt +---x---+
    //      |       |
    //      <---<---+
 
    Field_G Cup2;
    upper(Cup2, U, nu, mu);
 
    Field_G Umu;
    copy(Umu, 0, U, mu);
 
    Field_G Unu;
    copy(Unu, 0, U, nu);
 
    Field_G v;
    m_shift.backward(v, Cup2, mu);
 
    Field_G c;
    m_shift.backward(c, Umu, nu);
 
    for (int z = 0; z < Nz; z++) {
      for (int y = 0; y < Ny; y++) {
        for (int x = 0; x < Nx; x++) {
          int t = 0;
          if (Communicator::ipe(3) == 0) {
            int site = m_index.site(x, y, z, t);
 
            Mat_SU_N wmat(Nc);
            wmat = m_wk * v.mat(site);
 
            Mat_SU_N cmat(Nc);
            cmat = wmat * c.mat_dag(site);
 
            wmat    = cmat * Unu.mat_dag(site);
            rect01 += ReTr(m_i_omega0 * wmat);
 
            wmat    = m_i_omega0 * v.mat(site);
            cmat    = wmat * c.mat_dag(site);
            wmat    = cmat * Unu.mat_dag(site);
            rect02 += ReTr(m_wk * wmat);
          }
 
          t = Nt - 1;
          if (Communicator::ipe(3) == NPEt - 1) {
            int site = m_index.site(x, y, z, t);
 
            Mat_SU_N cmat(Nc);
            cmat = m_i_omega0 * m_wkpr;
 
            Mat_SU_N wmat(Nc);
            wmat = cmat * v.mat_dag(site);
 
            cmat    = Unu.mat(site) * wmat;
            rectt1 += ReTr(cmat * U.mat_dag(site, mu));
 
            cmat    = m_i_omega0 * v.mat_dag(site);
            wmat    = m_wkpr * cmat;
            cmat    = Unu.mat(site) * wmat;
            rectt2 += ReTr(cmat * U.mat_dag(site, mu));
          }
        }
      }
    }
 
    // rect03
    //      +---+
    //      |   |
    //      v   ^
    //      |   |
    //  t=0 x--->
    //   omega0
 
    Field_G Cup1;
    upper(Cup1, U, mu, nu);
 
    m_shift.backward(v, Unu, mu);
    m_shift.backward(c, Cup1, nu);
 
    for (int z = 0; z < Nz; z++) {
      for (int y = 0; y < Ny; y++) {
        for (int x = 0; x < Nx; x++) {
          int t = 0;
          if (Communicator::ipe(3) == 0) {
            int site = m_index.site(x, y, z, t);
 
            Mat_SU_N wmat(Nc);
            wmat = c.mat(site) * v.mat_dag(site);
 
            Mat_SU_N cmat(Nc);
            cmat = Unu.mat(site) * wmat;
 
            wmat    = m_wk * cmat.dag();
            rect03 += ReTr(m_i_omega0 * wmat);
          }
        }
      }
    }
 
    // rectt3
    //   omega0
    // t=Nt x--->
    //      |   |
    //      ^   v
    //      |   |
    //      +---+
 
    Field_G Cdn1;
    lower(Cdn1, U, mu, nu);
 
    m_shift.backward(v, Unu, mu);
 
    for (int z = 0; z < Nz; z++) {
      for (int y = 0; y < Ny; y++) {
        for (int x = 0; x < Nx; x++) {
          int t = Nt - 1;
          if (Communicator::ipe(3) == NPEt - 1) {
            int site = m_index.site(x, y, z, t);
 
            Mat_SU_N wmat(Nc);
            wmat = m_i_omega0 * m_wkpr;
 
            Mat_SU_N cmat(Nc);
            cmat = wmat * v.mat_dag(site);
 
            wmat    = cmat * Cdn1.mat_dag(site);
            rectt3 += ReTr(Unu.mat(site) * wmat);
          }
        }
      }
    }
  }
  rect01 = Communicator::reduce_sum(rect01);
  rect02 = Communicator::reduce_sum(rect02);
  rect03 = Communicator::reduce_sum(rect03);
  rectt1 = Communicator::reduce_sum(rectt1);
  rectt2 = Communicator::reduce_sum(rectt2);
  rectt3 = Communicator::reduce_sum(rectt3);
 
  rect01 *= ctr / Lx;
  rect02 *= ctr / Lx;
  rect03 /= Lx;
  rectt1 *= ctr / Lx;
  rectt2 *= ctr / Lx;
  rectt3 /= Lx;
 
  double rect = -rect01 - rect02 - rect03 + rectt1 + rectt2 + rectt3;
 
  vout.general(m_vl, "SF_delSg_rect, at 01, 02, 03, at T1, T2, T3 = %.8f %.8f %.8f %.8f %.8f %.8f %.8f\n",
               rect, rect01, rect02, rect03, rectt1, rectt2, rectt3);
 
  return rect;
}
 
 
//====================================================================
 
double Staple_SF::plaquette(const Field_G& U)
{
  if (!m_initialized) {
    vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
  //  return (plaq_s(U) +plaq_t(U))/2;
  return(plaq_s(U) + plaq_t(U));
}
 
 
//====================================================================
 
double Staple_SF::plaquette_ct(const Field_G& U, const double ct)
{
  if (!m_initialized) {
    vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
  return(plaq_s(U) + plaq_t_ct(U, ct));
}
 
 
//====================================================================
 
double Staple_SF::plaq_s(const Field_G& U)
{
  const int Nvol = CommonParameters::Nvol();
 
  Field_G staple;
  double  plaq = 0.0;
 
  upper(staple, U, 0, 1);
  for (int site = 0; site < Nvol; site++) {
    plaq += ReTr(U.mat(site, 0) * staple.mat_dag(site));   // P_xy
  }
 
  upper(staple, U, 1, 2);
  for (int site = 0; site < Nvol; site++) {
    plaq += ReTr(U.mat(site, 1) * staple.mat_dag(site));   // P_yz
  }
 
  upper(staple, U, 2, 0);
  for (int site = 0; site < Nvol; site++) {
    plaq += ReTr(U.mat(site, 2) * staple.mat_dag(site));   // P_zx
  }
 
  plaq = Communicator::reduce_sum(plaq);
 
  return plaq;
}
 
 
//====================================================================
 
double Staple_SF::plaq_t(const Field_G& U)
{
  if (!m_initialized) {
    vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
 
  const int Ndim = CommonParameters::Ndim();
  const int Nvol = CommonParameters::Nvol();
 
  Field_G staple;
  double  plaq = 0.0;
 
  for (int nu = 0; nu < Ndim - 1; nu++) {
    lower(staple, U, 3, nu);
    //    staple = upper(U,3,nu);
    for (int site = 0; site < Nvol; site++) {
      plaq += ReTr(U.mat(site, 3) * staple.mat_dag(site));   // P_tk
    }
  }
 
  plaq = Communicator::reduce_sum(plaq);
 
  return plaq;
}
 
 
//====================================================================
 
double Staple_SF::plaq_t_ct(const Field_G& U, const double ct)
{
  if (!m_initialized) {
    vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
 
  const int Ndim = CommonParameters::Ndim();
  const int Nt   = CommonParameters::Nt();
  const int Nvol = CommonParameters::Nvol();
  const int NPEt = CommonParameters::NPEt();
 
  Field_G staple;
  double  plaq = 0.0;
 
  for (int nu = 0; nu < Ndim - 1; nu++) {
    lower(staple, U, 3, nu);
    // If the node is at the boundary the temporal plaquette is multiplied with ct.
    if (Communicator::ipe(3) == 0) {
      Field_SF::mult_ct_boundary(staple, 0, ct);
    }
    if (Communicator::ipe(3) == NPEt - 1) {
      Field_SF::mult_ct_boundary(staple, Nt - 1, ct);
    }
    for (int site = 0; site < Nvol; site++) {
      plaq += ReTr(U.mat(site, 3) * staple.mat_dag(site));   // P_tk
    }
  }
 
  plaq = Communicator::reduce_sum(plaq);
 
  return plaq;
}
 
 
//====================================================================
 
void Staple_SF::staple(Field_G& W, const Field_G& U, const int mu)
{
  if (!m_initialized) {
    vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
 
  const int Ndim = CommonParameters::Ndim();
 
  W.set(0.0);
 
  for (int nu = 0; nu < Ndim; nu++) {
    if (nu != mu) {
      Field_G c_tmp;
 
      upper(c_tmp, U, mu, nu);
      axpy(W, 1.0, c_tmp);
 
      lower(c_tmp, U, mu, nu);
      axpy(W, 1.0, c_tmp);
    }
  }
}
 
 
//====================================================================
 
void Staple_SF::staple_ct(Field_G& W, const Field_G& U, const int mu, const double ct)
{
  const int Ndim = CommonParameters::Ndim();
  const int Nt   = CommonParameters::Nt();
  const int NPEt = CommonParameters::NPEt();
 
  W.set(0.0);
 
  for (int nu = 0; nu < Ndim; nu++) {
    if (nu != mu) {
      Field_G staple_upper;
      Field_G staple_lower;
 
      upper(staple_upper, U, mu, nu);
      lower(staple_lower, U, mu, nu);
 
      if (Communicator::ipe(3) == 0) {
        if (mu == 3) {
          Field_SF::mult_ct_boundary(staple_upper, 0, ct);
          Field_SF::mult_ct_boundary(staple_lower, 0, ct);
        }
        if (nu == 3) {
          Field_SF::mult_ct_boundary(staple_lower, 1, ct);
        }
      }
 
      if (Communicator::ipe(3) == NPEt - 1) {
        if (mu == 3) {
          Field_SF::mult_ct_boundary(staple_upper, Nt - 1, ct);
          Field_SF::mult_ct_boundary(staple_lower, Nt - 1, ct);
        }
        if (nu == 3) {
          Field_SF::mult_ct_boundary(staple_upper, Nt - 1, ct);
        }
      }
 
      axpy(W, 1.0, staple_upper);
      axpy(W, 1.0, staple_lower);
    }
  }
}
 
 
//====================================================================
void Staple_SF::upper(Field_G& c, const Field_G& U, const int mu, const int nu)
{
  if (!m_initialized) {
    vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
 
  const int Nvol = CommonParameters::Nvol();
 
  // (1)  mu (2)
  //    +-->--+
  // nu |     |
  //   i+     +
  Field_G Umu;
  copy(Umu, 0, U, mu);
 
  Field_G Unu;
  copy(Unu, 0, U, nu);
  if (mu != 3) Field_SF::set_boundary_wk(Umu, m_wk);
  if (nu != 3) Field_SF::set_boundary_wk(Unu, m_wk);
 
  Field_G v;
  m_shift.backward(v, Unu, mu);
  m_shift.backward(c, Umu, nu);
  if (mu == 3) Field_SF::set_boundary_wkpr(v, m_wkpr);
  if (nu == 3) Field_SF::set_boundary_wkpr(c, m_wkpr);
 
  Field_G w;
  mult_Field_Gnd(w, 0, c, 0, v, 0);
  mult_Field_Gnn(c, 0, Unu, 0, w, 0);
  if (mu != 3) Field_SF::set_boundary_zero(c);
}
 
 
//====================================================================
void Staple_SF::lower(Field_G& c, const Field_G& U, const int mu, const int nu)
{
  if (!m_initialized) {
    vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n",
                 class_name.c_str());
    exit(EXIT_FAILURE);
  }
 
  const int Nvol = CommonParameters::Nvol();
 
  //    +     +
  // nu |     |
  //   i+-->--+
  //  (1)  mu (2)
  Field_G Umu;
  copy(Umu, 0, U, mu);
 
  Field_G Unu;
  copy(Unu, 0, U, nu);
  if (mu != 3) Field_SF::set_boundary_wk(Umu, m_wk);
  if (nu != 3) Field_SF::set_boundary_wk(Unu, m_wk);
 
  Field_G w;
  m_shift.backward(w, Unu, mu);
  if (mu == 3) Field_SF::set_boundary_wkpr(w, m_wkpr);
 
  Field_G v;
  mult_Field_Gnn(v, 0, Umu, 0, w, 0);
  mult_Field_Gdn(w, 0, Unu, 0, v, 0);
 
  m_shift.forward(c, w, nu);
  if (mu != 3) Field_SF::set_boundary_zero(c);
}
 
 
//====================================================================
 
void Staple_SF::print_plaquette(const Field_G& U)
{
  const int Lx = CommonParameters::Lx();
  const int Ly = CommonParameters::Ly();
  const int Lz = CommonParameters::Lz();
  const int Lt = CommonParameters::Lt();
 
  const double plaq  = plaquette(U);
  const double plaq2 = plaq + 3 * 3 * Lx * Ly * Lz;
 
  vout.general(m_vl, "plaq_SF without boundary spatial plaq = %.8f\n",
               plaq / (3 * Lx * Ly * Lz * (6 * Lt - 3)));
  vout.general(m_vl, "plaq_SF with boundary spatial plaq = %.8f\n",
               plaq2 / (3 * 6 * Lx * Ly * Lz * Lt));
}
 
 
//============================================================END=====