Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
test_Spectrum_Clover_2ptFunction_withFieldIO.cpp
Go to the documentation of this file.
1 
13 #include <sstream>
14 #include "parameterManager_YAML.h"
15 #include "parameters_factory.h"
16 
17 #include "bridgeIO.h"
18 using Bridge::vout;
19 
20 #include "gaugeConfig.h"
21 #include "staples.h"
22 
23 #include "randomNumbers_Mseries.h"
24 
25 #include "fopr_Clover.h"
26 #include "fprop_Standard_lex.h"
27 
28 #include "director_Smear.h"
29 #include "fopr_Smeared.h"
30 
31 #include "gammaMatrixSet.h"
32 #include "gaugeFixing.h"
33 #include "projection.h"
34 #include "smear.h"
35 #include "solver.h"
36 #include "source.h"
37 
38 #include "corr2pt_4spinor.h"
39 
40 #include "fieldIO_Binary.h"
41 
42 #ifdef USE_TEST
43 #include "test.h"
44 #endif
45 
46 #ifdef USE_TESTMANAGER_AUTOREGISTER
47 #include "testManager.h"
48 #endif
49 
50 //====================================================================
52 
67 namespace Test_Spectrum_Clover {
68  //- test-private parameters
69  namespace {
70  const std::string filename_input = "test_Spectrum_Clover_Hadron2ptFunction.yaml";
71  const std::string filename_output = "stdout";
72 
73  class Parameters_Test_Spectrum_Clover : public Parameters {
74  public:
75  Parameters_Test_Spectrum_Clover()
76  {
77  Register_string("gauge_config_type_input", "NULL");
78  Register_string("config_filename_input", "NULL");
79 
80  Register_string("gauge_config_type_output", "NULL");
81  Register_string("config_filename_output", "NULL");
82 
83  Register_int("number_of_valence_quarks", 0);
84 
85  Register_string("verbose_level", "NULL");
86 
87  Register_double("expected_result", 0.0);
88  }
89  };
90 
91  class Parameters_Quark : virtual public Parameters {
92  public:
93  Parameters_Quark()
94  {
95  Register_string("temporary_filename_base", "NULL");
96 
97  Register_Parameters("Fopr_Clover", ParametersFactory::New("Fopr.Clover"));
98  Register_Parameters("Source", ParametersFactory::New("Source"));
99  }
100  };
101  }
102 
103  //- prototype declaration
104  int hadron_2ptFunction_withFileIO(void);
105 
106  int calculate_quark_propagator();
107  int calculate_hadron_correlator();
108 
109  string generate_filename(const string&, const int);
110  string generate_filename(const string&, const int, const int);
111 
112 #ifdef USE_TESTMANAGER_AUTOREGISTER
113  namespace {
114  static const bool is_registered = TestManager::RegisterTest(
115  "Spectrum.Clover.Hadron2ptFunction_withFileIO",
116  hadron_2ptFunction_withFileIO
117  );
118  }
119 #endif
120 
121  //====================================================================
122  int hadron_2ptFunction_withFileIO(void)
123  {
124  int result = 0;
125 
126  result += calculate_quark_propagator();
127  result += calculate_hadron_correlator();
128 
129  return result;
130  }
131 
132 
133  //====================================================================
134  string generate_filename(const string& base, const int idx)
135  {
136  char buf[1024];
137 
138  sprintf(buf, "%s_%02d.dat", base.c_str(), idx);
139 
140  return string(buf);
141  }
142 
143 
144  //====================================================================
145  string generate_filename(const string& base, const int icolor, const int ispin)
146  {
147  char buf[1024];
148 
149  sprintf(buf, "%s_c%d_s%d.dat", base.c_str(), (icolor + 1), (ispin + 1));
150 
151  return string(buf);
152  }
153 
154 
155  //====================================================================
156  int calculate_quark_propagator()
157  {
158  // #### parameter setup ####
159  int Nc = CommonParameters::Nc();
160  int Nd = CommonParameters::Nd();
161  int Ndim = CommonParameters::Ndim();
162  int Nvol = CommonParameters::Nvol();
163 
164  ParameterManager_YAML params_manager;
165 
166  Parameters *params_test = new Parameters_Test_Spectrum_Clover;
167 
168  Parameters *params_pre = new Parameters;
169 
170  params_pre->Register_Parameters("Test_Spectrum_Clover", params_test);
171 
172  params_manager.read_params(filename_input, params_pre);
173 
174 
175  const int N_quark = params_test->get_int("number_of_valence_quarks");
176  std::vector<Parameters *> params_quark(N_quark);
177 
178  for (int iq = 0; iq < N_quark; ++iq) {
179  params_quark[iq] = new Parameters_Quark;
180  }
181 
182  Parameters *params_gfix = ParametersFactory::New("GaugeFixing");
183  Parameters *params_proj = ParametersFactory::New("Projection");
184  Parameters *params_smear = ParametersFactory::New("Smear");
185  Parameters *params_dr_smear = ParametersFactory::New("Director_Smear");
186  Parameters *params_solver = ParametersFactory::New("Solver");
187 
188  Parameters *params_all = new Parameters;
189 
190 
191  std::stringstream oss;
192 
193  for (int iq = 0; iq < N_quark; ++iq) {
194  oss.str("");
195  oss << "Quark_" << iq + 1;
196  string str = oss.str();
197  params_all->Register_Parameters(str.c_str(), params_quark[iq]);
198  }
199 
200  params_all->Register_Parameters("GaugeFixing", params_gfix);
201  params_all->Register_Parameters("Projection", params_proj);
202  params_all->Register_Parameters("Smear", params_smear);
203  params_all->Register_Parameters("Director_Smear", params_dr_smear);
204  params_all->Register_Parameters("Solver", params_solver);
205 
206  params_manager.read_params(filename_input, params_all);
207 
208  const string str_gconf_read = params_test->get_string("gauge_config_type_input");
209  const string readfile = params_test->get_string("config_filename_input");
210  const string str_vlevel = params_test->get_string("verbose_level");
211 
212  const string str_gfix_type = params_gfix->get_string("gauge_fixing_type");
213 
214  std::vector<std::string> savefile_base(N_quark);
215  std::vector<Parameters *> params_clover(N_quark);
216  std::vector<Parameters *> params_source(N_quark);
217 
218  for (int iq = 0; iq < N_quark; ++iq) {
219  savefile_base[iq] = params_quark[iq]->get_string("temporary_filename_base");
220  params_clover[iq] = params_quark[iq]->get_Parameters("Fopr_Clover");
221  params_source[iq] = params_quark[iq]->get_Parameters("Source");
222  }
223 
224  const string str_proj_type = params_proj->get_string("projection_type");
225  const string str_smear_type = params_smear->get_string("smear_type");
226  const string str_solver_type = params_solver->get_string("solver_type");
227 
228  Bridge::VerboseLevel vl = vout.set_verbose_level(str_vlevel);
229 
230  //- print input parameters
231  vout.general(vl, " gconf_read = %s\n", str_gconf_read.c_str());
232  vout.general(vl, " readfile = %s\n", readfile.c_str());
233  vout.general(vl, " vlevel = %s\n", str_vlevel.c_str());
234  vout.general(vl, " gfix_type = %s\n", str_gfix_type.c_str());
235  vout.general(vl, " proj_type = %s\n", str_proj_type.c_str());
236  vout.general(vl, " smear_type = %s\n", str_smear_type.c_str());
237  vout.general(vl, " solver_type = %s\n", str_solver_type.c_str());
238 
239  // NB. all str_gmset_type are supposed to be the same.
240  string str_gmset_type = params_clover[0]->get_string("gamma_matrix_type");
241  vout.general(vl, " gmset_type = %s\n", str_gmset_type.c_str());
242 
243  std::vector<std::string> str_source_type(N_quark);
244 
245  for (int iq = 0; iq < N_quark; ++iq) {
246  vout.general(vl, " Quark_%d:\n", iq + 1);
247 
248  str_source_type[iq] = params_source[iq]->get_string("source_type");
249  vout.general(vl, " source_type = %s\n", str_source_type[iq].c_str());
250  }
251  vout.general(vl, "\n");
252 
253  //- input parameter check
254  int err = 0;
255  err += ParameterCheck::non_NULL(str_gconf_read);
256  err += ParameterCheck::non_NULL(readfile);
257 
258  if (err) {
259  vout.crucial(vl, "Test_Spectrum_Clover: Input parameters have not been set.\n");
260  abort();
261  }
262 
263 
264  // #### Set up a gauge configuration ####
265  Field_G *U = new Field_G(Nvol, Ndim);
266  GaugeConfig *gconf_read = new GaugeConfig(str_gconf_read);
267  gconf_read->read_file((Field *)U, readfile);
268  // gconf_read->set_cold((Field*)U);
269 
270  Projection *proj = Projection::New(str_proj_type);
271  Smear *smear = Smear::New(str_smear_type, proj);
272  smear->set_parameters(*params_smear);
273 
274  Director_Smear *dr_smear = new Director_Smear((Smear *)smear);
275  dr_smear->set_parameters(*params_dr_smear);
276  dr_smear->set_config(U);
277 
278  int Nsmear = dr_smear->get_Nsmear();
279  Field_G *Usmear = (Field_G *)dr_smear->getptr_smearedConfig(Nsmear);
280 
281 
282  // #### Gauge fixing ####
283  Staples *staple = new Staples;
284  Field_G *Ufix = new Field_G(Nvol, Ndim);
285 
286  int ndelay = 1000;
287  RandomNumbers *rand = new RandomNumbers_Mseries(ndelay);
288 
289  GaugeFixing *gfix = GaugeFixing::New(str_gfix_type, rand);
290  gfix->set_parameters(*params_gfix);
291 
292  double plaq = staple->plaquette(*Usmear);
293  vout.general(vl, "plaq(original) = %18.14f\n", plaq);
294 
295  gfix->fix(*Ufix, *Usmear);
296 
297  double plaq2 = staple->plaquette(*Ufix);
298  vout.general(vl, "plaq(fixed) = %18.14f\n", plaq2);
299  vout.general(vl, "plaq(diff) = %18.10e\n", plaq - plaq2);
300 
301 
302  // #### object setup #####
303  std::vector<Fopr_Clover *> fopr_c(N_quark);
304  std::vector<Solver *> solver(N_quark);
305  std::vector<Fprop *> fprop_lex(N_quark);
306  std::vector<Source *> source(N_quark);
307 
308  for (int iq = 0; iq < N_quark; ++iq) {
309  fopr_c[iq] = new Fopr_Clover(str_gmset_type);
310  fopr_c[iq]->set_parameters(*params_clover[iq]);
311  fopr_c[iq]->set_config(Ufix);
312 
313  solver[iq] = Solver::New(str_solver_type, fopr_c[iq]);
314  solver[iq]->set_parameters(*params_solver);
315 
316  fprop_lex[iq] = new Fprop_Standard_lex(solver[iq]);
317 
318  source[iq] = Source::New(str_source_type[iq]);
319  source[iq]->set_parameters(*params_source[iq]);
320  }
321  vout.general(vl, "\n");
322 
323  FieldIO *field_io = new FieldIO_Binary(IO_Format::Trivial);
324 
325 
326  // #### Execution main part ####
327  Field_F xq, b;
328  b = 0.0;
329 
330  int Nconv;
331  double diff, diff2;
332 
333  for (int iq = 0; iq < N_quark; ++iq) {
334  vout.general(vl, "Solving quark propagator, flavor = %d:\n", iq + 1);
335  vout.general(vl, " color spin Nconv diff diff2\n");
336 
337  for (int ispin = 0; ispin < Nd; ++ispin) {
338  for (int icolor = 0; icolor < Nc; ++icolor) {
339  int idx = icolor + Nc * ispin;
340  source[iq]->set(b, idx);
341 
342  fprop_lex[iq]->invert_D(xq, b, Nconv, diff);
343 
344  Field_F y(b);
345  fopr_c[iq]->set_mode("D");
346  y -= (Field_F)fopr_c[iq]->mult(xq);
347  diff2 = y.norm2();
348 
349  vout.general(vl, " %2d %2d %6d %12.4e %12.4e\n",
350  icolor, ispin, Nconv, diff, diff2);
351 
352  // fieldio->write_file(&xq, generate_filename(savefile_base[iq], idx));
353  field_io->write_file(&xq, generate_filename(savefile_base[iq], icolor, ispin));
354  }
355  }
356 
357  vout.general(vl, "\n");
358  }
359 
360 
361  // #### tidy up ####
362  delete gconf_read;
363  delete U;
364  delete Ufix;
365 
366  delete dr_smear;
367 
368  delete rand;
369  delete staple;
370 
371  for (int iq = 0; iq < N_quark; ++iq) {
372  delete fopr_c[iq];
373  delete solver[iq];
374  delete fprop_lex[iq];
375  delete source[iq];
376  }
377 
378  delete field_io;
379 
380  delete gfix;
381  delete proj;
382  delete smear;
383 
384  delete params_pre;
385  delete params_test;
386  delete params_gfix;
387  delete params_proj;
388  delete params_smear;
389  delete params_dr_smear;
390  delete params_solver;
391 
392  for (int iq = 0; iq < N_quark; ++iq) {
393  delete params_clover[iq];
394  delete params_source[iq];
395  delete params_quark[iq];
396  }
397 
398  delete params_all;
399 
400 
401  return EXIT_SUCCESS;
402  }
403 
404 
405  //====================================================================
406  int calculate_hadron_correlator()
407  {
408  // #### parameter setup ####
409  int Nc = CommonParameters::Nc();
410  int Nd = CommonParameters::Nd();
411  int Ndim = CommonParameters::Ndim();
412  int Nvol = CommonParameters::Nvol();
413 
414  ParameterManager_YAML params_manager;
415 
416  Parameters *params_test = new Parameters_Test_Spectrum_Clover;
417 
418  Parameters *params_pre = new Parameters;
419 
420  params_pre->Register_Parameters("Test_Spectrum_Clover", params_test);
421 
422  params_manager.read_params(filename_input, params_pre);
423 
424 
425  const int N_quark = params_test->get_int("number_of_valence_quarks");
426  std::vector<Parameters *> params_quark(N_quark);
427 
428  for (int iq = 0; iq < N_quark; ++iq) {
429  params_quark[iq] = new Parameters_Quark;
430  }
431 
432  Parameters *params_all = new Parameters;
433 
434 
435  std::ostringstream oss;
436 
437  for (int iq = 0; iq < N_quark; ++iq) {
438  oss.str("");
439  oss << "Quark_" << iq + 1;
440  string str = oss.str();
441  params_all->Register_Parameters(str.c_str(), params_quark[iq]);
442  }
443 
444  params_manager.read_params(filename_input, params_all);
445 
446  const string str_vlevel = params_test->get_string("verbose_level");
447 #ifdef USE_TEST
448  const double expected_result = params_test->get_double("expected_result");
449 #endif
450 
451  std::vector<std::string> savefile_base(N_quark);
452  std::vector<Parameters *> params_clover(N_quark);
453 
454  for (int iq = 0; iq < N_quark; ++iq) {
455  savefile_base[iq] = params_quark[iq]->get_string("temporary_filename_base");
456  params_clover[iq] = params_quark[iq]->get_Parameters("Fopr_Clover");
457  }
458 
459  Bridge::VerboseLevel vl = vout.set_verbose_level(str_vlevel);
460 
461  //- print input parameters
462  vout.general(vl, " vlevel = %s\n", str_vlevel.c_str());
463 
464  // NB. all str_gmset_type are supposed to be the same.
465  string str_gmset_type = params_clover[0]->get_string("gamma_matrix_type");
466  vout.general(vl, " gmset_type = %s\n", str_gmset_type.c_str());
467 
468 
469  // #### Set up objects #####
470  GammaMatrixSet *gmset = GammaMatrixSet::New(str_gmset_type);
471 
472  FieldIO *field_io = new FieldIO_Binary(IO_Format::Trivial);
473 
474 
475  // #### Execution main part ####
476  typedef std::valarray<Field_F> PropagatorSet;
477 
478  std::vector<PropagatorSet> sq(N_quark);
479  for (int iq = 0; iq < N_quark; ++iq) {
480  sq[iq].resize(Nc * Nd);
481 
482  for (int i = 0; i < Nc * Nd; ++i) {
483  sq[iq][i] = 0.0;
484  }
485  }
486 
487 
488  for (int iq = 0; iq < N_quark; ++iq) {
489  for (int ispin = 0; ispin < Nd; ++ispin) {
490  for (int icolor = 0; icolor < Nc; ++icolor) {
491  int idx = icolor + Nc * ispin;
492 
493  field_io->read_file(&sq[iq][idx], generate_filename(savefile_base[iq], icolor, ispin));
494  }
495  }
496  }
497 
498 
499  //- meson correlators
500  std::ofstream log_file;
501  if (filename_output != "stdout") {
502  log_file.open(filename_output.c_str());
503  vout.init(log_file);
504  }
505 
506  vout.general(vl, "2-point correlator:\n");
507  Corr2pt_4spinor corr(gmset);
508  valarray<double> result(N_quark);
509 
510  //- case(iq_1 == iq_2)
511  for (int iq = 0; iq < N_quark; ++iq) {
512  vout.general(vl, "Flavor combination = %d, %d\n", iq + 1, iq + 1);
513  result[iq] = corr.meson_all(sq[iq], sq[iq]);
514  vout.general(vl, "\n");
515  }
516 
517 
518  //- case(iq_1 < iq_2)
519  for (int iq = 0; iq < N_quark; ++iq) {
520  for (int jq = iq + 1; jq < N_quark; ++jq) {
521  vout.general(vl, "Flavor combination = %d, %d\n", iq + 1, jq + 1);
522  double result_2 = corr.meson_all(sq[iq], sq[jq]);
523  vout.general(vl, "\n");
524  }
525  }
526 
527  if (filename_output != "stdout") {
528  log_file.close();
529  vout.init(std::cout);
530  }
531 
532 
533  // #### tidy up ####
534  delete field_io;
535 
536  delete gmset;
537 
538  delete params_pre;
539  delete params_test;
540 
541  for (int iq = 0; iq < N_quark; ++iq) {
542  delete params_clover[iq];
543  delete params_quark[iq];
544  }
545 
546  delete params_all;
547 
548 
549 #ifdef USE_TEST
550  return Test::verify(expected_result, result[0]);
551 
552 #else
553  return EXIT_SUCCESS;
554 #endif
555  }
556 } // namespace Test_Spectrum_Clover