Bridge++  Ver. 2.0.2
test_Spectrum_2ptFunction_eo.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 
22 #include "Measurements/Fermion/corr2pt_4spinor.h"
23 #include "Measurements/Fermion/fprop_Standard_eo.h"
24 #include "Measurements/Fermion/source.h"
25 #include "Measurements/Gauge/gaugeFixing.h"
26 
27 #include "Tools/filename.h"
28 #include "Tools/gammaMatrixSet.h"
29 #include "Tools/randomNumberManager.h"
30 
31 //====================================================================
33 
45 namespace Test_Spectrum {
46  const std::string test_name = "Spectrum.Hadron2ptFunction_eo";
47 
48  //- test-private parameters
49  namespace {
50  // const std::string filename_input = "test_Spectrum_Clover_Hadron2ptFunction.yaml";
51  }
52 
53  //- prototype declaration
54  int hadron_2ptFunction_eo(const std::string&);
55 
56  //- hadron_2ptFunction for various fermions
57  int hadron_2ptFunction_eo_Clover()
58  {
59  return hadron_2ptFunction_eo("test_Spectrum_Clover_Hadron2ptFunction.yaml");
60  }
61 
62 
63  int hadron_2ptFunction_eo_Wilson()
64  {
65  return hadron_2ptFunction_eo("test_Spectrum_Wilson_Hadron2ptFunction_2.yaml");
66  }
67 
68 
69 #ifdef USE_TESTMANAGER_AUTOREGISTER
70  namespace {
71 #if defined(USE_GROUP_SU2)
72  // Nc=2 is not available.
73 #else
74  static const bool is_registered_Clover_eo = TestManager::RegisterTest(
75  "Spectrum.Clover.Hadron2ptFunction_eo",
76  hadron_2ptFunction_eo_Clover
77  );
78 
79  //- NB. test_Spectrum_Wilson is implemented separately for beginners
80  // static const bool is_registered_Wilson_eo = TestManager::RegisterTest(
81  // "Spectrum.Wilson.Hadron2ptFunction_eo",
82  // hadron_2ptFunction_eo_Wilson
83  // );
84 #endif
85  }
86 #endif
87 
88  //====================================================================
89  int hadron_2ptFunction_eo(const std::string& filename_input)
90  {
91  // #### parameter setup ####
92  const int Nc = CommonParameters::Nc();
93  const int Nd = CommonParameters::Nd();
94  const int Ndim = CommonParameters::Ndim();
95  const int Nvol = CommonParameters::Nvol();
96 
97  const Parameters params_all = ParameterManager::read(filename_input);
98 
99  const Parameters params_test = params_all.lookup("Test_Spectrum");
100 
101  const int N_quark = params_test.get_int("number_of_valence_quarks");
102 
103  const std::string str_gconf_status = params_test.get_string("gauge_config_status");
104  const std::string str_gconf_read = params_test.get_string("gauge_config_type_input");
105  const vector<int> Nshift_origin = params_test.get_int_vector("shift_origin");
106  const std::string readfile = params_test.get_string("config_filename_input");
107  const std::string str_rand_type = params_test.get_string("random_number_type");
108  const unsigned long seed = params_test.get_unsigned_long("seed_for_random_number");
109  const std::string str_vlevel = params_test.get_string("verbose_level");
110 
111  const bool do_check = params_test.is_set("expected_result");
112  const double expected_result = do_check ? params_test.get_double("expected_result") : 0.0;
113 
114  const std::string str_gfix_type = params_all.lookup("GaugeFixing").get_string("gauge_fixing_type");
115  const std::string str_proj_type = params_all.lookup("Projection").get_string("projection_type");
116  const std::string str_smear_type = params_all.lookup("Smear").get_string("smear_type");
117  const int Nsmear = params_all.lookup("Director_Smear").get_int("number_of_smearing");
118 
119  const Bridge::VerboseLevel vl = vout.set_verbose_level(str_vlevel);
120 
121  //- print input parameters
122  vout.general(vl, " gconf_status = %s\n", str_gconf_status.c_str());
123  vout.general(vl, " gconf_read = %s\n", str_gconf_read.c_str());
124  vout.general(vl, " readfile = %s\n", readfile.c_str());
125  for (int mu = 0; mu < Ndim; ++mu) {
126  vout.general(vl, " shift_origin[%d] = %d\n", mu, Nshift_origin[mu]);
127  }
128  vout.general(vl, " rand_type = %s\n", str_rand_type.c_str());
129  vout.general(vl, " seed = %lu\n", seed);
130  vout.general(vl, " vlevel = %s\n", str_vlevel.c_str());
131  vout.general(vl, " gfix_type = %s\n", str_gfix_type.c_str());
132  vout.general(vl, " proj_type = %s\n", str_proj_type.c_str());
133  vout.general(vl, " smear_type = %s\n", str_smear_type.c_str());
134 
135  vector<Parameters> params_quark;
136  for (int iq = 0; iq < N_quark; ++iq) {
137  std::string qlabel = Filename("Quark_{id}").format(iq + 1);
138  params_quark.push_back(params_all.lookup(qlabel));
139  }
140 
141  // NB. all str_gmset_type are supposed to be the same.
142  const std::string str_gmset_type = params_quark[0].lookup("Fopr").get_string("gamma_matrix_type");
143  vout.general(vl, " gmset_type = %s\n", str_gmset_type.c_str());
144 
145  for (int iq = 0; iq < N_quark; ++iq) {
146  vout.general(vl, " Quark_%d:\n", iq + 1);
147  vout.general(vl, " solver_type = %s\n",
148  params_quark[iq].lookup("Solver").get_string("solver_type").c_str());
149  vout.general(vl, " source_type = %s\n",
150  params_quark[iq].lookup("Source").get_string("source_type").c_str());
151 
152  std::string str_solver_type = params_quark[iq].lookup("Solver").get_string("solver_type");
153  vout.general(vl, " solver_type = %s\n", str_solver_type.c_str());
154 
155  vector<int> source_position = params_quark[iq].lookup("Source").get_int_vector("source_position");
156  for (int mu = 0; mu < Ndim; ++mu) {
157  vout.general(vl, " source_position[%d] = %d\n", mu, source_position[mu]);
158  }
159  }
160  vout.general(vl, "\n");
161 
162  //- input parameter check
163  int err = 0;
164  err += ParameterCheck::non_NULL(str_gconf_status);
165 
166  if (err) {
167  vout.crucial(vl, "Error at %s: input parameters have not been set\n", test_name.c_str());
168  exit(EXIT_FAILURE);
169  }
170 
171  for (int iq = 0; iq < N_quark; ++iq) {
172  std::string str_solver_type = params_quark[iq].lookup("Solver").get_string("solver_type");
173 
174  if (str_solver_type == "CG") {
175  vout.crucial(vl, "Error at %s: CG can not be adopted. Use CGNE,CGNR, instead.\n", test_name.c_str());
176  exit(EXIT_FAILURE);
177  }
178  }
179 
180  if ((Nsmear > 0) && (str_proj_type == "Stout_SU3")) {
181  if (CommonParameters::Nc() != 3) {
182  vout.crucial(vl, "check skipped: Nc = 3 is needed, but Nc = %d.\n\n", CommonParameters::Nc());
183  return EXIT_SKIP;
184  }
185  }
186 
187  if ((str_gfix_type == "Coulomb") || (str_gfix_type == "Landau")) {
188  if (CommonParameters::Nc() != 3) {
189  vout.crucial(vl, "check skipped: Nc = 3 is needed, but Nc = %d.\n\n", CommonParameters::Nc());
190  return EXIT_SKIP;
191  }
192  }
193 
194 
195  RandomNumberManager::initialize(str_rand_type, seed);
196 
197 
198  // #### Set up a gauge configuration ####
199  Field_G U(Nvol, Ndim);
200 
201  if (str_gconf_status == "Continue") {
202  GaugeConfig(str_gconf_read).read(U, readfile);
203  } else if (str_gconf_status == "Cold_start") {
204  GaugeConfig("Unit").read(U);
205  } else if (str_gconf_status == "Hot_start") {
206  GaugeConfig("Random").read(U);
207  } else {
208  vout.crucial(vl, "Error at %s: unsupported gconf status \"%s\"\n", test_name.c_str(), str_gconf_status.c_str());
209  exit(EXIT_FAILURE);
210  }
211 
212  {
213  ShiftField_lex shift;
214  for (int i_dir = 0; i_dir < Ndim; ++i_dir) {
215  for (int i_shift = 0; i_shift < Nshift_origin[i_dir]; ++i_shift) {
216  Field_G Ushift(Nvol, Ndim);
217  shift.backward(Ushift, U, i_dir);
218  copy(U, Ushift);
219  }
220  }
221  }
222 
223  if (Nsmear > 0) {
224  unique_ptr<Projection> proj(Projection::New(str_proj_type, params_all.lookup("Projection")));
225  unique_ptr<Smear> smear(Smear::New(str_smear_type, proj.get(), params_all.lookup("Smear")));
226  unique_ptr<Director_Smear> dr_smear(new Director_Smear(smear.get(), params_all.lookup("Director_Smear")));
227  dr_smear->set_config(&U);
228 
229  const Field_G *Usmear = (Field_G *)dr_smear->getptr_smearedConfig(Nsmear);
230  copy(U, *Usmear);
231  }
232 
233  // #### Gauge fixing ####
234  {
235  Field_G Ufix(Nvol, Ndim);
236  unique_ptr<GaugeFixing> gfix(GaugeFixing::New(str_gfix_type, params_all.lookup("GaugeFixing")));
237 
238  gfix->fix(Ufix, U);
239  copy(U, Ufix);
240  }
241 
242 
243  // #### object setup #####
244  unique_ptr<GammaMatrixSet> gmset(GammaMatrixSet::New(str_gmset_type));
245 
246  struct QuarkType
247  {
248  unique_ptr<Fopr> fopr;
249  unique_ptr<Fopr> fopr_eo;
250  unique_ptr<Solver> solver;
251  unique_ptr<Fprop> fprop_eo;
252  unique_ptr<Source> source;
253  };
254  // NB. Fopr is used only for check of diff2 below.
255 
256  std::vector<QuarkType> quark(N_quark);
257 
258  for (int iq = 0; iq < N_quark; ++iq) {
259  const Parameters& params_fopr = params_quark[iq].lookup("Fopr");
260  const Parameters& params_solver = params_quark[iq].lookup("Solver");
261  const Parameters& params_source = params_quark[iq].lookup("Source");
262 
263  quark[iq].fopr.reset(Fopr::New(params_fopr.get_string("fermion_type"),
264  params_fopr));
265  quark[iq].fopr->set_config(&U);
266 
267  quark[iq].fopr_eo.reset(Fopr::New(params_fopr.get_string("fermion_type") + "_eo",
268  params_fopr));
269  quark[iq].fopr_eo->set_config(&U);
270 
271  quark[iq].solver.reset(Solver::New(params_solver.get_string("solver_type"),
272  quark[iq].fopr_eo.get(),
273  params_solver));
274 
275  quark[iq].fprop_eo.reset(new Fprop_Standard_eo(quark[iq].solver.get()));
276 
277  quark[iq].source.reset(Source::New(params_source.get_string("source_type"),
278  params_source));
279  }
280 
281  Corr2pt_4spinor corr(gmset.get(), params_all.lookup("Corr2pt_4spinor"));
282 
283  Timer timer(test_name);
284 
285 
286  // #### Execution main part ####
287  timer.start();
288 
289  typedef std::vector<Field_F> PropagatorSet;
290  std::vector<PropagatorSet> sq(N_quark);
291  for (int iq = 0; iq < N_quark; ++iq) {
292  sq[iq].resize(Nc * Nd);
293 
294  for (int i_cd = 0; i_cd < Nc * Nd; ++i_cd) {
295  sq[iq][i_cd].set(0.0);
296  }
297  }
298 
299  for (int iq = 0; iq < N_quark; ++iq) {
300  vout.general(vl, "Solving quark propagator, flavor = %d:\n", iq + 1);
301  vout.general(vl, " color spin Nconv diff diff2\n");
302 
303  for (int ispin = 0; ispin < Nd; ++ispin) {
304  for (int icolor = 0; icolor < Nc; ++icolor) {
305  int i_cd = icolor + Nc * ispin;
306 
307  Field_F b; // b.set(0.0);
308  quark[iq].source->set(b, i_cd);
309 
310  int Nconv;
311  double diff;
312  quark[iq].fprop_eo->invert_D(sq[iq][i_cd], b, Nconv, diff);
313 
314  Field_F y(b);
315  quark[iq].fopr->set_mode("D");
316  quark[iq].fopr->mult(y, sq[iq][i_cd]); // y = fopr[iq]->mult(sq[iq][i_cd]);
317  axpy(y, -1.0, b); // y -= b;
318  double diff2 = y.norm2() / b.norm2();
319 
320  vout.general(vl, " %2d %2d %6d %12.4e %12.4e\n",
321  icolor, ispin, Nconv, diff, diff2);
322  }
323  }
324 
325  vout.general(vl, "\n");
326  }
327 
328 
329  //- meson correlators
330  std::vector<double> result(N_quark);
331 
332  vout.general(vl, "2-point correlator:\n");
333 
334  //- case(iq_1 == iq_2)
335  for (int iq = 0; iq < N_quark; ++iq) {
336  vout.general(vl, "Flavor combination = %d, %d\n", iq + 1, iq + 1);
337  result[iq] = corr.meson_all(sq[iq], sq[iq]);
338  vout.general(vl, "\n");
339  }
340 
341 
342  //- case(iq_1 < iq_2)
343  for (int iq = 0; iq < N_quark; ++iq) {
344  for (int jq = iq + 1; jq < N_quark; ++jq) {
345  vout.general(vl, "Flavor combination = %d, %d\n", iq + 1, jq + 1);
346  double result_2 = corr.meson_all(sq[iq], sq[jq]);
347  vout.general(vl, "\n");
348  }
349  }
350 
351 
352  timer.report();
353 
354  RandomNumberManager::finalize();
355 
356 
357  if (do_check) {
358  return Test::verify(result[0], expected_result);
359  } else {
360  vout.detailed(vl, "check skipped: expected_result not set.\n\n");
361  return EXIT_SKIP;
362  }
363  }
364 } // 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
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
Parameters::get_int_vector
vector< int > get_int_vector(const string &key) const
Definition: parameters.cpp:267
Test_Spectrum::hadron_2ptFunction_eo_Clover
int hadron_2ptFunction_eo_Clover()
Definition: test_Spectrum_2ptFunction_eo.cpp:57
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
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
Timer::start
void start()
Definition: timer.cpp:44
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
Test_Spectrum::hadron_2ptFunction_eo_Wilson
int hadron_2ptFunction_eo_Wilson()
Definition: test_Spectrum_2ptFunction_eo.cpp:63
ParameterCheck::vl
Bridge::VerboseLevel vl
Definition: parameterCheck.cpp:18
test.h
ParameterManager::read
static void read(const std::string &params_file, Parameters &params)
Definition: parameterManager.cpp:33
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
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
Test_Spectrum::hadron_2ptFunction_eo
int hadron_2ptFunction_eo(const std::string &)
Definition: test_Spectrum_2ptFunction_eo.cpp:89
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
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