Bridge++  Ver. 2.0.2
force_G_Rectangle_SF.cpp
Go to the documentation of this file.
1 
14 #include "force_G_Rectangle_SF.h"
15 
16 #ifdef USE_FACTORY_AUTOREGISTER
17 namespace {
18  bool init = Force_G_Rectangle_SF::register_factory();
19 }
20 #endif
21 
22 const std::string Force_G_Rectangle_SF::class_name = "Force_G_Rectangle_SF";
23 
24 //====================================================================
25 void Force_G_Rectangle_SF::set_parameters(const Parameters& params)
26 {
27  std::string vlevel;
28  if (!params.fetch_string("verbose_level", vlevel)) {
29  m_vl = vout.set_verbose_level(vlevel);
30  }
31 
32  //- fetch and check input parameters
33  double beta, c_plaq, c_rect;
34  double ct0, ct1, ct2, ctr0, ctr1, ctr2;
35  std::vector<double> phi, phipr;
36 
37  int err = 0;
38  err += params.fetch_double("beta", beta);
39  err += params.fetch_double("c_plaq", c_plaq);
40  err += params.fetch_double("c_rect", c_rect);
41 
42  err += params.fetch_double("ct0", ct0);
43  err += params.fetch_double("ct1", ct1);
44  err += params.fetch_double("ct2", ct2);
45 
46  err += params.fetch_double("ctr0", ctr0);
47  err += params.fetch_double("ctr1", ctr1);
48  err += params.fetch_double("ctr2", ctr2);
49 
50  err += params.fetch_double_vector("phi", phi);
51  err += params.fetch_double_vector("phipr", phipr);
52 
53  if (err) {
54  vout.crucial(m_vl, "Error at %s: input parameter not found.\n", class_name.c_str());
55  exit(EXIT_FAILURE);
56  }
57 
58 
59  const double gg = 6.0 / beta;
60 
61  const double ct = ct0 + ct1 * gg + ct2 * gg * gg;
62  const double ctr = ctr0 + ctr1 * gg + ctr2 * gg * gg;
63 
64  set_parameters(beta, c_plaq, c_rect, &phi[0], &phipr[0], ct, ctr);
65 }
66 
67 
68 //====================================================================
69 void Force_G_Rectangle_SF::get_parameters(Parameters& params) const
70 {
71  params.set_double("beta", m_beta);
72  params.set_double("c_plaq", m_c_plaq);
73  params.set_double("c_rect", m_c_rect);
74 
75  // params.set_double("ct0", m_ct0);
76  // params.set_double("ct1", m_ct1);
77  // params.set_double("ct2", m_ct2);
78 
79  // params.set_double("ctr0", m_ctr0);
80  // params.set_double("ctr1", m_ctr1);
81  // params.set_double("ctr2", m_ctr2);
82 
83  params.set_double_vector("phi", m_phi);
84  params.set_double_vector("phipr", m_phipr);
85 
86  params.set_double("ct", m_ct);
87  params.set_double("ctr", m_ctr);
88 
89  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
90 }
91 
92 
93 //====================================================================
94 
108 void Force_G_Rectangle_SF::set_parameters(const double beta, const double c_plaq, const double c_rect,
109  double *phi, double *phipr, const double ct, const double ctr)
110 {
111  //- print input parameters
112  vout.general(m_vl, "%s:\n", class_name.c_str());
113  vout.general(m_vl, " beta = %12.6f\n", beta);
114  vout.general(m_vl, " c_plaq = %12.6f\n", c_plaq);
115  vout.general(m_vl, " c_rect = %12.6f\n", c_rect);
116  vout.general(m_vl, " phi1 = %12.6f\n", phi[0]);
117  vout.general(m_vl, " phi2 = %12.6f\n", phi[1]);
118  vout.general(m_vl, " phi3 = %12.6f\n", phi[2]);
119  vout.general(m_vl, " phipr1 = %12.6f\n", phipr[0]);
120  vout.general(m_vl, " phipr2 = %12.6f\n", phipr[1]);
121  vout.general(m_vl, " phipr3 = %12.6f\n", phipr[2]);
122  vout.general(m_vl, " ct = %12.6f\n", ct);
123  vout.general(m_vl, " ctr = %12.6f\n", ctr);
124 
125  //- range check
126  // NB. beta,c_plaq,c_rect,phi,phipr,ct,ctr = 0 is allowed.
127 
128  //- store values
129  m_beta = beta;
130  m_c_plaq = c_plaq;
131  m_c_rect = c_rect;
132 
133  m_ct = ct;
134  m_ctr = ctr;
135  // m_phi = phi ;
136  // m_phipr = phipr;
137 
138  m_phi.resize(3);
139  m_phipr.resize(3);
140  for (int i = 0; i < 3; ++i) {
141  m_phi[i] = phi[i];
142  m_phipr[i] = phipr[i];
143  }
144 
145  //- post-process
146  m_staple.set_parameters(m_phi, m_phipr);
147 
148  const int Lx = CommonParameters::Lx();
149 
150  double c0r = cos(phi[0] / Lx);
151  double c0i = sin(phi[0] / Lx);
152  double c1r = cos(phi[1] / Lx);
153  double c1i = sin(phi[1] / Lx);
154  double c2r = cos(phi[2] / Lx);
155  double c2i = sin(phi[2] / Lx);
156  m_wk.zero();
157  m_wk.set(0, 0, c0r, c0i);
158  m_wk.set(1, 1, c1r, c1i);
159  m_wk.set(2, 2, c2r, c2i);
160 
161  c0r = cos(phipr[0] / Lx);
162  c0i = sin(phipr[0] / Lx);
163  c1r = cos(phipr[1] / Lx);
164  c1i = sin(phipr[1] / Lx);
165  c2r = cos(phipr[2] / Lx);
166  c2i = sin(phipr[2] / Lx);
167  m_wkpr.zero();
168  m_wkpr.set(0, 0, c0r, c0i);
169  m_wkpr.set(1, 1, c1r, c1i);
170  m_wkpr.set(2, 2, c2r, c2i);
171 }
172 
173 
174 //====================================================================
175 
230 void Force_G_Rectangle_SF::force_core(Field& force)
231 {
232  const int Nvol = CommonParameters::Nvol();
233  const int Ndim = CommonParameters::Ndim();
234 
235  const int Nt = CommonParameters::Nt();
236  const int NPEt = CommonParameters::NPEt();
237 
238  const int IPE3 = Communicator::ipe(3);
239 
240  assert(m_U->nin() == m_Nc * m_Nc * 2);
241  assert(m_U->nvol() == Nvol);
242  assert(m_U->nex() == Ndim);
243 
244  assert(force.nin() == m_Nc * m_Nc * 2);
245  assert(force.nvol() == Nvol);
246  assert(force.nex() == Ndim);
247 
248  Mat_SU_N wzero(m_Nc);
249  wzero.zero();
250 
251  for (int mu = 0; mu < Ndim; ++mu) {
252  Field_G force1;
253  force1.set(0.0);
254 
255  for (int nu = 0; nu < Ndim; ++nu) {
256  if (nu == mu) continue;
257 
258  Field_G Cup1;
259  m_staple.upper(Cup1, *m_U, mu, nu);
260 
261  Field_G Cup2;
262  m_staple.upper(Cup2, *m_U, nu, mu);
263 
264  Field_G Cdn1;
265  m_staple.lower(Cdn1, *m_U, mu, nu);
266 
267  Field_G Cdn2;
268  m_staple.lower(Cdn2, *m_U, nu, mu);
269 
270  // plaquette term
271  Field_G Umu = Cup1;
272  Field_G Unu = Cdn1;
273  if (IPE3 == 0) {
274  if (mu == 3) {
275  Field_SF::mult_ct_boundary(Umu, 0, m_ct);
276  Field_SF::mult_ct_boundary(Unu, 0, m_ct);
277  }
278  if (nu == 3) {
279  Field_SF::mult_ct_boundary(Unu, 1, m_ct);
280  }
281  }
282  if (IPE3 == NPEt - 1) {
283  if (mu == 3) {
284  Field_SF::mult_ct_boundary(Umu, Nt - 1, m_ct);
285  Field_SF::mult_ct_boundary(Unu, Nt - 1, m_ct);
286  }
287  if (nu == 3) {
288  Field_SF::mult_ct_boundary(Umu, Nt - 1, m_ct);
289  }
290  }
291  axpy(force1, 0, m_c_plaq, Umu, 0);
292  axpy(force1, 0, m_c_plaq, Unu, 0);
293 
294  // rectangular term
295  copy(Umu, 0, *m_U, mu);
296  copy(Unu, 0, *m_U, nu);
297  // For the boundary spatial link, use Wk.
298  if ((IPE3 == 0) && (nu == 3)) Field_SF::set_boundary_wk(Umu, m_wk);
299  if ((IPE3 == 0) && (mu == 3)) Field_SF::set_boundary_wk(Unu, m_wk);
300 
301  // +---+---+
302  // | | term
303  // x <---+
304  Field_G v;
305  m_shift.backward(v, Cup2, mu);
306 
307  Field_G c;
308  m_shift.backward(c, Umu, nu);
309  // The force for the spatial link near the boundary. Multiplied with ctr.
310  if ((IPE3 == NPEt - 1) && (nu == 3)) {
311  Field_SF::set_boundary_wkpr(c, m_wkpr);
312  Field_SF::mult_ct_boundary(c, Nt - 1, m_ctr);
313  }
314 
315  Field_G w;
316  mult_Field_Gnd(w, 0, c, 0, v, 0);
317  mult_Field_Gnn(c, 0, Unu, 0, w, 0);
318  // The force for the boundary spatial link is set to zero.
319  if ((IPE3 == 0) && (nu == 3)) Field_SF::set_boundary_zero(c);
320  axpy(force1, 0, m_c_rect, c, 0);
321 
322  // +---+
323  // | |
324  // + + term
325  // | |
326  // x v
327 
328  m_shift.backward(v, Unu, mu);
329  m_shift.backward(c, Cup1, nu);
330  // The force for the temporal link near the boundary. Multiplied with ctr.
331  if ((IPE3 == NPEt - 1) && (mu == 3)) {
332  Field_SF::set_boundary_wkpr(v, m_wkpr);
333  Field_SF::mult_ct_boundary(v, Nt - 1, m_ctr);
334  }
335  mult_Field_Gnd(w, 0, c, 0, v, 0);
336  mult_Field_Gnn(c, 0, Unu, 0, w, 0);
337  // The force for the boundary spatial link is set to zero.
338  if ((IPE3 == 0) && (nu == 3)) Field_SF::set_boundary_zero(c);
339  // The force for the boundary temporal link.
340  if ((IPE3 == 0) && (mu == 3)) Field_SF::mult_ct_boundary(c, 0, m_ctr);
341  axpy(force1, 0, m_c_rect, c, 0);
342 
343  // +---+---+
344  // | | term
345  // +---x v
346 
347  m_shift.backward(v, Unu, mu);
348  m_shift.backward(c, Umu, nu);
349  if ((IPE3 == NPEt - 1) && (nu == 3)) {
350  Field_SF::set_boundary_wkpr(c, m_wkpr);
351  Field_SF::mult_ct_boundary(c, Nt - 1, m_ctr);
352  }
353  if ((IPE3 == NPEt - 1) && (mu == 3)) Field_SF::set_boundary_wkpr(v, m_wkpr);
354  mult_Field_Gnd(w, 0, c, 0, v, 0);
355  mult_Field_Gnn(c, 0, Cdn2, 0, w, 0);
356  if ((IPE3 == 0) && (nu == 3)) Field_SF::set_boundary_zero(c);
357  axpy(force1, 0, m_c_rect, c, 0);
358 
359  // x <---+
360  // | | term
361  // +---+---+
362 
363  m_shift.backward(v, Cup2, mu);
364  mult_Field_Gnn(w, 0, Umu, 0, v, 0);
365  mult_Field_Gdn(v, 0, Unu, 0, w, 0);
366  if ((IPE3 == 0) && (nu == 3)) Field_SF::mult_ct_boundary(v, 0, m_ctr);
367  m_shift.forward(c, v, nu);
368  if ((IPE3 == 0) && (nu == 3)) Field_SF::set_boundary_zero(c);
369  axpy(force1, 0, m_c_rect, c, 0);
370 
371  // x ^
372  // | |
373  // + + term
374  // | |
375  // +---+
376 
377  m_shift.backward(v, Unu, mu);
378  if ((IPE3 == NPEt - 1) && (mu == 3)) {
379  Field_SF::set_boundary_wkpr(v, m_wkpr);
380  Field_SF::mult_ct_boundary(v, Nt - 1, m_ctr);
381  }
382  mult_Field_Gnn(w, 0, Cdn1, 0, v, 0);
383  mult_Field_Gdn(v, 0, Unu, 0, w, 0);
384  if ((IPE3 == 0) && (mu == 3)) Field_SF::mult_ct_boundary(v, 0, m_ctr);
385  m_shift.forward(c, v, nu);
386  if ((IPE3 == 0) && (nu == 3)) Field_SF::set_boundary_zero(c);
387  axpy(force1, 0, m_c_rect, c, 0);
388 
389  // +---x ^
390  // | | term
391  // +---+---+
392 
393  m_shift.backward(v, Unu, mu);
394  if ((IPE3 == NPEt - 1) && (mu == 3)) Field_SF::set_boundary_wkpr(v, m_wkpr);
395  mult_Field_Gnn(w, 0, Umu, 0, v, 0);
396  mult_Field_Gdn(v, 0, Cdn2, 0, w, 0);
397  if ((IPE3 == 0) && (nu == 3)) Field_SF::mult_ct_boundary(v, 0, m_ctr);
398  m_shift.forward(c, v, nu);
399  if ((IPE3 == 0) && (nu == 3)) Field_SF::set_boundary_zero(c);
400  axpy(force1, 0, m_c_rect, c, 0);
401  }
402 
403  Field_G force2;
404  mult_Field_Gnd(force2, 0, *m_U, mu, force1, 0);
405  at_Field_G(force2, 0);
406 
407  axpy(force, mu, -(m_beta / m_Nc), force2, 0);
408  }
409 
410  double Fave, Fmax, Fdev;
411  force.stat(Fave, Fmax, Fdev);
412  vout.general(m_vl, " Fave = %12.6f Fmax = %12.6f Fdev = %12.6f\n",
413  Fave, Fmax, Fdev);
414 }
415 
416 
417 //============================================================END=====
Field_SF::set_boundary_wkpr
void set_boundary_wkpr(Field_G &u, const Mat_SU_N &wkpr)
Definition: field_SF.cpp:63
Force_G_Rectangle_SF::get_parameters
void get_parameters(Parameters &params) const
Definition: force_G_Rectangle_SF.cpp:69
mult_Field_Gdn
void mult_Field_Gdn(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)
Definition: field_G_imp.cpp:134
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
Staple_SF::set_parameters
void set_parameters(const Parameters &params)
Definition: staple_SF.cpp:35
ShiftField_lex::forward
void forward(Field &, const Field &, const int mu)
Definition: shiftField_lex.cpp:79
Force_G_Rectangle_SF::m_phipr
std::vector< double > m_phipr
Definition: force_G_Rectangle_SF.h:135
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
Force_G_Rectangle_SF::m_c_plaq
double m_c_plaq
Definition: force_G_Rectangle_SF.h:117
Force_G_Rectangle_SF::m_vl
Bridge::VerboseLevel m_vl
Definition: force_G_Rectangle_SF.h:109
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
force_G_Rectangle_SF.h
Parameters
Class for parameters.
Definition: parameters.h:46
Force_G_Rectangle_SF::m_shift
ShiftField_lex m_shift
Definition: force_G_Rectangle_SF.h:123
Force_G_Rectangle_SF::m_wkpr
Mat_SU_N m_wkpr
Definition: force_G_Rectangle_SF.h:126
Parameters::set_double
void set_double(const string &key, const double value)
Definition: parameters.cpp:33
Force_G_Rectangle_SF::m_beta
double m_beta
Definition: force_G_Rectangle_SF.h:116
Field::nex
int nex() const
Definition: field.h:128
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
Field::stat
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:169
Field_SF::set_boundary_wk
void set_boundary_wk(Field_G &u, const Mat_SU_N &wk)
Definition: field_SF.cpp:32
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:380
Parameters::set_double_vector
void set_double_vector(const string &key, const vector< double > &value)
Definition: parameters.cpp:42
Field::nin
int nin() const
Definition: field.h:126
Staple_SF::upper
void upper(Field_G &, const Field_G &, const int, const int)
Definition: staple_SF.cpp:877
at_Field_G
void at_Field_G(Field_G &W, const int ex)
Definition: field_G_imp.cpp:419
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:212
SU_N::Mat_SU_N::set
void set(int c, const double &re, const double &im)
Definition: mat_SU_N.h:137
Field_SF::mult_ct_boundary
void mult_ct_boundary(Field_G &u, const int t, const double ct)
Definition: field_SF.cpp:185
Field_SF::set_boundary_zero
void set_boundary_zero(Field_G &u)
Definition: field_SF.cpp:96
CommonParameters::Lx
static int Lx()
Definition: commonParameters.h:91
Force_G_Rectangle_SF::m_ctr
double m_ctr
SF boundary improvement coefficient for the rectangle action.
Definition: force_G_Rectangle_SF.h:132
SU_N::Mat_SU_N::zero
Mat_SU_N & zero()
Definition: mat_SU_N.h:429
CommonParameters::Nt
static int Nt()
Definition: commonParameters.h:108
Force_G_Rectangle_SF::m_Nc
int m_Nc
Definition: force_G_Rectangle_SF.h:114
SU_N::Mat_SU_N
Definition: mat_SU_N.h:36
Field::nvol
int nvol() const
Definition: field.h:127
Force_G_Rectangle_SF::m_ct
double m_ct
SF boundary improvement coefficient for the plaquatte action.
Definition: force_G_Rectangle_SF.h:130
Force_G_Rectangle_SF::m_staple
Staple_SF m_staple
Definition: force_G_Rectangle_SF.h:122
Staple_SF::lower
void lower(Field_G &, const Field_G &, const int, const int)
Definition: staple_SF.cpp:912
Force_G_Rectangle_SF::set_parameters
void set_parameters(const Parameters &params)
Definition: force_G_Rectangle_SF.cpp:25
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
mult_Field_Gnn
void mult_Field_Gnn(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)
Definition: field_G_imp.cpp:95
Force_G::m_U
Field_G * m_U
Definition: force_G.h:34
ShiftField_lex::backward
void backward(Field &, const Field &, const int mu)
Definition: shiftField_lex.cpp:59
Force_G_Rectangle_SF::class_name
static const std::string class_name
Definition: force_G_Rectangle_SF.h:106
Communicator::ipe
static int ipe(const int dir)
logical coordinate of current proc.
Definition: communicator.cpp:105
Parameters::fetch_string
int fetch_string(const string &key, string &value) const
Definition: parameters.cpp:378
Parameters::fetch_double
int fetch_double(const string &key, double &value) const
Definition: parameters.cpp:327
Force_G_Rectangle_SF::m_c_rect
double m_c_rect
Definition: force_G_Rectangle_SF.h:118
CommonParameters::NPEt
static int NPEt()
Definition: commonParameters.h:100
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:180
Field
Container of Field-type object.
Definition: field.h:46
Force_G_Rectangle_SF::m_wk
Mat_SU_N m_wk
SF boundary condition.
Definition: force_G_Rectangle_SF.h:126
Field_G
SU(N) gauge field.
Definition: field_G.h:38
Parameters::fetch_double_vector
int fetch_double_vector(const string &key, vector< double > &value) const
Definition: parameters.cpp:410
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
Force_G_Rectangle_SF::force_core
void force_core(Field &)
Definition: force_G_Rectangle_SF.cpp:230
mult_Field_Gnd
void mult_Field_Gnd(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)
Definition: field_G_imp.cpp:173
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:512
Force_G_Rectangle_SF::m_phi
std::vector< double > m_phi
Definition: force_G_Rectangle_SF.h:135
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:154