Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
test_HMC_Wilson_Leapfrog_Nf2_TopologyFixing.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 "force_F_Wilson_Nf2.h"
26 #include "fopr_TMWilson.h"
27 #include "force_F_TMWilson_Nf2.h"
28 
29 #include "randomNumbers_Mseries.h"
30 
31 #include "director_Smear.h"
32 #include "fopr_Smeared.h"
33 #include "force_F_Smeared.h"
34 
35 #include "action_F_Standard_lex.h"
36 #include "fprop_Standard_lex.h"
37 #include "action_F_Ratio_lex.h"
38 
39 #include "forceSmear.h"
40 #include "projection.h"
41 #include "smear.h"
42 #include "solver.h"
43 
44 #include "hmc_Leapfrog.h"
45 
46 #ifdef USE_TEST
47 #include "test.h"
48 #endif
49 
50 #ifdef USE_TESTMANAGER_AUTOREGISTER
51 #include "testManager.h"
52 #endif
53 
54 //====================================================================
56 
68 namespace Test_HMC_Wilson {
69  //- test-private parameters
70  namespace {
71  const std::string filename_input = "test_HMC_Wilson_Leapfrog_Nf2_TopologyFixing.yaml";
72  const std::string filename_output = "stdout";
73 
74  class Parameters_Test_HMC_Wilson : public Parameters {
75  public:
76  Parameters_Test_HMC_Wilson()
77  {
78  Register_string("gauge_config_status", "NULL");
79 
80  Register_string("gauge_config_type_input", "NULL");
81  Register_string("config_filename_input", "NULL");
82 
83  Register_string("gauge_config_type_output", "NULL");
84  Register_string("config_filename_output", "NULL");
85 
86  Register_int("trajectory_number", 0);
87  Register_int("trajectory_number_step", 0);
88  Register_int("save_config_interval", 0);
89 
90  Register_string("verbose_level", "NULL");
91 
92  Register_double("expected_result", 0.0);
93  }
94  };
95  }
96 
97  //- prototype declaration
99 
100 #ifdef USE_TESTMANAGER_AUTOREGISTER
101  namespace {
102  static const bool is_registered = TestManager::RegisterTest(
103  "HMC.Wilson.Leapfrog_Nf2_TopologyFixing",
105  );
106  }
107 #endif
108 
109  //====================================================================
111  {
112  // ##### parameter setup #####
113  int Nc = CommonParameters::Nc();
114  int Nd = CommonParameters::Nd();
115  int Nvol = CommonParameters::Nvol();
116  int Ndim = CommonParameters::Ndim();
117  int NinG = 2 * Nc * Nc;
118 
119  Parameters *params_test = new Parameters_Test_HMC_Wilson;
120  Parameters *params_action_G = ParametersFactory::New("Action.G_Rectangle");
121  Parameters *params_wilson = ParametersFactory::New("Fopr.Wilson");
122  Parameters *params_tmwilson = ParametersFactory::New("Fopr.TMWilson");
123  Parameters *params_proj = ParametersFactory::New("Projection");
124  Parameters *params_smear = ParametersFactory::New("Smear");
125  Parameters *params_dr_smear = ParametersFactory::New("Director_Smear");
126  Parameters *params_solver_MD = ParametersFactory::New("Solver");
127  Parameters *params_solver_H = ParametersFactory::New("Solver");
128  Parameters *params_hmc = ParametersFactory::New("HMC.Leapfrog");
129 
130  Parameters *params_all = new Parameters;
131 
132  params_all->Register_Parameters("Test_HMC_Wilson", params_test);
133  params_all->Register_Parameters("Action_G_Rectangle", params_action_G);
134  params_all->Register_Parameters("Fopr_Wilson", params_wilson);
135  params_all->Register_Parameters("Fopr_TMWilson", params_tmwilson);
136  params_all->Register_Parameters("Projection", params_proj);
137  params_all->Register_Parameters("Smear", params_smear);
138  params_all->Register_Parameters("Director_Smear", params_dr_smear);
139  params_all->Register_Parameters("Solver_MD", params_solver_MD);
140  params_all->Register_Parameters("Solver_H", params_solver_H);
141  params_all->Register_Parameters("HMC_Leapfrog", params_hmc);
142 
143  ParameterManager_YAML params_manager;
144  params_manager.read_params(filename_input, params_all);
145 
146  const string str_gconf_status = params_test->get_string("gauge_config_status");
147  const string str_gconf_read = params_test->get_string("gauge_config_type_input");
148  const string readfile = params_test->get_string("config_filename_input");
149  const string str_gconf_write = params_test->get_string("gauge_config_type_output");
150  const string writefile = params_test->get_string("config_filename_output");
151  int iconf = params_test->get_int("trajectory_number");
152  int Ntraj = params_test->get_int("trajectory_number_step");
153  const int i_save_conf = params_test->get_int("save_config_interval");
154  const string str_vlevel = params_test->get_string("verbose_level");
155 #ifdef USE_TEST
156  const double expected_result = params_test->get_double("expected_result");
157 #endif
158 
159  //- NB. gamma_matrix_type of Wilson and TMWilson must be the same.
160  const string str_gmset_type = params_wilson->get_string("gamma_matrix_type");
161  const string str_proj_type = params_proj->get_string("projection_type");
162  const string str_smear_type = params_smear->get_string("smear_type");
163  const string str_solver_MD_type = params_solver_MD->get_string("solver_type");
164  const string str_solver_H_type = params_solver_H->get_string("solver_type");
165 
167 
168  //- print input parameters
169  vout.general(vl, " gconf_status = %s\n", str_gconf_status.c_str());
170  vout.general(vl, " gconf_read = %s\n", str_gconf_read.c_str());
171  vout.general(vl, " readfile = %s\n", readfile.c_str());
172  vout.general(vl, " gconf_write = %s\n", str_gconf_write.c_str());
173  vout.general(vl, " writefile = %s\n", writefile.c_str());
174  vout.general(vl, " iconf = %d\n", iconf);
175  vout.general(vl, " Ntraj = %d\n", Ntraj);
176  vout.general(vl, " i_save_conf = %d\n", i_save_conf);
177  vout.general(vl, " vlevel = %s\n", str_vlevel.c_str());
178 
179  vout.general(vl, " gmset_type = %s\n", str_gmset_type.c_str());
180  vout.general(vl, " proj_type = %s\n", str_proj_type.c_str());
181  vout.general(vl, " smear_type = %s\n", str_smear_type.c_str());
182  vout.general(vl, " solver_MD_type = %s\n", str_solver_MD_type.c_str());
183  vout.general(vl, " solver_H_type = %s\n", str_solver_H_type.c_str());
184  vout.general(vl, "\n");
185 
186  //- input parameter check
187  int err = 0;
188  err += ParameterCheck::non_NULL(str_gconf_read);
189  err += ParameterCheck::non_NULL(readfile);
190  err += ParameterCheck::non_NULL(str_gconf_write);
191  err += ParameterCheck::non_NULL(writefile);
192  err += ParameterCheck::non_zero(iconf);
193  err += ParameterCheck::non_zero(Ntraj);
194  err += ParameterCheck::non_zero(i_save_conf);
195 
196  if (err) {
197  vout.crucial(vl, "Test_HMC_Wilson: Input parameters have not been set.\n");
198  abort();
199  }
200 
201 
202  // ##### object setup #####
203  Field_G *U = new Field_G(Nvol, Ndim);
204  GaugeConfig *gconf_read = new GaugeConfig(str_gconf_read);
205 
206  if (str_gconf_status == "Continue") {
207  gconf_read->read_file((Field *)U, readfile);
208  } else if (str_gconf_status == "Start_cold") {
209  gconf_read->set_cold((Field *)U);
210  } else {
211  vout.crucial(vl, "Test_HMC_Wilson: unsupported gconf status \"%s\".\n", str_gconf_status.c_str());
212  abort();
213  }
214 
215 
216  Action_G_Rectangle *action_G = new Action_G_Rectangle;
217  action_G->set_parameters(*params_action_G);
218 
219 
220  Projection *proj = Projection::New(str_proj_type);
221  Smear *smear = Smear::New(str_smear_type, proj);
222  smear->set_parameters(*params_smear);
223 
224  ForceSmear *force_smear = ForceSmear::New(str_smear_type, proj);
225  force_smear->set_parameters(*params_smear);
226 
227  Director_Smear *dr_smear = new Director_Smear((Smear *)smear);
228  dr_smear->set_parameters(*params_dr_smear);
229 
230 
231  Fopr_Wilson *fopr_w = new Fopr_Wilson(str_gmset_type);
232  fopr_w->set_parameters(*params_wilson);
233 
234  Force_F_Wilson_Nf2 *force_fopr_w = new Force_F_Wilson_Nf2(str_gmset_type);
235  force_fopr_w->set_parameters(*params_wilson);
236 
237  Fopr_Smeared *fopr_smear = new Fopr_Smeared((Fopr *)fopr_w, dr_smear);
238  Force_F_Smeared *force_fopr_smear
239  = new Force_F_Smeared((Force *)force_fopr_w, (ForceSmear *)force_smear, dr_smear);
240 
241 
242  Solver *solver_MD = Solver::New(str_solver_MD_type, fopr_smear);
243  solver_MD->set_parameters(*params_solver_MD);
244  Fprop *fprop_MD = new Fprop_Standard_lex(solver_MD);
245 
246  Solver *solver_H = Solver::New(str_solver_H_type, fopr_smear);
247  solver_H->set_parameters(*params_solver_H);
248  Fprop *fprop_H = new Fprop_Standard_lex(solver_H);
249 
250 
251  //- preconditioner
252  Fopr_TMWilson *fopr_prec = new Fopr_TMWilson(str_gmset_type);
253  fopr_prec->set_parameters(*params_tmwilson);
254 
255  Force_F_TMWilson_Nf2 *force_fopr_prec = new Force_F_TMWilson_Nf2(str_gmset_type);
256  force_fopr_prec->set_parameters(*params_tmwilson);
257 
258  Fopr_Smeared *fopr_smear_prec = new Fopr_Smeared((Fopr *)fopr_prec, dr_smear);
259  Force_F_Smeared *force_fopr_smear_prec
260  = new Force_F_Smeared((Force *)force_fopr_prec, (ForceSmear *)force_smear, dr_smear);
261 
262 
263  Solver *solver_H_prec = Solver::New(str_solver_H_type, fopr_smear_prec);
264  solver_H_prec->set_parameters(*params_solver_H);
265  Fprop *fprop_H_prec = new Fprop_Standard_lex(solver_H_prec);
266 
267 
268  Action_F_Ratio_lex *action_F
269  = new Action_F_Ratio_lex((Fopr *)fopr_smear_prec, (Force *)force_fopr_smear_prec,
270  (Fopr *)fopr_smear, (Force *)force_fopr_smear,
271  fprop_H_prec,
272  fprop_MD, fprop_H);
273 
274 
275  valarray<Action *> actions(2);
276  actions[0] = (Action *)action_G;
277  actions[1] = (Action *)action_F;
278 
279  valarray<Director *> directors(1);
280  directors[0] = (Director *)dr_smear;
281 
282  RandomNumbers *rand = new RandomNumbers_Mseries(iconf);
283 
284  HMC_Leapfrog hmc(actions, directors, rand);
285  hmc.set_parameters(*params_hmc);
286 
287 
288  // #### Execution main part ####
289  vout.general(vl, "HMC: Ntraj = %d\n", Ntraj);
290 
291  double result = 0.0;
292  for (int traj = 0; traj < Ntraj; ++traj) {
293  vout.general(vl, "\n");
294  vout.general(vl, "traj = %d\n", traj);
295 
296  result = hmc.update(*U);
297  }
298 
299 
300  // ##### tidy up #####
301  delete params_test;
302  delete params_action_G;
303  delete params_wilson;
304  delete params_tmwilson;
305  delete params_smear;
306  delete params_dr_smear;
307  delete params_solver_MD;
308  delete params_solver_H;
309  delete params_hmc;
310  delete params_all;
311 
312  delete U;
313  delete gconf_read;
314 
315  delete rand;
316 
317  delete action_G;
318  delete action_F;
319 
320  delete proj;
321  delete smear;
322  delete force_smear;
323  delete dr_smear;
324 
325  delete fopr_w;
326  delete force_fopr_w;
327  delete fopr_smear;
328  delete force_fopr_smear;
329  delete solver_MD;
330  delete fprop_MD;
331  delete solver_H;
332  delete fprop_H;
333 
334  delete fopr_prec;
335  delete force_fopr_prec;
336  delete fopr_smear_prec;
337  delete force_fopr_smear_prec;
338  delete solver_H_prec;
339  delete fprop_H_prec;
340 
341 
342 #ifdef USE_TEST
343  return Test::verify(expected_result, result);
344 
345 #else
346  return EXIT_SUCCESS;
347 #endif
348  }
349 } // namespace Test_HMC_Wilson