HILA
Loading...
Searching...
No Matches
test_forces.cpp
1// Representations introduces a datatype and must be before field.h
2#include "datatypes/representations.h"
3// Test.h includes field.h
4#include "test.h"
5#include "hmc/hmc.h"
6#include "hmc/gauge_field.h"
7#include "hmc/smearing.h"
8#include "dirac/wilson.h"
9#include "dirac/staggered.h"
10#include "hmc/fermion_field.h"
11
12constexpr int N = 2;
13
14template <typename fermion_action, typename dirac, typename gauge_field_type>
15void check_forces(fermion_action &fa, dirac &D, gauge_field_type &gauge) {
16 using sun = typename gauge_field_type::fund_type;
17 using forcetype = SquareMatrix<gauge_field_type::N, Complex<double>>;
18 Field<forcetype> force[NDIM];
19 double eps = 1e-5;
20
21 for (int ng = 0; ng < 1; ng++) { // matrix::generator_count(); ng++){
22 // gauge.set_unity();
23 fa.draw_gaussian_fields();
24 gauge.zero_momentum();
25 foralldir(dir) { force[dir][ALL] = 0; }
26
27 CoordinateVector coord(0);
28
29 sun g1 = gauge.get_gauge(0).get_element(coord);
30 sun h = sun(1) + Complex<double>(0, eps) * sun::generator(ng);
31 sun g12 = h * g1;
32
33 Field<typename dirac::vector_type> psi, chi, tmp, tmp2;
34#if NDIM > 3
35 psi.set_boundary_condition(e_t, hila::bc::ANTIPERIODIC);
36 psi.set_boundary_condition(-e_t, hila::bc::ANTIPERIODIC);
37#endif
38 chi.copy_boundary_condition(psi);
39 tmp.copy_boundary_condition(psi);
40 tmp2.copy_boundary_condition(psi);
41
42 onsites(D.par) {
43 psi[X].gaussian_random();
44 chi[X].gaussian_random();
45 }
46
47 double s1 = 0;
48 D.apply(psi, tmp);
49 onsites(D.par) { s1 += chi[X].rdot(tmp[X]); }
50
51 gauge.get_gauge(0).set_element(g12, coord);
52 gauge.refresh();
53
54 double s2 = 0;
55 D.apply(psi, tmp);
56 onsites(D.par) { s2 += chi[X].rdot(tmp[X]); }
57
58 gauge.get_gauge(0).set_element(g1, coord);
59 gauge.refresh();
60
61 D.force(chi, psi, force, 1);
62
63 gauge.add_momentum(force);
64 sun f = gauge.get_momentum(0).get_element(coord);
65 double f1 = (s2 - s1) / eps;
66 double f2 = (f * Complex<double>(0, 1) * sun::generator(ng)).trace().re;
67 double diff = f2 - f1;
68
69 if (hila::myrank() == 0) {
70 // hila::out0 << "Action 1 " << s1 << "\n";
71 // hila::out0 << "Action 2 " << s2 << "\n";
72 // hila::out0 << "Calculated deriv " << f2 << "\n";
73 // hila::out0 << "Actual deriv " << (s2-s1)/eps << "\n";
74 // hila::out0 << "Fermion deriv " << ng << " diff " << diff << "\n";
75 assert(diff * diff < eps && "Fermion deriv");
76 }
77
78 onsites(D.par) {
79 psi[X].gaussian_random();
80 chi[X].gaussian_random();
81 }
82 gauge.zero_momentum();
83
84 s1 = 0;
85 D.dagger(psi, tmp);
86 onsites(D.par) { s1 += chi[X].rdot(tmp[X]); }
87
88 gauge.get_gauge(0).set_element(g12, coord);
89 gauge.refresh();
90
91 s2 = 0;
92 D.dagger(psi, tmp);
93 onsites(D.par) { s2 += chi[X].rdot(tmp[X]); }
94
95 gauge.get_gauge(0).set_element(g1, coord);
96 gauge.refresh();
97
98 D.force(chi, psi, force, -1);
99 gauge.add_momentum(force);
100 f = gauge.get_momentum(0).get_element(coord);
101 f1 = (s2 - s1) / eps;
102 f2 = (f * Complex<double>(0, 1) * sun::generator(ng)).trace().re;
103 diff = f2 - f1;
104
105 if (hila::myrank() == 0) {
106 // hila::out0 << "Action 1 " << s1 << "\n";
107 // hila::out0 << "Action 2 " << s2 << "\n";
108 // hila::out0 << "Calculated deriv " << f2 << "\n";
109 // hila::out0 << "Actual deriv " << (s2-s1)/eps << "\n";
110 // hila::out0 << "Fermion dg deriv " << ng << " diff " << diff << "\n";
111 assert(diff * diff < eps && "Fermion dg deriv");
112 }
113
114 gauge.zero_momentum();
115 Field<double> sf1, sf2;
116 sf1[ALL] = 0;
117 sf2[ALL] = 0;
118
119 gauge.get_gauge(0).set_element(g12, coord);
120 gauge.refresh();
121 fa.action(sf2);
122
123 gauge.get_gauge(0).set_element(g1, coord);
124 gauge.refresh();
125 fa.action(sf1);
126
127 double dS = 0;
128 s1 = s2 = 0;
129 onsites(ALL) {
130 dS += sf2[X] - sf1[X];
131 s1 += sf1[X];
132 s2 += sf2[X];
133 }
134
135 fa.force_step(1.0);
136 f = gauge.get_momentum(0).get_element(coord);
137 f1 = dS / eps;
138 f2 = (f * Complex<double>(0, 1) * sun::generator(ng)).trace().re;
139 diff = f2 - f1;
140
141 if (hila::myrank() == 0) {
142 // hila::out0 << "Action 1 " << s1 << "\n";
143 // hila::out0 << "Action 2 " << s2 << " " << dS << " " << dS/s2 << "\n";
144 // hila::out0 << "Calculated force " << f2 << "\n";
145 // hila::out0 << "Actual force " << f1 << "\n";
146 // hila::out0 << "Fermion force " << ng << " diff " << diff << "\n";
147 assert(diff * diff < eps && "Fermion force");
148 }
149 }
150}
151
152int main(int argc, char **argv) {
153
154/* Use a smaller lattice size since
155 * the inversion takes a while */
156#if NDIM == 1
157 const CoordinateVector nd = {64};
158#elif NDIM == 2
159 const CoordinateVector nd = {32, 8};
160#elif NDIM == 3
161 const CoordinateVector nd = {16, 8, 8};
162#elif NDIM == 4
163 const CoordinateVector nd = {16, 8, 8, 8};
164#endif
165 hila::initialize(argc, argv);
166 lattice.setup(nd);
167 hila::seed_random(2);
168
169 // Test the force calculation by varying one gauge link
170 // (Needs to be moved to tests)
171 gauge_field<SU<N, double>> gauge;
172 gauge.random();
173 double eps = 1e-5;
174
175 gauge_action ga(gauge, 1.0);
176
177 for (int ng = 0; ng < SU<N>::generator_count(); ng++) {
178 foralldir(dir) {
179 onsites(ALL) { gauge.momentum[dir][X] = 0; }
180 }
181 SU<N> g1 = gauge.gauge[0].get_value_at(50);
182 SU<N> h = SU<N>(1) + Complex<double>(0, eps) * SU<N>::generator(ng);
183 SU<N> g12 = h * g1;
184
185 double s1 = ga.action();
186
187 if (hila::myrank() == 0)
188 gauge.gauge[0].set_value_at(g12, 50);
189 gauge.gauge[0].mark_changed(ALL);
190
191 double s2 = ga.action();
192
193 if (hila::myrank() == 0)
194 gauge.gauge[0].set_value_at(g1, 50);
195 gauge.gauge[0].mark_changed(ALL);
196
197 ga.force_step(1.0);
198 SU<N> f = gauge.momentum[0].get_value_at(50);
199 double diff =
200 (f * Complex<double>(0, 1) * SU<N>::generator(ng)).trace().re - (s2 - s1) / eps;
201
202 if (hila::myrank() == 0) {
203 // hila::out << "Force " <<
204 // (f*Complex<double>(0,1)*SU<N>::generator(ng)).trace().re << "\n"; hila::out
205 // << "Force " << (f*SU<N>::generator(ng)).trace().re << "\n"; hila::out <<
206 // "Deriv " << (s2-s1)/eps << "\n"; hila::out << "Force " << ng << " diff "
207 // << diff << "\n";
208 h = Complex<double>(0, 1) * SU<N>::generator(ng);
209 assert(diff * diff < eps * 10 && "Gauge force");
210 }
211 }
212
213 using SUN = SU<N, double>;
214 using adj = adjointRep<N, double>;
215 using sym = symmetric<N, double>;
216 using asym = antisymmetric<N, double>;
217
218 // Define represented gauge fields for the Dirac force tests
219 gauge.random();
220 represented_gauge_field<adj> adj_gauge(gauge);
221 represented_gauge_field<sym> sym_gauge(gauge);
222 represented_gauge_field<asym> asym_gauge(gauge);
223 stout_smeared_field<SUN> stout_gauge(gauge, 0.1, 4, 10);
224 HEX_smeared_field<SUN> hex_gauge(gauge, 10);
225
226 // Staggered Forces
227 hila::out0 << "Checking staggered forces:\n";
228 {
229 dirac_staggered D(5.0, gauge);
230 fermion_action fa(D, gauge);
231 check_forces(fa, D, gauge);
232 }
233
234 {
235 hila::out0 << "Checking evenodd preconditioned staggered forces:\n";
236 dirac_staggered_evenodd D(5.0, gauge);
237 fermion_action fa(D, gauge);
238 check_forces(fa, D, gauge);
239 }
240
241 {
242 hila::out0 << "Checking stout smeared forces:\n";
243 dirac_staggered_evenodd D(5.0, stout_gauge);
244 fermion_action fa(D, stout_gauge);
245 check_forces(fa, D, stout_gauge);
246 }
247
248 {
249 hila::out0 << "Checking HEX smeared forces:\n";
250 dirac_staggered_evenodd D(5.0, hex_gauge);
251 fermion_action fa(D, hex_gauge);
252 check_forces(fa, D, hex_gauge);
253 }
254
255 {
256 hila::out0 << "Checking adjoint staggered forces:\n";
257 dirac_staggered_evenodd D(5.0, adj_gauge);
258 fermion_action fa(D, adj_gauge);
259 check_forces(fa, D, adj_gauge);
260 }
261
262 {
263 hila::out0 << "Checking symmetric staggered forces:\n";
264 dirac_staggered_evenodd D(5.0, sym_gauge);
265 fermion_action fa(D, sym_gauge);
266 check_forces(fa, D, sym_gauge);
267 }
268
269 {
270 hila::out0 << "Checking antisymmetric staggered forces:\n";
271 dirac_staggered_evenodd D(5.0, asym_gauge);
272 fermion_action fa(D, asym_gauge);
273 check_forces(fa, D, asym_gauge);
274 }
275
276 // Wilson forces
277 {
278 hila::out0 << "Checking Wilson forces:\n";
279 Dirac_Wilson D(0.05, gauge);
280 fermion_action fa(D, gauge);
281 check_forces(fa, D, gauge);
282 }
283
284 {
285 hila::out0 << "Checking evenodd Wilson forces:\n";
286 Dirac_Wilson_evenodd D(0.12, gauge);
287 fermion_action fa(D, gauge);
288 check_forces(fa, D, gauge);
289
290 hila::out0 << "Checking hasenbusch 1:\n";
291 Hasenbusch_action_1 fa1(D, gauge, 0);
292 check_forces(fa1, D, gauge);
293
294 hila::out0 << "Checking hasenbusch 2:\n";
295 Hasenbusch_action_2 fa2(D, gauge, 0);
296 check_forces(fa2, D, gauge);
297 }
298
299 {
300 hila::out0 << "Checking adjoint Wilson forces:\n";
301 Dirac_Wilson_evenodd D(0.05, adj_gauge);
302 fermion_action fa(D, adj_gauge);
303 check_forces(fa, D, adj_gauge);
304 }
305
306 {
307 hila::out0 << "Checking symmetric Wilson forces:\n";
308 Dirac_Wilson_evenodd D(0.05, sym_gauge);
309 fermion_action fa(D, sym_gauge);
310 check_forces(fa, D, sym_gauge);
311 }
312
313 {
314 hila::out0 << "Checking antisymmetric Wilson forces:\n";
315 Dirac_Wilson_evenodd D(0.05, asym_gauge);
316 fermion_action fa(D, asym_gauge);
317 check_forces(fa, D, asym_gauge);
318 }
319
320 // Check also the momentum action and derivative
321 for (int ng = 0; ng < SU<N>::generator_count(); ng++) {
322 ga.draw_gaussian_fields();
323
324 double s1 = ga.action();
325 SU<N> h = gauge.momentum[0].get_value_at(0);
326 h += eps * Complex<double>(0, 1) * SU<N>::generator(ng);
327 if (hila::myrank() == 0)
328 gauge.momentum[0].set_value_at(h, 0);
329 double s2 = ga.action();
330
331 double diff = (h * SU<N>::generator(ng)).trace().re + (s2 - s1) / eps;
332 if (hila::myrank() == 0) {
333 // hila::out << "Mom 1 " << (h*SU<N>::generator(ng)).trace().re << "\n";
334 // hila::out << "Mom 2 " << (s2-s1)/eps << "\n";
335 // hila::out << "Mom " << ng << " diff " << diff << "\n";
336 h = Complex<double>(0, 1) * SU<N>::generator(ng);
337 assert(diff * diff < eps * 10 && "Momentum derivative");
338 }
339 }
340
341 hila::finishrun();
342}
Complex
Complex definition.
Definition cmplx.h:50
CoordinateVector_t< int >
Dirac_Wilson_evenodd
Definition wilson.h:234
Dirac_Wilson
Definition wilson.h:146
Field
The field class implements the standard methods for accessing Fields. Hilapp replaces the parity acce...
Definition field.h:61
Field::set_boundary_condition
void set_boundary_condition(Direction dir, hila::bc bc)
Set the boundary condition in a given Direction (periodic or antiperiodic)
Definition field.h:580
Field::set_element
const T set_element(const CoordinateVector &coord, const A &value)
Get singular element which will be broadcast to all nodes.
Definition field.h:1241
Field::get_value_at
auto get_value_at(const unsigned i) const
Get an individual element outside a loop. This is also used as a getter in the vanilla code.
Definition field.h:684
Field::set_value_at
void set_value_at(const A &value, unsigned i)
Set an individual element outside a loop. This is also used as a getter in the vanilla code.
Definition field.h:709
Field::copy_boundary_condition
void copy_boundary_condition(const Field< A > &rhs)
Copy the boundary condition from another field.
Definition field.h:650
Hasenbusch_action_1
Definition fermion_field.h:182
Hasenbusch_action_2
Definition fermion_field.h:207
Matrix
Matrix class which defines matrix operations.
Definition matrix.h:1679
SU
Class for SU(N) matrix.
Definition sun_matrix.h:37
dirac_staggered_evenodd
Definition libraries/dirac/staggered.h:193
dirac_staggered
Definition libraries/dirac/staggered.h:92
fermion_action
Definition fermion_field.h:22
fermion_action::draw_gaussian_fields
void draw_gaussian_fields()
Definition fermion_field.h:123
fermion_action::action
double action()
Definition fermion_field.h:90
fermion_action::force_step
void force_step(double eps)
Definition fermion_field.h:146
gauge_action
Definition gauge_field.h:448
gauge_field_base::momentum
Field< sun > momentum[4]
Definition gauge_field.h:24
gauge_field_base::gauge
Field< sun > gauge[4]
A matrix field for each Direction.
Definition gauge_field.h:21
gauge_field
Definition gauge_field.h:193
gauge_field::random
void random()
Draw a random gauge field.
Definition gauge_field.h:215
lattice_struct::setup
void setup(const CoordinateVector &siz)
General lattice setup.
Definition lattice.cpp:33
represented_gauge_field
Definition gauge_field.h:297
stout_smeared_field
Definition smearing.h:76
foralldir
#define foralldir(d)
Macro to loop over (all) Direction(s)
Definition coordinates.h:78
ALL
constexpr Parity ALL
bit pattern: 011
Definition coordinates.h:180
hila::myrank
int myrank()
rank of this node
Definition com_mpi.cpp:235
hila::out0
std::ostream out0
This writes output only from main process (node 0)
hila::initialize
void initialize(int argc, char **argv)
Read in command line arguments. Initialise default stream and MPI communication.
Definition initialize.cpp:65
hila::seed_random
void seed_random(uint64_t seed, bool device_rng=true)
Seed random generators with 64-bit unsigned value. On MPI shuffles the seed so that different MPI ran...
Definition random.cpp:86
hila::finishrun
void finishrun()
Normal, controlled exit - all nodes must call this. Prints timing information and information about c...
Definition initialize.cpp:316