action_F_Overlap_Nf2.cpp

Go to the documentation of this file.

#include "action_F_Overlap_Nf2.h"

//====================================================================
void Action_F_Overlap_Nf2::set_parameters(const Parameters& params)
{
  Bridge::VerboseLevel vl = params.get_VerboseLevel();

  set_parameter_verboselevel(vl);
}


//====================================================================
void Action_F_Overlap_Nf2::set_parameters()
{
  // ####  parameter setup  ####
  int Ndim = CommonParameters::Ndim();
  int Nc   = CommonParameters::Nc();
  int Nd   = CommonParameters::Nd();
  int Nvol = CommonParameters::Nvol();
  int NinG = 2 * Nc * Nc;

  m_index = new Index_lex;

  // setting Wilson fermion operator
  m_fopr_w = new Fopr_Wilson;
  m_M0     = 1.6;
  m_kappa  = 0.5 / (4.0 - m_M0);

  m_boundary.resize(Ndim);
  m_boundary[0] = 1;
  m_boundary[1] = 1;
  m_boundary[2] = 1;
  m_boundary[3] = 1;

  m_fopr_w->set_parameters(m_kappa, m_boundary);
  m_fopr_w->set_config(m_U);


  // setting overlap fermion operator
  m_fopr_ov      = new Fopr_Overlap(m_fopr_w);
  m_mq           = 0.2;
  m_Np           = 16;
  m_x_min        = 0.01;
  m_x_max        = 8.0;
  m_Niter_ms     = 1000;
  m_Stop_cond_ms = 1.0e-20;

  m_fopr_ov->set_parameters(m_M0, m_mq, m_Np, m_x_min, m_x_max,
                            m_Niter_ms, m_Stop_cond_ms,
                            m_boundary);

  // Zolotarev coefficients
  m_sigma.resize(m_Np);
  m_cl.resize(2 * m_Np);
  m_bl.resize(m_Np);

  // low-modes
  int Nmm = 100;
  m_TDa.resize(Nmm);
  m_vk.resize(Nmm);

  m_Nsbt = 0;

  int NinF = 2 * Nc * Nd;
  for (int k = 0; k < Nmm; ++k) {
    m_vk[k].reset(NinF, Nvol, 1);
  }

  m_sort  = new SortField_AbsLow;
  m_eigen = new Eigensolver_IRLanczos(m_fopr_w, m_sort);
  int    Nk          = 20;
  int    Np          = 50;
  double Enorm_eigen = 1.e-22;
  int    Niter_eigen = 500;
  double Vthreshold  = 0.15;
  m_eigen->set_parameters(Nk, Np, Enorm_eigen, Niter_eigen, Vthreshold);

  // force saved
  m_force.reset(NinG, Nvol, Ndim);

  m_status_linkv = 0;
}


//====================================================================
double Action_F_Overlap_Nf2::langevin(RandomNumbers *rand)
{
  int Nc       = CommonParameters::Nc();
  int Nd       = CommonParameters::Nd();
  int Lvol     = CommonParameters::Lvol();
  int size_psf = Lvol * Nc * Nd;

  Field_F xi;

  rand->gauss_lex_global(xi);

  //  m_psf = (Field_F) m_fopr_ov->Ddag(xi);
  //m_psf = (Field_F) m_fopr_ov->H(xi);

  m_fopr_ov->set_mode("H");
  m_psf = (Field_F)m_fopr_ov->mult(xi);

  double xi2   = xi.norm();
  double H_psf = xi2 * xi2;

  vout.general(m_vl, "H_Foverlap = %18.8f\n", H_psf);
  vout.general(m_vl, "H_F/dof    = %18.8f\n", H_psf / size_psf);

  m_status_linkv = 0;

  return H_psf;
}


//====================================================================
double Action_F_Overlap_Nf2::calcH()
{
  int Nc       = CommonParameters::Nc();
  int Nd       = CommonParameters::Nd();
  int Lvol     = CommonParameters::Lvol();
  int size_psf = Lvol * Nc * Nd;

  int    Niter     = 500;
  double Stop_cond = 1.0e-24;

  int    Nconv;
  double diff;

  Field_F v1;

  Solver *solver = new Solver_CG(m_fopr_ov);

  solver->set_parameters(Niter, Stop_cond);

  m_fopr_ov->set_mode("DdagD");
  solver->solve(v1, m_psf, Nconv, diff);

  vout.general(m_vl, "  Nconv = %d\n", Nconv);
  vout.general(m_vl, "  diff  = %.8e\n", diff);

  double H_psf = v1 * m_psf;

  vout.general(m_vl, "H_F_overlap = %18.8f\n", H_psf);
  vout.general(m_vl, "H_F/dof     = %18.8f\n", H_psf / size_psf);

  delete solver;

  return H_psf;
}


//====================================================================
const Field Action_F_Overlap_Nf2::force()
{
  if (m_status_linkv == 0) {
    int Nc       = CommonParameters::Nc();
    int Nd       = CommonParameters::Nd();
    int Nvol     = CommonParameters::Nvol();
    int Lvol     = CommonParameters::Lvol();
    int Ndim     = CommonParameters::Ndim();
    int NinG     = 2 * Nc * Nc;
    int size_psf = Lvol * Nc * Nd;

    //- solving psi = (H^2)^(-1) phi
    Field_F psi(Nvol, 1);

    int    Niter     = 500;
    double Stop_cond = 1.0e-24;

    int    Nconv;
    double diff;

    Solver *solver = new Solver_CG(m_fopr_ov);
    solver->set_parameters(Niter, Stop_cond);

    m_fopr_ov->set_mode("DdagD");
    solver->solve(psi, m_psf, Nconv, diff);

    vout.general(m_vl, "  Solver(ov):  Nconv = %d  diff  = %.8e\n", Nconv, diff);
    delete solver;

    Force_F_Overlap_Nf2 force_ov;
    force_ov.set_parameters(m_M0, m_mq,
                            m_Np, m_x_min, m_x_max,
                            m_Niter_ms, m_Stop_cond_ms,
                            m_boundary);
    force_ov.set_config(m_U);

    Field force(NinG, Nvol, Ndim);
    force = force_ov.force_core(psi);

    m_force = force;
    ++m_status_linkv;

    return m_force;
  } else {
    vout.general(m_vl, "  F_overlap_Nf2 returns previous force.\n");
    return m_force;
  }
}


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