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