Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
staples_SF.cpp
Go to the documentation of this file.
1 
14 #include "staples_SF.h"
15 
16 #ifdef USE_PARAMETERS_FACTORY
17 #include "parameters_factory.h"
18 #endif
19 
20 using std::valarray;
21 
38 //- parameter entries
39 namespace {
40  void append_entry(Parameters& param)
41  {
42  param.Register_double_vector("phi", std::valarray<double>());
43  param.Register_double_vector("phipr", std::valarray<double>());
44  param.Register_double_vector("pomega", std::valarray<double>());
45 
46  param.Register_string("verbose_level", "NULL");
47  }
48 
49 
50 #ifdef USE_PARAMETERS_FACTORY
51  bool init_param = ParametersFactory::Register("Staples_SF", append_entry);
52 #endif
53 }
54 //- end
55 
56 //- parameters class
57 Parameters_Staples_SF::Parameters_Staples_SF() { append_entry(*this); }
58 //- end
59 
60 //********************************************************************
61 void Staples_SF::set_parameters(const Parameters& params)
62 {
63  const string str_vlevel = params.get_string("verbose_level");
64 
65  m_vl = vout.set_verbose_level(str_vlevel);
66 
67  //- fetch and check input parameters
68  valarray<double> phi, phipr, pomega;
69 
70  int err = 0;
71  err += params.fetch_double_vector("phi", phi);
72  err += params.fetch_double_vector("phipr", phipr);
73  err += params.fetch_double_vector("pomega", pomega);
74 
75  if (err) {
76  vout.crucial(m_vl, "Staples_SF: fetch error, input parameter not found.\n");
77  abort();
78  }
79 
80 
81  set_parameters(phi, phipr, pomega); // call valarray version
82 }
83 
84 
85 //====================================================================
86 void Staples_SF::set_parameters(valarray<double>& phi, valarray<double>& phipr,
87  valarray<double>& pomega)
88 {
89  //- print input parameters
90  vout.general(m_vl, "Parameters of Staples_SF:\n");
91  vout.general(m_vl, " phi1 = %12.6f\n", phi[0]);
92  vout.general(m_vl, " phi2 = %12.6f\n", phi[1]);
93  vout.general(m_vl, " phi3 = %12.6f\n", phi[2]);
94  vout.general(m_vl, " phipr1 = %12.6f\n", phipr[0]);
95  vout.general(m_vl, " phipr2 = %12.6f\n", phipr[1]);
96  vout.general(m_vl, " phipr3 = %12.6f\n", phipr[2]);
97  vout.general(m_vl, " pomega1 = %12.6f\n", pomega[0]);
98  vout.general(m_vl, " pomega2 = %12.6f\n", pomega[1]);
99  vout.general(m_vl, " pomega3 = %12.6f\n", pomega[2]);
100 
101  //- range check
102  // NB. phi,phipr,pomega == 0 is allowed.
103  assert(phi.size() == 3);
104  assert(phipr.size() == 3);
105  assert(pomega.size() == 3);
106 
107  //- store values
108  set_parameters(&phi[0], &phipr[0], &pomega[0]); // call double[] version
109 }
110 
111 
112 //====================================================================
113 void Staples_SF::set_parameters(double *phi, double *phipr)
114 {
115  double c0r, c0i, c1r, c1i, c2r, c2i;
116 
117  c0r = cos(phi[0] / Lx);
118  c0i = sin(phi[0] / Lx);
119  c1r = cos(phi[1] / Lx);
120  c1i = sin(phi[1] / Lx);
121  c2r = cos(phi[2] / Lx);
122  c2i = sin(phi[2] / Lx);
123 
124  wk.zero();
125  wk.set(0, 0, c0r, c0i);
126  wk.set(1, 1, c1r, c1i);
127  wk.set(2, 2, c2r, c2i);
128 
129  c0r = cos(phipr[0] / Lx);
130  c0i = sin(phipr[0] / Lx);
131  c1r = cos(phipr[1] / Lx);
132  c1i = sin(phipr[1] / Lx);
133  c2r = cos(phipr[2] / Lx);
134  c2i = sin(phipr[2] / Lx);
135 
136  wkpr.zero();
137  wkpr.set(0, 0, c0r, c0i);
138  wkpr.set(1, 1, c1r, c1i);
139  wkpr.set(2, 2, c2r, c2i);
140 
141  iomega0.zero();
142  iomega0.set(0, 0, 0.0, 1.0);
143  iomega0.set(1, 1, 0.0, -0.5);
144  iomega0.set(2, 2, 0.0, -0.5);
145 
146  initialized = 1;
147 }
148 
149 
150 //====================================================================
151 
167 void Staples_SF::set_parameters(const double *phi, const double *phipr, const double *pomega)
168 {
169  double c0r, c0i, c1r, c1i, c2r, c2i;
170 
171  c0r = cos(phi[0] / Lx);
172  c0i = sin(phi[0] / Lx);
173  c1r = cos(phi[1] / Lx);
174  c1i = sin(phi[1] / Lx);
175  c2r = cos(phi[2] / Lx);
176  c2i = sin(phi[2] / Lx);
177  wk.zero();
178  wk.set(0, 0, c0r, c0i);
179  wk.set(1, 1, c1r, c1i);
180  wk.set(2, 2, c2r, c2i);
181 
182  c0r = cos(phipr[0] / Lx);
183  c0i = sin(phipr[0] / Lx);
184  c1r = cos(phipr[1] / Lx);
185  c1i = sin(phipr[1] / Lx);
186  c2r = cos(phipr[2] / Lx);
187  c2i = sin(phipr[2] / Lx);
188  wkpr.zero();
189  wkpr.set(0, 0, c0r, c0i);
190  wkpr.set(1, 1, c1r, c1i);
191  wkpr.set(2, 2, c2r, c2i);
192 
193  iomega0.zero();
194  iomega0.set(0, 0, 0.0, pomega[0]);
195  iomega0.set(1, 1, 0.0, pomega[1]);
196  iomega0.set(2, 2, 0.0, pomega[2]);
197 
198  initialized = 1;
199 }
200 
201 
202 //====================================================================
203 
247 double Staples_SF::sf_coupling_plaq(const Field_G& U, double ct)
248 {
249  if (!initialized) {
250  vout.general(m_vl, "Staples_SF: Paramer is not initialized.\n");
251  abort();
252  }
253 
254  double plaq = 0.0;
255  double plaqt0 = 0.0;
256  double plaqtT = 0.0;
257  double scr;
258  Mat_SU_N up(Nc);
259  static Field_G staple;
260 
261  int site;
262  for (int nu = 0; nu < Ndim - 1; nu++) {
263  staple = upper(U, 3, nu);
264 
265  for (int z = 0; z < Nz; z++) {
266  for (int y = 0; y < Ny; y++) {
267  for (int x = 0; x < Nx; x++) {
268  int t;
269  // boundary
270  if (Communicator::ipe(3) == 0) {
271  t = 0;
272  site = index.site(x, y, z, t);
273  up = staple.mat(site) * U.mat_dag(site, 3);
274  scr = ReTr(iomega0 * up);
275 
276  /*
277  up.unit();
278  up *= staple.mat(site);
279  up *= U->mat_dag(site,3);
280  up *= iomega0;
281  scr = ReTr( up );
282  */
283  plaq -= scr;
284  plaqt0 += scr;
285  }
286  // boundary
287  if (Communicator::ipe(3) == NPEt - 1) {
288  t = Nt - 1;
289  site = index.site(x, y, z, t);
290  up = staple.mat_dag(site) * U.mat(site, 3);
291  scr = ReTr(iomega0 * up);
292 
293  /*
294  up.unit();
295  up *= staple.mat_dag(site);
296  up *= U->mat(site,3);
297  up *= iomega0;
298  scr = ReTr( up );
299  */
300  plaq += scr;
301  plaqtT += scr;
302  }
303  }
304  }
305  }
306  }
307 
308  plaq = Communicator::reduce_sum(plaq);
309  plaqt0 = Communicator::reduce_sum(plaqt0);
310  plaqtT = Communicator::reduce_sum(plaqtT);
311  plaq *= ct / Lx;
312  plaqt0 *= ct / Lx;
313  plaqtT *= ct / Lx;
314 
315  vout.general(m_vl, "SF_delSg_plaq, from 0, from T = %.8f %.8f %.8f\n", plaq, plaqt0, plaqtT);
316 
317  return plaq;
318 }
319 
320 
321 //====================================================================
322 
421 double Staples_SF::sf_coupling_rect(const Field_G& m_U, double ctr)
422 {
423  if (!initialized) {
424  vout.general(m_vl, "Staples_SF: Paramer is not initialized.\n");
425  abort();
426  }
427 
428  double rect;
429  double rect01 = 0.0;
430  double rect02 = 0.0;
431  double rect03 = 0.0;
432  double rectt1 = 0.0;
433  double rectt2 = 0.0;
434  double rectt3 = 0.0;
435  double scr;
436 
437  Field_G Cup1(Nvol, 1), Cup2(Nvol, 1);
438  Field_G Cdn1(Nvol, 1), Cdn2(Nvol, 1);
439  Field_G Umu(Nvol, 1), Unu(Nvol, 1);
440  Field_G v(Nvol, 1), c(Nvol, 1);
441 
442  Mat_SU_N wmat(Nc), cmat(Nc);
443 
444  int site;
445  int nu = 3;
446  for (int mu = 0; mu < Ndim - 1; mu++) {
447  // rect01
448  // <---<---+
449  // | |
450  // t=0 x--->---+
451  // omega0
452 
453  // rect02
454  // <---<---+
455  // | |
456  // t=0 +---x---+
457  // omega0
458 
459  // rectt1
460  // omega0
461  // t=Nt x--->---+
462  // | |
463  // +---<---+
464 
465  // rectt2
466  // omega0
467  // t=Nt +---x---+
468  // | |
469  // <---<---+
470  Cup2 = upper(m_U, nu, mu);
471 
472  Umu.setpart_ex(0, m_U, mu);
473  Unu.setpart_ex(0, m_U, nu);
474  shift.backward(v, Cup2, mu);
475  shift.backward(c, Umu, nu);
476  for (int z = 0; z < Nz; z++) {
477  for (int y = 0; y < Ny; y++) {
478  for (int x = 0; x < Nx; x++) {
479  int t;
480  t = 0;
481  if (Communicator::ipe(3) == 0) {
482  site = index.site(x, y, z, t);
483  wmat = wk * v.mat(site);
484  cmat = wmat * c.mat_dag(site);
485  wmat = cmat * Unu.mat_dag(site);
486  rect01 += ReTr(iomega0 * wmat);
487  wmat = iomega0 * v.mat(site);
488  cmat = wmat * c.mat_dag(site);
489  wmat = cmat * Unu.mat_dag(site);
490  rect02 += ReTr(wk * wmat);
491  }
492  t = Nt - 1;
493  if (Communicator::ipe(3) == NPEt - 1) {
494  site = index.site(x, y, z, t);
495  cmat = iomega0 * wkpr;
496  wmat = cmat * v.mat_dag(site);
497  cmat = Unu.mat(site) * wmat;
498  rectt1 += ReTr(cmat * m_U.mat_dag(site, mu));
499  cmat = iomega0 * v.mat_dag(site);
500  wmat = wkpr * cmat;
501  cmat = Unu.mat(site) * wmat;
502  rectt2 += ReTr(cmat * m_U.mat_dag(site, mu));
503  }
504  }
505  }
506  }
507 
508  // rect03
509  // +---+
510  // | |
511  // v ^
512  // | |
513  // t=0 x--->
514  // omega0
515  Cup1 = upper(m_U, mu, nu);
516 
517  shift.backward(v, Unu, mu);
518  shift.backward(c, Cup1, nu);
519  for (int z = 0; z < Nz; z++) {
520  for (int y = 0; y < Ny; y++) {
521  for (int x = 0; x < Nx; x++) {
522  int t;
523  t = 0;
524  if (Communicator::ipe(3) == 0) {
525  site = index.site(x, y, z, t);
526  wmat = c.mat(site) * v.mat_dag(site);
527  cmat = Unu.mat(site) * wmat;
528  wmat = wk * cmat.dag();
529  rect03 += ReTr(iomega0 * wmat);
530  }
531  }
532  }
533  }
534 
535  // rectt3
536  // omega0
537  // t=Nt x--->
538  // | |
539  // ^ v
540  // | |
541  // +---+
542  Cdn1 = lower(m_U, mu, nu);
543 
544  shift.backward(v, Unu, mu);
545  for (int z = 0; z < Nz; z++) {
546  for (int y = 0; y < Ny; y++) {
547  for (int x = 0; x < Nx; x++) {
548  int t;
549  t = Nt - 1;
550  site = index.site(x, y, z, t);
551  if (Communicator::ipe(3) == NPEt - 1) {
552  wmat = iomega0 * wkpr;
553  cmat = wmat * v.mat_dag(site);
554  wmat = cmat * Cdn1.mat_dag(site);
555  rectt3 += ReTr(Unu.mat(site) * wmat);
556  }
557  }
558  }
559  }
560  }
561 
562  rect01 = Communicator::reduce_sum(rect01);
563  rect02 = Communicator::reduce_sum(rect02);
564  rect03 = Communicator::reduce_sum(rect03);
565  rectt1 = Communicator::reduce_sum(rectt1);
566  rectt2 = Communicator::reduce_sum(rectt2);
567  rectt3 = Communicator::reduce_sum(rectt3);
568  rect01 *= ctr / Lx;
569  rect02 *= ctr / Lx;
570  rect03 /= Lx;
571  rectt1 *= ctr / Lx;
572  rectt2 *= ctr / Lx;
573  rectt3 /= Lx;
574  rect = -rect01 - rect02 - rect03 + rectt1 + rectt2 + rectt3;
575 
576  vout.general(m_vl, "SF_delSg_rect, at 01, 02, 03, at T1, T2, T3 = %.8f %.8f %.8f %.8f %.8f %.8f %.8f\n",
577  rect, rect01, rect02, rect03, rectt1, rectt2, rectt3);
578 
579  return rect;
580 }
581 
582 
583 //====================================================================
584 
592 double Staples_SF::plaquette(const Field_G& U)
593 {
594  if (!initialized) {
595  vout.general(m_vl, "Staples_SF: Paramer is not initialized.\n");
596  abort();
597  }
598  // return (plaq_s(U) +plaq_t(U))/2;
599  return(plaq_s(U) + plaq_t(U));
600 }
601 
602 
603 //====================================================================
604 
618 double Staples_SF::plaquette_ct(const Field_G& U, double ct)
619 {
620  if (!initialized) {
621  vout.general(m_vl, "Staples_SF: Paramer is not initialized.\n");
622  abort();
623  }
624  return(plaq_s(U) + plaq_t_ct(U, ct));
625 }
626 
627 
628 //====================================================================
629 
638 double Staples_SF::plaq_s(const Field_G& U)
639 {
640  double plaq = 0.0;
641  static Field_G staple;
642 
643  staple = upper(U, 0, 1);
644  // staple = lower(U,0,1);
645  for (int site = 0; site < Nvol; site++) {
646  plaq += ReTr(U.mat(site, 0) * staple.mat_dag(site)); // P_xy
647  }
648 
649  staple = upper(U, 1, 2);
650  // staple = lower(U,1,2);
651  for (int site = 0; site < Nvol; site++) {
652  plaq += ReTr(U.mat(site, 1) * staple.mat_dag(site)); // P_yz
653  }
654 
655  staple = upper(U, 2, 0);
656  // staple = lower(U,2,0);
657  for (int site = 0; site < Nvol; site++) {
658  plaq += ReTr(U.mat(site, 2) * staple.mat_dag(site)); // P_zx
659  }
660 
661  plaq = Communicator::reduce_sum(plaq);
662 
663  // return plaq/(Lvol*Nc*3.0);
664  return plaq;
665 }
666 
667 
668 //====================================================================
669 
678 double Staples_SF::plaq_t(const Field_G& U)
679 {
680  if (!initialized) {
681  vout.general(m_vl, "Staples_SF: Paramer is not initialized.\n");
682  abort();
683  }
684 
685  double plaq = 0.0;
686  static Field_G staple;
687 
688  for (int nu = 0; nu < Ndim - 1; nu++) {
689  staple = lower(U, 3, nu);
690  // staple = upper(U,3,nu);
691  for (int site = 0; site < Nvol; site++) {
692  plaq += ReTr(U.mat(site, 3) * staple.mat_dag(site)); // P_tk
693  }
694  }
695 
696  plaq = Communicator::reduce_sum(plaq);
697 
698  // return plaq/(Lvol*Nc*3.0);
699  return plaq;
700 }
701 
702 
703 //====================================================================
704 
719 double Staples_SF::plaq_t_ct(const Field_G& U, double ct)
720 {
721  if (!initialized) {
722  vout.general(m_vl, "Staples_SF: Paramer is not initialized.\n");
723  abort();
724  }
725 
726  double plaq = 0.0;
727  static Field_G_SF staple;
728 
729  for (int nu = 0; nu < Ndim - 1; nu++) {
730  staple = lower(U, 3, nu);
731  // If the node is at the boundary the temporal plaquette is multiplied with ct.
732  if (Communicator::ipe(3) == 0) {
733  staple.mult_ct_boundary(0, ct);
734  }
735  if (Communicator::ipe(3) == NPEt - 1) {
736  staple.mult_ct_boundary(Nt - 1, ct);
737  }
738  for (int site = 0; site < Nvol; site++) {
739  plaq += ReTr(U.mat(site, 3) * staple.mat_dag(site)); // P_tk
740  }
741  }
742 
743  plaq = Communicator::reduce_sum(plaq);
744 
745  // return plaq/(Lvol*Nc*3.0);
746  return plaq;
747 }
748 
749 
750 //====================================================================
751 
770 void Staples_SF::staple(Field_G& W, const Field_G& U, const int mu)
771 {
772  if (!initialized) {
773  vout.general(m_vl, "Staples_SF: Paramer is not initialized.\n");
774  abort();
775  }
776 
777  W = 0.0;
778  for (int nu = 0; nu < Ndim; nu++) {
779  if (nu != mu) {
780  W += upper(U, mu, nu);
781  W += lower(U, mu, nu);
782  }
783  }
784 }
785 
786 
787 //====================================================================
788 
807 void Staples_SF::staple_ct(Field_G& W, const Field_G& U, const int mu, double ct)
808 {
809  W = 0.0;
810  Field_G_SF staple_upper;
811  Field_G_SF staple_lower;
812 
813  for (int nu = 0; nu < Ndim; nu++) {
814  if (nu != mu) {
815  staple_upper = upper(U, mu, nu);
816  staple_lower = lower(U, mu, nu);
817  if (Communicator::ipe(3) == 0) {
818  if (mu == 3) {
819  staple_upper.mult_ct_boundary(0, ct);
820  staple_lower.mult_ct_boundary(0, ct);
821  }
822  if (nu == 3) {
823  staple_lower.mult_ct_boundary(1, ct);
824  }
825  }
826  if (Communicator::ipe(3) == NPEt - 1) {
827  if (mu == 3) {
828  staple_upper.mult_ct_boundary(Nt - 1, ct);
829  staple_lower.mult_ct_boundary(Nt - 1, ct);
830  }
831  if (nu == 3) {
832  staple_upper.mult_ct_boundary(Nt - 1, ct);
833  }
834  }
835  W += (Field_G)staple_upper;
836  W += (Field_G)staple_lower;
837  }
838  }
839 }
840 
841 
842 //====================================================================
843 
862 Field_G_SF Staples_SF::upper(const Field_G& U, const int mu, const int nu)
863 {
864  if (!initialized) {
865  vout.general(m_vl, "Staples_SF: Paramer is not initialized.\n");
866  abort();
867  }
868  // (1) mu (2)
869  // +-->--+
870  // nu | |
871  // i+ +
872 
873  Field_G_SF c;
874  Umu.setpart_ex(0, U, mu);
875  Unu.setpart_ex(0, U, nu);
876  if (mu != 3) Umu.set_boundary_wk(wk);
877  if (nu != 3) Unu.set_boundary_wk(wk);
878 
879  shift.backward(v, Unu, mu);
880  shift.backward(c, Umu, nu);
881  if (mu == 3) v.set_boundary_wkpr(wkpr);
882  if (nu == 3) c.set_boundary_wkpr(wkpr);
883 
884  w.mult_Field_Gnd(0, c, 0, v, 0);
885  c.mult_Field_Gnn(0, Unu, 0, w, 0);
886  if (mu != 3) c.set_boundary_zero();
887 
888  return c;
889 }
890 
891 
892 //====================================================================
893 
908 Field_G_SF Staples_SF::lower(const Field_G& U, const int mu, const int nu)
909 {
910  if (!initialized) {
911  vout.general(m_vl, "Staples_SF: Paramer is not initialized.\n");
912  abort();
913  }
914  // + +
915  // nu | |
916  // i+-->--+
917  // (1) mu (2)
918 
919  Field_G_SF c;
920  Umu.setpart_ex(0, U, mu);
921  Unu.setpart_ex(0, U, nu);
922  if (mu != 3) Umu.set_boundary_wk(wk);
923  if (nu != 3) Unu.set_boundary_wk(wk);
924 
925  shift.backward(w, Unu, mu);
926  if (mu == 3) w.set_boundary_wkpr(wkpr);
927 
928  v.mult_Field_Gnn(0, Umu, 0, w, 0);
929  w.mult_Field_Gdn(0, Unu, 0, v, 0);
930 
931  shift.forward(c, w, nu);
932  if (mu != 3) c.set_boundary_zero();
933 
934  return c;
935 }
936 
937 
938 //====================================================================
939 
957 void Staples_SF::print_plaquette(const Field_G& U)
958 {
959  double plaq = plaquette(U);
960  double plaq2 = plaq + 3 * 3 * Lx * Ly * Lz;
961 
962  vout.general(m_vl, "plaq_SF without boundary spatial plaq = %.8f\n",
963  plaq / (3 * Lx * Ly * Lz * (6 * Lt - 3)));
964  vout.general(m_vl, "plaq_SF with boundary spatial plaq = %.8f\n",
965  plaq2 / (3 * 6 * Lx * Ly * Lz * Lt));
966 }
967 
968 
969 //============================================================END=====