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