source_MomentumWall.cpp

Go to the documentation of this file.

#include "source_MomentumWall.h"

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

using std::valarray;

#ifdef USE_FACTORY
namespace {
  Source *create_object()
  {
    return new Source_MomentumWall();
  }


  bool init = Source::Factory::Register("MomentumWall", create_object);
}
#endif

//- parameter entry
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_int_vector("source_position", valarray<int>());
    param.Register_double_vector("source_momentum", valarray<double>());

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


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

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

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

//====================================================================
void Source_MomentumWall::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;
  valarray<double> source_momentum;

  int err = 0;
  err += params.fetch_int_vector("source_position", source_position);
  err += params.fetch_double_vector("source_momentum", source_momentum);

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

  set_parameters(source_position, source_momentum);
}


//====================================================================
void Source_MomentumWall::set_parameters(const valarray<int>& source_position,
                                         const valarray<double>& source_momentum)
{
  // ####  parameter setup  ####
  int Ndim = CommonParameters::Ndim();

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

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

  const int t_dir = Ndim - 1;

  //- print input parameters
  vout.general(m_vl, "Source for spinor field - momentum wall smeared:\n");
  for (int mu = 0; mu < Ndim; ++mu) {
    vout.general(m_vl, "  source_position[%d] = %d\n",
                 mu, source_position[mu]);
  }
  for (int mu = 0; mu < Ndim - 1; ++mu) {
    vout.general(m_vl, "  source_momentum[%d] = %g\n",
                 mu, source_momentum[mu]);
  }


  //- 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]));
  }

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

  //- 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_source_momentum = source_momentum;

  //- post-process

  //- PE location in t-direction.
  int tpe = m_source_position[t_dir] / Nsize[t_dir];

  m_in_node = false;

  if (tpe == Communicator::ipe(t_dir)) {
    m_in_node = true;
  }
}


//====================================================================
void Source_MomentumWall::set(Field& src, int j)
{
  int Ndim = CommonParameters::Ndim();

  //- global lattice size
  int Lx = CommonParameters::Lx();
  int Ly = CommonParameters::Ly();
  int Lz = CommonParameters::Lz();

  const int Lvol3 = Lx * Ly * Lz;

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

  //- clear field
  src = 0.0;

  if (m_in_node) {
    int t = m_source_position[3] % Nsize[3];

    for (int z = 0; z < Nsize[2]; ++z) {
      for (int y = 0; y < Nsize[1]; ++y) {
        for (int x = 0; x < Nsize[0]; ++x) {
          int isite = m_index.site(x, y, z, t);

          int gX = x + Nsize[0] * Communicator::ipe(0);
          int gY = y + Nsize[1] * Communicator::ipe(1);
          int gZ = z + Nsize[2] * Communicator::ipe(2);

          double theta
            = double(gX - m_source_position[0]) * m_source_momentum[0]
              + double(gY - m_source_position[1]) * m_source_momentum[1]
              + double(gZ - m_source_position[2]) * m_source_momentum[2];

          //XXX field layout: complex as two doubles
          src.set(2 * j + 0, isite, 0, cos(theta) / Lvol3);
          src.set(2 * j + 1, isite, 0, sin(theta) / Lvol3);
        }
      }
    }
  }
}


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