integrator_UpdateU.cpp

Go to the documentation of this file.

#include "integrator_UpdateU.h"

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

//- parameter entry
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_double("step_size", 0.0);
    param.Register_int("order_of_exp_iP", 0);

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


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

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

//====================================================================
void Integrator_UpdateU::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
  double Estep;
  int    Nprec;

  int err = 0;
  err += params.fetch_double("step_size", Estep);
  err += params.fetch_int("order_of_exp_iP", Nprec);

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


  set_parameters(Estep, Nprec);
}


//====================================================================
void Integrator_UpdateU::set_parameters(double Estep, int Nprec)
{
  //- print input parameters
  vout.general(m_vl, "Integrator (updateU) setup:\n");
  vout.general(m_vl, "  Estep    = %10.6f\n", m_Estep);
  vout.general(m_vl, "  Nprec     = %d\n", m_Nprec);

  //- range check
  int err = 0;
  // NB. Estep == 0 is allowed.
  err += ParameterCheck::non_zero(Nprec);

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

  //- store values
  m_Estep = Estep;
  m_Nprec = Nprec;
}


//====================================================================
void Integrator_UpdateU::evolve(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 = m_Estep / (m_Nprec - iprec);
        u2  = h1 * u1;
        u2 *= exf;
        u1  = u2;
        u1 += u0;
      }

      u1.reunit();
      U.set_mat(site, ex, u1);
    }
  }

  for (int i = 0; i < m_action.size(); ++i) {
    m_action[i]->notify_linkv();
  }

  for (int i = 0; i < m_director.size(); ++i) {
    m_director[i]->notify_linkv();
  }
}


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