action_F_Rational.cpp

Go to the documentation of this file.

#include "action_F_Rational.h"

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

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_int("Np", 0);
    param.Register_int("n_exp", 0);
    param.Register_int("d_exp", 0);
    param.Register_double("x_min", 0.0);
    param.Register_double("x_max", 0.0);
    param.Register_int("Niter", 0);
    param.Register_double("Stop_cond", 0.0);

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


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

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

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

//====================================================================
void Action_F_Rational::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    Np, n_exp, d_exp;
  double x_min, x_max;
  int    Niter;
  double Stop_cond;

  int err = 0;
  err += params.fetch_int("Np", Np);
  err += params.fetch_int("n_exp", n_exp);
  err += params.fetch_int("d_exp", d_exp);
  err += params.fetch_double("x_min", x_min);
  err += params.fetch_double("x_max", x_max);
  err += params.fetch_int("Niter", Niter);
  err += params.fetch_double("Stop_cond", Stop_cond);

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


  set_parameters(Np, n_exp, d_exp, x_min, x_max, Niter, Stop_cond);
}


//====================================================================
void Action_F_Rational::set_parameters(int Np, int n_exp, int d_exp,
                                       double x_min, double x_max,
                                       int Niter, double Stop_cond)
{
  //- print input parameters
  vout.general(m_vl, "%s:\n", class_name.c_str());
  vout.general(m_vl, "  Np        = %4d\n", Np);
  vout.general(m_vl, "  n_exp     = %4d\n", n_exp);
  vout.general(m_vl, "  d_exp     = %4d\n", d_exp);
  vout.general(m_vl, "  x_min     = %8.4f\n", x_min);
  vout.general(m_vl, "  x_max     = %8.4f\n", x_max);
  vout.general(m_vl, "  Niter     = %6d\n", Niter);
  vout.general(m_vl, "  Stop_cond = %8.2e\n", Stop_cond);

  //- range check
  int err = 0;
  err += ParameterCheck::non_zero(Np);
  err += ParameterCheck::non_zero(n_exp);
  err += ParameterCheck::non_zero(d_exp);
  // NB. x_min,x_max == 0 is allowed.
  err += ParameterCheck::non_negative(Niter);
  err += ParameterCheck::square_non_zero(Stop_cond);

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

  //- store values
  m_Np        = Np;
  m_n_exp     = n_exp;
  m_d_exp     = d_exp;
  m_x_min     = x_min;
  m_x_max     = x_max;
  m_Niter     = Niter;
  m_Stop_cond = Stop_cond;

  //- post-process
  setup();
}


//====================================================================
void Action_F_Rational::setup()
{
  int Nc   = CommonParameters::Nc();
  int Nvol = CommonParameters::Nvol();
  int Ndim = CommonParameters::Ndim();
  int NinG = 2 * Nc * Nc;

  m_force.reset(NinG, Nvol, Ndim);

  // rational operator for Langevin step.
  m_fopr_langev = new Fopr_Rational(m_fopr);
  int d_exp2 = 2 * m_d_exp;
  m_fopr_langev->set_parameters(m_Np, m_n_exp, d_exp2,
                                m_x_min, m_x_max, m_Niter, m_Stop_cond);

  // rational operator for Hamiltonian calculation.
  m_fopr_H = new Fopr_Rational(m_fopr);
  int n_expm = -m_n_exp;
  m_fopr_H->set_parameters(m_Np, n_expm, m_d_exp,
                           m_x_min, m_x_max, m_Niter, m_Stop_cond);

  // force of rational operator
  m_force_rational = new Force_F_Rational(m_fopr, m_fopr_force);
  m_force_rational->set_parameters(m_Np, n_expm, m_d_exp,
                                   m_x_min, m_x_max, m_Niter, m_Stop_cond);

  // link variable update flag
  m_status_linkv = 0;
}


//====================================================================
void Action_F_Rational::tidyup()
{
  delete m_force_rational;
  delete m_fopr_H;
  delete m_fopr_langev;
}


//====================================================================
double Action_F_Rational::langevin(RandomNumbers *rand)
{
  int NinF     = m_fopr->field_nin();
  int NvolF    = m_fopr->field_nvol();
  int NexF     = m_fopr->field_nex();
  int size_psf = NinF * NvolF * NexF * CommonParameters::NPE();

  m_psf.reset(NinF, NvolF, NexF);

  vout.general(m_vl, "  %s: %s\n", class_name.c_str(), m_label.c_str());

  Field xi(NinF, NvolF, NexF);
  rand->gauss_lex_global(xi);

  m_fopr_langev->set_config(m_U);
  m_psf = m_fopr_langev->mult(xi);

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

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

  return H_psf;
}


//====================================================================
double Action_F_Rational::calcH()
{
  int NinF     = m_fopr->field_nin();
  int NvolF    = m_fopr->field_nvol();
  int NexF     = m_fopr->field_nex();
  int size_psf = NinF * NvolF * NexF * CommonParameters::NPE();

  Field v1(NinF, NvolF, NexF);

  vout.general(m_vl, "  %s: %s\n", class_name.c_str(), m_label.c_str());

  m_fopr_H->set_config(m_U);
  v1 = m_fopr_H->mult(m_psf);

  double H_psf = v1 * m_psf;

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

  return H_psf;
}


//====================================================================
const Field Action_F_Rational::force()
{
  if (m_status_linkv == 0) {
    int Nvol = m_U->nvol();
    int Nex  = m_U->nex();

    Field_G force(Nvol, Nex), force1(Nvol, Nex);

    vout.general(m_vl, "  %s: %s\n", class_name.c_str(), m_label.c_str());

    m_force_rational->set_config(m_U);
    force = m_force_rational->force_core(m_psf);

    m_force = (Field)force;
    ++m_status_linkv;

    double Fave, Fmax, Fdev;
    m_force.stat(Fave, Fmax, Fdev);
    vout.general(m_vl, "    Frational_ave = %12.6f  Frational_max = %12.6f  Frational_dev = %12.6f\n",
                 Fave, Fmax, Fdev);

    return m_force;
  } else {
    vout.general(m_vl, "  %s returns previous force.\n", class_name.c_str());
    return m_force;
  }
}


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