Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
test_HMC_Clover_SF_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 "staples_SF.h"
23 #include "action_G_Plaq_SF.h"
24 #include "action_G_Rectangle_SF.h"
25 
26 #include "fopr_Clover_SF.h"
27 #include "force_F_Clover_SF.h"
28 
29 #include "randomNumbers_Mseries.h"
30 #include "randomNumbers_MT19937.h"
31 
32 #include "director_Smear.h"
33 #include "fopr_Smeared.h"
34 #include "force_F_Smeared.h"
35 
36 #include "action_F_Standard_SF.h"
37 
38 #include "fopr_Rational_SF.h"
39 #include "action_F_Rational_Frame_SF.h"
40 
41 #include "forceSmear.h"
42 #include "projection.h"
43 #include "smear.h"
44 #include "solver.h"
45 
46 #include "hmc_General.h"
47 #include "builder_Integrator.h"
48 #include "integrator.h"
49 
50 #ifdef USE_TEST
51 #include "test.h"
52 #endif
53 
54 #ifdef USE_TESTMANAGER_AUTOREGISTER
55 #include "testManager.h"
56 #endif
57 
58 //====================================================================
60 
86 namespace Test_HMC_Clover_SF {
87  //- test-private parameters
88  namespace {
89  const std::string filename_input = "test_HMC_Clover_SF_RHMC_Nf2p1.yaml";
90  const std::string filename_output = "stdout";
91 
92  class Parameters_Test_HMC_Clover_SF : public Parameters {
93  public:
94  Parameters_Test_HMC_Clover_SF()
95  {
96  Register_string("gauge_config_status", "NULL");
97 
98  Register_string("gauge_config_type_input", "NULL");
99  Register_string("config_filename_input", "NULL");
100 
101  Register_string("gauge_config_type_output", "NULL");
102  Register_string("config_filename_output", "NULL");
103 
104  Register_string("verbose_level", "NULL");
105 
106  Register_int("trajectory_number", 0);
107  Register_int("trajectory_number_step", 0);
108  Register_int("save_config_interval", 0);
109 
110  Register_double("expected_result", 0.0);
111  }
112  };
113  }
114 
115  //- prototype declaration
116  int RHMC_Nf2p1(void);
117 
118 #ifdef USE_TESTMANAGER_AUTOREGISTER
119  namespace {
120  static const bool is_registered = TestManager::RegisterTest(
121  "HMC.Clover_SF.RHMC_Nf2p1",
122  RHMC_Nf2p1
123  );
124  }
125 #endif
126 
127  //====================================================================
128  int RHMC_Nf2p1(void)
129  {
130  // ##### parameter setup #####
131  int Nc = CommonParameters::Nc();
132  int Nvol = CommonParameters::Nvol();
133  int Ndim = CommonParameters::Ndim();
134  int NinG = 2 * Nc * Nc;
135 
136  Parameters *params_test = new Parameters_Test_HMC_Clover_SF;
137  Parameters *params_action_G = ParametersFactory::New("Action.G_Rectangle_SF");
138  Parameters *params_clover = ParametersFactory::New("Fopr.Clover_SF");
139  Parameters *params_clover1 = ParametersFactory::New("Fopr.Clover_SF");
140  Parameters *params_proj = ParametersFactory::New("Projection");
141  Parameters *params_smear = ParametersFactory::New("Smear");
142  Parameters *params_dr_smear = ParametersFactory::New("Director_Smear");
143  Parameters *params_rational_H = ParametersFactory::New("Fopr.Rational");
144  Parameters *params_rational_MD = ParametersFactory::New("Fopr.Rational");
145  Parameters *params_integrator = ParametersFactory::New("Builder_Integrator");
146  Parameters *params_hmc = ParametersFactory::New("HMC.General");
147 
148  Parameters *params_all = new Parameters;
149 
150  params_all->Register_Parameters("Test_HMC_Clover_SF", params_test);
151  params_all->Register_Parameters("Action_G_Rectangle_SF", params_action_G);
152  params_all->Register_Parameters("Fopr_Clover_Nf2_SF", params_clover);
153  params_all->Register_Parameters("Fopr_Clover_Nf1_SF", params_clover1);
154  params_all->Register_Parameters("Projection", params_proj);
155  params_all->Register_Parameters("Smear", params_smear);
156  params_all->Register_Parameters("Director_Smear", params_dr_smear);
157  params_all->Register_Parameters("Fopr_Rational_H", params_rational_H);
158  params_all->Register_Parameters("Fopr_Rational_MD", params_rational_MD);
159  params_all->Register_Parameters("Builder_Integrator", params_integrator);
160  params_all->Register_Parameters("HMC_General", params_hmc);
161 
162  ParameterManager_YAML params_manager;
163  params_manager.read_params(filename_input, params_all);
164 
165  const string str_gconf_status = params_test->get_string("gauge_config_status");
166  const string str_gconf_read = params_test->get_string("gauge_config_type_input");
167  const string readfile = params_test->get_string("config_filename_input");
168  const string str_gconf_write = params_test->get_string("gauge_config_type_output");
169  const string writefile = params_test->get_string("config_filename_output");
170  int iconf = params_test->get_int("trajectory_number");
171  int Ntraj = params_test->get_int("trajectory_number_step");
172  const int i_save_conf = params_test->get_int("save_config_interval");
173  const string str_vlevel = params_test->get_string("verbose_level");
174 #ifdef USE_TEST
175  const double expected_result = params_test->get_double("expected_result");
176 #endif
177 
178  const string str_gmset_type = params_clover->get_string("gamma_matrix_type");
179  const string str_proj_type = params_proj->get_string("projection_type");
180  const string str_smear_type = params_smear->get_string("smear_type");
181 
182  Bridge::VerboseLevel vl = vout.set_verbose_level(str_vlevel);
183 
184  //- print input parameters
185  vout.general(vl, " gconf_status = %s\n", str_gconf_status.c_str());
186  vout.general(vl, " gconf_read = %s\n", str_gconf_read.c_str());
187  vout.general(vl, " readfile = %s\n", readfile.c_str());
188  vout.general(vl, " gconf_write = %s\n", str_gconf_write.c_str());
189  vout.general(vl, " writefile = %s\n", writefile.c_str());
190  vout.general(vl, " iconf = %d\n", iconf);
191  vout.general(vl, " Ntraj = %d\n", Ntraj);
192  vout.general(vl, " i_save_conf = %d\n", i_save_conf);
193  vout.general(vl, " vlevel = %s\n", str_vlevel.c_str());
194 
195  vout.general(vl, " gmset_type = %s\n", str_gmset_type.c_str());
196  vout.general(vl, " proj_type = %s\n", str_proj_type.c_str());
197  vout.general(vl, " smear_type = %s\n", str_smear_type.c_str());
198  vout.general(vl, "\n");
199 
200  //- input parameter check
201  int err = 0;
202  err += ParameterCheck::non_NULL(str_gconf_read);
203  err += ParameterCheck::non_NULL(readfile);
204  err += ParameterCheck::non_NULL(str_gconf_write);
205  err += ParameterCheck::non_NULL(writefile);
206  err += ParameterCheck::non_zero(iconf);
207  err += ParameterCheck::non_zero(Ntraj);
208  err += ParameterCheck::non_zero(i_save_conf);
209 
210  if (err) {
211  vout.crucial(vl, "Test_HMC_Clover_SF: Input parameters have not been set.\n");
212  abort();
213  }
214 
215 
216  // ##### object setup #####
217  Field_G *U = new Field_G(Nvol, Ndim);
218  GaugeConfig *gconf_read = new GaugeConfig(str_gconf_read);
219 
220  if (str_gconf_status == "Continue") {
221  gconf_read->read_file((Field *)U, readfile);
222  } else if (str_gconf_status == "Start_cold") {
223  gconf_read->set_cold((Field *)U);
224  } else {
225  vout.crucial(vl, "Test_HMC_Clover_SF: unsupported gconf status \"%s\".\n", str_gconf_status.c_str());
226  abort();
227  }
228 
229 
230  Action_G_Rectangle_SF *action_G = new Action_G_Rectangle_SF;
231  action_G->set_parameters(*params_action_G);
232 
233 
234  //-- N_f=2 part
235  Fopr_Clover_SF *fopr_w = new Fopr_Clover_SF();
236  fopr_w->set_parameters(*params_clover);
237  Force_F_Clover_SF *force_fopr_w = new Force_F_Clover_SF();
238  force_fopr_w->set_parameters(*params_clover);
239 
240 
241  Projection *proj = Projection::New(str_proj_type);
242  Smear *smear = Smear::New(str_smear_type, proj);
243  smear->set_parameters(*params_smear);
244 
245  ForceSmear *force_smear = ForceSmear::New(str_smear_type, proj);
246  force_smear->set_parameters(*params_smear);
247 
248  Director_Smear *dr_smear = new Director_Smear((Smear *)smear);
249  dr_smear->set_parameters(*params_dr_smear);
250 
251  Fopr_Smeared *fopr_smear = new Fopr_Smeared((Fopr *)fopr_w, dr_smear);
252  Force_F_Smeared *force_fopr_smear
253  = new Force_F_Smeared((Force *)force_fopr_w, (ForceSmear *)force_smear, dr_smear);
254 
255 
256  Action_F_Standard_SF *action_F_Nf2
257  = new Action_F_Standard_SF((Fopr *)fopr_smear, (Force *)force_fopr_smear);
258 
259 
260  //-- N_f=1 part
261  Fopr_Clover_SF *fopr_w1 = new Fopr_Clover_SF();
262  fopr_w1->set_parameters(*params_clover1);
263  Force_F_Clover_SF *force_fopr_w1 = new Force_F_Clover_SF();
264  force_fopr_w1->set_parameters(*params_clover1);
265 
266  Fopr_Rational_SF *fopr_r1 = new Fopr_Rational_SF(fopr_w1);
267  fopr_r1->set_parameters(*params_rational_H);
268  Fopr_Smeared *fopr_langev = new Fopr_Smeared((Fopr *)fopr_r1, dr_smear);
269 
270  Fopr_Rational_SF *fopr_r2 = new Fopr_Rational_SF(fopr_w1);
271  fopr_r2->set_parameters(*params_rational_MD);
272  Fopr_Smeared *fopr_H = new Fopr_Smeared((Fopr *)fopr_r2, dr_smear);
273 
274  Force_F_Rational *force_fopr_r2 = new Force_F_Rational(fopr_w1, force_fopr_w1);
275  force_fopr_r2->set_parameters(*params_rational_MD);
276  Force_F_Smeared *force_fopr_MD
277  = new Force_F_Smeared((Force *)force_fopr_r2,
278  (ForceSmear *)force_smear, dr_smear);
279 
280  Action_F_Rational_frame_SF *action_F_Nf1
281  = new Action_F_Rational_frame_SF((Fopr *)fopr_langev,
282  (Fopr *)fopr_H,
283  (Force *)force_fopr_MD);
284 
285 
286  valarray<Action *> actions(3);
287  actions[0] = (Action *)action_F_Nf1;
288  actions[1] = (Action *)action_F_Nf2;
289  actions[2] = (Action *)action_G;
290 
291  valarray<Director *> directors(1);
292  directors[0] = (Director *)dr_smear;
293 
294  Builder_Integrator *builder = new Builder_Integrator(actions, directors);
295  builder->set_parameters(*params_integrator);
296  Integrator *integrator = builder->build();
297 
298 
299  //- Mersenne Twister for random number generator.
300  RandomNumbers_MT19937 *rand;
301  if (iconf == 0) {
302  std::vector<unsigned long> seed(4);
303  seed[0] = 0x6a92;
304  seed[1] = 0x3708;
305  seed[2] = 0xab41;
306  seed[3] = 0x5c52;
307  rand = new RandomNumbers_MT19937(seed);
308  } else {
309  rand = new RandomNumbers_MT19937(iconf);
310  }
311 
312 
313  HMC_General hmc(actions, directors, integrator, (RandomNumbers *)rand);
314  hmc.set_parameters(*params_hmc);
315 
316 
317  // #### Execution main part ####
318  vout.general(vl, "RHMC start: Ntraj = %d\n", Ntraj);
319 
320  double result = 0.0;
321  for (int traj = 0; traj < Ntraj; ++traj) {
322  vout.general(vl, "\n");
323  vout.general(vl, "---------------------------------------------------\n");
324  vout.general(vl, "traj = %d\n", traj);
325 
326  result = hmc.update(*U);
327  }
328 
329  if (Communicator::nodeid() == 0) {
330  rand->writefile("rand_MT");
331  }
332 
333 
334  // ##### tidy up #####
335  delete params_test;
336  delete params_action_G;
337  delete params_clover;
338  delete params_clover1;
339  delete params_smear;
340  delete params_dr_smear;
341  delete params_rational_H;
342  delete params_rational_MD;
343  delete params_integrator;
344  delete params_hmc;
345  delete params_all;
346 
347  delete U;
348  delete gconf_read;
349 
350  delete rand;
351 
352  delete action_G;
353 
354  delete fopr_w;
355  delete force_fopr_w;
356  delete action_F_Nf2;
357 
358  delete proj;
359  delete smear;
360  delete force_smear;
361  delete dr_smear;
362  delete fopr_smear;
363  delete force_fopr_smear;
364 
365  delete fopr_w1;
366  delete force_fopr_w1;
367  delete fopr_r1;
368  delete fopr_r2;
369  delete force_fopr_r2;
370  delete action_F_Nf1;
371 
372  delete fopr_langev;
373  delete fopr_H;
374  delete force_fopr_MD;
375 
376  delete builder;
377 
378 
379 #ifdef USE_TEST
380  return Test::verify(expected_result, result);
381 
382 #else
383  return EXIT_SUCCESS;
384 #endif
385  }
386 } // namespace Test_HMC_Clover_SF