Bridge++  Version 1.5.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
source_Random.cpp
Go to the documentation of this file.
1 
14 #include "source_Random.h"
15 
16 #ifdef USE_FACTORY_AUTOREGISTER
17 namespace {
18  bool init = Source_Random::register_factory();
19 }
20 #endif
21 
22 const std::string Source_Random::class_name = "Source_Random";
23 
24 //====================================================================
25 void Source_Random::set_parameters(const Parameters& params)
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  string str_noise_type;
35 
36  int err = 0;
37  err += params.fetch_int_vector("source_position", source_position);
38  err += params.fetch_int_vector("source_momentum", source_momentum);
39  err += params.fetch_string("noise_type", str_noise_type);
40 
41  if (err) {
42  vout.crucial(m_vl, "Error at %s: input parameter not found.\n", class_name.c_str());
43  exit(EXIT_FAILURE);
44  }
45 
46 
47  set_parameters(source_position, source_momentum, str_noise_type);
48 }
49 
50 
51 //====================================================================
52 void Source_Random::set_parameters(const std::vector<int>& source_position,
53  const std::vector<int>& source_momentum,
54  const string str_noise_type)
55 {
56  // #### parameter setup ####
57  const int Ndim = CommonParameters::Ndim();
58  const int t_dir = Ndim - 1;
59 
60  //- global lattice size
61  std::vector<int> Lsize(Ndim);
62 
63  Lsize[0] = CommonParameters::Lx();
64  Lsize[1] = CommonParameters::Ly();
65  Lsize[2] = CommonParameters::Lz();
66  Lsize[3] = CommonParameters::Lt();
67 
68  //- local size
69  const int Nt = CommonParameters::Nt();
70 
71  //- print input parameters
72  vout.general(m_vl, "Source for spinor field - Random:\n");
73  for (int mu = 0; mu < Ndim; ++mu) {
74  vout.general(m_vl, " source_position[%d] = %d\n",
75  mu, source_position[mu]);
76  }
77 
78  for (int mu = 0; mu < Ndim - 1; ++mu) {
79  vout.general(m_vl, " source_momentum[%d] = %d\n",
80  mu, source_momentum[mu]);
81  }
82 
83  vout.general(m_vl, " noise_type = %s\n", str_noise_type.c_str());
84 
85 
86  //- range check
87  int err = 0;
88  for (int mu = 0; mu < Ndim; ++mu) {
89  // NB. Lsize[mu] > abs(source_position[mu])
90  err += ParameterCheck::non_negative(Lsize[mu] - abs(source_position[mu]));
91  }
92 
93  if (err) {
94  vout.crucial(m_vl, "Error at %s: parameter range check failed.\n", class_name.c_str());
95  exit(EXIT_FAILURE);
96  }
97 
98 
99  //- store values
100  m_source_position.resize(Ndim);
101  for (int mu = 0; mu < Ndim; ++mu) {
102  m_source_position[mu] = (source_position[mu] + Lsize[mu]) % Lsize[mu];
103  }
104  m_source_momentum = source_momentum;
105  m_str_noise_type = str_noise_type;
106 
107 
108  //-- post-process
109  // m_rand->reset(static_cast<unsigned long>(i_seed_noise));
110 
111  //- PE location in t-direction.
112  const int tpe = m_source_position[t_dir] / Nt;
113 
114  m_in_node = false;
115 
116  if (tpe == Communicator::ipe(t_dir)) {
117  m_in_node = true;
118  }
119 }
120 
121 
122 //====================================================================
123 void Source_Random::set(Field& src, const int idx)
124 {
125  const int Ndim = CommonParameters::Ndim();
126 
127  //- local size
128  std::vector<int> Nsize(Ndim);
129 
130  Nsize[0] = CommonParameters::Nx();
131  Nsize[1] = CommonParameters::Ny();
132  Nsize[2] = CommonParameters::Nz();
133  Nsize[3] = CommonParameters::Nt();
134 
135  Field v(src);
136 
137  if (m_str_noise_type == "Gaussian") {
138  m_rand->gauss_lex_global(v);
139  } else if (m_str_noise_type == "U1") {
140  m_rand->U1_lex_global(v);
141  } else if (m_str_noise_type == "Z2") {
142  m_rand->Z2_lex_global(v);
143  } else {
144  vout.crucial(m_vl, "Error at %s: Noise type of %s is not supported.\n", class_name.c_str(), m_str_noise_type.c_str());
145  exit(EXIT_FAILURE);
146  }
147 
148  //- clear field
149  src.set(0.0);
150 
151  const int idx_r = 2 * idx;
152  const int idx_i = 2 * idx + 1;
153 
154  if (m_in_node) {
155  int t = m_source_position[3] % Nsize[3];
156 
157  for (int z = 0; z < Nsize[2]; ++z) {
158  for (int y = 0; y < Nsize[1]; ++y) {
159  for (int x = 0; x < Nsize[0]; ++x) {
160  int site = m_index.site(x, y, z, t);
161 
162  double rn_r = v.cmp(idx_r, site, 0);
163  double rn_i = v.cmp(idx_i, site, 0);
164 
165  //XXX field layout: complex as two doubles
166  src.set(idx_r, site, 0, rn_r);
167  src.set(idx_i, site, 0, rn_i);
168 
169  //- for debug
170  vout.paranoiac(m_vl, "%s: src(x,y,z) = %e %e (%d %d %d)\n",
171  class_name.c_str(), rn_r, rn_i, x, y, z);
172  }
173  }
174  }
175  }
176 
177 #ifdef DEBUG
178  //- global lattice size
179  const int Lx = CommonParameters::Lx();
180  const int Ly = CommonParameters::Ly();
181  const int Lz = CommonParameters::Lz();
182  const int Lvol3 = Lx * Ly * Lz;
183 
184  const int Nvol = CommonParameters::Nvol();
185 
186  double src1_r = 0.0;
187  double src1_i = 0.0;
188  for (int site = 0; site < Nvol; ++site) {
189  src1_r += src.cmp(idx_r, site, 0);
190  src1_i += src.cmp(idx_i, site, 0);
191  }
192  const double src_r = Communicator::reduce_sum(src1_r) / Lvol3;
193  const double src_i = Communicator::reduce_sum(src1_i) / Lvol3;
194  vout.general(m_vl, "%s: sum_x src = %e %e\n", class_name.c_str(), src_r, src_i);
195 
196  const double src2 = src.norm2() / Lvol3;
197  vout.general(m_vl, "%s: src.norm2 = %e\n", class_name.c_str(), src2);
198 #endif
199 }
200 
201 
202 //====================================================================
203 void Source_Random::set(Field& src, const int i_color, const int i_spin)
204 {
205  const int Nc = CommonParameters::Nc();
206  const int idx = i_color + Nc * i_spin;
207 
208  set(src, idx);
209 }
210 
211 
212 //====================================================================
213 void Source_Random::set_all_color(Field& src, const int i_spin)
214 {
215  const int Nc = CommonParameters::Nc();
216  const int Ndim = CommonParameters::Ndim();
217 
218  //- local size
219  std::vector<int> Nsize(Ndim);
220 
221  Nsize[0] = CommonParameters::Nx();
222  Nsize[1] = CommonParameters::Ny();
223  Nsize[2] = CommonParameters::Nz();
224  Nsize[3] = CommonParameters::Nt();
225 
226  Field v(src);
227 
228  if (m_str_noise_type == "Gaussian") {
229  m_rand->gauss_lex_global(v);
230  } else if (m_str_noise_type == "U1") {
231  m_rand->U1_lex_global(v);
232  } else if (m_str_noise_type == "Z2") {
233  m_rand->Z2_lex_global(v);
234  } else {
235  vout.crucial(m_vl, "Error at %s: Noise type of %s is not supported.\n", class_name.c_str(), m_str_noise_type.c_str());
236  exit(EXIT_FAILURE);
237  }
238 
239  //- clear field
240  src.set(0.0);
241 
242  if (m_in_node) {
243  int t = m_source_position[3] % Nsize[3];
244 
245  for (int z = 0; z < Nsize[2]; ++z) {
246  for (int y = 0; y < Nsize[1]; ++y) {
247  for (int x = 0; x < Nsize[0]; ++x) {
248  int site = m_index.site(x, y, z, t);
249 
250  for (int i_color = 0; i_color < Nc; ++i_color) {
251  int idx = i_color + Nc * i_spin;
252  int idx_r = 2 * idx;
253  int idx_i = 2 * idx + 1;
254 
255  double rn_r = v.cmp(idx_r, site, 0);
256  double rn_i = v.cmp(idx_i, site, 0);
257 
258  //XXX field layout: complex as two doubles
259  src.set(idx_r, site, 0, rn_r);
260  src.set(idx_i, site, 0, rn_i);
261 
262  //- for debug
263  vout.paranoiac(m_vl, "%s: src(idx, x,y,z) = %e %e (%d %d %d %d)\n",
264  class_name.c_str(), rn_r, rn_i, idx, x, y, z);
265  }
266  }
267  }
268  }
269  }
270 
271 #ifdef DEBUG
272  //- global lattice size
273  const int Lx = CommonParameters::Lx();
274  const int Ly = CommonParameters::Ly();
275  const int Lz = CommonParameters::Lz();
276  const int Lvol3 = Lx * Ly * Lz;
277 
278  const int Nvol = CommonParameters::Nvol();
279 
280  const int Nd = CommonParameters::Nd();
281 
282  double src1_r = 0.0;
283  double src1_i = 0.0;
284  for (int site = 0; site < Nvol; ++site) {
285  for (int idx = 0; idx < Nc * Nd; ++idx) {
286  int idx_r = 2 * idx;
287  int idx_i = 2 * idx + 1;
288 
289  src1_r += src.cmp(idx_r, site, 0);
290  src1_i += src.cmp(idx_i, site, 0);
291  }
292  }
293  const double src_r = Communicator::reduce_sum(src1_r) / (Lvol3 * Nc * Nd);
294  const double src_i = Communicator::reduce_sum(src1_i) / (Lvol3 * Nc * Nd);
295  vout.general(m_vl, "%s: sum_x src = %e %e\n", class_name.c_str(), src_r, src_i);
296 
297  const double src2 = src.norm2() / (Lvol3 * Nc * Nd);
298  vout.general(m_vl, "%s: src.norm2 = %e\n", class_name.c_str(), src2);
299 #endif
300 }
301 
302 
303 //====================================================================
304 void Source_Random::set_all_color_spin(Field& src)
305 {
306  const int Nc = CommonParameters::Nc();
307  const int Nd = CommonParameters::Nd();
308  const int Ndim = CommonParameters::Ndim();
309 
310  //- local size
311  std::vector<int> Nsize(Ndim);
312 
313  Nsize[0] = CommonParameters::Nx();
314  Nsize[1] = CommonParameters::Ny();
315  Nsize[2] = CommonParameters::Nz();
316  Nsize[3] = CommonParameters::Nt();
317 
318  Field v(src);
319 
320  if (m_str_noise_type == "Gaussian") {
321  m_rand->gauss_lex_global(v);
322  } else if (m_str_noise_type == "U1") {
323  m_rand->U1_lex_global(v);
324  } else if (m_str_noise_type == "Z2") {
325  m_rand->Z2_lex_global(v);
326  } else {
327  vout.crucial(m_vl, "Error at %s: Noise type of %s is not supported.\n", class_name.c_str(), m_str_noise_type.c_str());
328  exit(EXIT_FAILURE);
329  }
330 
331  //- clear field
332  src.set(0.0);
333 
334  if (m_in_node) {
335  int t = m_source_position[3] % Nsize[3];
336 
337  for (int z = 0; z < Nsize[2]; ++z) {
338  for (int y = 0; y < Nsize[1]; ++y) {
339  for (int x = 0; x < Nsize[0]; ++x) {
340  int site = m_index.site(x, y, z, t);
341 
342  for (int idx = 0; idx < Nc * Nd; ++idx) {
343  int idx_r = 2 * idx;
344  int idx_i = 2 * idx + 1;
345 
346  double rn_r = v.cmp(idx_r, site, 0);
347  double rn_i = v.cmp(idx_i, site, 0);
348 
349  //XXX field layout: complex as two doubles
350  src.set(idx_r, site, 0, rn_r);
351  src.set(idx_i, site, 0, rn_i);
352 
353  //- for debug
354  vout.paranoiac(m_vl, "%s: src(idx, x,y,z) = %e %e (%d %d %d %d)\n",
355  class_name.c_str(), rn_r, rn_i, idx, x, y, z);
356  }
357  }
358  }
359  }
360  }
361 
362 #ifdef DEBUG
363  //- global lattice size
364  const int Lx = CommonParameters::Lx();
365  const int Ly = CommonParameters::Ly();
366  const int Lz = CommonParameters::Lz();
367  const int Lvol3 = Lx * Ly * Lz;
368 
369  const int Nvol = CommonParameters::Nvol();
370 
371  double src1_r = 0.0;
372  double src1_i = 0.0;
373  for (int site = 0; site < Nvol; ++site) {
374  for (int idx = 0; idx < Nc * Nd; ++idx) {
375  int idx_r = 2 * idx;
376  int idx_i = 2 * idx + 1;
377 
378  src1_r += src.cmp(idx_r, site, 0);
379  src1_i += src.cmp(idx_i, site, 0);
380  }
381  }
382  const double src_r = Communicator::reduce_sum(src1_r) / (Lvol3 * Nc * Nd);
383  const double src_i = Communicator::reduce_sum(src1_i) / (Lvol3 * Nc * Nd);
384  vout.general(m_vl, "%s: sum_x src = %e %e\n", class_name.c_str(), src_r, src_i);
385 
386  const double src2 = src.norm2() / (Lvol3 * Nc * Nd);
387  vout.general(m_vl, "%s: src.norm2 = %e\n", class_name.c_str(), src2);
388 #endif
389 }
390 
391 
392 //====================================================================
393 void Source_Random::set_all_space_time(Field& src, const int ic)
394 {
395  // This implementation assumes the given field src is complex field.
396  int nc = CommonParameters::Nc();
397  int Nex = src.nex();
398  int Nvol = src.nvol();
399  int Nin = src.nin();
400 
401  assert(ic < nc);
402  assert(Nin >= 2 * nc);
403  int nd = Nin / nc / 2;
404 // double rn, rn2, rn3;
405  double rn;
406 // double RF2 = 1.0 / sqrt(2.0);
407 
408  Field v(2, Nvol, Nex);
409  if (m_str_noise_type == "Gaussian") {
410  m_rand->gauss_lex_global(v);
411  } else if (m_str_noise_type == "U1") {
412  m_rand->U1_lex_global(v);
413  } else if (m_str_noise_type == "Z2") {
414  m_rand->Z2_lex_global(v);
415  } else {
416  vout.crucial(m_vl, "Error at %s: Noise type of %s is not supported.\n", class_name.c_str(), m_str_noise_type.c_str());
417  exit(EXIT_FAILURE);
418  }
419  // m_rand->uniform_lex_global(v);
420  src.set(0.0);
421 
422  for (int ex = 0; ex < Nex; ++ex) {
423  for (int site = 0; site < Nvol; ++site) {
424  for (int ri = 0; ri < 2; ++ri) {
425  rn = v.cmp(ri, site, ex);
426  for (int id = 0; id < nd; ++id) {
427  int in = ri + 2 * (ic + nc * id);
428 // rn2 = floor(2.0 * rn);
429 // rn3 = (2.0 * rn2 - 1.0) * RF2;
430 // src.set(in, site, ex, rn3);
431  src.set(in, site, ex, rn);
432  }
433  }
434  }
435  }
436 }
437 
438 
439 //====================================================================
440 void Source_Random::set_all_space_time(Field& src, const int ic, const int is)
441 {
442  // This implementation assumes the given field src is complex field.
443  int nc = CommonParameters::Nc();
444  int nd = CommonParameters::Nd();
445  int Nex = src.nex();
446  int Nvol = src.nvol();
447  int Nin = src.nin();
448 
449  assert(ic < nc);
450  assert(is < nd);
451  assert(Nin == 2 * nc * nd);
452 
453  Field v(2, Nvol, Nex);
454  if (m_str_noise_type == "Gaussian") {
455  m_rand->gauss_lex_global(v);
456  } else if (m_str_noise_type == "U1") {
457  m_rand->U1_lex_global(v);
458  } else if (m_str_noise_type == "Z2") {
459  m_rand->Z2_lex_global(v);
460  } else {
461  vout.crucial(m_vl, "Error at %s: Noise type of %s is not supported.\n", class_name.c_str(), m_str_noise_type.c_str());
462  exit(EXIT_FAILURE);
463  }
464  // m_rand->uniform_lex_global(v);
465  src.set(0.0);
466 
467 // double rn, rn2, rn3;
468  double rn;
469 // double RF2 = 1.0 / sqrt(2.0);
470 
471  for (int ex = 0; ex < Nex; ++ex) {
472  for (int site = 0; site < Nvol; ++site) {
473  for (int ri = 0; ri < 2; ++ri) {
474  int in = ri + 2 * (ic + nc * is);
475  rn = v.cmp(ri, site, ex);
476 // rn2 = floor(2.0 * rn);
477 // rn3 = (2.0 * rn2 - 1.0) * RF2;
478 // src.set(in, site, ex, rn3);
479  src.set(in, site, ex, rn);
480  }
481  }
482  }
483 }
484 
485 
486 //====================================================================
487 //============================================================END=====
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:503
Source_Random::class_name
static const std::string class_name
Definition: source_Random.h:38
Field::norm2
double norm2() const
Definition: field.cpp:592
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
Index_lex::site
int site(const int &x, const int &y, const int &z, const int &t) const
Definition: index_lex.h:53
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:197
Source_Random::m_in_node
bool m_in_node
Definition: source_Random.h:46
Field
Container of Field-type object.
Definition: field.h:45
Field::nvol
int nvol() const
Definition: field.h:127
Source_Random::set_all_color_spin
void set_all_color_spin(Field &src)
Definition: source_Random.cpp:304
RandomNumbers::U1_lex_global
virtual void U1_lex_global(Field &)
U(1) random number defined on global lattice.
Definition: randomNumbers.cpp:113
Field::cmp
double cmp(const int jin, const int site, const int jex) const
Definition: field.h:143
Parameters
Class for parameters.
Definition: parameters.h:46
Communicator::ipe
static int ipe(const int dir)
logical coordinate of current proc.
Definition: communicator.cpp:105
Parameters::fetch_string
int fetch_string(const string &key, string &value) const
Definition: parameters.cpp:378
RandomNumbers::gauss_lex_global
virtual void gauss_lex_global(Field &)
gaussian random number defined on global lattice.
Definition: randomNumbers.cpp:95
Source_Random::set_all_space_time
void set_all_space_time(Field &src, const int ic)
Setting a noise vector. Filling all the sites and spin indices for color index "ic". The same random number is set for all the spin indices at one site.
Definition: source_Random.cpp:393
Source_Random::m_source_position
std::vector< int > m_source_position
Definition: source_Random.h:43
source_Random.h
Field::nin
int nin() const
Definition: field.h:126
Source_Random::m_str_noise_type
std::string m_str_noise_type
Definition: source_Random.h:45
Source_Random::set_all_color
void set_all_color(Field &src, const int i_spin)
Definition: source_Random.cpp:213
Source_Random::m_rand
RandomNumbers * m_rand
Definition: source_Random.h:41
Field::nex
int nex() const
Definition: field.h:128
Bridge::BridgeIO::paranoiac
void paranoiac(const char *format,...)
Definition: bridgeIO.cpp:235
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:178
Source_Random::m_index
Index_lex m_index
Definition: source_Random.h:42
Source_Random::m_source_momentum
std::vector< int > m_source_momentum
Definition: source_Random.h:44
RandomNumbers::Z2_lex_global
virtual void Z2_lex_global(Field &)
Z(2) random number defined on global lattice.
Definition: randomNumbers.cpp:142
Source::m_vl
Bridge::VerboseLevel m_vl
Definition: source.h:40
ParameterCheck::non_negative
int non_negative(const int v)
Definition: parameterCheck.cpp:21
Communicator::reduce_sum
static int reduce_sum(int count, double *recv_buf, double *send_buf, int pattern=0)
make a global sum of an array of double over the communicator. pattern specifies the dimensions to be...
Definition: communicator.cpp:244
Source_Random::set
void set(Field &src, const int idx)
Definition: source_Random.cpp:123
Source_Random::set_parameters
void set_parameters(const Parameters &params)
Definition: source_Random.cpp:25
Parameters::get_string
string get_string(const string &key) const
Definition: parameters.cpp:221
Parameters::fetch_int_vector
int fetch_int_vector(const string &key, vector< int > &value) const
Definition: parameters.cpp:429
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:131