Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
test_HMC_Clover_Isochemical_RHMC_Nf2p1.cpp
Go to the documentation of this file.
1 
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 
22 #include "action_G_Plaq.h"
23 #include "action_G_Rectangle.h"
24 
25 #include "fopr_Clover_Isochemical.h"
26 #include "force_F_Clover_Nf2_Isochemical.h"
27 
28 #include "randomNumbers_Mseries.h"
29 
30 #include "director_Smear.h"
31 #include "fopr_Smeared.h"
32 #include "force_F_Smeared.h"
33 
34 #include "action_F_Standard_lex.h"
35 #include "fprop_Standard_lex.h"
36 #include "action_F_Ratio_lex.h"
37 
38 #include "fopr_Rational.h"
39 #include "force_F_Rational.h"
40 #include "math_Rational.h"
41 #include "action_F_Rational_Frame.h"
42 
43 #include "forceSmear.h"
44 #include "projection.h"
45 #include "smear.h"
46 #include "solver.h"
47 
48 #include "hmc_General.h"
49 #include "builder_Integrator.h"
50 #include "integrator.h"
51 
52 #include "polyakovLoop.h"
53 
54 #include "noiseVector_Z2.h"
55 #include "quarkNumberSusceptibility_Wilson.h"
56 
57 #ifdef USE_TEST
58 #include "test.h"
59 #endif
60 
61 #ifdef USE_TESTMANAGER_AUTOREGISTER
62 #include "testManager.h"
63 #endif
64 
65 //====================================================================
67 
78 namespace Test_HMC_Clover_Isochemical {
79  //- test-private parameters
80  namespace {
81  const std::string filename_input = "test_HMC_Clover_Isochemical_RHMC_Nf2p1.yaml";
82  const std::string filename_output = "stdout";
83 
84  class Parameters_Test_HMC_Clover_Isochemical : public Parameters {
85  public:
86  Parameters_Test_HMC_Clover_Isochemical()
87  {
88  Register_string("gauge_config_status", "NULL");
89 
90  Register_string("gauge_config_type_input", "NULL");
91  Register_string("config_filename_input", "NULL");
92 
93  Register_string("gauge_config_type_output", "NULL");
94  Register_string("config_filename_output", "NULL");
95 
96  Register_int("trajectory_number", 0);
97  Register_int("trajectory_number_step", 0);
98  Register_int("save_config_interval", 0);
99 
100  Register_int("seed_for_noise", 0);
101  Register_int("number_of_noises", 0);
102 
103  Register_string("verbose_level", "NULL");
104 
105  Register_double("expected_result", 0.0);
106  }
107  };
108  }
109 
110  //- prototype declaration
111  int RHMC_Nf2p1(void);
112 
113 #ifdef USE_TESTMANAGER_AUTOREGISTER
114  namespace {
115  static const bool is_registered = TestManager::RegisterTest(
116  "HMC.Clover_Isochemical.RHMC_Nf2p1",
117  RHMC_Nf2p1
118  );
119  }
120 #endif
121 
122  //====================================================================
123  int RHMC_Nf2p1(void)
124  {
125  // ##### parameter setup #####
126  int Nc = CommonParameters::Nc();
127  int Nvol = CommonParameters::Nvol();
128  int Ndim = CommonParameters::Ndim();
129  int NinG = 2 * Nc * Nc;
130 
131  Parameters *params_test = new Parameters_Test_HMC_Clover_Isochemical;
132  Parameters *params_action_G = ParametersFactory::New("Action.G_Rectangle");
133  Parameters *params_clover = ParametersFactory::New("Fopr.Clover_Isochemical");
134  Parameters *params_clover_prec = ParametersFactory::New("Fopr.Clover_Isochemical");
135  Parameters *params_clover1 = ParametersFactory::New("Fopr.Clover_Isochemical");
136  Parameters *params_proj = ParametersFactory::New("Projection");
137  Parameters *params_smear = ParametersFactory::New("Smear");
138  Parameters *params_dr_smear = ParametersFactory::New("Director_Smear");
139  Parameters *params_rational_MD = ParametersFactory::New("Fopr.Rational");
140  Parameters *params_rational_H = ParametersFactory::New("Fopr.Rational");
141  Parameters *params_solver_MD = ParametersFactory::New("Solver");
142  Parameters *params_solver_H = ParametersFactory::New("Solver");
143  Parameters *params_integrator = ParametersFactory::New("Builder_Integrator");
144  Parameters *params_hmc = ParametersFactory::New("HMC.General");
145 
146  Parameters *params_all = new Parameters;
147 
148  params_all->Register_Parameters("Test_HMC_Clover_Isochemical", params_test);
149  params_all->Register_Parameters("Action_G_Rectangle", params_action_G);
150  params_all->Register_Parameters("Fopr_Clover_Isochemical_Nf2", params_clover);
151  params_all->Register_Parameters("Fopr_Clover_Isochemical_Nf2_prec", params_clover_prec);
152  params_all->Register_Parameters("Fopr_Clover_Isochemical_Nf1", params_clover1);
153  params_all->Register_Parameters("Projection", params_proj);
154  params_all->Register_Parameters("Smear", params_smear);
155  params_all->Register_Parameters("Director_Smear", params_dr_smear);
156  params_all->Register_Parameters("Fopr_Rational_MD", params_rational_MD);
157  params_all->Register_Parameters("Fopr_Rational_H", params_rational_H);
158  params_all->Register_Parameters("Solver_MD", params_solver_MD);
159  params_all->Register_Parameters("Solver_H", params_solver_H);
160  params_all->Register_Parameters("Builder_Integrator", params_integrator);
161  params_all->Register_Parameters("HMC_General", params_hmc);
162 
163  ParameterManager_YAML params_manager;
164  params_manager.read_params(filename_input, params_all);
165 
166  const string str_gconf_status = params_test->get_string("gauge_config_status");
167  const string str_gconf_read = params_test->get_string("gauge_config_type_input");
168  const string readfile = params_test->get_string("config_filename_input");
169  const string str_gconf_write = params_test->get_string("gauge_config_type_output");
170  const string writefile = params_test->get_string("config_filename_output");
171  int iconf = params_test->get_int("trajectory_number");
172  const int Ntraj = params_test->get_int("trajectory_number_step");
173  const int i_save_conf = params_test->get_int("save_config_interval");
174  int i_seed_noise = params_test->get_int("seed_for_noise");
175  const string str_vlevel = params_test->get_string("verbose_level");
176 #ifdef USE_TEST
177  const double expected_result = params_test->get_double("expected_result");
178 #endif
179 
180  const string str_gmset_type = params_clover->get_string("gamma_matrix_type");
181  const string str_proj_type = params_proj->get_string("projection_type");
182  const string str_smear_type = params_smear->get_string("smear_type");
183  const string str_solver_MD_type = params_solver_MD->get_string("solver_type");
184  const string str_solver_H_type = params_solver_H->get_string("solver_type");
185 
186  Bridge::VerboseLevel vl = vout.set_verbose_level(str_vlevel);
187 
188  //- print input parameters
189  vout.general(vl, " gconf_status = %s\n", str_gconf_status.c_str());
190  vout.general(vl, " gconf_read = %s\n", str_gconf_read.c_str());
191  vout.general(vl, " readfile = %s\n", readfile.c_str());
192  vout.general(vl, " gconf_write = %s\n", str_gconf_write.c_str());
193  vout.general(vl, " writefile = %s\n", writefile.c_str());
194  vout.general(vl, " iconf = %d\n", iconf);
195  vout.general(vl, " Ntraj = %d\n", Ntraj);
196  vout.general(vl, " i_save_conf = %d\n", i_save_conf);
197  vout.general(vl, " vlevel = %s\n", str_vlevel.c_str());
198 
199  vout.general(vl, " gmset_type = %s\n", str_gmset_type.c_str());
200  vout.general(vl, " proj_type = %s\n", str_proj_type.c_str());
201  vout.general(vl, " smear_type = %s\n", str_smear_type.c_str());
202  vout.general(vl, " solver_MD_type = %s\n", str_solver_MD_type.c_str());
203  vout.general(vl, " solver_H_type = %s\n", str_solver_H_type.c_str());
204 
205  //- input parameter check
206  int err = 0;
207  err += ParameterCheck::non_NULL(str_gconf_read);
208  err += ParameterCheck::non_NULL(readfile);
209  err += ParameterCheck::non_NULL(str_gconf_write);
210  err += ParameterCheck::non_NULL(writefile);
211  err += ParameterCheck::non_zero(iconf);
212  err += ParameterCheck::non_zero(Ntraj);
213  err += ParameterCheck::non_zero(i_save_conf);
214 
215  if (err) {
216  vout.crucial(vl, "test_HMC_Clover_Isochemical: Input parameters have not been set.\n");
217  abort();
218  }
219 
220 
221  // ##### object setup #####
222  Field_G *U = new Field_G(Nvol, Ndim);
223  GaugeConfig *gconf_read = new GaugeConfig(str_gconf_read);
224 
225  if (str_gconf_status == "Continue") {
226  gconf_read->read_file((Field *)U, readfile);
227  } else if (str_gconf_status == "Start_cold") {
228  gconf_read->set_cold((Field *)U);
229  } else {
230  vout.crucial(vl, "Test_HMC_Clover_Isochemical: unsupported gconf status \"%s\".\n", str_gconf_status.c_str());
231  abort();
232  }
233 
234 
235  Action_G_Rectangle *action_G = new Action_G_Rectangle;
236  action_G->set_parameters(*params_action_G);
237 
238 
239  //-- N_f=2 part
240  Fopr_Clover_Isochemical *fopr_c = new Fopr_Clover_Isochemical(str_gmset_type);
241  fopr_c->set_parameters(*params_clover);
242 
243  Force_F_Clover_Nf2_Isochemical *force_fopr_c = new Force_F_Clover_Nf2_Isochemical(str_gmset_type);
244  force_fopr_c->set_parameters(*params_clover);
245 
246  Fopr_Clover_Isochemical *fopr_c_prec = new Fopr_Clover_Isochemical(str_gmset_type);
247  fopr_c_prec->set_parameters(*params_clover_prec);
248 
249  Force_F_Clover_Nf2_Isochemical *force_fopr_c_prec = new Force_F_Clover_Nf2_Isochemical(str_gmset_type);
250  force_fopr_c_prec->set_parameters(*params_clover_prec);
251 
252 
253  Projection *proj = Projection::New(str_proj_type);
254  Smear *smear = Smear::New(str_smear_type, proj);
255  smear->set_parameters(*params_smear);
256 
257  ForceSmear *force_smear = ForceSmear::New(str_smear_type, proj);
258  force_smear->set_parameters(*params_smear);
259 
260  Director_Smear *dr_smear = new Director_Smear((Smear *)smear);
261  dr_smear->set_parameters(*params_dr_smear);
262 
263  Fopr_Smeared *fopr_smear = new Fopr_Smeared((Fopr *)fopr_c, dr_smear);
264  Force_F_Smeared *force_fopr_smear
265  = new Force_F_Smeared((Force *)force_fopr_c, (ForceSmear *)force_smear, dr_smear);
266 
267  Fopr_Smeared *fopr_smear_prec = new Fopr_Smeared((Fopr *)fopr_c_prec, dr_smear);
268  Force_F_Smeared *force_fopr_smear_prec
269  = new Force_F_Smeared((Force *)force_fopr_c_prec, (ForceSmear *)force_smear, dr_smear);
270 
271 
272  Solver *solver_MD = Solver::New(str_solver_MD_type, fopr_smear);
273  solver_MD->set_parameters(*params_solver_MD);
274  Fprop *fprop_MD = new Fprop_Standard_lex(solver_MD);
275 
276  Solver *solver_H = Solver::New(str_solver_H_type, fopr_smear);
277  solver_H->set_parameters(*params_solver_H);
278  Fprop *fprop_H = new Fprop_Standard_lex(solver_H);
279 
280  Solver *solver_MD_prec = Solver::New(str_solver_MD_type, fopr_smear_prec);
281  solver_MD_prec->set_parameters(*params_solver_MD);
282  Fprop *fprop_MD_prec = new Fprop_Standard_lex(solver_MD_prec);
283 
284  Solver *solver_H_prec = Solver::New(str_solver_H_type, fopr_smear_prec);
285  solver_H_prec->set_parameters(*params_solver_H);
286  Fprop *fprop_H_prec = new Fprop_Standard_lex(solver_H_prec);
287 
288 
289  Action_F_Standard_lex *action_F_Nf2_prec
290  = new Action_F_Standard_lex((Fopr *)fopr_smear_prec, (Force *)force_fopr_smear_prec,
291  fprop_MD_prec, fprop_H_prec);
292 
293  Action_F_Ratio_lex *action_F_Nf2_d
294  = new Action_F_Ratio_lex((Fopr *)fopr_smear_prec, (Force *)force_fopr_smear_prec,
295  (Fopr *)fopr_smear, (Force *)force_fopr_smear,
296  fprop_H_prec,
297  fprop_MD, fprop_H);
298 
299 
300  //-- N_f=1 part
301  Fopr_Clover_Isochemical *fopr_c1 = new Fopr_Clover_Isochemical(str_gmset_type);
302  fopr_c1->set_parameters(*params_clover1);
303 
304  Force_F_Clover_Nf2_Isochemical *force_fopr_c1 = new Force_F_Clover_Nf2_Isochemical(str_gmset_type);
305  force_fopr_c1->set_parameters(*params_clover1);
306 
307  Fopr_Rational *fopr_r1 = new Fopr_Rational(fopr_c1);
308  fopr_r1->set_parameters(*params_rational_H);
309  Fopr_Smeared *fopr_langev = new Fopr_Smeared((Fopr *)fopr_r1, dr_smear);
310 
311  Fopr_Rational *fopr_r2 = new Fopr_Rational(fopr_c1);
312  fopr_r2->set_parameters(*params_rational_MD);
313  Fopr_Smeared *fopr_H = new Fopr_Smeared((Fopr *)fopr_r2, dr_smear);
314 
315  Force_F_Rational *force_fopr_r2 = new Force_F_Rational(fopr_c1, force_fopr_c1);
316  force_fopr_r2->set_parameters(*params_rational_MD);
317  Force_F_Smeared *force_fopr_MD
318  = new Force_F_Smeared((Force *)force_fopr_r2, (ForceSmear *)force_smear, dr_smear);
319 
320 
321  Action_F_Rational_frame *action_F_Nf1
322  = new Action_F_Rational_frame((Fopr *)fopr_langev, (Fopr *)fopr_H, (Force *)force_fopr_MD);
323 
324 
325  valarray<Action *> actions(4);
326  actions[0] = (Action *)action_F_Nf2_d;
327  actions[1] = (Action *)action_F_Nf2_prec;
328  actions[2] = (Action *)action_F_Nf1;
329  actions[3] = (Action *)action_G;
330 
331  valarray<Director *> directors(1);
332  directors[0] = (Director *)dr_smear;
333 
334  Builder_Integrator *builder = new Builder_Integrator(actions, directors);
335  builder->set_parameters(*params_integrator);
336  Integrator *integrator = builder->build();
337 
338 
339  //- Random number is initialized with a parameter specified by iconf
340  RandomNumbers *rand = new RandomNumbers_Mseries(iconf);
341 
342 
343  HMC_General hmc(actions, directors, integrator, rand);
344  hmc.set_parameters(*params_hmc);
345 
346 
347  PolyakovLoop *ploop = new PolyakovLoop();
348 
349  RandomNumbers *rand_nv = new RandomNumbers_Mseries(i_seed_noise);
350  NoiseVector *nv = new NoiseVector_Z2(rand_nv);
351  QuarkNumberSusceptibility_Wilson *quark_suscept
352  = new QuarkNumberSusceptibility_Wilson(fopr_smear, fprop_H, nv);
353  quark_suscept->set_parameters(*params_test);
354 
355 
356  // #### Execution main part ####
357  vout.general(vl, "HMC start: Ntraj = %d\n", Ntraj);
358 
359  double ploop0 = ploop->measure_ploop(*U);
360  vout.general(vl, "Polyakov loop = %f\n", ploop0);
361 
362  fopr_smear->set_config(U);
363  double result_quark_suscept0 = quark_suscept->measure();
364 
365  double result = 0.0;
366  for (int traj = 0; traj < Ntraj; ++traj) {
367  vout.general(vl, "\n");
368  vout.general(vl, "---------------------------------------------------\n");
369  vout.general(vl, "traj = %d\n", traj);
370 
371  result = hmc.update(*U);
372 
373  double ploop1 = ploop->measure_ploop(*U);
374  vout.general(vl, "Polyakov loop = %f\n", ploop1);
375 
376  fopr_smear->set_config(U);
377  double result_quark_suscept1 = quark_suscept->measure();
378  }
379 
380 
381  // ##### tidy up #####
382  delete U;
383  delete gconf_read;
384 
385  delete action_G;
386 
387  delete fopr_c;
388  delete force_fopr_c;
389  delete fopr_c_prec;
390  delete force_fopr_c_prec;
391 
392  delete proj;
393  delete smear;
394  delete force_smear;
395  delete dr_smear;
396  delete fopr_smear;
397  delete force_fopr_smear;
398  delete fopr_smear_prec;
399  delete force_fopr_smear_prec;
400 
401  delete solver_MD, fprop_MD;
402  delete solver_H, fprop_H;
403  delete solver_H_prec, fprop_H_prec;
404 
405  delete action_F_Nf2_prec;
406  delete action_F_Nf2_d;
407 
408  delete fopr_c1;
409  delete force_fopr_c1;
410  delete fopr_r1;
411  delete fopr_langev;
412  delete fopr_r2;
413  delete fopr_H;
414  delete force_fopr_r2;
415  delete force_fopr_MD;
416  delete action_F_Nf1;
417 
418  delete builder;
419 
420  delete rand;
421 
422  delete ploop;
423  delete rand_nv;
424  delete nv;
425  delete quark_suscept;
426 
427  delete params_test;
428  delete params_action_G;
429  delete params_clover;
430  delete params_clover_prec;
431  delete params_clover1;
432  delete params_proj;
433  delete params_smear;
434  delete params_dr_smear;
435  delete params_rational_MD;
436  delete params_rational_H;
437  delete params_solver_MD;
438  delete params_solver_H;
439  delete params_integrator;
440  delete params_hmc;
441  delete params_all;
442 
443 
444 #ifdef USE_TEST
445  return Test::verify(expected_result, result);
446 
447 #else
448  return EXIT_SUCCESS;
449 #endif
450  }
451 } // namespace Test_HMC_Clover_Isochemical