integrator_Omelyan.cpp

Go to the documentation of this file.

#include "integrator_Omelyan.h"

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

using Bridge::vout;

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_int("level", 0);
    param.Register_double("step_size", 0.0);
    param.Register_int("number_of_steps", 0);
    param.Register_double("lambda_Omelyan", 0.0);

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


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

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

//====================================================================
void Integrator_Omelyan::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    level;
  double Estep;
  int    Nstep;
  double lambda_Omelyan;

  int err = 0;
  err += params.fetch_int("level", level);
  err += params.fetch_double("step_size", Estep);
  err += params.fetch_int("number_of_steps", Nstep);
  err += params.fetch_double("lambda_Omelyan", lambda_Omelyan);

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


  set_parameters(level, Estep, Nstep, lambda_Omelyan);
}


//====================================================================
void Integrator_Omelyan::set_parameters(int level, double Estep, int Nstep,
                                        double lambda_Omelyan)
{
  //- print input parameters
  vout.general(m_vl, "Integrator (Omelyan) setup:\n");
  vout.general(m_vl, "  Level: %4d\n", m_level);
  vout.general(m_vl, "  Estep    = %10.6f\n", m_Estep);
  vout.general(m_vl, "  Nstep    = %4d\n", m_Nstep);
  vout.general(m_vl, "  Number of actions: %4d\n", m_action.size());
  vout.general(m_vl, "  lambda_Omelyan = %10.6f\n", m_lambda);

  //- range check
  // NB. level,Estep,Nstep,lambda_Omelyan == 0 is allowed.

  //- store values
  m_Estep  = Estep;
  m_Nstep  = Nstep;
  m_level  = level;
  m_lambda = lambda_Omelyan;
}


//====================================================================
void Integrator_Omelyan::evolve(Field_G& iP, Field_G& U)
{
  int Nin  = U.nin();
  int Nvol = U.nvol();
  int Nex  = U.nex();
  int Nc   = CommonParameters::Nc();

  double estep = m_Estep;
  //  double esteph = 0.5 * estep;
  double estep2;

  Field force(Nin, Nvol, Nex), force1(Nin, Nvol, Nex);


  vout.general(m_vl, "Integration level-%d start.\n", m_level);

  // Initial half step of update of iP
  estep2 = estep * m_lambda;
  if (m_Nstep > 0) {
    int istep = 0;
    vout.general(m_vl, "istep = %d\n", istep);
    force = 0.0;
    for (int i = 0; i < m_action.size(); ++i) {
      force1 = m_action[i]->force();
      force += estep2 * force1;
    }
    iP += (Field_G)force;
  }

  // Molecular dynamics step
  for (int istep = 1; istep < m_Nstep + 1; istep++) {
    vout.general(m_vl, "istep = %d\n", istep);

    m_integ_next->evolve(iP, U);

    estep2 = estep * (1.0 - 2.0 * m_lambda);
    force  = 0.0;
    for (int i = 0; i < m_action.size(); ++i) {
      force1 = m_action[i]->force();
      force += estep2 * force1;
    }
    iP += (Field_G)force;

    m_integ_next->evolve(iP, U);

    estep2 = estep * 2.0 * m_lambda;
    if (istep == m_Nstep) estep2 = estep * m_lambda;

    force = 0.0;
    for (int i = 0; i < m_action.size(); ++i) {
      force1 = m_action[i]->force();
      force += estep2 * force1;
    }
    iP += (Field_G)force;
  }  // here istep loop ends

  vout.general(m_vl, "Integration level-%d finished.\n", m_level);
}


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