Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
test_HMC_Clover_Isochemical_Leapfrog_Nf2.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 
37 #include "forceSmear.h"
38 #include "projection.h"
39 #include "smear.h"
40 #include "solver.h"
41 
42 #include "hmc_Leapfrog.h"
43 
44 #include "polyakovLoop.h"
45 
46 #include "noiseVector_Z2.h"
47 #include "quarkNumberSusceptibility_Wilson.h"
48 
49 #ifdef USE_TEST
50 #include "test.h"
51 #endif
52 
53 #ifdef USE_TESTMANAGER_AUTOREGISTER
54 #include "testManager.h"
55 #endif
56 
57 //====================================================================
59 
70 namespace Test_HMC_Clover_Isochemical {
71  //- test-private parameters
72  namespace {
73  const std::string filename_input = "test_HMC_Clover_Isochemical_Leapfrog_Nf2.yaml";
74  const std::string filename_output = "stdout";
75 
76  class Parameters_Test_HMC_Clover_Isochemical : public Parameters {
77  public:
78  Parameters_Test_HMC_Clover_Isochemical()
79  {
80  Register_string("gauge_config_status", "NULL");
81 
82  Register_string("gauge_config_type_input", "NULL");
83  Register_string("config_filename_input", "NULL");
84 
85  Register_string("gauge_config_type_output", "NULL");
86  Register_string("config_filename_output", "NULL");
87 
88  Register_int("trajectory_number", 0);
89  Register_int("trajectory_number_step", 0);
90  Register_int("save_config_interval", 0);
91 
92  Register_int("seed_for_noise", 0);
93  Register_int("number_of_noises", 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 leapfrog_Nf2(void);
104 
105 #ifdef USE_TESTMANAGER_AUTOREGISTER
106  namespace {
107  static const bool is_registered = TestManager::RegisterTest(
108  "HMC.Clover_Isochemical.Leapfrog_Nf2",
109  leapfrog_Nf2
110  );
111  }
112 #endif
113 
114  //====================================================================
115  int leapfrog_Nf2(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_Isochemical;
124  Parameters *params_action_G = ParametersFactory::New("Action.G_Rectangle");
125  Parameters *params_clover = ParametersFactory::New("Fopr.Clover_Isochemical");
126  Parameters *params_proj = ParametersFactory::New("Projection");
127  Parameters *params_smear = ParametersFactory::New("Smear");
128  Parameters *params_dr_smear = ParametersFactory::New("Director_Smear");
129  Parameters *params_solver_MD = ParametersFactory::New("Solver");
130  Parameters *params_solver_H = ParametersFactory::New("Solver");
131  Parameters *params_hmc = ParametersFactory::New("HMC.Leapfrog");
132 
133  Parameters *params_all = new Parameters;
134 
135  params_all->Register_Parameters("Test_HMC_Clover_Isochemical", params_test);
136  params_all->Register_Parameters("Action_G_Rectangle", params_action_G);
137  params_all->Register_Parameters("Fopr_Clover_Isochemical", params_clover);
138  params_all->Register_Parameters("Projection", params_proj);
139  params_all->Register_Parameters("Smear", params_smear);
140  params_all->Register_Parameters("Director_Smear", params_dr_smear);
141  params_all->Register_Parameters("Solver_MD", params_solver_MD);
142  params_all->Register_Parameters("Solver_H", params_solver_H);
143  params_all->Register_Parameters("HMC_Leapfrog", params_hmc);
144 
145  ParameterManager_YAML params_manager;
146  params_manager.read_params(filename_input, params_all);
147 
148  const string str_gconf_status = params_test->get_string("gauge_config_status");
149  const string str_gconf_read = params_test->get_string("gauge_config_type_input");
150  const string readfile = params_test->get_string("config_filename_input");
151  const string str_gconf_write = params_test->get_string("gauge_config_type_output");
152  const string writefile = params_test->get_string("config_filename_output");
153  int iconf = params_test->get_int("trajectory_number");
154  const int Ntraj = params_test->get_int("trajectory_number_step");
155  const int i_save_conf = params_test->get_int("save_config_interval");
156  int i_seed_noise = params_test->get_int("seed_for_noise");
157  const string str_vlevel = params_test->get_string("verbose_level");
158 #ifdef USE_TEST
159  const double expected_result = params_test->get_double("expected_result");
160 #endif
161 
162  const string str_gmset_type = params_clover->get_string("gamma_matrix_type");
163  const string str_proj_type = params_proj->get_string("projection_type");
164  const string str_smear_type = params_smear->get_string("smear_type");
165  const string str_solver_MD_type = params_solver_MD->get_string("solver_type");
166  const string str_solver_H_type = params_solver_H->get_string("solver_type");
167 
168  Bridge::VerboseLevel vl = vout.set_verbose_level(str_vlevel);
169 
170  //- print input parameters
171  vout.general(vl, " gconf_status = %s\n", str_gconf_status.c_str());
172  vout.general(vl, " gconf_read = %s\n", str_gconf_read.c_str());
173  vout.general(vl, " readfile = %s\n", readfile.c_str());
174  vout.general(vl, " gconf_write = %s\n", str_gconf_write.c_str());
175  vout.general(vl, " writefile = %s\n", writefile.c_str());
176  vout.general(vl, " iconf = %d\n", iconf);
177  vout.general(vl, " Ntraj = %d\n", Ntraj);
178  vout.general(vl, " i_save_conf = %d\n", i_save_conf);
179  vout.general(vl, " i_seed_noise = %d\n", i_seed_noise);
180  vout.general(vl, " vlevel = %s\n", str_vlevel.c_str());
181 
182  vout.general(vl, " gmset_type = %s\n", str_gmset_type.c_str());
183  vout.general(vl, " proj_type = %s\n", str_proj_type.c_str());
184  vout.general(vl, " smear_type = %s\n", str_smear_type.c_str());
185  vout.general(vl, " solver_MD_type = %s\n", str_solver_MD_type.c_str());
186  vout.general(vl, " solver_H_type = %s\n", str_solver_H_type.c_str());
187  vout.general(vl, "\n");
188 
189  //- input parameter check
190  int err = 0;
191  err += ParameterCheck::non_NULL(str_gconf_read);
192  err += ParameterCheck::non_NULL(readfile);
193  err += ParameterCheck::non_NULL(str_gconf_write);
194  err += ParameterCheck::non_NULL(writefile);
195  err += ParameterCheck::non_zero(iconf);
196  err += ParameterCheck::non_zero(Ntraj);
197  err += ParameterCheck::non_zero(i_save_conf);
198 
199  if (err) {
200  vout.crucial(vl, "test_HMC_Clover_Isochemical: Input parameters have not been set.\n");
201  abort();
202  }
203 
204 
205  // ##### object setup #####
206  Field_G *U = new Field_G(Nvol, Ndim);
207  GaugeConfig *gconf_read = new GaugeConfig(str_gconf_read);
208 
209  if (str_gconf_status == "Continue") {
210  gconf_read->read_file((Field *)U, readfile);
211  } else if (str_gconf_status == "Start_cold") {
212  gconf_read->set_cold((Field *)U);
213  } else {
214  vout.crucial(vl, "Test_HMC_Clover_Isochemical: unsupported gconf status \"%s\".\n", str_gconf_status.c_str());
215  abort();
216  }
217 
218 
219  Action_G_Rectangle *action_G = new Action_G_Rectangle; // plaq.+rectangle (SA)
220  action_G->set_parameters(*params_action_G);
221 
222 
223  Fopr_Clover_Isochemical *fopr_c = new Fopr_Clover_Isochemical(str_gmset_type); // define fermion operator (SA)
224  fopr_c->set_parameters(*params_clover);
225 
226  Force_F_Clover_Nf2_Isochemical *force_fopr_c = new Force_F_Clover_Nf2_Isochemical(str_gmset_type);
227  // define fermion force (SA)
228  force_fopr_c->set_parameters(*params_clover);
229 
230  // define smearing method (SA)
231  Projection *proj = Projection::New(str_proj_type);
232  Smear *smear = Smear::New(str_smear_type, proj);
233  smear->set_parameters(*params_smear);
234 
235  // define force smearing method (SA)
236  ForceSmear *force_smear = ForceSmear::New(str_smear_type, proj);
237  force_smear->set_parameters(*params_smear);
238 
239  Director_Smear *dr_smear = new Director_Smear((Smear *)smear);
240  dr_smear->set_parameters(*params_dr_smear);
241 
242  Fopr_Smeared *fopr_smear = new Fopr_Smeared((Fopr *)fopr_c, dr_smear);
243  // define smeared fermion operator (SA)
244  Force_F_Smeared *force_fopr_smear
245  = new Force_F_Smeared((Force *)force_fopr_c, (ForceSmear *)force_smear, dr_smear);
246  // define smeared fermion force (SA)
247 
248 
249  Solver *solver_MD = Solver::New(str_solver_MD_type, fopr_smear);
250  solver_MD->set_parameters(*params_solver_MD);
251  Fprop *fprop_MD = new Fprop_Standard_lex(solver_MD);
252 
253  Solver *solver_H = Solver::New(str_solver_H_type, fopr_smear);
254  solver_H->set_parameters(*params_solver_H);
255  Fprop *fprop_H = new Fprop_Standard_lex(solver_H);
256 
257  Action_F_Standard_lex *action_F
258  = new Action_F_Standard_lex((Fopr *)fopr_smear, (Force *)force_fopr_smear,
259  fprop_MD, fprop_H);
260  // define fermion action (SA)
261 
262 
263  valarray<Action *> actions(2);
264  actions[0] = (Action *)action_G; // register action[0] (SA)
265  actions[1] = (Action *)action_F; // register action[1] (SA)
266 
267  valarray<Director *> directors(1);
268  directors[0] = (Director *)dr_smear; // register director[0] (SA)
269 
270 
271  //- Random number is initialized with a parameter specified by iconf
272  RandomNumbers *rand = new RandomNumbers_Mseries(iconf);
273 
274 
275  HMC_Leapfrog hmc(actions, directors, rand); // define hmc_leapfrog (SA)
276  hmc.set_parameters(*params_hmc);
277 
278 
279  PolyakovLoop *ploop = new PolyakovLoop();
280 
281  RandomNumbers *rand_nv = new RandomNumbers_Mseries(i_seed_noise);
282  NoiseVector *nv = new NoiseVector_Z2(rand_nv);
283  QuarkNumberSusceptibility_Wilson *quark_suscept
284  = new QuarkNumberSusceptibility_Wilson(fopr_smear, fprop_H, nv);
285  quark_suscept->set_parameters(*params_test);
286 
287 
288  // #### Execution main part ####
289  vout.general(vl, "HMC: Ntraj = %d\n", Ntraj);
290 
291  double ploop0 = ploop->measure_ploop(*U);
292  vout.general(vl, "Polyakov loop = %f\n", ploop0);
293 
294  double result = 0.0;
295  for (int traj = 0; traj < Ntraj; ++traj) {
296  vout.general(vl, "\n");
297  vout.general(vl, "traj = %d\n", traj);
298 
299  result = hmc.update(*U); // hmc update (SA)
300 
301  double ploop1 = ploop->measure_ploop(*U);
302  vout.general(vl, "Polyakov loop = %f\n", ploop1);
303 
304  double result_quark_suscept1 = quark_suscept->measure();
305  }
306 
307 
308  // ##### tidy up #####
309  delete U;
310  delete gconf_read;
311 
312  delete action_G;
313 
314  delete fopr_c;
315  delete force_fopr_c;
316 
317  delete proj;
318  delete smear;
319  delete force_smear;
320  delete dr_smear;
321  delete fopr_smear;
322  delete force_fopr_smear;
323 
324  delete solver_MD, fprop_MD;
325  delete solver_H, fprop_H;
326 
327  delete action_F;
328 
329  delete rand;
330 
331  delete ploop;
332  delete rand_nv;
333  delete nv;
334  delete quark_suscept;
335 
336  delete params_test;
337  delete params_action_G;
338  delete params_clover;
339  delete params_proj;
340  delete params_smear;
341  delete params_dr_smear;
342  delete params_solver_MD;
343  delete params_solver_H;
344  delete params_hmc;
345  delete params_all;
346 
347 
348 #ifdef USE_TEST
349  return Test::verify(expected_result, result);
350 
351 #else
352  return EXIT_SUCCESS;
353 #endif
354  }
355 } // namespace Test_HMC_Clover_Isochemical