integrator_Leapfrog.cpp

Go to the documentation of this file.

#include "integrator_Leapfrog.h"

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

//- 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_string("verbose_level", "NULL");
  }


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

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

//====================================================================
void Integrator_Leapfrog::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;

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

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


  set_parameters(level, Estep, Nstep);
}


//====================================================================
void Integrator_Leapfrog::set_parameters(int level, double Estep, int Nstep)
{
  //- print input parameters
  vout.general(m_vl, "Integrator (leapfrog) 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());

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

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


//====================================================================
void Integrator_Leapfrog::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
  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 += esteph * 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;
    if (istep == m_Nstep) estep2 = esteph;

    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=====