builder_Integrator.cpp

Go to the documentation of this file.

#include "builder_Integrator.h"

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

using std::valarray;
using std::string;

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_string("integrator", "NULL");
    param.Register_double("step_size", 0.0);
    param.Register_int("number_of_levels", 0);
    param.Register_int_vector("number_of_actions", valarray<int>());
    param.Register_int_vector("number_of_steps", valarray<int>());
    param.Register_int("order_of_exp_iP", 0);
    param.Register_double("lambda_Omelyan", 0.0);

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


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

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

//====================================================================
void Builder_Integrator::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
  string        str_integrator_type;
  double        Estep;
  int           Nlevel;
  valarray<int> Naction;
  valarray<int> Nstep;
  int           Nprec;
  double        lambda_Omelyan;

  int err = 0;
  err += params.fetch_string("integrator", str_integrator_type);
  err += params.fetch_double("step_size", Estep);
  err += params.fetch_int("number_of_levels", Nlevel);
  err += params.fetch_int_vector("number_of_actions", Naction);
  err += params.fetch_int_vector("number_of_steps", Nstep);
  err += params.fetch_int("order_of_exp_iP", Nprec);
  err += params.fetch_double("lambda_Omelyan", lambda_Omelyan);

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


  set_parameters(str_integrator_type, Estep, Nlevel, Naction, Nstep, Nprec, lambda_Omelyan);
}


//====================================================================
void Builder_Integrator::set_parameters(string str_integrator_type,
                                        double Estep, int Nlevel,
                                        const valarray<int>& Naction,
                                        const valarray<int>& Nstep,
                                        const int Nprec,
                                        const double lambda_Omelyan)
{
  //- print input parameters
  vout.general(m_vl, "MD integrator builder:\n");
  vout.general(m_vl, "  Integrator = %s\n", str_integrator_type.c_str());
  vout.general(m_vl, "  Estep      = %10.6f\n", Estep);
  vout.general(m_vl, "  Nlevel     = %2d\n", Nlevel);
  for (int lv = 0; lv < Nlevel; ++lv) {
    vout.general(m_vl, "   level     = %2d:  Naction = %2d  Nstep = %4d\n",
                 lv, Naction[lv], Nstep[lv]);
  }
  vout.general(m_vl, "  Nprec      = %4d\n", Nprec);

  //- range check
  int err = 0;
  err += ParameterCheck::non_NULL(str_integrator_type);
  // NB. Estep == 0 is allowed.
  err += ParameterCheck::non_zero(Nlevel);
  for (int lv = 0; lv < Nlevel; ++lv) {
    err += ParameterCheck::non_zero(Naction[lv]);
  }
  // NB. Nstep == 0 is allowed.
  err += ParameterCheck::non_zero(Nprec);
  // NB. lambda_Omelyan == 0 is allowed.

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

  //- store values
  m_str_integrator_type = str_integrator_type;

  m_Estep  = Estep;
  m_Nlevel = Nlevel;

  m_Naction.resize(Nlevel);
  m_Nstep.resize(Nlevel);

  for (int lv = 0; lv < m_Nlevel; ++lv) {
    m_Naction[lv] = Naction[lv];
    m_Nstep[lv]   = Nstep[lv];
  }

  m_Nprec          = Nprec;
  m_lambda_Omelyan = lambda_Omelyan;
}


//====================================================================
Integrator *Builder_Integrator::build()
{
  //- select leapfrog or omelyan
  Integrator *integrator;

  if (m_str_integrator_type == "Leapfrog") {
    integrator = build_leapfrog();
  } else if (m_str_integrator_type == "Omelyan") {
    integrator = build_omelyan();
  } else {
    vout.crucial("Builder_Integrator::build : unsupported smear type \"%s\".\n",
                 m_str_integrator_type.c_str());
    abort();
  }

  return integrator;
}


//====================================================================
Integrator *Builder_Integrator::build_leapfrog()
{
  // #####  Following setup is for PUP order.
  m_integs.resize(m_Nlevel + 1);

  double estep_lowest = m_Estep;
  for (int lv = 1; lv < m_Nlevel; ++lv) {
    estep_lowest *= 1.0 / ((double)m_Nstep[lv]);
  }
  Integrator_UpdateU *updateU
    = new Integrator_UpdateU(m_action, m_director);
  updateU->set_parameters(estep_lowest, m_Nprec);
  m_integs[m_Nlevel] = (Integrator *)updateU;

  double estep   = estep_lowest;
  int    jaction = m_action.size();

  for (int lv = m_Nlevel - 1; lv >= 0; --lv) {
    valarray<Action *> actions(m_Naction[lv]);

    jaction -= m_Naction[lv];
    for (int lv2 = 0; lv2 < m_Naction[lv]; ++lv2) {
      actions[lv2] = m_action[jaction + lv2];
    }

    Integrator_Leapfrog *leapfrog
      = new Integrator_Leapfrog(actions, m_integs[lv + 1]);
    leapfrog->set_parameters(lv, estep, m_Nstep[lv]);
    m_integs[lv] = (Integrator *)leapfrog;

    estep *= ((double)m_Nstep[lv]);
  }

  return (Integrator *)m_integs[0];
}


//====================================================================
Integrator *Builder_Integrator::build_omelyan()
{
  // #####  Following setup is for PUP order.
  m_integs.resize(m_Nlevel + 1);

  double estep_lowest = m_Estep;
  for (int lv = 1; lv < m_Nlevel; ++lv) {
    estep_lowest *= 0.5 / ((double)m_Nstep[lv]);
  }
  Integrator_UpdateU *updateU
    = new Integrator_UpdateU(m_action, m_director);
  updateU->set_parameters(0.5 * estep_lowest, m_Nprec);
  m_integs[m_Nlevel] = (Integrator *)updateU;

  double estep   = estep_lowest;
  int    jaction = m_action.size();

  for (int lv = m_Nlevel - 1; lv >= 0; --lv) {
    valarray<Action *> actions(m_Naction[lv]);

    jaction -= m_Naction[lv];
    for (int lv2 = 0; lv2 < m_Naction[lv]; ++lv2) {
      actions[lv2] = m_action[jaction + lv2];
    }

    Integrator_Omelyan *omelyan
      = new Integrator_Omelyan(actions, m_integs[lv + 1]);
    omelyan->set_parameters(lv, estep, m_Nstep[lv], m_lambda_Omelyan);
    m_integs[lv] = (Integrator *)omelyan;

    estep *= 2.0 * ((double)m_Nstep[lv]);
  }

  return (Integrator *)m_integs[0];
}


//====================================================================
void Builder_Integrator::tidyup()
{
  for (int lv = 0; lv < m_Nlevel + 1; ++lv) {
    delete m_integs[lv];
  }
}


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