Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
force_F_Overlap_Nf2.cpp
Go to the documentation of this file.
1 
14 #include "force_F_Overlap_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("quark_mass", 0.0);
27  param.Register_double("domain_wall_height", 0.0);
28  param.Register_int("number_of_poles", 0);
29  param.Register_double("lower_bound", 0.0);
30  param.Register_double("upper_bound", 0.0);
31  param.Register_int("maximum_number_of_iteration", 0);
32  param.Register_double("convergence_criterion_squared", 0.0);
33  param.Register_int_vector("boundary_condition", std::valarray<int>());
34 
35  param.Register_string("verbose_level", "NULL");
36  }
37 
38 
39 #ifdef USE_PARAMETERS_FACTORY
40  bool init_param = ParametersFactory::Register("Force.F_Overlap_Nf2", append_entry);
41 #endif
42 }
43 //- end
44 
45 //- parameters class
46 Parameters_Force_F_Overlap_Nf2::Parameters_Force_F_Overlap_Nf2() { append_entry(*this); }
47 //- end
48 
49 //====================================================================
50 void Force_F_Overlap_Nf2::set_parameters(const Parameters& params)
51 {
52  const string str_vlevel = params.get_string("verbose_level");
53 
54  m_vl = vout.set_verbose_level(str_vlevel);
55 
56  //- fetch and check input parameters
57  double mq, M0;
58  int Np;
59  double x_min, x_max;
60  int Niter_ms;
61  double Stop_cond_ms;
62  valarray<int> bc;
63 
64  int err = 0;
65  err += params.fetch_double("quark_mass", mq);
66  err += params.fetch_double("domain_wall_height", M0);
67  err += params.fetch_int("number_of_poles", Np);
68  err += params.fetch_double("lower_bound", x_min);
69  err += params.fetch_double("upper_bound", x_max);
70  err += params.fetch_int("maximum_number_of_iteration", Niter_ms);
71  err += params.fetch_double("convergence_criterion_squared", Stop_cond_ms);
72  err += params.fetch_int_vector("boundary_condition", bc);
73 
74  if (err) {
75  vout.crucial(m_vl, "Force_F_Overlap_Nf2: fetch error, input parameter not found.\n");
76  abort();
77  }
78 
79 
80  set_parameters(mq, M0, Np, x_min, x_max, Niter_ms, Stop_cond_ms, bc);
81 }
82 
83 
84 //====================================================================
85 void Force_F_Overlap_Nf2::set_parameters(const double M0, const double mq,
86  const int Np, const double x_min, const double x_max,
87  const int Niter_ms, const double Stop_cond_ms,
88  const valarray<int> bc)
89 {
90  int Ndim = CommonParameters::Ndim();
91 
92  //- print input parameters
93  vout.general(m_vl, "Force_F_Overlap_Nf2 paramters:\n");
94  vout.general(m_vl, " M0 = %10.6f\n", M0);
95  vout.general(m_vl, " mq = %10.6f\n", mq);
96  vout.general(m_vl, " Np = %4d\n", Np);
97  vout.general(m_vl, " x_min = %12.8f\n", x_min);
98  vout.general(m_vl, " x_max = %12.6f\n", x_max);
99  vout.general(m_vl, " Niter_ms = %6d\n", Niter_ms);
100  vout.general(m_vl, " Stop_cond_ms = %12.4e\n", Stop_cond_ms);
101  for (int mu = 0; mu < Ndim; ++mu) {
102  vout.general(m_vl, " boundary[%d] = %2d\n", mu, bc[mu]);
103  }
104 
105  //- range check
106  int err = 0;
107  err += ParameterCheck::non_zero(M0);
108  err += ParameterCheck::non_zero(mq);
109  err += ParameterCheck::non_zero(Np);
110  // NB. x_min,x_max == 0 is allowed.
111  err += ParameterCheck::non_zero(Niter_ms);
112  err += ParameterCheck::square_non_zero(Stop_cond_ms);
113 
114  if (err) {
115  vout.crucial(m_vl, "Force_F_Overlap_Nf2: parameter range check failed.\n");
116  abort();
117  }
118 
119  assert(bc.size() == Ndim);
120 
121  //- store values
122  m_M0 = M0;
123  m_mq = mq;
124  m_Np = Np;
125  m_x_min = x_min;
126  m_x_max = x_max;
127  m_Niter_ms = Niter_ms;
128  m_Stop_cond_ms = Stop_cond_ms;
129 
130  m_boundary.resize(Ndim);
131  for (int mu = 0; mu < Ndim; ++mu) {
132  m_boundary[mu] = bc[mu];
133  }
134 
135  //- post-process
136  m_kappa = 0.5 / (4.0 - m_M0);
137 
138  m_fopr_w = new Fopr_Wilson();
139  m_fopr_w->set_parameters(m_kappa, m_boundary);
140  m_fopr_w->set_config(m_U);
141 
142  m_force_w = new Force_F_Wilson_Nf2;
143  m_force_w->set_parameters(m_kappa, m_boundary);
144  m_force_w->set_config(m_U);
145 
146  //- Zolotarev coefficients
147  m_sigma.resize(m_Np);
148  m_cl.resize(2 * m_Np);
149  m_bl.resize(m_Np);
150 }
151 
152 
153 //====================================================================
154 Field Force_F_Overlap_Nf2::force_core(const Field& psi)
155 {
156  int Nc = CommonParameters::Nc();
157  int Nd = CommonParameters::Nd();
158  int Nvol = CommonParameters::Nvol();
159  int Lvol = CommonParameters::Lvol();
160  int Ndim = CommonParameters::Ndim();
161  int size_psf = Lvol * Nc * Nd;
162 
163  //- solving psi = (H^2)^(-1) phi
164  Field_F psi5(Nvol, 1);
165 
166  psi5 = (Field_F)m_fopr_w->mult_gm5(psi);
167 
168  double ff = m_M0 * m_M0 - 0.25 * m_mq * m_mq;
169  psi5 *= ff;
170 
171  // force determination
172 
173  int Nin = Nc * Nc * 2;
174  Field force(Nin, Nvol, Ndim);
175 
176  force = force_core1(psi, psi5);
177  force *= -1.0;
178 
179  return force;
180 }
181 
182 
183 //====================================================================
184 Field Force_F_Overlap_Nf2::force_core1(const Field_F& psi,
185  const Field_F& psi5)
186 {
187  int Nc = CommonParameters::Nc();
188  int Nd = CommonParameters::Nd();
189  int Nvol = CommonParameters::Nvol();
190  int Ndim = CommonParameters::Ndim();
191  int NinG = 2 * Nc * Nc;
192  int NinF = 2 * Nc * Nd;
193 
194  //- Zolotarev coefficient defined
195  double bmax = m_x_max / m_x_min;
196  Math_Sign_Zolotarev sign_func(m_Np, bmax);
197 
198  sign_func.get_sign_parameters(m_cl, m_bl);
199 
200  for (int i = 0; i < m_Np; ++i) {
201  m_sigma[i] = m_cl[2 * i] * m_x_min * m_x_min;
202  }
203 
204  for (int i = 0; i < m_Np; ++i) {
205  vout.general(m_vl, " %3d %12.4e %12.4e %12.4e\n",
206  i, m_cl[i], m_cl[i + m_Np], m_bl[i]);
207  }
208 
209  //- Shiftsolver
210  int Nshift = m_Np;
211 
212  valarray<Field> psi_l(Nshift);
213  for (int i = 0; i < Nshift; ++i) {
214  psi_l[i].reset(NinF, Nvol, 1);
215  }
216 
217  valarray<Field> psi5_l(Nshift);
218  for (int i = 0; i < Nshift; ++i) {
219  psi5_l[i].reset(NinF, Nvol, 1);
220  }
221 
222  int Nconv;
223  double diff;
224 
225  vout.general(m_vl, "Shift solver in MD\n");
226  vout.general(m_vl, " Number of shift values = %d\n", m_sigma.size());
227 
228  m_fopr_w->set_mode("DdagD");
229 
230  Shiftsolver_CG *solver = new Shiftsolver_CG(m_fopr_w,
231  m_Niter_ms, m_Stop_cond_ms);
232 
233  solver->solve(psi_l, m_sigma, psi, Nconv, diff);
234  solver->solve(psi5_l, m_sigma, psi5, Nconv, diff);
235 
236  delete solver;
237 
238  Field force(NinG, Nvol, Ndim);
239  force = 0.0;
240 
241  Field_F psid(Nvol, 1);
242  Field_F psi5d(Nvol, 1);
243  psid = 0.0;
244  psi5d = 0.0;
245 
246  Field force1(NinG, Nvol, Ndim);
247  Field_F w1(Nvol, 1);
248  Field_F w2(Nvol, 1);
249  Field_F w3(Nvol, 1);
250  Field_F w4(Nvol, 1);
251 
252  for (int ip = 0; ip < m_Np; ++ip) {
253  psid.addpart_ex(0, psi_l[ip], 0, m_bl[ip]);
254  psi5d.addpart_ex(0, psi5_l[ip], 0, m_bl[ip]);
255 
256  // w1 = m_fopr_w->H(psi_l[ip]);
257  // w2 = m_fopr_w->H(psi5_l[ip]);
258  m_fopr_w->set_mode("H");
259  w1 = (Field_F)m_fopr_w->mult(psi_l[ip]);
260  w2 = (Field_F)m_fopr_w->mult(psi5_l[ip]);
261 
262  double coeff_l = (m_cl[2 * ip] - m_cl[2 * m_Np - 1])
263  * m_bl[ip] * m_x_min;
264 
265  //- first term
266  // w3 = m_fopr_w->H(w2);
267  w3 = (Field_F)m_fopr_w->mult(w2);
268  // force1 = m_force_w->force_core1(w3,(Field_F)psi_l[ip]);
269  w4 = (Field_F)psi_l[ip];
270 
271  force1 = m_force_w->force_core1(w3, w4);
272  force += coeff_l * force1;
273  force1 = m_force_w->force_core1(w2, w1);
274  force += coeff_l * force1;
275 
276  //- second term
277  // w3 = m_fopr_w->H(w1);
278  w3 = (Field_F)m_fopr_w->mult(w1);
279  // force1 = m_force_w->force_core1(w3,(Field_F)psi5_l[ip]);
280  w4 = (Field_F)psi5_l[ip];
281 
282  force1 = m_force_w->force_core1(w3, w4);
283  force += coeff_l * force1;
284  force1 = m_force_w->force_core1(w1, w2);
285  force += coeff_l * force1;
286  }
287 
288  double coeff1 = m_cl[2 * m_Np - 1] * m_x_min * m_x_min;
289  double coeff2 = 1.0 / m_x_min;
290 
291  //- first term
292  // w1 = m_fopr_w->H(psid);
293  // w2 = m_fopr_w->H(w1);
294  w1 = (Field_F)m_fopr_w->mult(psid);
295  w2 = (Field_F)m_fopr_w->mult(w1);
296  w2.addpart_ex(0, psid, 0, coeff1);
297 
298  force1 = m_force_w->force_core1(psi5, w2);
299  force += coeff2 * force1;
300 
301  //- second term
302  // w1 = m_fopr_w->H(psi5d);
303  // w2 = m_fopr_w->H(w1);
304  w1 = (Field_F)m_fopr_w->mult(psi5d);
305  w2 = (Field_F)m_fopr_w->mult(w1);
306  w2.addpart_ex(0, psi5d, 0, coeff1);
307 
308  force1 = m_force_w->force_core1(psi, w2);
309  force += coeff2 * force1;
310 
311 
312  return force;
313 }
314 
315 
316 //====================================================================
317 //===========================================================END======