HILA
Loading...
Searching...
No Matches
suN_hmc.cpp
1#include "hila.h"
2#include "gauge/polyakov.h"
4
5#ifndef NCOLOR
6#define NCOLOR 3
7#endif
8
9#ifdef HMCACTION
10#define HMCS HMCACTION
11#else
12#define HMCS 1
13#endif
14
15#if HMCS == 1
17#elif HMCS > 1
19#if HMCS < 5
21#elif HMCS == 5
23#endif
24#endif
25
26
27#include "gauge/gradient_flow.h"
28#include "tools/string_format.h"
30
31#ifdef STOUTSMEAR
32#include "gauge/stout_smear.h"
33#define STOUTSTEPS STOUTSMEAR
34#else
35#define STOUTSTEPS 0
36#endif
37
38
39using ftype = double;
41
42ftype stoutc = 0.15;
43int stout_nsteps = STOUTSTEPS;
44
45
46// define a struct to hold the input parameters: this
47// makes it simpler to pass the values around
48struct parameters {
49 ftype beta; // inverse gauge coupling
50 ftype dt; // HMC time step
51 int trajlen; // number of HMC time steps per trajectory
52 int n_traj; // number of trajectories to generate
53 int n_therm; // number of thermalization trajectories (counts only accepted traj.)
54 int gflow_freq; // number of trajectories between gflow measurements
55 ftype gflow_max_l; // flow scale at which gradient flow stops
56 ftype gflow_l_step; // flow scale interval between flow measurements
57 ftype gflow_a_accu; // desired absolute accuracy of gradient flow integration steps
58 ftype gflow_r_accu; // desired relative accuracy of gradient flow integration steps
59 int n_save; // number of trajectories between config. check point
60 std::string config_file;
61 ftype time_offset;
62};
63
64
65///////////////////////////////////////////////////////////////////////////////////
66// HMC functions
67
68template <typename group, typename atype = hila::arithmetic_type<group>>
69double measure_s(const GaugeField<group> &U) {
70 // compute the value of the chosen gauge action (without beta/N factor)
71#if STOUTSTEPS > 0
72
74 stout_smear(U, tU, stoutc, stout_nsteps);
75
76#else
77
78 const GaugeField<group> &tU = U;
79
80#endif
81
82#if HMCS == 1 // BP
83 double res = measure_s_bp(tU);
84#elif HMCS == 2 // LW
85 atype c12 = -1.0 / 12.0; // rectangle weight
86 atype c11 = 1.0 - 8.0 * c12; // plaquette weight
87 double res = measure_s_impr(tU, c11, c12);
88#elif HMCS == 3 // IWASAKI
89 atype c12 = -0.331; // rectangle weight
90 atype c11 = 1.0 - 8.0 * c12; // plaquette weight
91 double res = measure_s_impr(tU, c11, c12);
92#elif HMCS == 4 // DBW2
93 atype c12 = -1.4088; // rectangle weight
94 atype c11 = 1.0 - 8.0 * c12; // plaquette weight
95 double res = measure_s_impr(tU, c11, c12);
96#elif HMCS == 5 // LOG-PLAQUETTE
97 double res = measure_s_log_plaq(tU);
98#else // WILSON
99 double res = measure_s_wplaq(tU);
100#endif // END HMCS
101
102 return res;
103}
104
105template <typename group, typename atype = hila::arithmetic_type<group>>
106void update_E(const GaugeField<group> &U, VectorField<Algebra<group>> &E, atype delta) {
107 // compute the force for the chosen action and use it to evolve E
108 atype eps = delta / group::size();
109
110#if STOUTSTEPS > 0
111
112 std::vector<GaugeField<group>> tUl(stout_nsteps + 1);
113 std::vector<VectorField<group>> tUKl(stout_nsteps);
114 std::vector<VectorField<group>> tstapl(stout_nsteps);
115 stout_smeark(U, tUl, tstapl, tUKl, stoutc);
116
118
119 foralldir(d1) onsites(ALL) tE[d1][X] = 0;
120
121 GaugeField<group> &tU = tUl[stout_nsteps];
122
123#else // STOUTSTEPS==0
124
125 const GaugeField<group> &tU = U;
127
128#endif // END STOUTSTEPS
129
130#if HMCS == 1 // BP
131 get_force_bp_add(tU, tE, eps);
132#elif HMCS == 2 // LW
133 atype c12 = -1.0 / 12.0; // rectangle weight
134 atype c11 = 1.0 - 8.0 * c12; // plaquette weight
135 get_force_impr_add(tU, tE, eps * c11, eps * c12);
136#elif HMCS == 3 // IWASAKI
137 atype c12 = -0.331; // rectangle weight
138 atype c11 = 1.0 - 8.0 * c12; // plaquette weight
139 get_force_impr_add(tU, tE, eps * c11, eps * c12);
140#elif HMCS == 4 // DBW2
141 atype c12 = -1.4088; // rectangle weight
142 atype c11 = 1.0 - 8.0 * c12; // plaquette weight
143 get_force_impr_add(tU, tE, eps * c11, eps * c12);
144#elif HMCS == 5 // LOG-PLAQUETTE
145 get_force_log_plaq_add(tU, tE, eps);
146#else // WILSON
147 get_force_wplaq_add(tU, tE, eps);
148#endif // END HMCS
149
150#if STOUTSTEPS > 0
151
153 stout_smeark_force(tUl, tstapl, tUKl, tE, KS, stoutc);
154
155 foralldir(d1) {
156 onsites(ALL) E[d1][X] += KS[d1][X];
157 }
158
159#endif // END STOUTSTEPS
160
161}
162
163template <typename group, typename atype = hila::arithmetic_type<group>>
164void update_U(GaugeField<group> &U, const VectorField<Algebra<group>> &E, atype delta) {
165 // evolve U with momentum E over time step delta
166 foralldir(d) {
167 onsites(ALL) U[d][X] = chexp(E[d][X].expand_scaled(delta)) * U[d][X];
168 }
169}
170
171template <typename group>
172double measure_e2(const VectorField<Algebra<group>> &E) {
173 // compute gauge kinetic energy from momentum field E
174 Reduction<double> e2 = 0;
175 e2.allreduce(false).delayed(true);
176 foralldir(d) {
177 onsites(ALL) e2 += E[d][X].squarenorm();
178 }
179 return e2.value() / 2;
180}
181
182template <typename group>
183double measure_action(const GaugeField<group> &U, const VectorField<Algebra<group>> &E,
184 const parameters &p, double &plaq) {
185 // measure the total action, consisting of plaquette and momentum term
186 plaq = p.beta * measure_s(U);
187 double e2 = measure_e2(E);
188 return plaq + e2 / 2;
189}
190
191template <typename group>
192double measure_action(const GaugeField<group> &U, const VectorField<Algebra<group>> &E,
193 const parameters &p) {
194 // measure the total action, consisting of plaquette and momentum term
195 double plaq = 0;
196 return measure_action(U, E, p, plaq);
197}
198
199
200template <typename group, typename atype = hila::arithmetic_type<group>>
201void do_trajectory(GaugeField<group> &U, VectorField<Algebra<group>> &E, const parameters &p) {
202 // leap frog integration for normal action
203
204 // start trajectory: advance U by half a time step
205 update_U(U, E, p.dt / 2);
206 // main trajectory integration:
207 for (int n = 0; n < p.trajlen - 1; ++n) {
208 update_E(U, E, p.beta * p.dt);
209 update_U(U, E, p.dt);
210 }
211 // end trajectory: bring U and E to the same time
212 update_E(U, E, p.beta * p.dt);
213 update_U(U, E, p.dt / 2);
214
216}
217
218// end HMC functions
219///////////////////////////////////////////////////////////////////////////////////
220// measurement functions
221
222template <typename group>
224 // perform measurements on current gauge and momentum pair (U, E) and
225 // print results in formatted form to standard output
226 static bool first = true;
227 if (first) {
228 // print legend for measurement output
230 << "LMEAS: s-local plaq E^2 P.real P.imag\n";
231 first = false;
232 }
233 auto slocal = measure_s(U) / (lattice.volume() * NDIM * (NDIM - 1) / 2);
234 auto plaq = measure_s_wplaq(U) / (lattice.volume() * NDIM * (NDIM - 1) / 2);
235 auto e2 = measure_e2(E) / (lattice.volume() * NDIM);
236 auto poly = measure_polyakov(U, e_t);
237 hila::out0 << string_format("MEAS % 0.6e % 0.6e % 0.6e % 0.6e % 0.6e", slocal, plaq, e2,
238 poly.real(), poly.imag())
239 << '\n';
240}
241
242// end measurement functions
243///////////////////////////////////////////////////////////////////////////////////
244// load/save config functions
245
246template <typename group>
247void checkpoint(const GaugeField<group> &U, int trajectory, const parameters &p) {
248 double t = hila::gettime();
249 // name of config with extra suffix
250 std::string config_file =
251 p.config_file + "_" + std::to_string(abs((trajectory + 1) / p.n_save) % 2);
252 // save config
253 U.config_write(config_file);
254 // write run_status file
255 if (hila::myrank() == 0) {
256 std::ofstream outf;
257 outf.open("run_status", std::ios::out | std::ios::trunc);
258 outf << "trajectory " << trajectory + 1 << '\n';
259 outf << "seed " << static_cast<uint64_t>(hila::random() * (1UL << 61)) << '\n';
260 outf << "time " << hila::gettime() << '\n';
261 // write config name to status file:
262 outf << "config name " << config_file << '\n';
263 outf.close();
264 }
265 std::stringstream msg;
266 msg << "Checkpointing, time " << hila::gettime() - t;
267 hila::timestamp(msg.str().c_str());
268}
269
270template <typename group>
271bool restore_checkpoint(GaugeField<group> &U, int &trajectory, parameters &p) {
272 uint64_t seed;
273 bool ok = true;
274 p.time_offset = 0;
275 hila::input status;
276 if (status.open("run_status", false, false)) {
277 hila::out0 << "RESTORING FROM CHECKPOINT:\n";
278 trajectory = status.get("trajectory");
279 seed = status.get("seed");
280 p.time_offset = status.get("time");
281 // get config name with suffix from status file:
282 std::string config_file = status.get("config name");
283 status.close();
284 hila::seed_random(seed);
285 U.config_read(config_file);
286 ok = true;
287 } else {
288 std::ifstream in;
289 in.open(p.config_file, std::ios::in | std::ios::binary);
290 if (in.is_open()) {
291 in.close();
292 hila::out0 << "READING initial config\n";
293 U.config_read(p.config_file);
294 ok = true;
295 } else {
296 ok = false;
297 }
298 }
299 return ok;
300}
301
302// end load/save config functions
303///////////////////////////////////////////////////////////////////////////////////
304
305
306int main(int argc, char **argv) {
307
308 // hila::initialize should be called as early as possible
309 hila::initialize(argc, argv);
310
311 // hila provides an input class hila::input, which is
312 // a convenient way to read in parameters from input files.
313 // parameters are presented as key - value pairs, as an example
314 // " lattice size 64, 64, 64, 64"
315 // is read below.
316 //
317 // Values are broadcast to all MPI nodes.
318 //
319 // .get() -method can read many different input types,
320 // see file "input.h" for documentation
321
322#if HMCS == 1
323 hila::out0 << "SU(" << mygroup::size() << ") HMC with bulk-prevention gauge action\n";
324#elif HMCS == 2
325 hila::out0 << "SU(" << mygroup::size() << ") HMC with Luscher-Weisz gauge action\n";
326#elif HMCS == 3
327 hila::out0 << "SU(" << mygroup::size() << ") HMC with Iwasaki gauge action\n";
328#elif HMCS == 4
329 hila::out0 << "SU(" << mygroup::size() << ") HMC with DBW2 gauge action\n";
330#elif HMCS == 5
331 hila::out0 << "SU(" << mygroup::size() << ") HMC with log-plaquette gauge action\n";
332#elif HMCS == 6
333 hila::out0 << "SU(" << mygroup::size() << ") HMC with topology suppressing gauge action\n";
334#else
335 hila::out0 << "SU(" << mygroup::size() << ") HMC with Wilson's plaquette gauge action\n";
336#endif
337
338 hila::out0 << "Using floating point epsilon: " << fp<ftype>::epsilon << "\n";
339
340 parameters p;
341
342 hila::input par("parameters");
343
344 CoordinateVector lsize;
345 // reads NDIM numbers
346 lsize = par.get("lattice size");
347 // inverse gauge coupling
348 p.beta = par.get("beta");
349 // HMC step size
350 p.dt = par.get("dt");
351 // trajectory length in steps
352 p.trajlen = par.get("trajectory length");
353 // number of trajectories
354 p.n_traj = par.get("number of trajectories");
355 // number of thermalization trajectories
356 p.n_therm = par.get("thermalization trajs");
357 // wilson flow frequency (number of traj. between w. flow measurement)
358 p.gflow_freq = par.get("gflow freq");
359 // wilson flow max. flow-distance
360 p.gflow_max_l = par.get("gflow max lambda");
361 // wilson flow flow-distance step size
362 p.gflow_l_step = par.get("gflow lambda step");
363 // wilson flow absolute accuracy (per integration step)
364 p.gflow_a_accu = par.get("gflow abs. accuracy");
365 // wilson flow relative accuracy (per integration step)
366 p.gflow_r_accu = par.get("gflow rel. accuracy");
367 // random seed = 0 -> get seed from time
368 long seed = par.get("random seed");
369 // save config and checkpoint
370 p.n_save = par.get("trajs/saved");
371 // measure surface properties and print "profile"
372 p.config_file = par.get("config name");
373
374 par.close(); // file is closed also when par goes out of scope
375
376 // set up the lattice
377 lattice.setup(lsize);
378
379 // We need random number here
380 hila::seed_random(seed);
381
382 // Alloc gauge field and momenta (E)
385
386 if (0 && hila::myrank() == 0) {
387 // test matrix exponential and corresponding differential computations
388 Alg_gen<NCOLOR, ftype> genlist[NCOLOR * NCOLOR - 1];
390 if(false) {
391 for (int i = 0; i < NCOLOR * NCOLOR - 1; ++i) {
392 hila::out0 << hila::prettyprint(genlist[i].to_matrix()) << "\n\n";
393 }
394 }
395 Matrix<NCOLOR * NCOLOR - 1, NCOLOR * NCOLOR - 1, ftype> omat;
397 Matrix<NCOLOR, NCOLOR, Complex<ftype>> dV[NCOLOR][NCOLOR];
398 chexp(2.0*genlist[NCOLOR+1].to_matrix()+1.5*genlist[1].to_matrix(), V, dV);
399 hila::out0 << "exp(A):\n" << hila::prettyprint(V) << "\n";
400 for (int i = 0; i < NCOLOR; ++i) {
401 for (int j = 0; j < NCOLOR; ++j) {
402 dV[i][j] = dV[i][j] * V.dagger();
403 hila::out0 << "dexp(A)/dA[" << i << "][" << j << "]*exp(-A):\n"
404 << hila::prettyprint(dV[i][j]) << "\n";
405 }
406 }
407 project_to_algebra_bilinear(dV, omat, genlist);
408 hila::out0 << "dexp(A)^i/dA^j*exp(-A):\n" << hila::prettyprint(omat) << "\n";
409
410 if (false) {
411 Alg_gen<NCOLOR, ftype> genprodlist[NCOLOR * NCOLOR - 1][NCOLOR * NCOLOR - 1];
412 Algebra<mygroup>::generator_product_list(genlist, genprodlist);
413
414 for (int i = 0; i < NCOLOR * NCOLOR - 1; ++i) {
415 for (int j = 0; j < NCOLOR * NCOLOR - 1; ++j) {
416 hila::out0 << hila::prettyprint(genprodlist[i][j].to_matrix()) << "\n\n";
417 }
418 }
419 project_to_algebra_bilinear(V, omat, genprodlist);
420 hila::out0 << "algebra_bilinear omat[][]=2*ReTr(t[].dagger()*V*t[]):\n"
421 << hila::prettyprint(omat) << "\n";
422 }
423 }
424
425 int start_traj = 0;
426 if (!restore_checkpoint(U, start_traj, p)) {
427 U = 1;
428 hila::out0 << "cold start: initial link variables set to identity\n";
429 }
430
431 hila::timer update_timer("Updates");
432 hila::timer measure_timer("Measurements");
433 hila::timer gf_timer("Gradient Flow");
434
435 auto orig_dt = p.dt;
436 auto orig_trajlen = p.trajlen;
438 int nreject = 0;
439 ftype t_step0 = 0.0;
440 double g_act_old, act_old, g_act_new, act_new;
441 g_act_old = p.beta * measure_s(U);
442
443 for (int trajectory = start_traj; trajectory < p.n_traj; ++trajectory) {
444 if (trajectory < p.n_therm) {
445 // during thermalization: start with 10% of normal step size (and trajectory length)
446 // and increse linearly with number of accepted thermalization trajectories. normal
447 // step size is reached after 3/4*p.n_therm accepted thermalization trajectories.
448 if (trajectory < p.n_therm * 3.0 / 4.0) {
449 p.dt = orig_dt * (0.1 + 0.9 * 4.0 / 3.0 * trajectory / p.n_therm);
450 p.trajlen = orig_trajlen;
451 } else {
452 p.dt = orig_dt;
453 p.trajlen = orig_trajlen;
454 }
455 if (nreject > 1) {
456 // if two consecutive thermalization trajectories are rejected, decrese the
457 // step size by factor 0.5 (but keeping trajectory length constant, since
458 // thermalization becomes inefficient if trajectory length decreases)
459 for (int irej = 0; irej < nreject - 1; ++irej) {
460 p.dt *= 0.5;
461 p.trajlen *= 2;
462 }
463 hila::out0 << " thermalization step size(reduzed due to multiple reject) dt="
464 << std::setprecision(8) << p.dt << '\n';
465 } else {
466 hila::out0 << " thermalization step size dt=" << std::setprecision(8) << p.dt
467 << '\n';
468 }
469 } else if (trajectory == p.n_therm) {
470 p.dt = orig_dt;
471 p.trajlen = orig_trajlen;
472 hila::out0 << " normal stepsize dt=" << std::setprecision(8) << p.dt << '\n';
473 }
474
475 update_timer.start();
476
477 U_old = U;
478 ftype ttime = hila::gettime();
479
480 foralldir(d) onsites(ALL) E[d][X].gaussian_random();
481
482 act_old = g_act_old + measure_e2(E) / 2;
483
484 do_trajectory(U, E, p);
485
486 act_new = measure_action(U, E, p, g_act_new);
487
488
489 bool reject = hila::broadcast(exp(act_old - act_new) < hila::random());
490
491 if (trajectory < p.n_therm) {
492 // during thermalization: keep track of number of rejected trajectories
493 if (reject) {
494 ++nreject;
495 --trajectory;
496 } else {
497 if (nreject > 0) {
498 --nreject;
499 }
500 }
501 }
502
503 hila::out0 << std::setprecision(12) << "HMC " << trajectory << " S_TOT_start " << act_old
504 << " dS_TOT " << std::setprecision(6) << act_new - act_old
505 << std::setprecision(12);
506 if (reject) {
507 hila::out0 << " REJECT" << " --> S_GAUGE " << g_act_old;
508 U = U_old;
509 } else {
510 hila::out0 << " ACCEPT" << " --> S_GAUGE " << g_act_new;
511 g_act_old = g_act_new;
512 }
513 update_timer.stop();
514
515 hila::out0 << " time " << std::setprecision(3) << hila::gettime() - ttime << '\n';
516
517
518 hila::out0 << "Measure_start " << trajectory << '\n';
519
520 measure_timer.start();
521
522 measure_stuff(U, E);
523
524 measure_timer.stop();
525
526 hila::out0 << "Measure_end " << trajectory << '\n';
527
528 if (trajectory >= p.n_therm) {
529
530 if (p.gflow_freq > 0 && trajectory % p.gflow_freq == 0) {
531 int gtrajectory = trajectory / p.gflow_freq;
532 if (p.gflow_l_step > 0) {
533
534 gf_timer.start();
535
536 int nflow_steps = (int)(p.gflow_max_l / p.gflow_l_step);
537
538 ftype gftime = hila::gettime();
539 hila::out0 << "Gflow_start " << gtrajectory << '\n';
540
542
543 ftype t_step = t_step0;
544 measure_gradient_flow_stuff(V, (ftype)0.0, t_step);
545 t_step = do_gradient_flow_adapt(V, (ftype)0.0, p.gflow_l_step, p.gflow_a_accu,
546 p.gflow_r_accu, t_step);
547 measure_gradient_flow_stuff(V, p.gflow_l_step, t_step);
548 t_step0 = t_step;
549
550 for (int i = 1; i < nflow_steps; ++i) {
551
552 t_step =
553 do_gradient_flow_adapt(V, i * p.gflow_l_step, (i + 1) * p.gflow_l_step,
554 p.gflow_a_accu, p.gflow_r_accu, t_step);
555
556 measure_gradient_flow_stuff(V, (i + 1) * p.gflow_l_step, t_step);
557 }
558
559 gf_timer.stop();
560
561 hila::out0 << "Gflow_end " << gtrajectory << " time " << std::setprecision(3)
562 << hila::gettime() - gftime << '\n';
563 }
564 }
565 }
566
567 if (p.n_save > 0 && trajectory>=0 && (trajectory + 1) % p.n_save == 0) {
568 checkpoint(U, trajectory, p);
569 }
570 }
571
573}
Array< n, m, T > exp(Array< n, m, T > a)
Exponential.
Definition array.h:1059
Definition su2.h:7
Gauge field class.
Definition gaugefield.h:22
void reunitarize_gauge()
Reunitarize Gauge Field consisting of matrices.
Definition gaugefield.h:94
void config_write(const std::string &filename) const
config_write writes the gauge field to file, with additional "verifying" header
Definition gaugefield.h:155
static constexpr int size()
Returns size of Vector or square Matrix.
Definition matrix.h:248
Rtype dagger() const
Hermitian conjugate of matrix.
Definition matrix.h:1002
Matrix class which defines matrix operations.
Definition matrix.h:1620
Special reduction class: enables delayed and non-blocking reductions, which are not possible with the...
Definition reduction.h:14
const T value()
Return value of the reduction variable. Wait for the comms if needed.
Definition reduction.h:157
Reduction & allreduce(bool b=true)
allreduce(bool) turns allreduce on or off. By default on.
Definition reduction.h:130
Reduction & delayed(bool b=true)
deferred(bool) turns deferred on or off. By default turns on.
Definition reduction.h:148
Class for SU(N) matrix.
Definition sun_matrix.h:110
hila::input - Class for parsing runtime parameter files.
Definition input.h:52
void close()
Closes input parameter file.
Definition input.cpp:79
returntype get(const std::string &key)
Get next value in stack of read in input string from parameters file.
Definition input.h:269
bool open(const std::string &fname, bool use_cmdline=true, bool exit_on_error=true)
Open file that parameters are read from.
Definition input.cpp:22
void setup(const CoordinateVector &siz)
General lattice setup.
Definition lattice.cpp:33
T abs(const Complex< T > &a)
Return absolute value of Complex number.
Definition cmplx.h:1187
#define foralldir(d)
Macro to loop over (all) Direction(s)
Definition coordinates.h:78
constexpr Parity ALL
bit pattern: 011
int chexp(const Matrix_t< n, m, T, MT > &mat, Matrix_t< n, m, T, MT > &omat, Mt(&pl)[n])
Calculate exp of a square matrix.
Definition matrix.h:2933
double random()
Real valued uniform random number generator.
Definition hila_gpu.cpp:121
int myrank()
rank of this node
Definition com_mpi.cpp:234
std::ostream out0
This writes output only from main process (node 0)
void initialize(int argc, char **argv)
Read in command line arguments. Initialise default stream and MPI communication.
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
T broadcast(T &var, int rank=0)
Broadcast the value of var to all MPI ranks from rank (default=0).
Definition com_mpi.h:152
double gettime()
Definition timing.cpp:163
void finishrun()
Normal, controlled exit - all nodes must call this. Prints timing information and information about c...
Complex< double > measure_polyakov(const GaugeField< T > &U, Direction dir=Direction(NDIM - 1))
Measure Polyakov lines to direction dir.
Definition polyakov.h:17
void do_trajectory(GaugeField< group > &U, const parameters &p)
Evolve gauge field.
void measure_stuff(const GaugeField< group > &U, const parameters &p)
Measures Polyakov lines and Wilson action.
Definition suN_gauge.cpp:41