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