Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
force_F_Clover_SF.cpp
Go to the documentation of this file.
1 
14 #include "force_F_Clover_SF.h"
15 
16 #ifdef USE_PARAMETERS_FACTORY
17 #include "parameters_factory.h"
18 #endif
19 
20 using std::valarray;
21 
22 //- parameter entries
23 namespace {
24  void append_entry(Parameters& param)
25  {
26  param.Register_double("hopping_parameter", 0.0);
27  param.Register_double("clover_coefficient", 0.0);
28  param.Register_int_vector("boundary_condition", std::valarray<int>());
29  param.Register_double_vector("phi", std::valarray<double>());
30  param.Register_double_vector("phipr", std::valarray<double>());
31 
32  param.Register_string("verbose_level", "NULL");
33  }
34 
35 
36 #ifdef USE_PARAMETERS_FACTORY
37  bool init_param = ParametersFactory::Register("Force.F_Clover_SF", append_entry);
38 #endif
39 }
40 //- end
41 
42 //- parameters class
43 Parameters_Force_F_Clover_SF::Parameters_Force_F_Clover_SF() { append_entry(*this); }
44 //- end
45 
46 //====================================================================
47 void Force_F_Clover_SF::set_parameters(const Parameters& params)
48 {
49  const string str_vlevel = params.get_string("verbose_level");
50 
51  m_vl = vout.set_verbose_level(str_vlevel);
52 
53  //- fetch and check input parameters
54  double kappa, cSW;
55  valarray<int> bc;
56  valarray<double> phi, phipr;
57 
58  int err = 0;
59  err += params.fetch_double("hopping_parameter", kappa);
60  err += params.fetch_double("clover_coefficient", cSW);
61  err += params.fetch_int_vector("boundary_condition", bc);
62  err += params.fetch_double_vector("phi", phi);
63  err += params.fetch_double_vector("phipr", phipr);
64 
65  if (err) {
66  vout.crucial(m_vl, "Force_F_Clover_SF: fetch error, input parameter not found.\n");
67  abort();
68  }
69 
70 
71  set_parameters(kappa, cSW, bc, &phi[0], &phipr[0]);
72 }
73 
74 
75 //====================================================================
76 void Force_F_Clover_SF::set_parameters(double kappa, double cSW, const valarray<int> bc,
77  double *phi, double *phipr)
78 {
79  int Ndim = CommonParameters::Ndim();
80 
81  //- print input parameters
82  vout.general(m_vl, "Parameters of Force_F_Clover_SF:\n");
83  vout.general(m_vl, " kappa = %8.4f\n", kappa);
84  vout.general(m_vl, " cSW = %8.4f\n", cSW);
85  for (int mu = 0; mu < Ndim; ++mu) {
86  vout.general(m_vl, " boundary[%d] = %2d\n", mu, bc[mu]);
87  }
88 
89  vout.general(m_vl, " phi1 = %12.6f\n", phi[0]);
90  vout.general(m_vl, " phi2 = %12.6f\n", phi[1]);
91  vout.general(m_vl, " phi3 = %12.6f\n", phi[2]);
92  vout.general(m_vl, " phipr1= %12.6f\n", phipr[0]);
93  vout.general(m_vl, " phipr2= %12.6f\n", phipr[1]);
94  vout.general(m_vl, " phipr3= %12.6f\n", phipr[2]);
95 
96  //- range check
97  // NB. kappa,cSW == 0 is allowed.
98  assert(bc.size() == Ndim);
99  // NB. phi,phipr == 0 is allowed.
100 
101  //- store values
102  m_kappa = kappa;
103  m_cSW = cSW;
104 
105  for (int i = 0; i < 3; ++i) {
106  m_phi[i] = phi[i];
107  m_phipr[i] = phipr[i];
108  }
109 
110  m_boundary.resize(Ndim);
111  for (int mu = 0; mu < Ndim; ++mu) {
112  m_boundary[mu] = bc[mu];
113  }
114 
115  //- propagate parameters
116  m_fopr_c->set_parameters(m_kappa, m_cSW, m_boundary, m_phi, m_phipr);
117 
118  m_force_w->set_parameters(m_kappa, m_boundary);
119 
120  set_wk.set_parameters(m_phi, m_phipr);
121 }
122 
123 
124 //====================================================================
125 Field Force_F_Clover_SF::force_core(const Field& eta)
126 {
127  int Nc = CommonParameters::Nc();
128  int Nvol = CommonParameters::Nvol();
129  int Ndim = CommonParameters::Ndim();
130  int NinG = 2 * Nc * Nc;
131 
132  Field_G force(Nvol, Ndim), force1(Nvol, Ndim);
133  Mat_SU_N ut(Nc);
134 
135  force1 = force_udiv(eta);
136 
137  for (int mu = 0; mu < Ndim; ++mu) {
138  force.mult_Field_Gnn(mu, *m_U, mu, force1, mu);
139  force.at_Field_G(mu);
140  }
141  force *= -2.0;
142 
143  return (Field)force;
144 }
145 
146 
147 //====================================================================
148 Field Force_F_Clover_SF::force_core1(const Field& zeta,
149  const Field& eta)
150 {
151  int Nc = CommonParameters::Nc();
152  int Nvol = CommonParameters::Nvol();
153  int Ndim = CommonParameters::Ndim();
154  int NinG = 2 * Nc * Nc;
155 
156  Field_G force(Nvol, Ndim), force1(Nvol, Ndim);
157  Mat_SU_N ut(Nc);
158 
159  force1 = force_udiv1(zeta, eta);
160 
161  for (int mu = 0; mu < Ndim; ++mu) {
162  force.mult_Field_Gnn(mu, *m_U, mu, force1, mu);
163  force.at_Field_G(mu);
164  }
165  force *= -2.0;
166 
167  return (Field)force;
168 }
169 
170 
171 //====================================================================
172 Field Force_F_Clover_SF::force_udiv(const Field& eta)
173 {
174  int Nc = CommonParameters::Nc();
175  int Nvol = CommonParameters::Nvol();
176  int Ndim = CommonParameters::Ndim();
177  int NinG = 2 * Nc * Nc;
178 
179  Field force(NinG, Nvol, Ndim);
180 
181  Field_F zeta(Nvol, 1);
182 
183  zeta = m_fopr_c->H(eta);
184 
185  force = force_udiv1((Field_F)eta, zeta);
186  force += force_udiv1(zeta, (Field_F)eta);
187 
188  set_wk.set_boundary_spatial_link_zero((Field_G&)force);
189 
190  return force;
191 }
192 
193 
194 //====================================================================
195 Field Force_F_Clover_SF::force_udiv1(const Field& zeta,
196  const Field& eta)
197 {
198  int Nc = CommonParameters::Nc();
199  int Nvol = CommonParameters::Nvol();
200  int Ndim = CommonParameters::Ndim();
201  int NinG = 2 * Nc * Nc;
202 
203  Field force(NinG, Nvol, Ndim);
204 
205  force = force_udiv1_impl((Field_F)zeta, (Field_F)eta);
206 
207  return force;
208 }
209 
210 
211 //====================================================================
212 Field Force_F_Clover_SF::force_udiv1_impl(const Field_F& zeta,
213  const Field_F& eta)
214 {
215  int Nc = CommonParameters::Nc();
216  int Nd = CommonParameters::Nd();
217  int Nvol = CommonParameters::Nvol();
218  int Ndim = CommonParameters::Ndim();
219  int NinG = 2 * Nc * Nc;
220 
221  ShiftField_lex shift;
222 
223  Field force(NinG, Nvol, Ndim);
224  Field_G force1(Nvol, 1), force2(Nvol, 1);
225  Field_G Umu(Nvol, 1), Unu(Nvol, 1), Utmp(Nvol, 1), Utmp2(Nvol, 1);
226  Field_F vt1(Nvol, 1), vt2(Nvol, 1), vt3(Nvol, 1), vt4(Nvol, 1);
227  Field_F zeta_mu(Nvol, 1);
228 
229  Mat_SU_N ut(Nc);
230  Vec_SU_N vec1(Nc), vec2(Nc);
231 
232  Field_F zeta1(zeta);
233 
234  set_zero.set_boundary_zero(zeta1);
235 
236  force = m_force_w->force_udiv1(zeta, eta);
237 
238  Field_F eta2(Nvol, 1), eta3(Nvol, 1);
239  eta2 = (Field_F)m_fopr_c->mult_gm5(eta);
240  set_zero.set_boundary_zero(eta2);
241 
242  for (int mu = 0; mu < Ndim; ++mu) {
243  for (int nu = 0; nu < Ndim; ++nu) {
244  if (nu == mu) continue;
245 
246  m_fopr_c->mult_isigma(eta3, eta2, mu, nu);
247 
248  Umu.setpart_ex(0, *m_U, mu);
249  Unu.setpart_ex(0, *m_U, nu);
250  if (mu != 3) set_wk.set_boundary_wk(Umu);
251  if (nu != 3) set_wk.set_boundary_wk(Unu);
252 
253  int ex = 0;
254  // R(1) and R(5)
255  vt1.mult_Field_Gd(0, *m_Cud, index_dir(mu, nu), eta3, ex);
256  tensorProd_Field_F(force1, zeta1, vt1);
257  force2 = force1;
258 
259  // R(2)
260  vt3.mult_Field_Gd(0, Umu, 0, eta3, ex);
261  shift.backward(vt1, vt3, nu);
262  shift.backward(vt2, zeta1, nu);
263  shift.backward(Utmp, Unu, mu);
264  if (mu == 3) set_wk.set_boundary_wkpr(Utmp);
265  vt3.mult_Field_Gn(0, Utmp, 0, vt1, ex);
266  vt4.mult_Field_Gn(0, Unu, 0, vt2, ex);
267  tensorProd_Field_F(force1, vt4, vt3);
268  force2 += force1;
269 
270  // R(4) and R(8)
271  shift.backward(vt1, eta3, mu);
272  shift.backward(zeta_mu, zeta1, mu);
273  vt4.mult_Field_Gn(0, *m_Cud, index_dir(mu, nu), zeta_mu, ex);
274  tensorProd_Field_F(force1, vt4, vt1);
275  force2 += force1;
276 
277  // R(3)
278  shift.backward(vt1, eta3, nu);
279  vt3.mult_Field_Gn(0, Unu, 0, vt1, ex);
280  vt4.mult_Field_Gn(0, Umu, 0, zeta_mu, ex);
281  shift.backward(vt1, vt3, mu);
282  shift.backward(vt2, vt4, nu);
283  vt4.mult_Field_Gn(0, Unu, 0, vt2, ex);
284  tensorProd_Field_F(force1, vt4, vt1);
285  force2 += force1;
286 
287  // R(6)
288  shift.backward(Utmp, Unu, mu);
289  if (mu == 3) set_wk.set_boundary_wkpr(Utmp);
290  Utmp2.mult_Field_Gdd(0, Utmp, 0, Umu, 0);
291  vt1.mult_Field_Gn(0, Utmp2, 0, eta3, ex);
292  vt2.mult_Field_Gd(0, Unu, 0, zeta1, ex);
293  shift.forward(vt3, vt1, nu);
294  shift.forward(vt4, vt2, nu);
295  tensorProd_Field_F(force1, vt4, vt3);
296  force2 -= force1;
297 
298  // R(7)
299  vt1.mult_Field_Gd(0, Unu, 0, eta3, ex);
300  vt2.mult_Field_Gn(0, Umu, 0, zeta_mu, ex);
301  shift.backward(vt3, vt1, mu);
302  shift.forward(vt1, vt3, nu);
303  vt4.mult_Field_Gd(0, Unu, 0, vt2, ex);
304  shift.forward(vt2, vt4, nu);
305  tensorProd_Field_F(force1, vt2, vt1);
306  force2 -= force1;
307 
308  force2 *= -m_kappa * m_cSW / 8.0;
309  force.addpart_ex(mu, force2, 0);
310  }
311  }
312 
313  return force;
314 }
315 
316 
317 //====================================================================
318 void Force_F_Clover_SF::set_component()
319 {
320  int Nc = CommonParameters::Nc();
321  int Nd = CommonParameters::Nd();
322  int Nvol = CommonParameters::Nvol();
323 
324  Staples_SF staple;
325 
326  staple.set_parameters(m_phi, m_phipr);
327 
328  Field_G Cmu_ud(Nvol, 1);
329 
330  for (int mu = 0; mu < m_Ndim; ++mu) {
331  for (int nu = 0; nu < m_Ndim; ++nu) {
332  if (nu == mu) continue;
333 
334  Cmu_ud = staple.upper(*m_U, mu, nu);
335  Cmu_ud -= staple.lower(*m_U, mu, nu);
336  m_Cud->setpart_ex(index_dir(mu, nu), Cmu_ud, 0);
337  }
338  }
339 }
340 
341 
342 //====================================================================
343 //============================================================END=====