Bridge++  Ver. 2.0.2
test_Spectrum_2ptFunction_eo_withFileIO.cpp
Go to the documentation of this file.
1 
14 #include "test.h"
15 
16 #include "Field/shiftField_lex.h"
17 
18 #include "Fopr/fopr_Smeared.h"
19 
20 #include "IO/gaugeConfig.h"
21 
26 
27 #include "Tools/filename.h"
28 #include "Tools/file_utils.h"
29 #include "Tools/gammaMatrixSet.h"
31 
32 //====================================================================
34 
55 namespace Test_Spectrum {
56  const std::string test_name = "Spectrum.Hadron2ptFunction_eo_withFileIO";
57 
58  //- test-private parameters
59  namespace {
60  // const std::string filename_input = "test_Spectrum_Clover_Hadron2ptFunction.yaml";
61  }
62 
63  //- prototype declaration
64  int hadron_2ptFunction_eo_withFileIO(const std::string&);
65 
66  int calculate_quark_propagator_eo_withFileIO(const std::string&);
67  int calculate_hadron_correlator_eo_withFileIO(const std::string&);
68 
69  //- hadron_2ptFunction for various fermions
71  {
72  return hadron_2ptFunction_eo_withFileIO("test_Spectrum_Clover_Hadron2ptFunction.yaml");
73  }
74 
75 
77  {
78  return hadron_2ptFunction_eo_withFileIO("test_Spectrum_Wilson_Hadron2ptFunction.yaml");
79  }
80 
81 
82 #ifdef USE_TESTMANAGER_AUTOREGISTER
83  namespace {
84 #if defined(USE_GROUP_SU2)
85  // Nc=2 is not available.
86 #else
87  static const bool is_registered_Clover = TestManager::RegisterTest(
88  "Spectrum.Clover.Hadron2ptFunction_eo_withFileIO",
90  );
91 
92  //- NB. test_Spectrum_Wilson is implemented separately for beginners
93  // static const bool is_registered_Wilson = TestManager::RegisterTest(
94  // "Spectrum.Wilson.Hadron2ptFunction_eo_withFileIO",
95  // hadron_2ptFunction_eo_withFileIO_Wilson
96  // );
97 #endif
98  }
99 #endif
100 
101  //====================================================================
102  int hadron_2ptFunction_eo_withFileIO(const std::string& filename_input)
103  {
104  int result = 0;
105 
106  Timer timer(test_name);
107 
108 
109  // #### Execution main part ####
110  timer.start();
111 
112  result += calculate_quark_propagator_eo_withFileIO(filename_input);
113  result += calculate_hadron_correlator_eo_withFileIO(filename_input);
114 
115  timer.report();
116 
117 
118  return result;
119  }
120 
121 
122  //====================================================================
123  int calculate_quark_propagator_eo_withFileIO(const std::string& filename_input)
124  {
125  // #### parameter setup ####
126  const int Nc = CommonParameters::Nc();
127  const int Nd = CommonParameters::Nd();
128  const int Ndim = CommonParameters::Ndim();
129  const int Nvol = CommonParameters::Nvol();
130 
131  const Parameters params_all = ParameterManager::read(filename_input);
132 
133  const Parameters params_test = params_all.lookup("Test_Spectrum");
134 
135  const int N_quark = params_test.get_int("number_of_valence_quarks");
136 
137  const std::string str_gconf_status = params_test.get_string("gauge_config_status");
138  const std::string str_gconf_read = params_test.get_string("gauge_config_type_input");
139  const std::string readfile = params_test.get_string("config_filename_input");
140  const vector<int> Nshift_origin = params_test.get_int_vector("shift_origin");
141  const std::string str_rand_type = params_test.get_string("random_number_type");
142  const unsigned long seed = params_test.get_unsigned_long("seed_for_random_number");
143  const std::string str_vlevel = params_test.get_string("verbose_level");
144 
145  const std::string str_gfix_type = params_all.lookup("GaugeFixing").get_string("gauge_fixing_type");
146  const std::string str_proj_type = params_all.lookup("Projection").get_string("projection_type");
147  const std::string str_smear_type = params_all.lookup("Smear").get_string("smear_type");
148  const int Nsmear = params_all.lookup("Director_Smear").get_int("number_of_smearing");
149 
150  const Bridge::VerboseLevel vl = vout.set_verbose_level(str_vlevel);
151 
152 
153  //- print input parameters
154  vout.general(vl, " gconf_status = %s\n", str_gconf_status.c_str());
155  vout.general(vl, " gconf_read = %s\n", str_gconf_read.c_str());
156  vout.general(vl, " readfile = %s\n", readfile.c_str());
157  for (int mu = 0; mu < Ndim; ++mu) {
158  vout.general(vl, " shift_origin[%d] = %d\n", mu, Nshift_origin[mu]);
159  }
160  vout.general(vl, " rand_type = %s\n", str_rand_type.c_str());
161  vout.general(vl, " seed = %lu\n", seed);
162  vout.general(vl, " vlevel = %s\n", str_vlevel.c_str());
163  vout.general(vl, " gfix_type = %s\n", str_gfix_type.c_str());
164  vout.general(vl, " proj_type = %s\n", str_proj_type.c_str());
165  vout.general(vl, " smear_type = %s\n", str_smear_type.c_str());
166 
167  vector<Parameters> params_quark;
168  for (int iq = 0; iq < N_quark; ++iq) {
169  std::string qlabel = Filename("Quark_{id}").format(iq + 1);
170  params_quark.push_back(params_all.lookup(qlabel));
171  }
172 
173  // NB. all str_gmset_type are supposed to be the same.
174  const string str_gmset_type = params_quark[0].lookup("Fopr").get_string("gamma_matrix_type");
175  vout.general(vl, " gmset_type = %s\n", str_gmset_type.c_str());
176 
177  std::vector<std::string> savefile_base(N_quark);
178  for (int iq = 0; iq < N_quark; ++iq) {
179  savefile_base[iq] = params_quark[iq].get_string("temporary_filename_base");
180  }
181 
182  for (int iq = 0; iq < N_quark; ++iq) {
183  std::string str_solver_type = params_quark[iq].lookup("Solver").get_string("solver_type");
184  std::string str_source_type = params_quark[iq].lookup("Source").get_string("source_type");
185  vector<int> source_position = params_quark[iq].lookup("Source").get_int_vector("source_position");
186 
187  vout.general(vl, " Quark_%d:\n", iq + 1);
188  vout.general(vl, " savefile_base[iq] = %s[%d]\n", savefile_base[iq].c_str(), iq);
189  vout.general(vl, " solver_type = %s\n", str_solver_type.c_str());
190  vout.general(vl, " source_type = %s\n", str_source_type.c_str());
191  for (int mu = 0; mu < Ndim; ++mu) {
192  vout.general(vl, " source_position[%d] = %d\n", mu, source_position[mu]);
193  }
194  }
195  vout.general(vl, "\n");
196 
197 
198  //- input parameter check
199  int err = 0;
200  err += ParameterCheck::non_NULL(str_gconf_status);
201 
202  if (err) {
203  vout.crucial(vl, "Error at %s: input parameters have not been set\n", test_name.c_str());
204  exit(EXIT_FAILURE);
205  }
206 
207  for (int iq = 0; iq < N_quark; ++iq) {
208  std::string str_solver_type = params_quark[iq].lookup("Solver").get_string("solver_type");
209 
210  if (str_solver_type == "CG") {
211  vout.crucial(vl, "Error at %s: CG can not be adopted. Use CGNE,CGNR, instead.\n", test_name.c_str());
212  exit(EXIT_FAILURE);
213  }
214  }
215 
216  if ((Nsmear > 0) && (str_proj_type == "Stout_SU3")) {
217  if (CommonParameters::Nc() != 3) {
218  vout.crucial(vl, "check skipped: Nc = 3 is needed, but Nc = %d.\n\n", CommonParameters::Nc());
219  //- NB. EXIT_SKIP is set in calculate_hadron_correlator_withFileIO()
220  // return EXIT_SKIP;
221  return EXIT_SUCCESS;
222  }
223  }
224 
225  if ((str_gfix_type == "Coulomb") || (str_gfix_type == "Landau")) {
226  if (CommonParameters::Nc() != 3) {
227  vout.crucial(vl, "check skipped: Nc = 3 is needed, but Nc = %d.\n\n", CommonParameters::Nc());
228  return EXIT_SKIP;
229  }
230  }
231 
232 
233  RandomNumberManager::initialize(str_rand_type, seed);
234 
235 
236  // #### Set up a gauge configuration ####
237  Field_G U(Nvol, Ndim);
238 
239  if (str_gconf_status == "Continue") {
240  GaugeConfig(str_gconf_read).read(U, readfile);
241  } else if (str_gconf_status == "Cold_start") {
242  GaugeConfig("Unit").read(U);
243  } else if (str_gconf_status == "Hot_start") {
244  GaugeConfig("Random").read(U);
245  } else {
246  vout.crucial(vl, "Error at %s: unsupported gconf status \"%s\"\n", test_name.c_str(), str_gconf_status.c_str());
247  exit(EXIT_FAILURE);
248  }
249 
250  {
252  for (int i_dir = 0; i_dir < Ndim; ++i_dir) {
253  for (int i_shift = 0; i_shift < Nshift_origin[i_dir]; ++i_shift) {
254  Field_G Ushift(Nvol, Ndim);
255  shift.backward(Ushift, U, i_dir);
256  copy(U, Ushift);
257  }
258  }
259  }
260 
261  if (Nsmear > 0) {
262  unique_ptr<Projection> proj(Projection::New(str_proj_type, params_all.lookup("Projection")));
263  unique_ptr<Smear> smear(Smear::New(str_smear_type, proj.get(), params_all.lookup("Smear")));
264  unique_ptr<Director_Smear> dr_smear(new Director_Smear(smear.get(), params_all.lookup("Director_Smear")));
265  dr_smear->set_config(&U);
266 
267  const Field_G *Usmear = (Field_G *)dr_smear->getptr_smearedConfig(Nsmear);
268  copy(U, *Usmear);
269  }
270 
271  // #### Gauge fixing ####
272  {
273  Field_G Ufix(Nvol, Ndim);
274 
275  unique_ptr<GaugeFixing> gfix(GaugeFixing::New(str_gfix_type, params_all.lookup("GaugeFixing")));
276 
277  gfix->fix(Ufix, U);
278  copy(U, Ufix);
279  }
280 
281 
282  // #### object setup #####
283  struct QuarkType
284  {
285  unique_ptr<Fopr> fopr;
286  unique_ptr<Fopr> fopr_eo;
287  unique_ptr<Solver> solver;
288  unique_ptr<Fprop> fprop_eo;
289  unique_ptr<Source> source;
290  };
291  // NB. Fopr is used only for check of diff2 below.
292 
293  std::vector<QuarkType> quark(N_quark);
294 
295  for (int iq = 0; iq < N_quark; ++iq) {
296  const Parameters& params_fopr = params_quark[iq].lookup("Fopr");
297  const Parameters& params_solver = params_quark[iq].lookup("Solver");
298  const Parameters& params_source = params_quark[iq].lookup("Source");
299 
300  quark[iq].fopr.reset(Fopr::New(params_fopr.get_string("fermion_type"),
301  params_fopr));
302  quark[iq].fopr->set_config(&U);
303 
304  quark[iq].fopr_eo.reset(Fopr::New(params_fopr.get_string("fermion_type") + "_eo",
305  params_fopr));
306  quark[iq].fopr_eo->set_config(&U);
307 
308  quark[iq].solver.reset(Solver::New(params_solver.get_string("solver_type"),
309  quark[iq].fopr_eo.get(),
310  params_solver));
311 
312  quark[iq].fprop_eo.reset(new Fprop_Standard_eo(quark[iq].solver.get()));
313 
314  quark[iq].source.reset(Source::New(params_source.get_string("source_type"),
315  params_source));
316  }
317  vout.general(vl, "\n");
318 
319  unique_ptr<FieldIO> field_io(new FieldIO_Binary(IO_Format::Trivial));
320 
321 
322 
323  // #### Execution main part ####
324  for (int iq = 0; iq < N_quark; ++iq) {
325  vout.general(vl, "Solving quark propagator, flavor = %d:\n", iq + 1);
326  vout.general(vl, " color spin Nconv diff diff2\n");
327 
328  for (int ispin = 0; ispin < Nd; ++ispin) {
329  for (int icolor = 0; icolor < Nc; ++icolor) {
330  int i_cd = icolor + Nc * ispin;
331 
332  Field_F b; // b.set(0.0);
333  quark[iq].source->set(b, i_cd);
334 
335  Field_F xq;
336  int Nconv;
337  double diff;
338  quark[iq].fprop_eo->invert_D(xq, b, Nconv, diff);
339 
340  Field_F y(b);
341  quark[iq].fopr->set_mode("D");
342  quark[iq].fopr->mult(y, xq); // y = fopr[iq]->mult(xq);
343  axpy(y, -1.0, b); // y -= b;
344  double diff2 = y.norm2() / b.norm2();
345 
346  vout.general(vl, " %2d %2d %6d %12.4e %12.4e\n",
347  icolor, ispin, Nconv, diff, diff2);
348 
349  std::string filename = FileUtils::generate_filename("%s_c%d_s%d.dat", savefile_base[iq].c_str(), (icolor + 1), (ispin + 1));
350  field_io->write_file(xq, filename);
351  }
352  }
353 
354  vout.general(vl, "\n");
355  }
356 
357  return EXIT_SUCCESS;
358  }
359 
360 
361  //====================================================================
362  int calculate_hadron_correlator_eo_withFileIO(const std::string& filename_input)
363  {
364  // #### parameter setup ####
365  const int Nc = CommonParameters::Nc();
366  const int Nd = CommonParameters::Nd();
367  const int Ndim = CommonParameters::Ndim();
368  const int Nvol = CommonParameters::Nvol();
369 
370  const Parameters params_all = ParameterManager::read(filename_input);
371 
372  const Parameters params_test = params_all.lookup("Test_Spectrum");
373 
374  const int N_quark = params_test.get_int("number_of_valence_quarks");
375 
376  const std::string str_vlevel = params_test.get_string("verbose_level");
377 
378  const bool do_check = params_test.is_set("expected_result");
379  const double expected_result = do_check ? params_test.get_double("expected_result") : 0.0;
380 
381  const std::string str_proj_type = params_all.lookup("Projection").get_string("projection_type");
382  const int Nsmear = params_all.lookup("Director_Smear").get_int("number_of_smearing");
383 
384  std::vector<std::string> savefile_base(N_quark);
385 
386  for (int iq = 0; iq < N_quark; ++iq) {
387  savefile_base[iq] = params_all.lookup(Filename("Quark_{id}").format(iq + 1)).get_string("temporary_filename_base");
388  }
389 
390  const Bridge::VerboseLevel vl = vout.set_verbose_level(str_vlevel);
391 
392  //- print input parameters
393  vout.general(vl, " vlevel = %s\n", str_vlevel.c_str());
394 
395  // NB. all str_gmset_type are supposed to be the same.
396  const std::string str_gmset_type = params_all.lookup("Quark_1").lookup("Fopr").get_string("gamma_matrix_type");
397  vout.general(vl, " gmset_type = %s\n", str_gmset_type.c_str());
398 
399  std::vector<std::string> str_source_type(N_quark);
400 
401  for (int iq = 0; iq < N_quark; ++iq) {
402  vout.general(vl, " Quark_%d:\n", iq + 1);
403  vout.general(vl, " savefile_base[iq] = %s[%d]\n", savefile_base[iq].c_str(), iq);
404  }
405  vout.general(vl, "\n");
406 
407  if ((Nsmear > 0) && (str_proj_type == "Stout_SU3")) {
408  if (CommonParameters::Nc() != 3) {
409  vout.crucial(vl, "check skipped: Nc = 3 is needed, but Nc = %d.\n\n", CommonParameters::Nc());
410  return EXIT_SKIP;
411  }
412  }
413 
414 
415  // #### Set up objects #####
416  unique_ptr<GammaMatrixSet> gmset(GammaMatrixSet::New(str_gmset_type));
417 
418  unique_ptr<FieldIO> field_io(new FieldIO_Binary(IO_Format::Trivial));
419 
420  Corr2pt_4spinor corr(gmset.get(), params_all.lookup("Corr2pt_4spinor"));
421 
422 
423  // #### Execution main part ####
424  typedef std::vector<Field_F> PropagatorSet;
425  std::vector<PropagatorSet> sq(N_quark);
426  for (int iq = 0; iq < N_quark; ++iq) {
427  sq[iq].resize(Nc * Nd);
428 
429  for (int i_cd = 0; i_cd < Nc * Nd; ++i_cd) {
430  sq[iq][i_cd].set(0.0);
431  }
432  }
433 
434  for (int iq = 0; iq < N_quark; ++iq) {
435  for (int ispin = 0; ispin < Nd; ++ispin) {
436  for (int icolor = 0; icolor < Nc; ++icolor) {
437  int i_cd = icolor + Nc * ispin;
438 
439  std::string filename = FileUtils::generate_filename("%s_c%d_s%d.dat", savefile_base[iq].c_str(), (icolor + 1), (ispin + 1));
440  field_io->read_file(sq[iq][i_cd], filename);
441  }
442  }
443  }
444 
445 
446  //- meson correlators
447  std::vector<double> result(N_quark);
448 
449  vout.general(vl, "2-point correlator:\n");
450 
451  //- case(iq_1 == iq_2)
452  for (int iq = 0; iq < N_quark; ++iq) {
453  vout.general(vl, "Flavor combination = %d, %d\n", iq + 1, iq + 1);
454  result[iq] = corr.meson_all(sq[iq], sq[iq]);
455  vout.general(vl, "\n");
456  }
457 
458 
459  //- case(iq_1 < iq_2)
460  for (int iq = 0; iq < N_quark; ++iq) {
461  for (int jq = iq + 1; jq < N_quark; ++jq) {
462  vout.general(vl, "Flavor combination = %d, %d\n", iq + 1, jq + 1);
463  double result_2 = corr.meson_all(sq[iq], sq[jq]);
464  vout.general(vl, "\n");
465  }
466  }
467 
469 
470 
471  if (do_check) {
472  return Test::verify(result[0], expected_result);
473  } else {
474  vout.detailed(vl, "check skipped: expected_result not set.\n\n");
475  return EXIT_SKIP;
476  }
477  }
478 } // namespace Test_Spectrum
Test::verify
int verify(const double result, const double expected, double eps)
Definition: test.cpp:27
ParameterCheck::non_NULL
int non_NULL(const std::string v)
Definition: parameterCheck.cpp:65
Test_Spectrum::test_name
const std::string test_name
Definition: test_Spectrum_2ptFunction.cpp:72
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
Parameters
Class for parameters.
Definition: parameters.h:46
Test_Spectrum::calculate_quark_propagator_eo_withFileIO
int calculate_quark_propagator_eo_withFileIO(const std::string &)
Definition: test_Spectrum_2ptFunction_eo_withFileIO.cpp:123
Corr2pt_4spinor::meson_all
double meson_all(const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2)
Definition: corr2pt_4spinor.cpp:53
shiftField_lex.h
Parameters::get_int
int get_int(const string &key) const
Definition: parameters.cpp:192
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:219
GaugeConfig::read
void read(Field_G &U, const string &filename=string())
Definition: gaugeConfig.cpp:121
gaugeConfig.h
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:380
filename.h
Test_Spectrum::hadron_2ptFunction_eo_withFileIO_Clover
int hadron_2ptFunction_eo_withFileIO_Clover()
Definition: test_Spectrum_2ptFunction_eo_withFileIO.cpp:70
Parameters::get_int_vector
vector< int > get_int_vector(const string &key) const
Definition: parameters.cpp:267
source.h
Timer
Definition: timer.h:31
RandomNumberManager::finalize
static void finalize()
Definition: randomNumberManager.cpp:80
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:212
IO_Format::Trivial
const Format * Trivial
Definition: io_format.cpp:27
fprop_Standard_eo.h
RandomNumberManager::initialize
static bool initialize(const std::string &rng_type, unsigned long seed)
Definition: randomNumberManager.cpp:57
Field::norm2
double norm2() const
Definition: field.cpp:113
Test_Spectrum::hadron_2ptFunction_eo_withFileIO_Wilson
int hadron_2ptFunction_eo_withFileIO_Wilson()
Definition: test_Spectrum_2ptFunction_eo_withFileIO.cpp:76
Timer::start
void start()
Definition: timer.cpp:44
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
ParameterCheck::vl
Bridge::VerboseLevel vl
Definition: parameterCheck.cpp:18
file_utils.h
test.h
ParameterManager::read
static void read(const std::string &params_file, Parameters &params)
Definition: parameterManager.cpp:33
FieldIO_Binary
FieldIO_Binary class for file I/O of Field data in binary format.
Definition: fieldIO_Binary.h:47
Parameters::get_double
double get_double(const string &key) const
Definition: parameters.cpp:175
Corr2pt_4spinor
Two-point correlator for Wilson-type fermions.
Definition: corr2pt_4spinor.h:42
Director_Smear
Manager of smeared configurations.
Definition: director_Smear.h:39
Test_Solver_Wilson::solver
int solver(const std::string &)
Definition: test_Solver_Wilson.cpp:134
ShiftField_lex
Methods to shift a field in the lexical site index.
Definition: shiftField_lex.h:39
Parameters::get_unsigned_long
unsigned long get_unsigned_long(const string &key) const
Definition: parameters.cpp:209
EXIT_SKIP
#define EXIT_SKIP
Definition: test.h:17
FileUtils::generate_filename
std::string generate_filename(const char *fmt,...)
Definition: file_utils.cpp:17
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
randomNumberManager.h
Test_Spectrum
Test of spectroscopy.
Definition: test_Spectrum_2ptFunction.cpp:71
GaugeConfig
GaugeConfig class for file I/O of gauge configuration.
Definition: gaugeConfig.h:80
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
Parameters::is_set
bool is_set(const string &key) const
Definition: parameters.cpp:525
Fprop_Standard_eo
Get quark propagator for Fopr with even-odd site index.
Definition: fprop_Standard_eo.h:38
Field_F
Wilson-type fermion field.
Definition: field_F.h:37
fopr_Smeared.h
Parameters::get_string
string get_string(const string &key) const
Definition: parameters.cpp:221
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:180
Test_Spectrum::hadron_2ptFunction_eo_withFileIO
int hadron_2ptFunction_eo_withFileIO(const std::string &)
Definition: test_Spectrum_2ptFunction_eo_withFileIO.cpp:102
Test_Spectrum::calculate_hadron_correlator_eo_withFileIO
int calculate_hadron_correlator_eo_withFileIO(const std::string &)
Definition: test_Spectrum_2ptFunction_eo_withFileIO.cpp:362
Bridge::VerboseLevel
VerboseLevel
Definition: bridgeIO.h:42
Field_G
SU(N) gauge field.
Definition: field_G.h:38
gaugeFixing.h
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
corr2pt_4spinor.h
Filename::format
std::string format()
Definition: filename.h:61
Filename
Filename utility.
Definition: filename.h:43
gammaMatrixSet.h
Parameters::lookup
Parameters lookup(const string &key) const
Definition: parameters.h:79
Timer::report
void report(const Bridge::VerboseLevel vl=Bridge::GENERAL)
Definition: timer.cpp:128
Test_Gauge::shift
int shift(void)
Definition: test_Gauge_Shift.cpp:58
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:512
TestManager::RegisterTest
static bool RegisterTest(const std::string &key, const Test_function func)
Definition: testManager.h:69