Bridge++  Version 1.5.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
source_MomentumWall.cpp
Go to the documentation of this file.
1 
14 #include "source_MomentumWall.h"
15 
16 #ifdef USE_FACTORY_AUTOREGISTER
17 namespace {
18  bool init = Source_MomentumWall::register_factory();
19 }
20 #endif
21 
22 const std::string Source_MomentumWall::class_name = "Source_MomentumWall";
23 
24 //====================================================================
26 {
27  const string str_vlevel = params.get_string("verbose_level");
28 
29  m_vl = vout.set_verbose_level(str_vlevel);
30 
31  //- fetch and check input parameters
32  std::vector<int> source_position;
33  std::vector<int> source_momentum;
34 
35  int err = 0;
36  err += params.fetch_int_vector("source_position", source_position);
37  err += params.fetch_int_vector("source_momentum", source_momentum);
38 
39  if (err) {
40  vout.crucial(m_vl, "Error at %s: input parameter not found.\n", class_name.c_str());
41  exit(EXIT_FAILURE);
42  }
43 
44  set_parameters(source_position, source_momentum);
45 }
46 
47 
48 //====================================================================
49 void Source_MomentumWall::set_parameters(const std::vector<int>& source_position,
50  const std::vector<int>& source_momentum)
51 {
52  // #### parameter setup ####
53  const int Ndim = CommonParameters::Ndim();
54  const int t_dir = Ndim - 1;
55 
56  //- global lattice size
57  std::vector<int> Lsize(Ndim);
58 
59  Lsize[0] = CommonParameters::Lx();
60  Lsize[1] = CommonParameters::Ly();
61  Lsize[2] = CommonParameters::Lz();
62  Lsize[3] = CommonParameters::Lt();
63 
64  //- local size
65  const int Nt = CommonParameters::Nt();
66 
67  //- print input parameters
68  vout.general(m_vl, "Source for spinor field - momentum wall smeared:\n");
69  for (int mu = 0; mu < Ndim; ++mu) {
70  vout.general(m_vl, " source_position[%d] = %d\n",
71  mu, source_position[mu]);
72  }
73  for (int mu = 0; mu < Ndim - 1; ++mu) {
74  vout.general(m_vl, " source_momentum[%d] = %d\n",
75  mu, source_momentum[mu]);
76  }
77 
78 
79  //- range check
80  int err = 0;
81  for (int mu = 0; mu < Ndim; ++mu) {
82  // NB. Lsize[mu] > abs(source_position[mu])
83  err += ParameterCheck::non_negative(Lsize[mu] - abs(source_position[mu]));
84  }
85 
86  if (err) {
87  vout.crucial(m_vl, "Error at %s: parameter range check failed.\n", class_name.c_str());
88  exit(EXIT_FAILURE);
89  }
90 
91  //- store values
92  m_source_position.resize(Ndim);
93  for (int mu = 0; mu < Ndim; ++mu) {
94  m_source_position[mu] = (source_position[mu] + Lsize[mu]) % Lsize[mu];
95  }
96 
97  m_source_momentum = source_momentum;
98 
99  //- post-process
100 
101  //- PE location in t-direction.
102  const int tpe = m_source_position[t_dir] / Nt;
103 
104  m_in_node = false;
105 
106  if (tpe == Communicator::ipe(t_dir)) {
107  m_in_node = true;
108  }
109 }
110 
111 
112 //====================================================================
113 void Source_MomentumWall::set(Field& src, const int idx)
114 {
115  const int Ndim = CommonParameters::Ndim();
116 
117  const int x_dir = 0;
118  const int y_dir = 1;
119  const int z_dir = 2;
120  const int t_dir = Ndim - 1;
121 
122  //- global lattice size
123  const int Lx = CommonParameters::Lx();
124  const int Ly = CommonParameters::Ly();
125  const int Lz = CommonParameters::Lz();
126  const int Lvol3 = Lx * Ly * Lz;
127 
128  //- local size
129  const int Nx = CommonParameters::Nx();
130  const int Ny = CommonParameters::Ny();
131  const int Nz = CommonParameters::Nz();
132  const int Nt = CommonParameters::Nt();
133 
134  static const double PI = 4.0 * atan(1.0);
135 
136  std::vector<double> p_unit(Ndim - 1);
137 
138  p_unit[x_dir] = (2.0 * PI / Lx) * m_source_momentum[x_dir];
139  p_unit[y_dir] = (2.0 * PI / Ly) * m_source_momentum[y_dir];
140  p_unit[z_dir] = (2.0 * PI / Lz) * m_source_momentum[z_dir];
141 
142  std::vector<int> ipe(Ndim - 1);
143  ipe[x_dir] = Communicator::ipe(x_dir);
144  ipe[y_dir] = Communicator::ipe(y_dir);
145  ipe[z_dir] = Communicator::ipe(z_dir);
146 
147  //- clear field
148  src.set(0.0);
149 
150  if (m_in_node) {
151  const int t = m_source_position[t_dir] % Nt;
152 
153  for (int z = 0; z < Nz; ++z) {
154  for (int y = 0; y < Ny; ++y) {
155  for (int x = 0; x < Nx; ++x) {
156  int site = m_index.site(x, y, z, t);
157 
158  int x_global = x + ipe[x_dir] * Nx;
159  int y_global = y + ipe[y_dir] * Ny;
160  int z_global = z + ipe[z_dir] * Nz;
161 
162  double p_x = p_unit[x_dir] * (x_global - m_source_position[x_dir]);
163  double p_y = p_unit[y_dir] * (y_global - m_source_position[y_dir]);
164  double p_z = p_unit[z_dir] * (z_global - m_source_position[z_dir]);
165 
166  double theta = p_x + p_y + p_z;
167 
168  int idx_r = 2 * idx;
169  int idx_i = 2 * idx + 1;
170 
171  //XXX field layout: complex as two doubles
172  src.set(idx_r, site, 0, cos(theta) / Lvol3);
173  src.set(idx_i, site, 0, sin(theta) / Lvol3);
174  }
175  }
176  }
177  }
178 }
179 
180 
181 //====================================================================
182 void Source_MomentumWall::set(Field& src, const int i_color, const int i_spin)
183 {
184  const int Nc = CommonParameters::Nc();
185  const int idx = i_color + Nc * i_spin;
186 
187  set(src, idx);
188 }
189 
190 
191 //====================================================================
192 void Source_MomentumWall::set_all_color(Field& src, const int i_spin)
193 {
194  const int Nc = CommonParameters::Nc();
195  const int Ndim = CommonParameters::Ndim();
196 
197  const int x_dir = 0;
198  const int y_dir = 1;
199  const int z_dir = 2;
200  const int t_dir = Ndim - 1;
201 
202  //- global lattice size
203  const int Lx = CommonParameters::Lx();
204  const int Ly = CommonParameters::Ly();
205  const int Lz = CommonParameters::Lz();
206  const int Lvol3 = Lx * Ly * Lz;
207 
208  //- local size
209  const int Nx = CommonParameters::Nx();
210  const int Ny = CommonParameters::Ny();
211  const int Nz = CommonParameters::Nz();
212  const int Nt = CommonParameters::Nt();
213 
214  static const double PI = 4.0 * atan(1.0);
215 
216  std::vector<double> p_unit(Ndim - 1);
217 
218  p_unit[x_dir] = (2.0 * PI / Lx) * m_source_momentum[x_dir];
219  p_unit[y_dir] = (2.0 * PI / Ly) * m_source_momentum[y_dir];
220  p_unit[z_dir] = (2.0 * PI / Lz) * m_source_momentum[z_dir];
221 
222  std::vector<int> ipe(Ndim - 1);
223  ipe[x_dir] = Communicator::ipe(x_dir);
224  ipe[y_dir] = Communicator::ipe(y_dir);
225  ipe[z_dir] = Communicator::ipe(z_dir);
226 
227  //- clear field
228  src.set(0.0);
229 
230  if (m_in_node) {
231  const int t = m_source_position[t_dir] % Nt;
232 
233  for (int z = 0; z < Nz; ++z) {
234  for (int y = 0; y < Ny; ++y) {
235  for (int x = 0; x < Nx; ++x) {
236  int site = m_index.site(x, y, z, t);
237 
238  int x_global = x + ipe[x_dir] * Nx;
239  int y_global = y + ipe[y_dir] * Ny;
240  int z_global = z + ipe[z_dir] * Nz;
241 
242  double p_x = p_unit[x_dir] * (x_global - m_source_position[x_dir]);
243  double p_y = p_unit[y_dir] * (y_global - m_source_position[y_dir]);
244  double p_z = p_unit[z_dir] * (z_global - m_source_position[z_dir]);
245 
246  double theta = p_x + p_y + p_z;
247 
248  for (int i_color = 0; i_color < Nc; ++i_color) {
249  int idx = i_color + Nc * i_spin;
250  int idx_r = 2 * idx;
251  int idx_i = 2 * idx + 1;
252 
253  //XXX field layout: complex as two doubles
254  src.set(idx_r, site, 0, cos(theta) / Lvol3);
255  src.set(idx_i, site, 0, sin(theta) / Lvol3);
256  }
257  }
258  }
259  }
260  }
261 }
262 
263 
264 //====================================================================
266 {
267  const int Nc = CommonParameters::Nc();
268  const int Nd = CommonParameters::Nd();
269  const int Ndim = CommonParameters::Ndim();
270 
271  const int x_dir = 0;
272  const int y_dir = 1;
273  const int z_dir = 2;
274  const int t_dir = Ndim - 1;
275 
276  //- global lattice size
277  const int Lx = CommonParameters::Lx();
278  const int Ly = CommonParameters::Ly();
279  const int Lz = CommonParameters::Lz();
280  const int Lvol3 = Lx * Ly * Lz;
281 
282  //- local size
283  const int Nx = CommonParameters::Nx();
284  const int Ny = CommonParameters::Ny();
285  const int Nz = CommonParameters::Nz();
286  const int Nt = CommonParameters::Nt();
287 
288  static const double PI = 4.0 * atan(1.0);
289 
290  std::vector<double> p_unit(Ndim - 1);
291 
292  p_unit[x_dir] = (2.0 * PI / Lx) * m_source_momentum[x_dir];
293  p_unit[y_dir] = (2.0 * PI / Ly) * m_source_momentum[y_dir];
294  p_unit[z_dir] = (2.0 * PI / Lz) * m_source_momentum[z_dir];
295 
296  std::vector<int> ipe(Ndim - 1);
297  ipe[x_dir] = Communicator::ipe(x_dir);
298  ipe[y_dir] = Communicator::ipe(y_dir);
299  ipe[z_dir] = Communicator::ipe(z_dir);
300 
301  //- clear field
302  src.set(0.0);
303 
304  if (m_in_node) {
305  const int t = m_source_position[t_dir] % Nt;
306 
307  for (int z = 0; z < Nz; ++z) {
308  for (int y = 0; y < Ny; ++y) {
309  for (int x = 0; x < Nx; ++x) {
310  int site = m_index.site(x, y, z, t);
311 
312  int x_global = x + ipe[x_dir] * Nx;
313  int y_global = y + ipe[y_dir] * Ny;
314  int z_global = z + ipe[z_dir] * Nz;
315 
316  double p_x = p_unit[x_dir] * (x_global - m_source_position[x_dir]);
317  double p_y = p_unit[y_dir] * (y_global - m_source_position[y_dir]);
318  double p_z = p_unit[z_dir] * (z_global - m_source_position[z_dir]);
319 
320  double theta = p_x + p_y + p_z;
321 
322  for (int idx = 0; idx < Nc * Nd; ++idx) {
323  int idx_r = 2 * idx;
324  int idx_i = 2 * idx + 1;
325 
326  //XXX field layout: complex as two doubles
327  src.set(idx_r, site, 0, cos(theta) / Lvol3);
328  src.set(idx_i, site, 0, sin(theta) / Lvol3);
329  }
330  }
331  }
332  }
333  }
334 }
335 
336 
337 //====================================================================
338 //============================================================END=====
BridgeIO vout
Definition: bridgeIO.cpp:503
void set_all_color_spin(Field &v)
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
int site(const int &x, const int &y, const int &z, const int &t) const
Definition: index_lex.h:53
std::vector< int > m_source_momentum
void general(const char *format,...)
Definition: bridgeIO.cpp:197
static const std::string class_name
Momentum wall source.
Container of Field-type object.
Definition: field.h:45
Class for parameters.
Definition: parameters.h:46
static int ipe(const int dir)
logical coordinate of current proc.
void set_all_color(Field &v, const int i_spin)
void crucial(const char *format,...)
Definition: bridgeIO.cpp:178
Bridge::VerboseLevel m_vl
Definition: source.h:40
int non_negative(const int v)
string get_string(const string &key) const
Definition: parameters.cpp:221
int fetch_int_vector(const string &key, vector< int > &value) const
Definition: parameters.cpp:429
std::vector< int > m_source_position
void set_parameters(const Parameters &params)
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:131
void set(Field &v, const int idx)