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