Bridge++  Ver. 1.3.x
action_G_Rectangle_SF.cpp
Go to the documentation of this file.
1 
14 #include "action_G_Rectangle_SF.h"
15 
16 
17 #ifdef USE_FACTORY
18 namespace {
19  Action *create_object()
20  {
21  return new Action_G_Rectangle_SF();
22  }
23 
24 
25  bool init = Action::Factory::Register("Action_G_Rectangle_SF", create_object);
26 }
27 #endif
28 
29 
30 //- parameter entries
31 namespace {
32  void append_entry(Parameters& param)
33  {
34  param.Register_double("beta", 0.0);
35  param.Register_double("c_plaq", 0.0);
36  param.Register_double("c_rect", 0.0);
37 
38  param.Register_double("ct0", 0.0);
39  param.Register_double("ct1", 0.0);
40  param.Register_double("ct2", 0.0);
41 
42  param.Register_double("ctr0", 0.0);
43  param.Register_double("ctr1", 0.0);
44  param.Register_double("ctr2", 0.0);
45 
46  param.Register_double_vector("phi", std::vector<double>());
47  param.Register_double_vector("phipr", std::vector<double>());
48 
49  // param.Register_double("ct" , 0.0);
50  // param.Register_double("ctr", 0.0);
51 
52  param.Register_string("verbose_level", "NULL");
53  }
54 
55 
56 #ifdef USE_PARAMETERS_FACTORY
57  bool init_param = ParametersFactory::Register("Action.G_Rectangle_SF", append_entry);
58 #endif
59 }
60 //- end
61 
62 //- parameters class
64 //- end
65 
66 const std::string Action_G_Rectangle_SF::class_name = "Action_G_Rectangle_SF";
67 
68 //====================================================================
70 {
71  const string str_vlevel = params.get_string("verbose_level");
72 
73  m_vl = vout.set_verbose_level(str_vlevel);
74 
75  //- fetch and check input parameters
76  double beta, c_plaq, c_rect;
77  double ct0, ct1, ct2, ctr0, ctr1, ctr2;
78  std::vector<double> phi, phipr;
79 
80  int err = 0;
81  err += params.fetch_double("beta", beta);
82  err += params.fetch_double("c_plaq", c_plaq);
83  err += params.fetch_double("c_rect", c_rect);
84 
85  err += params.fetch_double("ct0", ct0);
86  err += params.fetch_double("ct1", ct1);
87  err += params.fetch_double("ct2", ct2);
88 
89  err += params.fetch_double("ctr0", ctr0);
90  err += params.fetch_double("ctr1", ctr1);
91  err += params.fetch_double("ctr2", ctr2);
92 
93  err += params.fetch_double_vector("phi", phi);
94  err += params.fetch_double_vector("phipr", phipr);
95 
96  if (err) {
97  vout.crucial(m_vl, "%s: fetch error, input parameter not found.\n", class_name.c_str());
98  exit(EXIT_FAILURE);
99  }
100 
101 
102  double gg = 6.0 / beta;
103 
104  double ct = ct0 + ct1 * gg + ct2 * gg * gg;
105  double ctr = ctr0 + ctr1 * gg + ctr2 * gg * gg;
106 
107  set_parameters(beta, c_plaq, c_rect, &phi[0], &phipr[0], ct, ctr);
108 }
109 
110 
111 //====================================================================
112 
126 void Action_G_Rectangle_SF::set_parameters(double beta, double c_plaq, double c_rect,
127  double *phi, double *phipr, double ct, double ctr)
128 {
129  //- print input parameters
130  vout.general(m_vl, "%s:\n", class_name.c_str());
131  vout.general(m_vl, " beta = %12.6f\n", beta);
132  vout.general(m_vl, " c_plaq = %12.6f\n", c_plaq);
133  vout.general(m_vl, " c_rect = %12.6f\n", c_rect);
134  vout.general(m_vl, " phi1 = %12.6f\n", phi[0]);
135  vout.general(m_vl, " phi2 = %12.6f\n", phi[1]);
136  vout.general(m_vl, " phi3 = %12.6f\n", phi[2]);
137  vout.general(m_vl, " phipr1 = %12.6f\n", phipr[0]);
138  vout.general(m_vl, " phipr2 = %12.6f\n", phipr[1]);
139  vout.general(m_vl, " phipr3 = %12.6f\n", phipr[2]);
140  vout.general(m_vl, " ct = %12.6f\n", ct);
141  vout.general(m_vl, " ctr = %12.6f\n", ctr);
142 
143  //- range check
144  // NB. beta,c_plaq,c_rect,phi,phipr,ct,ctr = 0 is allowed.
145 
146  //- store values
147  m_beta = beta;
148  m_c_plaq = c_plaq;
149  m_c_rect = c_rect;
150 
151  m_ct = ct;
152  m_ctr = ctr;
153  // m_phi = phi ;
154  // m_phipr = phipr;
155 
156  //- post-process
157  m_staple.set_parameters(phi, phipr);
158 
159  int Lx = CommonParameters::Lx();
160  double c0r, c0i, c1r, c1i, c2r, c2i;
161  c0r = cos(phi[0] / Lx);
162  c0i = sin(phi[0] / Lx);
163  c1r = cos(phi[1] / Lx);
164  c1i = sin(phi[1] / Lx);
165  c2r = cos(phi[2] / Lx);
166  c2i = sin(phi[2] / Lx);
167  wk.zero();
168  wk.set(0, 0, c0r, c0i);
169  wk.set(1, 1, c1r, c1i);
170  wk.set(2, 2, c2r, c2i);
171 
172  c0r = cos(phipr[0] / Lx);
173  c0i = sin(phipr[0] / Lx);
174  c1r = cos(phipr[1] / Lx);
175  c1i = sin(phipr[1] / Lx);
176  c2r = cos(phipr[2] / Lx);
177  c2i = sin(phipr[2] / Lx);
178  wkpr.zero();
179  wkpr.set(0, 0, c0r, c0i);
180  wkpr.set(1, 1, c1r, c1i);
181  wkpr.set(2, 2, c2r, c2i);
182 
183 
184  int Nvol = CommonParameters::Nvol();
185  int Ndim = CommonParameters::Ndim();
186  int NinG = 2 * Nc * Nc;
187 }
188 
189 
190 //====================================================================
192 {
193  double H_U = calcH();
194 
195  return H_U;
196 }
197 
198 
199 //====================================================================
200 
243 {
244  int Ndim = CommonParameters::Ndim();
245  int Nvol = CommonParameters::Nvol();
246  int Lvol = CommonParameters::Lvol();
247 
248  int Nx = CommonParameters::Nx();
249  int Ny = CommonParameters::Ny();
250  int Nz = CommonParameters::Nz();
251  int Nt = CommonParameters::Nt();
252  int NPEt = CommonParameters::NPEt();
253 
254  // int Ndim2 = Ndim * (Ndim - 1) / 2;
255  // int size_U = Lvol * Ndim2;
256 
257  //Field_G Cup1(Nvol,1), Cup2(Nvol,1);
258  //Field_G Cdn1(Nvol,1), Cdn2(Nvol,1);
259  // Field_G Umu(Nvol,1), Unu(Nvol,1);
260  // Field_G v(Nvol,1), w(Nvol,1), c(Nvol,1);
261  Field_G_SF Cup1, Cup2;
262  Field_G_SF Cdn1, Cdn2;
263  Field_G_SF Umu, Unu;
264  Field_G_SF v, w, c;
265 
266  Mat_SU_N vmat(Nc), wmat(Nc), cmat(Nc);
267 
268  double plaqF = 0.0;
269  double rectF = 0.0;
270 
271  vout.general(m_vl, " %s: %s\n", class_name.c_str(), m_label.c_str());
272 
273  for (int mu = 0; mu < Ndim; ++mu) {
274  for (int nu = mu + 1; nu < Ndim; ++nu) {
275  m_staple.upper(Cup1, *m_U, mu, nu);
276  m_staple.upper(Cup2, *m_U, nu, mu);
277 
278  // plaquette term
279  Unu = Cup2;
280  // If the node is at the boundary the temporal plaquette is multiplied with ct.
281  if ((nu == 3) && (Communicator::ipe(3) == 0)) {
282  Unu.mult_ct_boundary(0, m_ct);
283  }
284  if ((nu == 3) && (Communicator::ipe(3) == NPEt - 1)) {
285  Unu.mult_ct_boundary(Nt - 1, m_ct);
286  }
287  for (int site = 0; site < Nvol; ++site) {
288  plaqF += ReTr(m_U->mat(site, nu) * Unu.mat_dag(site));
289  }
290 
291  // rectangular terms
292 
293  // +---+---+
294  // | | term
295  // x <---+
296 
297  Umu.setpart_ex(0, *m_U, mu);
298  Unu.setpart_ex(0, *m_U, nu);
299  if ((Communicator::ipe(3) == 0) && (nu == 3)) {
300  Umu.set_boundary_wk(wk);
301  }
302 
303  m_shift.backward(v, Cup2, mu);
304  m_shift.backward(c, Umu, nu);
305 
306  if ((Communicator::ipe(3) == 0) && (nu == 3)) {
307  c.mult_ct_boundary(0, m_ctr);
308  }
309  if ((Communicator::ipe(3) == NPEt - 1) && (nu == 3)) {
311  c.mult_ct_boundary(Nt - 1, m_ctr);
312  }
313 
314  mult_Field_Gnd(w, 0, c, 0, v, 0);
315  mult_Field_Gnn(c, 0, Unu, 0, w, 0);
316  for (int site = 0; site < Nvol; ++site) {
317  rectF += ReTr(Umu.mat(site) * c.mat_dag(site));
318  }
319 
320  // +---+
321  // | |
322  // + + term
323  // | |
324  // x v
325 
326  m_shift.backward(v, Unu, mu);
327  m_shift.backward(c, Cup1, nu);
328 
329  mult_Field_Gnd(w, 0, c, 0, v, 0);
330  mult_Field_Gnn(c, 0, Unu, 0, w, 0);
331  for (int site = 0; site < Nvol; ++site) {
332  rectF += ReTr(Umu.mat(site) * c.mat_dag(site));
333  }
334  }
335  }
336 
337  plaqF = Communicator::reduce_sum(plaqF);
338  rectF = Communicator::reduce_sum(rectF);
339 
340  // double plaq = plaqF/Nc;
341  // vout.general(m_vl," Plaquette = %18.8f\n",plaq/size_U);
342 
343  // double H_U = m_c_plaq * (Ndim2*Lvol - plaqF/Nc)
344  // + m_c_rect * (Ndim2*Lvol*2 - rectF/Nc);
345  double H_U = m_c_plaq * (-plaqF / Nc)
346  + m_c_rect * (-rectF / Nc);
347 
348  H_U = m_beta * H_U;
349 
350  vout.general(m_vl, " H_Grectangle = %18.8f\n", H_U);
351  // vout.general(m_vl," H_G/dof = %18.8f\n",H_U/size_U);
352 
353  return H_U;
354 }
355 
356 
357 //====================================================================
358 
414 {
415  int Nvol = CommonParameters::Nvol();
416  int Ndim = CommonParameters::Ndim();
417 
418  int Nx = CommonParameters::Nx();
419  int Ny = CommonParameters::Ny();
420  int Nz = CommonParameters::Nz();
421  int Nt = CommonParameters::Nt();
422  int NPEt = CommonParameters::NPEt();
423 
424  assert(m_U->nin() == Nc * Nc * 2);
425  assert(m_U->nvol() == Nvol);
426  assert(m_U->nex() == Ndim);
427 
428  assert(force.nin() == Nc * Nc * 2);
429  assert(force.nvol() == Nvol);
430  assert(force.nex() == Ndim);
431 
432  force.set(0.0);
433 
434  double betaNc = m_beta / Nc;
435  Mat_SU_N ut(Nc);
436 
437  Field_G force1(Nvol, 1);
438  Field_G force2(Nvol, 1);
439 
440  Field_G_SF Cup1, Cup2;
441  Field_G_SF Cdn1, Cdn2;
442  Field_G_SF Umu, Unu;
443  Field_G_SF v, w, c;
444 
445  Mat_SU_N vmat(Nc), wmat(Nc), cmat(Nc);
446  Mat_SU_N wzero(Nc);
447  wzero.zero();
448 
449  vout.general(m_vl, " %s: %s\n", class_name.c_str(), m_label.c_str());
450 
451  for (int mu = 0; mu < Ndim; ++mu) {
452  force1.set(0.0);
453  for (int nu = 0; nu < Ndim; ++nu) {
454  if (nu == mu) continue;
455 
456  m_staple.upper(Cup1, *m_U, mu, nu);
457  m_staple.upper(Cup2, *m_U, nu, mu);
458  m_staple.lower(Cdn1, *m_U, mu, nu);
459  m_staple.lower(Cdn2, *m_U, nu, mu);
460 
461  // plaquette term
462  Umu = Cup1;
463  Unu = Cdn1;
464  if (Communicator::ipe(3) == 0) {
465  if (mu == 3) {
466  Umu.mult_ct_boundary(0, m_ct);
467  Unu.mult_ct_boundary(0, m_ct);
468  }
469  if (nu == 3) {
470  Unu.mult_ct_boundary(1, m_ct);
471  }
472  }
473  if (Communicator::ipe(3) == NPEt - 1) {
474  if (mu == 3) {
475  Umu.mult_ct_boundary(Nt - 1, m_ct);
476  Unu.mult_ct_boundary(Nt - 1, m_ct);
477  }
478  if (nu == 3) {
479  Umu.mult_ct_boundary(Nt - 1, m_ct);
480  }
481  }
482  force1.addpart_ex(0, Umu, 0, m_c_plaq);
483  force1.addpart_ex(0, Unu, 0, m_c_plaq);
484 
485  // rectangular term
486  Umu.setpart_ex(0, *m_U, mu);
487  Unu.setpart_ex(0, *m_U, nu);
488  // For the boundary spatial link, use Wk.
489  if ((Communicator::ipe(3) == 0) && (nu == 3)) Umu.set_boundary_wk(wk);
490  if ((Communicator::ipe(3) == 0) && (mu == 3)) Unu.set_boundary_wk(wk);
491 
492  // +---+---+
493  // | | term
494  // x <---+
495 
496  m_shift.backward(v, Cup2, mu);
497  m_shift.backward(c, Umu, nu);
498  // The force for the spatial link near the boundary. Multiplied with ctr.
499  if ((Communicator::ipe(3) == NPEt - 1) && (nu == 3)) {
501  c.mult_ct_boundary(Nt - 1, m_ctr);
502  }
503  mult_Field_Gnd(w, 0, c, 0, v, 0);
504  mult_Field_Gnn(c, 0, Unu, 0, w, 0);
505  // The force for the boundary spatial link is set to zero.
506  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
507  force1.addpart_ex(0, c, 0, m_c_rect);
508 
509  // +---+
510  // | |
511  // + + term
512  // | |
513  // x v
514 
515  m_shift.backward(v, Unu, mu);
516  m_shift.backward(c, Cup1, nu);
517  // The force for the temporal link near the boundary. Multiplied with ctr.
518  if ((Communicator::ipe(3) == NPEt - 1) && (mu == 3)) {
520  v.mult_ct_boundary(Nt - 1, m_ctr);
521  }
522  mult_Field_Gnd(w, 0, c, 0, v, 0);
523  mult_Field_Gnn(c, 0, Unu, 0, w, 0);
524  // The force for the boundary spatial link is set to zero.
525  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
526  // The force for the boundary temporal link.
527  if ((Communicator::ipe(3) == 0) && (mu == 3)) c.mult_ct_boundary(0, m_ctr);
528  force1.addpart_ex(0, c, 0, m_c_rect);
529 
530  // +---+---+
531  // | | term
532  // +---x v
533 
534  m_shift.backward(v, Unu, mu);
535  m_shift.backward(c, Umu, nu);
536  if ((Communicator::ipe(3) == NPEt - 1) && (nu == 3)) {
538  c.mult_ct_boundary(Nt - 1, m_ctr);
539  }
540  if ((Communicator::ipe(3) == NPEt - 1) && (mu == 3)) v.set_boundary_wkpr(wkpr);
541  mult_Field_Gnd(w, 0, c, 0, v, 0);
542  mult_Field_Gnn(c, 0, Cdn2, 0, w, 0);
543  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
544  force1.addpart_ex(0, c, 0, m_c_rect);
545 
546  // x <---+
547  // | | term
548  // +---+---+
549 
550  m_shift.backward(v, Cup2, mu);
551  mult_Field_Gnn(w, 0, Umu, 0, v, 0);
552  mult_Field_Gdn(v, 0, Unu, 0, w, 0);
553  if ((Communicator::ipe(3) == 0) && (nu == 3)) v.mult_ct_boundary(0, m_ctr);
554  m_shift.forward(c, v, nu);
555  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
556  force1.addpart_ex(0, c, 0, m_c_rect);
557 
558  // x ^
559  // | |
560  // + + term
561  // | |
562  // +---+
563 
564  m_shift.backward(v, Unu, mu);
565  if ((Communicator::ipe(3) == NPEt - 1) && (mu == 3)) {
567  v.mult_ct_boundary(Nt - 1, m_ctr);
568  }
569  mult_Field_Gnn(w, 0, Cdn1, 0, v, 0);
570  mult_Field_Gdn(v, 0, Unu, 0, w, 0);
571  if ((Communicator::ipe(3) == 0) && (mu == 3)) v.mult_ct_boundary(0, m_ctr);
572  m_shift.forward(c, v, nu);
573  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
574  force1.addpart_ex(0, c, 0, m_c_rect);
575 
576  // +---x ^
577  // | | term
578  // +---+---+
579 
580  m_shift.backward(v, Unu, mu);
581  if ((Communicator::ipe(3) == NPEt - 1) && (mu == 3)) v.set_boundary_wkpr(wkpr);
582  mult_Field_Gnn(w, 0, Umu, 0, v, 0);
583  mult_Field_Gdn(v, 0, Cdn2, 0, w, 0);
584  if ((Communicator::ipe(3) == 0) && (nu == 3)) v.mult_ct_boundary(0, m_ctr);
585  m_shift.forward(c, v, nu);
586  if ((Communicator::ipe(3) == 0) && (nu == 3)) c.set_boundary_zero();
587  force1.addpart_ex(0, c, 0, m_c_rect);
588  }
589 
590  mult_Field_Gnd(force2, 0, *m_U, mu, force1, 0);
591  at_Field_G(force2, 0);
592 
593  axpy(force, mu, -betaNc, force2, 0);
594  }
595 
596  double Fave, Fmax, Fdev;
597  force.stat(Fave, Fmax, Fdev);
598  vout.general(m_vl, " Fave = %12.6f Fmax = %10.6f Fdev = %12.6f\n",
599  Fave, Fmax, Fdev);
600 }
601 
602 
603 //====================================================================
604 //============================================================END=====
void set_boundary_zero()
Set the boundary matrix to 0 for SF bc.
Definition: field_G_SF.cpp:101
BridgeIO vout
Definition: bridgeIO.cpp:278
void set_boundary_wkpr(const Mat_SU_N &U)
Set the boundary spatial link at t=Nt-1 for SF bc.
Definition: field_G_SF.cpp:76
void Register_string(const string &, const string &)
Definition: parameters.cpp:351
static int NPEt()
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:155
void upper(Field_G_SF &, const Field_G &, const int, const int)
Definition: staples_SF.cpp:860
void set_parameters(const Parameters &params)
Definition: staples_SF.cpp:62
void mult_Field_Gdn(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2)
void general(const char *format,...)
Definition: bridgeIO.cpp:65
Mat_SU_N & zero()
Definition: mat_SU_N.h:383
Container of Field-type object.
Definition: field.h:39
int nvol() const
Definition: field.h:116
Class for parameters.
Definition: parameters.h:38
static int ipe(const int dir)
logical coordinate of current proc.
static int Lvol()
void addpart_ex(int ex, const Field &w, int exw)
Definition: field.h:189
Base class of HMC action class family.
Definition: action.h:42
int nin() const
Definition: field.h:115
void mult_ct_boundary(int t, double ct)
Multiply the boundary improvement factor ct or ctr to an SU(N) matrix object which belongs to a site ...
Definition: field_G_SF.cpp:151
void set_parameters(const Parameters &params)
Mat_SU_N wk
SF boundary condition.
SU(N) gauge field.
Definition: field_G.h:38
void mult_Field_Gnd(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2)
void set_boundary_wk(const Mat_SU_N &U)
Set the boundary spatial link at t=0 for SF bc.
Definition: field_G_SF.cpp:50
void mult_Field_Gnn(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2)
double m_ct
SF boundary improvement coefficient for the plaquatte action.
void backward(Field &, const Field &, const int mu)
int nex() const
Definition: field.h:117
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:168
void at_Field_G(Field_G &w, const int ex)
void crucial(const char *format,...)
Definition: bridgeIO.cpp:48
static bool Register(const std::string &realm, const creator_callback &cb)
Bridge::VerboseLevel m_vl
Definition: action.h:81
Base class of random number generators.
Definition: randomNumbers.h:39
double m_ctr
SF boundary improvement coefficient for the rectangle action.
void Register_double_vector(const string &, const std::vector< double > &)
Definition: parameters.cpp:337
static const std::string class_name
void lower(Field_G_SF &, const Field_G &, const int, const int)
Definition: staples_SF.cpp:903
void stat(double &Fave, double &Fmax, double &Fdev) const
determines the statistics of the field. average, maximum value, and deviation is determined over glob...
Definition: field.cpp:516
Mat_SU_N mat_dag(const int site, const int mn=0) const
Definition: field_G.h:126
double langevin(RandomNumbers *)
Langevis step.
int fetch_double_vector(const string &key, std::vector< double > &val) const
Definition: parameters.cpp:158
static int reduce_sum(int count, double *recv_buf, double *send_buf, int pattern=0)
make a global sum of an array of double over the communicator. pattern specifies the dimensions to be...
void Register_double(const string &, const double)
Definition: parameters.cpp:323
void set(int c, double re, const double &im)
Definition: mat_SU_N.h:133
void setpart_ex(int ex, const Field &w, int exw)
Definition: field.h:177
int fetch_double(const string &key, double &val) const
Definition: parameters.cpp:124
string get_string(const string &key) const
Definition: parameters.cpp:87
Mat_SU_N mat(const int site, const int mn=0) const
Definition: field_G.h:113
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:28
double ReTr(const Mat_SU_N &m)
Definition: mat_SU_N.h:488
void forward(Field &, const Field &, const int mu)