Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
action_G_Rectangle_SF.cpp
Go to the documentation of this file.
1 
14 #include "action_G_Rectangle_SF.h"
15 
16 #ifdef USE_PARAMETERS_FACTORY
17 #include "parameters_factory.h"
18 #endif
19 
20 //- parameter entries
21 namespace {
22  void append_entry(Parameters& param)
23  {
24  param.Register_double("beta", 0.0);
25  param.Register_double("c_plaq", 0.0);
26  param.Register_double("c_rect", 0.0);
27 
28  param.Register_double("ct0", 0.0);
29  param.Register_double("ct1", 0.0);
30  param.Register_double("ct2", 0.0);
31 
32  param.Register_double("ctr0", 0.0);
33  param.Register_double("ctr1", 0.0);
34  param.Register_double("ctr2", 0.0);
35 
36  param.Register_double_vector("phi", std::valarray<double>());
37  param.Register_double_vector("phipr", std::valarray<double>());
38 
39  // param.Register_double("ct" , 0.0);
40  // param.Register_double("ctr", 0.0);
41 
42  param.Register_string("verbose_level", "NULL");
43  }
44 
45 
46 #ifdef USE_PARAMETERS_FACTORY
47  bool init_param = ParametersFactory::Register("Action.G_Rectangle_SF", append_entry);
48 #endif
49 }
50 //- end
51 
52 //- parameters class
53 Parameters_Action_G_Rectangle_SF::Parameters_Action_G_Rectangle_SF() { append_entry(*this); }
54 //- end
55 
56 //====================================================================
57 void Action_G_Rectangle_SF::set_parameters(const Parameters& params)
58 {
59  const string str_vlevel = params.get_string("verbose_level");
60 
61  m_vl = vout.set_verbose_level(str_vlevel);
62 
63  //- fetch and check input parameters
64  double beta, c_plaq, c_rect;
65  double ct0, ct1, ct2, ctr0, ctr1, ctr2;
66  std::valarray<double> phi, phipr;
67 
68  int err = 0;
69  err += params.fetch_double("beta", beta);
70  err += params.fetch_double("c_plaq", c_plaq);
71  err += params.fetch_double("c_rect", c_rect);
72 
73  err += params.fetch_double("ct0", ct0);
74  err += params.fetch_double("ct1", ct1);
75  err += params.fetch_double("ct2", ct2);
76 
77  err += params.fetch_double("ctr0", ctr0);
78  err += params.fetch_double("ctr1", ctr1);
79  err += params.fetch_double("ctr2", ctr2);
80 
81  err += params.fetch_double_vector("phi", phi);
82  err += params.fetch_double_vector("phipr", phipr);
83 
84  if (err) {
85  vout.crucial(m_vl, "Action_G_Rectangle_SF: fetch error, input parameter not found.\n");
86  abort();
87  }
88 
89 
90  double gg = 6.0 / beta;
91 
92  double ct = ct0 + ct1 * gg + ct2 * gg * gg;
93  double ctr = ctr0 + ctr1 * gg + ctr2 * gg * gg;
94 
95  set_parameters(beta, c_plaq, c_rect, &phi[0], &phipr[0], ct, ctr);
96 }
97 
98 
99 //====================================================================
100 
114 void Action_G_Rectangle_SF::set_parameters(double beta, double c_plaq, double c_rect,
115  double *phi, double *phipr, double ct, double ctr)
116 {
117  //- print input parameters
118  vout.general(m_vl, "Action_G_Rectangle_SF:\n");
119  vout.general(m_vl, " beta = %12.6f\n", beta);
120  vout.general(m_vl, " c_plaq = %12.6f\n", c_plaq);
121  vout.general(m_vl, " c_rect = %12.6f\n", c_rect);
122  vout.general(m_vl, " phi1 = %12.6f\n", phi[0]);
123  vout.general(m_vl, " phi2 = %12.6f\n", phi[1]);
124  vout.general(m_vl, " phi3 = %12.6f\n", phi[2]);
125  vout.general(m_vl, " phipr1 = %12.6f\n", phipr[0]);
126  vout.general(m_vl, " phipr2 = %12.6f\n", phipr[1]);
127  vout.general(m_vl, " phipr3 = %12.6f\n", phipr[2]);
128  vout.general(m_vl, " ct = %12.6f\n", ct);
129  vout.general(m_vl, " ctr = %12.6f\n", ctr);
130 
131  //- range check
132  // NB. beta,c_plaq,c_rect,phi,phipr,ct,ctr = 0 is allowed.
133 
134  //- store values
135  m_beta = beta;
136  m_c_plaq = c_plaq;
137  m_c_rect = c_rect;
138 
139  m_ct = ct;
140  m_ctr = ctr;
141  // m_phi = phi ;
142  // m_phipr = phipr;
143 
144  //- post-process
145  m_staple.set_parameters(phi, phipr);
146 
147  int Lx = CommonParameters::Lx();
148  double c0r, c0i, c1r, c1i, c2r, c2i;
149  c0r = cos(phi[0] / Lx);
150  c0i = sin(phi[0] / Lx);
151  c1r = cos(phi[1] / Lx);
152  c1i = sin(phi[1] / Lx);
153  c2r = cos(phi[2] / Lx);
154  c2i = sin(phi[2] / Lx);
155  wk.zero();
156  wk.set(0, 0, c0r, c0i);
157  wk.set(1, 1, c1r, c1i);
158  wk.set(2, 2, c2r, c2i);
159 
160  c0r = cos(phipr[0] / Lx);
161  c0i = sin(phipr[0] / Lx);
162  c1r = cos(phipr[1] / Lx);
163  c1i = sin(phipr[1] / Lx);
164  c2r = cos(phipr[2] / Lx);
165  c2i = sin(phipr[2] / Lx);
166  wkpr.zero();
167  wkpr.set(0, 0, c0r, c0i);
168  wkpr.set(1, 1, c1r, c1i);
169  wkpr.set(2, 2, c2r, c2i);
170 
171 
172  int Nvol = CommonParameters::Nvol();
173  int Ndim = CommonParameters::Ndim();
174  int NinG = 2 * Nc * Nc;
175  m_force.reset(NinG, Nvol, Ndim);
176 
177  m_status_linkv = 0;
178 }
179 
180 
181 //====================================================================
182 double Action_G_Rectangle_SF::langevin(RandomNumbers *rand)
183 {
184  m_status_linkv = 0;
185 
186  double H_U = calcH();
187 
188  return H_U;
189 }
190 
191 
192 //====================================================================
193 
235 double Action_G_Rectangle_SF::calcH()
236 {
237  int Ndim = CommonParameters::Ndim();
238  int Nvol = CommonParameters::Nvol();
239  int Lvol = CommonParameters::Lvol();
240 
241  int Nx = CommonParameters::Nx();
242  int Ny = CommonParameters::Ny();
243  int Nz = CommonParameters::Nz();
244  int Nt = CommonParameters::Nt();
245  int NPEt = CommonParameters::NPEt();
246 
247  int Ndim2 = Ndim * (Ndim - 1) / 2;
248  int size_U = Lvol * Ndim2;
249 
250  //Field_G Cup1(Nvol,1), Cup2(Nvol,1);
251  //Field_G Cdn1(Nvol,1), Cdn2(Nvol,1);
252  // Field_G Umu(Nvol,1), Unu(Nvol,1);
253  // Field_G v(Nvol,1), w(Nvol,1), c(Nvol,1);
254  Field_G_SF Cup1, Cup2;
255  Field_G_SF Cdn1, Cdn2;
256  Field_G_SF Umu, Unu;
257  Field_G_SF v, w, c;
258 
259  Mat_SU_N vmat(Nc), wmat(Nc), cmat(Nc);
260 
261  double plaqF = 0.0;
262  double rectF = 0.0;
263 
264  vout.general(m_vl, " Action_G_Rectangle_SF: %s\n", m_label.c_str());
265 
266  for (int mu = 0; mu < Ndim; ++mu) {
267  for (int nu = mu + 1; nu < Ndim; ++nu) {
268  Cup1 = m_staple.upper(*m_U, mu, nu);
269  Cup2 = m_staple.upper(*m_U, nu, mu);
270 
271  // plaquette term
272  Unu = Cup2;
273  // If the node is at the boundary the temporal plaquette is multiplied with ct.
274  if ((nu == 3) && (Communicator::ipe(3) == 0)) {
275  Unu.mult_ct_boundary(0, m_ct);
276  }
277  if ((nu == 3) && (Communicator::ipe(3) == NPEt - 1)) {
278  Unu.mult_ct_boundary(Nt - 1, m_ct);
279  }
280  for (int site = 0; site < Nvol; ++site) {
281  plaqF += ReTr(m_U->mat(site, nu) * Unu.mat_dag(site));
282  }
283 
284  // rectangular terms
285 
286  // +---+---+
287  // | | term
288  // x <---+
289 
290  Umu.setpart_ex(0, *m_U, mu);
291  Unu.setpart_ex(0, *m_U, nu);
292  if ((Communicator::ipe(3) == 0) && (nu == 3)) {
293  Umu.set_boundary_wk(wk);
294  }
295 
296  m_shift.backward(v, Cup2, mu);
297  m_shift.backward(c, Umu, nu);
298 
299  if ((Communicator::ipe(3) == 0) && (nu == 3)) {
300  c.mult_ct_boundary(0, m_ctr);
301  }
302  if ((Communicator::ipe(3) == NPEt - 1) && (nu == 3)) {
303  c.set_boundary_wkpr(wkpr);
304  c.mult_ct_boundary(Nt - 1, m_ctr);
305  }
306 
307  w.mult_Field_Gnd(0, c, 0, v, 0);
308  c.mult_Field_Gnn(0, Unu, 0, w, 0);
309  for (int site = 0; site < Nvol; ++site) {
310  rectF += ReTr(Umu.mat(site) * c.mat_dag(site));
311  }
312 
313  // +---+
314  // | |
315  // + + term
316  // | |
317  // x v
318 
319  m_shift.backward(v, Unu, mu);
320  m_shift.backward(c, Cup1, nu);
321 
322  w.mult_Field_Gnd(0, c, 0, v, 0);
323  c.mult_Field_Gnn(0, Unu, 0, w, 0);
324  for (int site = 0; site < Nvol; ++site) {
325  rectF += ReTr(Umu.mat(site) * c.mat_dag(site));
326  }
327  }
328  }
329 
330  plaqF = Communicator::reduce_sum(plaqF);
331  rectF = Communicator::reduce_sum(rectF);
332 
333  // double plaq = plaqF/Nc;
334  // vout.general(m_vl," Plaquette = %18.8f\n",plaq/size_U);
335 
336  // double H_U = m_c_plaq * (Ndim2*Lvol - plaqF/Nc)
337  // + m_c_rect * (Ndim2*Lvol*2 - rectF/Nc);
338  double H_U = m_c_plaq * (-plaqF / Nc)
339  + m_c_rect * (-rectF / Nc);
340 
341  H_U = m_beta * H_U;
342 
343  vout.general(m_vl, " H_Grectangle = %18.8f\n", H_U);
344  // vout.general(m_vl," H_G/dof = %18.8f\n",H_U/size_U);
345 
346  return H_U;
347 }
348 
349 
350 //====================================================================
351 
406 const Field Action_G_Rectangle_SF::force()
407 {
408  if (m_status_linkv == 0) {
409  int Nvol = CommonParameters::Nvol();
410  int Ndim = CommonParameters::Ndim();
411 
412  int Nx = CommonParameters::Nx();
413  int Ny = CommonParameters::Ny();
414  int Nz = CommonParameters::Nz();
415  int Nt = CommonParameters::Nt();
416  int NPEt = CommonParameters::NPEt();
417 
418  assert(m_U->nin() == Nc * Nc * 2);
419  assert(m_U->nvol() == Nvol);
420  assert(m_U->nex() == Ndim);
421 
422  double betaNc = m_beta / Nc;
423  Mat_SU_N ut(Nc);
424 
425  Field_G force(Nvol, Ndim);
426  Field_G force1(Nvol, 1);
427 
428  // Field_G Cup1(Nvol,1), Cup2(Nvol,1);
429  // Field_G Cdn1(Nvol,1), Cdn2(Nvol,1);
430  // Field_G Umu(Nvol,1), Unu(Nvol,1);
431  // Field_G v(Nvol,1), w(Nvol,1), c(Nvol,1);
432  Field_G_SF Cup1, Cup2;
433  Field_G_SF Cdn1, Cdn2;
434  Field_G_SF Umu, Unu;
435  Field_G_SF v, w, c;
436 
437  Mat_SU_N vmat(Nc), wmat(Nc), cmat(Nc);
438  Mat_SU_N wzero(Nc);
439  wzero.zero();
440 
441  vout.general(m_vl, " Action_G_Rectangle_SF: %s\n", m_label.c_str());
442 
443  for (int mu = 0; mu < Ndim; ++mu) {
444  force1 = 0.0;
445  for (int nu = 0; nu < Ndim; ++nu) {
446  if (nu == mu) continue;
447 
448  Cup1 = m_staple.upper(*m_U, mu, nu);
449  Cup2 = m_staple.upper(*m_U, nu, mu);
450  Cdn1 = m_staple.lower(*m_U, mu, nu);
451  Cdn2 = m_staple.lower(*m_U, nu, mu);
452 
453  // plaquette term
454  Umu = Cup1;
455  Unu = Cdn1;
456  if (Communicator::ipe(3) == 0) {
457  if (mu == 3) {
458  Umu.mult_ct_boundary(0, m_ct);
459  Unu.mult_ct_boundary(0, m_ct);
460  }
461  if (nu == 3) {
462  Unu.mult_ct_boundary(1, m_ct);
463  }
464  }
465  if (Communicator::ipe(3) == NPEt - 1) {
466  if (mu == 3) {
467  Umu.mult_ct_boundary(Nt - 1, m_ct);
468  Unu.mult_ct_boundary(Nt - 1, m_ct);
469  }
470  if (nu == 3) {
471  Umu.mult_ct_boundary(Nt - 1, m_ct);
472  }
473  }
474  force1.addpart_ex(0, Umu, 0, m_c_plaq);
475  force1.addpart_ex(0, Unu, 0, m_c_plaq);
476 
477  // rectangular term
478  Umu.setpart_ex(0, *m_U, mu);
479  Unu.setpart_ex(0, *m_U, nu);
480  // For the boundary spatial link, use Wk.
481  if ((Communicator::ipe(3) == 0) && (nu == 3)) Umu.set_boundary_wk(wk);
482  if ((Communicator::ipe(3) == 0) && (mu == 3)) Unu.set_boundary_wk(wk);
483 
484  // +---+---+
485  // | | term
486  // x <---+
487 
488  m_shift.backward(v, Cup2, mu);
489  m_shift.backward(c, Umu, nu);
490  // The force for the spatial link near the boundary. Multiplied with ctr.
491  if ((Communicator::ipe(3) == NPEt - 1) && (nu == 3)) {
492  c.set_boundary_wkpr(wkpr);
493  c.mult_ct_boundary(Nt - 1, m_ctr);
494  }
495  w.mult_Field_Gnd(0, c, 0, v, 0);
496  c.mult_Field_Gnn(0, Unu, 0, w, 0);
497  // The force for the boundary spatial link is set to zero.
498  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
499  force1.addpart_ex(0, c, 0, m_c_rect);
500 
501  // +---+
502  // | |
503  // + + term
504  // | |
505  // x v
506 
507  m_shift.backward(v, Unu, mu);
508  m_shift.backward(c, Cup1, nu);
509  // The force for the temporal link near the boundary. Multiplied with ctr.
510  if ((Communicator::ipe(3) == NPEt - 1) && (mu == 3)) {
511  v.set_boundary_wkpr(wkpr);
512  v.mult_ct_boundary(Nt - 1, m_ctr);
513  }
514  w.mult_Field_Gnd(0, c, 0, v, 0);
515  c.mult_Field_Gnn(0, Unu, 0, w, 0);
516  // The force for the boundary spatial link is set to zero.
517  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
518  // The force for the boundary temporal link.
519  if ((Communicator::ipe(3) == 0) && (mu == 3)) c.mult_ct_boundary(0, m_ctr);
520  force1.addpart_ex(0, c, 0, m_c_rect);
521 
522  // +---+---+
523  // | | term
524  // +---x v
525 
526  m_shift.backward(v, Unu, mu);
527  m_shift.backward(c, Umu, nu);
528  if ((Communicator::ipe(3) == NPEt - 1) && (nu == 3)) {
529  c.set_boundary_wkpr(wkpr);
530  c.mult_ct_boundary(Nt - 1, m_ctr);
531  }
532  if ((Communicator::ipe(3) == NPEt - 1) && (mu == 3)) v.set_boundary_wkpr(wkpr);
533  w.mult_Field_Gnd(0, c, 0, v, 0);
534  c.mult_Field_Gnn(0, Cdn2, 0, w, 0);
535  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
536  force1.addpart_ex(0, c, 0, m_c_rect);
537 
538  // x <---+
539  // | | term
540  // +---+---+
541 
542  m_shift.backward(v, Cup2, mu);
543  w.mult_Field_Gnn(0, Umu, 0, v, 0);
544  v.mult_Field_Gdn(0, Unu, 0, w, 0);
545  if ((Communicator::ipe(3) == 0) && (nu == 3)) v.mult_ct_boundary(0, m_ctr);
546  m_shift.forward(c, v, nu);
547  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
548  force1.addpart_ex(0, c, 0, m_c_rect);
549 
550  // x ^
551  // | |
552  // + + term
553  // | |
554  // +---+
555 
556  m_shift.backward(v, Unu, mu);
557  if ((Communicator::ipe(3) == NPEt - 1) && (mu == 3)) {
558  v.set_boundary_wkpr(wkpr);
559  v.mult_ct_boundary(Nt - 1, m_ctr);
560  }
561  w.mult_Field_Gnn(0, Cdn1, 0, v, 0);
562  v.mult_Field_Gdn(0, Unu, 0, w, 0);
563  if ((Communicator::ipe(3) == 0) && (mu == 3)) v.mult_ct_boundary(0, m_ctr);
564  m_shift.forward(c, v, nu);
565  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
566  force1.addpart_ex(0, c, 0, m_c_rect);
567 
568  // +---x ^
569  // | | term
570  // +---+---+
571 
572  m_shift.backward(v, Unu, mu);
573  if ((Communicator::ipe(3) == NPEt - 1) && (mu == 3)) v.set_boundary_wkpr(wkpr);
574  w.mult_Field_Gnn(0, Umu, 0, v, 0);
575  v.mult_Field_Gdn(0, Cdn2, 0, w, 0);
576  if ((Communicator::ipe(3) == 0) && (nu == 3)) v.mult_ct_boundary(0, m_ctr);
577  m_shift.forward(c, v, nu);
578  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
579  force1.addpart_ex(0, c, 0, m_c_rect);
580  }
581 
582  force.mult_Field_Gnd(mu, *m_U, mu, force1, 0);
583  force.at_Field_G(mu);
584  }
585 
586  force *= -betaNc;
587 
588  m_force = (Field)force;
589  ++m_status_linkv;
590 
591  double Fave, Fmax, Fdev;
592  m_force.stat(Fave, Fmax, Fdev);
593  vout.general(m_vl, " Fave = %12.6f Fmax = %10.6f Fdev = %12.6f\n",
594  Fave, Fmax, Fdev);
595 
596  return m_force;
597  } else {
598  vout.general(m_vl, " Action_G_Rectangle_SF returns previous force.\n");
599  return m_force;
600  }
601 }
602 
603 
604 //====================================================================
605 //============================================================END=====