source_4spinor_Exp.cpp

Go to the documentation of this file.

#include "source_4spinor_Exp.h"

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

using std::valarray;

//- parameter entry
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_int_vector("source_position", std::valarray<int>());
    param.Register_double("slope", 0.0);
    param.Register_double("power", 0.0);

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


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

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

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

//====================================================================
void Source_4spinor_Exp::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
  valarray<int> source_position;
  double        slope, power;

  int err = 0;
  err += params.fetch_int_vector("source_position", source_position);
  err += params.fetch_double("slope", slope);
  err += params.fetch_double("power", power);

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


  set_parameters(source_position, slope, power);
}


//====================================================================
void Source_4spinor_Exp::set_parameters(valarray<int>& source_position,
                                        double slope, double power)
{
  int Ndim = CommonParameters::Ndim();

  valarray<int> Lsize(Ndim);
  Lsize[0] = CommonParameters::Lx();
  Lsize[1] = CommonParameters::Ly();
  Lsize[2] = CommonParameters::Lz();
  Lsize[3] = CommonParameters::Lt();

  //- print input parameters
  vout.general(m_vl, "Source for 4-spinor field - exponential smeared:\n");
  for (int mu = 0; mu < Ndim; ++mu) {
    vout.general(m_vl, "  source_position[%d] = %d\n",
                 mu, source_position[mu]);
  }
  vout.general(m_vl, "  slope = %12.6f\n", slope);
  vout.general(m_vl, "  power = %12.6f\n", power);

  //- range check
  int err = 0;
  for (int mu = 0; mu < Ndim; ++mu) {
    // NB. Lsize[mu] > abs(source_position[mu])
    err += ParameterCheck::non_negative(Lsize[mu] - abs(source_position[mu]));
  }
  // NB. slope,power == 0 is allowed.

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

  assert(source_position.size() == Ndim);

  //- store values
  m_source_position.resize(Ndim);
  for (int mu = 0; mu < Ndim; ++mu) {
    m_source_position[mu] = (source_position[mu] + Lsize[mu]) % Lsize[mu];
  }

  m_slope = slope;
  m_power = power;


  //- post-process
  int Lvol3 = Lsize[0] * Lsize[1] * Lsize[2];
  m_src_func.reset(1, Lvol3, 1);
  m_src_func = 0.0;

  for (int z = -Lsize[2] + 1; z < Lsize[2]; ++z) {
    for (int y = -Lsize[1] + 1; y < Lsize[1]; ++y) {
      for (int x = -Lsize[0] + 1; x < Lsize[0]; ++x) {
        int    z2    = (m_source_position[2] + z + Lsize[2]) % Lsize[2];
        int    y2    = (m_source_position[1] + y + Lsize[1]) % Lsize[1];
        int    x2    = (m_source_position[0] + x + Lsize[0]) % Lsize[0];
        int    sites = x2 + Lsize[0] * (y2 + Lsize[1] * z2);
        double r     = sqrt((double)(x * x + y * y + z * z));
        double ex    = pow(r, m_power);
        double expf  = exp(-m_slope * ex);
        m_src_func.add(0, sites, 0, expf);
      }
    }
  }


  double Fnorm = 0.0;
  for (int i = 0; i < Lvol3; i++) {
    Fnorm += m_src_func.cmp(i) * m_src_func.cmp(i);
  }
  m_src_func *= 1.0 / sqrt(Fnorm);


  double epsilon_criterion = CommonParameters::epsilon_criterion();

  Fnorm = 0.0;
  for (int i = 0; i < Lvol3; i++) {
    Fnorm += m_src_func.cmp(i) * m_src_func.cmp(i);
  }

  if (!(abs(sqrt(Fnorm) - 1.0) < epsilon_criterion)) {
    vout.crucial(m_vl, "%s: norm check failed.\n", class_name.c_str());
    abort();
  }
}


//====================================================================
void Source_4spinor_Exp::set(Field_F& src, int ic, int id)
{
  int Nvol = CommonParameters::Nvol();
  int Ndim = CommonParameters::Ndim();

  valarray<int> Lsize(Ndim);
  Lsize[0] = CommonParameters::Lx();
  Lsize[1] = CommonParameters::Ly();
  Lsize[2] = CommonParameters::Lz();
  Lsize[3] = CommonParameters::Lt();

  valarray<int> Nsize(Ndim);
  Nsize[0] = CommonParameters::Nx();
  Nsize[1] = CommonParameters::Ny();
  Nsize[2] = CommonParameters::Nz();
  Nsize[3] = CommonParameters::Nt();

  int IPEx = Communicator::ipe(0);
  int IPEy = Communicator::ipe(1);
  int IPEz = Communicator::ipe(2);

  assert(ic < CommonParameters::Nc());
  assert(id < CommonParameters::Nd());
  assert(src.nvol() == Nvol);
  assert(src.nex() == 1);

  src = 0.0;

  valarray<int> site_src(Ndim);
  valarray<int> node_src(Ndim);
  for (int mu = 0; mu < Ndim; ++mu) {
    site_src[mu] = m_source_position[mu] % Nsize[mu];
    node_src[mu] = m_source_position[mu] / Nsize[mu];
  }

  if (node_src[3] == Communicator::ipe(3)) {
    int t = site_src[3];

    for (int z = 0; z < Nsize[2]; ++z) {
      int z2 = z + Nsize[2] * IPEz;
      for (int y = 0; y < Nsize[1]; ++y) {
        int y2 = y + Nsize[1] * IPEy;
        for (int x = 0; x < Nsize[0]; ++x) {
          int x2    = x + Nsize[0] * IPEx;
          int sites = x2 + Lsize[0] * (y2 + Lsize[1] * z2);
          int site  = m_index.site(x, y, z, t);
          src.set_ri(ic, id, site, 0, m_src_func.cmp(0, sites, 0), 0.0);
        }
      }
    }
  }
}


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