Bridge++  Ver.2.1.3
force_F_Domainwall.cpp
Go to the documentation of this file.
1 
12 
13 const std::string Force_F_Domainwall::class_name
14  = "Force_F_Domainwall";
15 
16 //====================================================================
18 {
20 
21  std::string vlevel;
22  if (!params.fetch_string("verbose_level", vlevel)) {
23  m_vl = vout.set_verbose_level(vlevel);
24  } else {
26  }
27 
28  vout.general(m_vl, "%s: construction\n", class_name.c_str());
30 
31  int Ndim = CommonParameters::Ndim();
32  int Nvol = CommonParameters::Nvol();
33 
34  m_boundary.resize(Ndim);
35  m_Ns = 0;
36 
37  set_parameters_impl(params);
38 
39  const double kappa = 1.0 / (8.0 - 2.0 * m_M0);
40 
41  Parameters params_kernel = params;
42  params_kernel.set_double("hopping_parameter", kappa);
43 
44  m_fopr_w = new Fopr_Wilson(params_kernel);
45 
46  m_fopr_dw = new Fopr_Domainwall(m_fopr_w, params);
47 
48  m_force_w = new Force_F_Wilson_Nf2(params_kernel);
49 
50  m_force1 = new Field_G(Nvol, Ndim);
51  m_force2 = new Field_G(Nvol, Ndim);
52 
53  m_eta4 = new Field_F(Nvol, 1);
54  m_zeta4 = new Field_F(Nvol, 1);
55 
56  m_Lw4 = new Field_F(Nvol, 1);
57  m_w4t1 = new Field_F(Nvol, 1);
58  m_w4t2 = new Field_F(Nvol, 1);
59 
61  vout.general(m_vl, "%s: construction finished.\n",
62  class_name.c_str());
63 
64 }
65 
66 //====================================================================
68 {
70 
72 
73  vout.general(m_vl, "%s: construction (obsolete)\n",
74  class_name.c_str());
76 
77  int Ndim = CommonParameters::Ndim();
78  int Nvol = CommonParameters::Nvol();
79 
80  m_boundary.resize(Ndim);
81  m_Ns = 0;
82 
83  m_fopr_w = new Fopr_Wilson();
86 
87  m_force1 = new Field_G(Nvol, Ndim);
88  m_force2 = new Field_G(Nvol, Ndim);
89 
90  m_eta4 = new Field_F(Nvol, 1);
91  m_zeta4 = new Field_F(Nvol, 1);
92 
93  m_Lw4 = new Field_F(Nvol, 1);
94  m_w4t1 = new Field_F(Nvol, 1);
95  m_w4t2 = new Field_F(Nvol, 1);
96 
98  vout.general(m_vl, "%s: construction finished.\n",
99  class_name.c_str());
100 
101 }
102 
103 //====================================================================
105 {
106  delete m_force_w;
107  delete m_fopr_dw;
108  delete m_fopr_w;
109 
110  delete m_force1;
111  delete m_force2;
112 
113  delete m_eta4;
114  delete m_zeta4;
115 
116  delete m_Lw4;
117  delete m_w4t1;
118  delete m_w4t2;
119 
120  delete m_zeta;
121  delete m_eta2;
122 
123 }
124 
125 //====================================================================
127 {
128  set_parameters_impl(params);
129 
130  //- Domain-wall operator
131  m_fopr_dw->set_parameters(params);
132 
133  //- Wilson opr/force
134  const double kappa = 1.0 / (8.0 - 2.0 * m_M0);
135 
136  Parameters params_kernel = params;
137  params_kernel.set_double("hopping_parameter", kappa);
138 
139  m_fopr_w->set_parameters(params_kernel);
140  m_force_w->set_parameters(params_kernel);
141 
142 }
143 
144 //====================================================================
146 {
147  const string str_vlevel = params.get_string("verbose_level");
148 
149  m_vl = vout.set_verbose_level(str_vlevel);
150 
151  //- fetch and check input parameters
152  double mq, M0, b, c;
153  int Ns;
154  std::vector<int> bc;
155  double alpha;
156 
157  int err = 0;
158  err += params.fetch_double("quark_mass", mq);
159  err += params.fetch_double("domain_wall_height", M0);
160  err += params.fetch_int("extent_of_5th_dimension", Ns);
161  err += params.fetch_int_vector("boundary_condition", bc);
162  err += params.fetch_double("coefficient_b", b);
163  err += params.fetch_double("coefficient_c", c);
164 
165  if (err) {
166  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
167  class_name.c_str());
168  exit(EXIT_FAILURE);
169  }
170 
171  std::vector<double> b_vec(Ns);
172  std::vector<double> c_vec(Ns);
173  for(int is = 0; is < Ns; ++is){
174  b_vec[is] = b;
175  c_vec[is] = c;
176  }
177 
178  int err3 = params.fetch_double("parameter_alpha", alpha);
179  if (err3) {
180  vout.general(m_vl, " parameter alpha is not provided: set to 1.0.\n");
181  alpha = 1.0;
182  }
183 
184  set_parameters_impl(mq, M0, Ns, bc, b_vec, c_vec, alpha);
185 
186 }
187 
188 //====================================================================
190  const double mq,
191  const double M0,
192  const int Ns,
193  const std::vector<int> bc,
194  const std::vector<double> b_vec,
195  const std::vector<double> c_vec,
196  const double alpha)
197 {
198  set_parameters_impl(mq, M0, Ns, bc, b_vec, c_vec, alpha);
199 
200  //- Domain-wall operator
202  m_b, m_c, m_alpha);
203  //- Wilson opr/force
204  const double kappa = 1.0 / (8.0 - 2.0 * m_M0);
207 
208 }
209 
210 //====================================================================
212  const double mq,
213  const double M0,
214  const int Ns,
215  const std::vector<int> bc,
216  const std::vector<double> b_vec,
217  const std::vector<double> c_vec,
218  const double alpha)
219 {
220 #pragma omp barrier
221 
222  int Ndim = CommonParameters::Ndim();
223 
224  //- range check
225  int err = 0;
226  err += ParameterCheck::non_zero(mq);
227  err += ParameterCheck::non_zero(M0);
228  err += ParameterCheck::non_zero(Ns);
229 
230  if (err) {
231  vout.crucial(m_vl, "Error at %s: parameter range check failed.\n",
232  class_name.c_str());
233  exit(EXIT_FAILURE);
234  }
235 
236  assert(bc.size() == Ndim);
237 
238  int ith = ThreadManager::get_thread_id();
239  if (ith == 0) {
240 
241  int Ns_prev = m_Ns;
242 
243  m_mq = mq;
244  m_M0 = M0;
245  m_Ns = Ns;
246  m_alpha = alpha;
247 
248  m_boundary.resize(Ndim);
249  for(int mu = 0; mu < Ndim; ++mu){
250  m_boundary[mu] = bc[mu];
251  }
252 
253  if (m_b.size() != m_Ns) {
254  m_b.resize(m_Ns);
255  m_c.resize(m_Ns);
256  }
257  for(int is = 0; is < m_Ns; ++is){
258  m_b[is] = b_vec[is];
259  m_c[is] = c_vec[is];
260  }
261 
262  if(m_Ns != Ns_prev){
263  if(Ns_prev != 0){
264  delete m_zeta;
265  delete m_eta2;
266  }
267 
268  int Nvol = CommonParameters::Nvol();
269  m_zeta = new Field_F(Nvol, m_Ns);
270  m_eta2 = new Field_F(Nvol, m_Ns);
271  }
272 
273  }
274 #pragma omp barrier
275 
276  //- print input parameters
277  vout.general(m_vl, "%s: parameters\n", class_name.c_str());
278  vout.general(m_vl, " mq = %12.8f\n", m_mq);
279  vout.general(m_vl, " M0 = %12.8f\n", m_M0);
280  vout.general(m_vl, " Ns = %4d\n", m_Ns);
281  for (int mu = 0; mu < Ndim; ++mu) {
282  vout.general(m_vl, " boundary[%d] = %2d\n", mu, m_boundary[mu]);
283  }
284  vout.general(m_vl, " coefficients:\n");
285  for (int is = 0; is < m_Ns; ++is) {
286  vout.general(m_vl, " b[%2d] = %16.10f c[%2d] = %16.10f\n",
287  is, m_b[is], is, m_c[is]);
288  }
289  vout.general(m_vl, " alpha = %8.4f\n", m_alpha);
290 
291 }
292 
293 //====================================================================
295 {
296  params.set_double("quark_mass", m_mq);
297  params.set_double("domain_wall_height", m_M0);
298  params.set_int("extent_of_5th_dimension", m_Ns);
299  params.set_int_vector("boundary_condition", m_boundary);
300  params.set_double_vector("coefficient_b", m_b);
301  params.set_double_vector("coefficient_c", m_c);
302  params.set_double("parameter_alpha", m_alpha);
303  params.set_string("gamma_matrix_type", m_repr);
304 
305  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
306 }
307 
308 //====================================================================
310 {
311 #pragma omp barrier
312 
313  int ith = ThreadManager::get_thread_id();
314  if (ith == 0) m_U = (Field_G *)U;
315 #pragma omp barrier
316 
317  m_fopr_w->set_config(U);
318  m_fopr_dw->set_config(U);
319  m_force_w->set_config(U);
320 }
321 
322 //====================================================================
323 void Force_F_Domainwall::set_mode(const std::string& mode)
324 {
325 #pragma omp barrier
326 
327  int ith = ThreadManager::get_thread_id();
328  if (ith == 0) m_mode = mode;
329 
330 #pragma omp barrier
331 
332 }
333 
334 //====================================================================
336  const Field& zeta,
337  const Field& eta)
338 {
339  vout.crucial(m_vl, "Error at %s: unimplemented.\n",
340  class_name.c_str());
341  exit(EXIT_FAILURE);
342 }
343 
344 //====================================================================
345 void Force_F_Domainwall::force_udiv(Field& force_, const Field& eta)
346 {
347  m_fopr_dw->set_mode("H");
348  m_fopr_dw->mult(*m_zeta, eta);
349 
350  force_udiv1_H(*m_force1, *m_zeta, eta);
351  copy(force_, *m_force1);
352 
354  axpy(force_, 1.0, *m_force1);
355 
356 }
357 
358 //====================================================================
360  const Field& zeta,
361  const Field& eta)
362 {
363  if(m_mode=="H"){
364  force_udiv1_H(force_, zeta, eta);
365  }else if(m_mode=="Hdag"){
366  force_udiv1_Hdag(force_, zeta, eta);
367  }else{
368  vout.crucial(m_vl, "Error at %s: irrelevant mult mode: %s.\n",
369  class_name.c_str(), m_mode.c_str());
370  exit(EXIT_FAILURE);
371  }
372 
373 }
374 
375 //====================================================================
377  const Field& zeta,
378  const Field& eta)
379 {
380  /*
381  int Nc = CommonParameters::Nc();
382  int Nd = CommonParameters::Nd();
383  int Nvol = CommonParameters::Nvol();
384  int Ndim = CommonParameters::Ndim();
385  int NinF = 2 * Nc * Nd;
386  int NinG = 2 * Nc * Nc;
387  */
388  //Field force1(NinG, Nvol, Ndim);
389 
390  force_.set(0.0);
391 
392  force_Himpl(*m_eta2, eta);
393 
394  for(int is = 0; is < m_Ns; ++is){
395 
396  copy(*m_zeta4, 0, zeta, is);
397  copy(*m_eta4, 0, *m_eta2, is);
398  scal(*m_eta4, 4.0-m_M0);
399 
401  axpy(force_, 1.0, *m_force2);
402 
403  }
404 
405 }
406 
407 //====================================================================
409  const Field& zeta,
410  const Field& eta)
411 {
412  /*
413  int Nc = CommonParameters::Nc();
414  int Nd = CommonParameters::Nd();
415  int Nvol = CommonParameters::Nvol();
416  int Ndim = CommonParameters::Ndim();
417  int NinF = 2 * Nc * Nd;
418  int NinG = 2 * Nc * Nc;
419  */
420  //Field force1(NinG, Nvol, Ndim);
421 
422  force_.set(0.0);
423 
424  force_Himpl(*m_eta2, zeta);
425 
426  for(int is = 0; is < m_Ns; ++is){
427 
428  copy(*m_zeta4, 0, *m_eta2, is);
429  copy(*m_eta4, 0, eta, is);
430  scal(*m_eta4, 4.0-m_M0);
431 
433  axpy(force_, 1.0, *m_force2);
434 
435  }
436 
437 }
438 
439 //====================================================================
441 {
442 #pragma omp barrier
443 
444  /*
445  int NinF = eta.nin();
446  int Nvol = eta.nvol();
447 
448  Field w4(NinF,Nvol,1), w4_up(NinF,Nvol,1), w4_dn(NinF,Nvol,1);
449  Field Lw4(NinF,Nvol,1);
450  */
451 
452  for(int is = 0; is < m_Ns; ++is){
453  int isR = m_Ns - 1 - is;
454 
455  m_Lw4->set(0.0);
456 
457  int is_up = (is+1) % m_Ns;
458  double Fup = -0.5;
459  if(is == m_Ns-1) Fup = m_mq/2.0;
460  copy(*m_w4t2, 0, eta, is_up);
462  axpy(*m_w4t2, -1.0, *m_w4t1);
463 
464  axpy(*m_Lw4, Fup, *m_w4t2);
465 
466  int is_dn = (is-1+m_Ns) % m_Ns;
467  double Fdn = -0.5;
468  if(is == 0) Fdn = m_mq/2.0;
469  copy(*m_w4t2, 0, eta, is_dn);
471  axpy(*m_w4t2, 1.0, *m_w4t1);
472  axpy(*m_Lw4, Fdn, *m_w4t2);
473 
474  copy(*m_w4t1, 0, eta, is);
475  scal(*m_w4t1, m_b[is]);
476  axpy(*m_w4t1, -m_c[is], *m_Lw4);
477 
478  copy(eta2, isR, *m_w4t1, 0);
479 #pragma omp barrier
480  }
481 
482 }
483 
484 //============================================================END=====
Force_F_Domainwall::class_name
static const std::string class_name
Definition: force_F_Domainwall.h:34
Force_F_Domainwall::force_udiv1_Hdag
void force_udiv1_Hdag(Field &force_, const Field &zeta, const Field &eta)
Definition: force_F_Domainwall.cpp:408
Force_F_Domainwall::m_Ns
int m_Ns
size of fifth-dimension
Definition: force_F_Domainwall.h:42
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
Force_F_Domainwall::m_c
std::vector< double > m_c
coefficient parameter
Definition: force_F_Domainwall.h:45
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
Force_F_Wilson_Nf2::force_udiv1
void force_udiv1(Field &force, const Field &zeta, const Field &eta)
Definition: force_F_Wilson_Nf2.cpp:241
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
Parameters
Class for parameters.
Definition: parameters.h:46
Force_F_Domainwall::m_alpha
double m_alpha
Neff's improving parameter.
Definition: force_F_Domainwall.h:46
Parameters::set_double
void set_double(const string &key, const double value)
Definition: parameters.cpp:33
Force_F_Domainwall::m_zeta
Field_F * m_zeta
5d working field
Definition: force_F_Domainwall.h:59
Bridge::BridgeIO::decrease_indent
void decrease_indent()
Definition: bridgeIO.cpp:518
Force_F_Domainwall::m_b
std::vector< double > m_b
coefficient parameter
Definition: force_F_Domainwall.h:44
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.cpp:508
Force_F_Domainwall::force_udiv1
void force_udiv1(Field &force, const Field &zeta, const Field &eta)
Definition: force_F_Domainwall.cpp:359
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
Force_F_Domainwall::m_M0
double m_M0
domain-wall height
Definition: force_F_Domainwall.h:41
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:381
Parameters::set_double_vector
void set_double_vector(const string &key, const vector< double > &value)
Definition: parameters.cpp:42
AFopr_Domainwall::mult
void mult(AFIELD &v, const AFIELD &w)
multiplies fermion operator to a given field.
Definition: afopr_Domainwall-tmpl.h:553
Force_F_Domainwall::set_mode
void set_mode(const std::string &mode)
in Force, setting the mode is optional when H is nonhermitian.
Definition: force_F_Domainwall.cpp:323
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:213
Force_F_Domainwall::init
void init()
initial setup.
Definition: force_F_Domainwall.cpp:67
Force_F_Domainwall::m_fopr_w
Fopr_Wilson * m_fopr_w
Definition: force_F_Domainwall.h:53
Force_F_Domainwall::m_w4t2
Field_F * m_w4t2
4d working field
Definition: force_F_Domainwall.h:66
Force_F_Domainwall::force_udiv
void force_udiv(Field &force, const Field &eta)
Definition: force_F_Domainwall.cpp:345
Force_F_Domainwall::force_Himpl
void force_Himpl(Field &eta2, const Field &eta)
Definition: force_F_Domainwall.cpp:440
Parameters::fetch_int_vector
int fetch_int_vector(const string &key, vector< int > &value) const
Definition: parameters.cpp:429
Force_F_Domainwall::m_zeta4
Field_F * m_zeta4
4d working field
Definition: force_F_Domainwall.h:63
Force_F_Domainwall::m_force1
Field_G * m_force1
Definition: force_F_Domainwall.h:56
Fopr_Domainwall
AFopr_Domainwall< Field > Fopr_Domainwall
Domain-wall fermion operator.
Definition: fopr_Domainwall.h:33
Org::Fopr_Wilson::set_parameters
void set_parameters(const Parameters &params)
sets parameters by a Parameter object: to be implemented in a subclass.
Definition: fopr_Wilson_impl.cpp:127
Force_F_Domainwall::m_fopr_dw
Fopr_Domainwall * m_fopr_dw
Definition: force_F_Domainwall.h:52
Force_F_Domainwall::m_force2
Field_G * m_force2
Definition: force_F_Domainwall.h:57
threadManager.h
Force_F_Domainwall::force_udiv1_H
void force_udiv1_H(Field &force_, const Field &zeta, const Field &eta)
Definition: force_F_Domainwall.cpp:376
Force_F_Domainwall::m_eta4
Field_F * m_eta4
4d working field
Definition: force_F_Domainwall.h:62
Force_F_Domainwall::m_vl
Bridge::VerboseLevel m_vl
verbose level
Definition: force_F_Domainwall.h:48
Force_F_Domainwall::set_parameters
void set_parameters(const Parameters &params)
sets parameters by a Parameter object: to be implemented in a subclass.
Definition: force_F_Domainwall.cpp:126
Force_F_Domainwall::m_mode
std::string m_mode
Definition: force_F_Domainwall.h:50
Parameters::set_int_vector
void set_int_vector(const string &key, const vector< int > &value)
Definition: parameters.cpp:45
Org::Fopr_Wilson::mult_gm5
void mult_gm5(Field &v, const Field &w)
multiplies gamma_5 matrix.
Definition: fopr_Wilson_impl.cpp:370
Force_F_Domainwall::tidyup
void tidyup()
final tidyup.
Definition: force_F_Domainwall.cpp:104
Force_F_Domainwall::m_repr
std::string m_repr
gamma matrix representation
Definition: force_F_Domainwall.h:47
Force_F_Domainwall::get_parameters
void get_parameters(Parameters &params) const
Definition: force_F_Domainwall.cpp:294
AFopr_Domainwall::set_config
void set_config(Field *U)
sets the gauge configuration.
Definition: afopr_Domainwall.h:128
Force_F_Domainwall::m_w4t1
Field_F * m_w4t1
4d working field
Definition: force_F_Domainwall.h:65
Force_F_Domainwall::force_core1
void force_core1(Field &force, const Field &zeta, const Field &eta)
force_core1 overrides default implementation in the base class.
Definition: force_F_Domainwall.cpp:335
Force_F_Domainwall::m_eta2
Field_F * m_eta2
5d working field
Definition: force_F_Domainwall.h:60
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
Force_F_Wilson_Nf2::set_parameters
void set_parameters(const Parameters &params)
sets parameters by a Parameter object: to be implemented in a subclass.
Definition: force_F_Wilson_Nf2.cpp:118
AForce_F< Field >::m_U
Field_G * m_U
Gauge configuration.
Definition: aforce_F.h:42
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:195
Force_F_Domainwall::m_Lw4
Field_F * m_Lw4
4d working field
Definition: force_F_Domainwall.h:67
AFopr_Domainwall::set_parameters
void set_parameters(const Parameters &params)
sets parameters by a Parameter object: to be implemented in a subclass.
Definition: afopr_Domainwall-tmpl.h:201
Force_F_Domainwall::set_parameters_impl
void set_parameters_impl(const Parameters &params)
Definition: force_F_Domainwall.cpp:145
ParameterCheck::non_zero
int non_zero(const double v)
Definition: parameterCheck.cpp:32
Force_F_Domainwall::m_mq
double m_mq
quark mass
Definition: force_F_Domainwall.h:40
Parameters::set_int
void set_int(const string &key, const int value)
Definition: parameters.cpp:36
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:262
Org::Fopr_Wilson::set_config
void set_config(Field *U)
sets the gauge configuration.
Definition: fopr_Wilson_impl.cpp:188
Parameters::fetch_string
int fetch_string(const string &key, string &value) const
Definition: parameters.cpp:378
Force_F_Wilson_Nf2
Force for the standard Wilson fermion operator.
Definition: force_F_Wilson_Nf2.h:33
Field_F
Wilson-type fermion field.
Definition: field_F.h:37
Parameters::fetch_double
int fetch_double(const string &key, double &value) const
Definition: parameters.cpp:327
Force_F_Domainwall::set_config
void set_config(Field *U)
sets verbose level.
Definition: force_F_Domainwall.cpp:309
Parameters::get_string
string get_string(const string &key) const
Definition: parameters.cpp:221
Force_F_Domainwall::m_boundary
std::vector< int > m_boundary
boundary conditions
Definition: force_F_Domainwall.h:43
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:242
AFopr_Domainwall::set_mode
void set_mode(std::string mode)
setting the mode of multiplication if necessary. Default implementation here is just to avoid irrelev...
Definition: afopr_Domainwall-tmpl.h:540
Field
Container of Field-type object.
Definition: field.h:46
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
Field_G
SU(N) gauge field.
Definition: field_G.h:38
Parameters::fetch_int
int fetch_int(const string &key, int &value) const
Definition: parameters.cpp:346
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:262
Force_F_Wilson_Nf2::set_config
void set_config(Field *U)
sets verbose level.
Definition: force_F_Wilson_Nf2.cpp:207
Force_F_Domainwall::m_force_w
Force_F_Wilson_Nf2 * m_force_w
Definition: force_F_Domainwall.h:54
force_F_Domainwall.h
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:572
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:216
Fopr_Wilson
Org::Fopr_Wilson Fopr_Wilson
Wilson fermion operator.
Definition: fopr_Wilson.h:50