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