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