From a43ff1f98e9b9814f858bccb11c174b418458491 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 10 Oct 2018 21:45:32 -0400 Subject: big rearrangement of files to make libraries and example (research) files clearer, and changed to c++ std lib random numbers --- lib/include/wolff.hpp | 35 ++++++++++++ lib/include/wolff/cluster.hpp | 111 ++++++++++++++++++++++++++++++++++++ lib/include/wolff/finite_states.hpp | 40 +++++++++++++ lib/include/wolff/graph.hpp | 21 +++++++ lib/include/wolff/state.hpp | 89 +++++++++++++++++++++++++++++ lib/include/wolff/types.h | 36 ++++++++++++ 6 files changed, 332 insertions(+) create mode 100644 lib/include/wolff.hpp create mode 100644 lib/include/wolff/cluster.hpp create mode 100644 lib/include/wolff/finite_states.hpp create mode 100644 lib/include/wolff/graph.hpp create mode 100644 lib/include/wolff/state.hpp create mode 100644 lib/include/wolff/types.h (limited to 'lib/include') diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp new file mode 100644 index 0000000..c10a211 --- /dev/null +++ b/lib/include/wolff.hpp @@ -0,0 +1,35 @@ + +#include "wolff/cluster.hpp" +#include "wolff/state.hpp" + +template +void wolff(count_t N, state_t & s, std::function gen_R, std::function &)> measurements, std::mt19937& r, bool silent) { + +#ifdef FINITE_STATES +#ifdef NOFIELD + initialize_probs(s.J, s.T); +#else + initialize_probs(s.J, s.H, s.T); +#endif +#endif + + std::uniform_int_distribution dist(0, s.nv); + + if (!silent) printf("\n"); + for (count_t steps = 0; steps < N; steps++) { + if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s.E, s.last_cluster_size); + + v_t v0 = dist(r); + R_t step = gen_R(r, s.spins[v0]); + flip_cluster (s, v0, step, r); + + measurements(s); + } + + if (!silent) { + printf("\033[F\033[J"); + } + printf("WOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, s.last_cluster_size); + +} + diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp new file mode 100644 index 0000000..104f3c2 --- /dev/null +++ b/lib/include/wolff/cluster.hpp @@ -0,0 +1,111 @@ + +#pragma once + +#include +#include +#include +#include + +#include "types.h" +#include "state.hpp" +#include "graph.hpp" + +template +void flip_cluster(state_t& s, v_t v0, const R_t& r, std::mt19937& rand) { + std::uniform_real_distribution dist(0.0,1.0); + v_t nv = 0; + + std::stack stack; + stack.push(v0); + + std::vector marks(s.g.nv, false); + + while (!stack.empty()) { + v_t v = stack.top(); + stack.pop(); + + if (!marks[v]) { // don't reprocess anyone we've already visited! + marks[v] = true; + + X_t si_new; +#ifndef NOFIELD + R_t R_new; + + bool v_is_ghost = (v == s.nv); // ghost site has the last index + + if (v_is_ghost) { + R_new = r.act(s.R); // if we are, then we're moving the transformation + } else +#endif + { + si_new = r.act(s.spins[v]); // otherwise, we're moving the spin at our site + } + + for (const v_t &vn : s.g.v_adj[v]) { + double dE, prob; + +#ifndef NOFIELD + bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost + + if (v_is_ghost || vn_is_ghost) { // this is a ghost-involved bond + X_t rs_old, rs_new; + v_t non_ghost; + + if (vn_is_ghost) { + // if our neighbor is the ghost, the current site is a normal + // spin - rotate it back! + rs_old = s.R.act_inverse(s.spins[v]); + rs_new = s.R.act_inverse(si_new); + non_ghost = v; + } else { + /* if we're the ghost, we need to rotate our neighbor back in + both the old and new ways */ + rs_old = s.R.act_inverse(s.spins[vn]); + rs_new = R_new.act_inverse(s.spins[vn]); + non_ghost = vn; + } + + dE = s.H(rs_old) - s.H(rs_new); + +#ifdef FINITE_STATES + prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)]; +#endif + + s.update_magnetization(rs_old, rs_new); + s.update_fourierZero(non_ghost, rs_old, rs_new); + } else // this is a perfectly normal bond! +#endif + { + dE = s.J(s.spins[v], s.spins[vn]) - s.J(si_new, s.spins[vn]); + +#ifdef FINITE_STATES + prob = J_probs[state_to_ind(s.spins[v])][state_to_ind(si_new)][state_to_ind(s.spins[vn])]; +#endif + } + + s.update_energy(dE); + +#ifndef FINITE_STATES + prob = 1.0 - exp(-dE / s.T); +#endif + + if (dist(rand) < prob) { + stack.push(vn); // push the neighboring vertex to the stack + } + } + +#ifndef NOFIELD + if (v_is_ghost) { + s.R = R_new; + } else +#endif + { + s.spins[v] = si_new; + nv++; + } + } + } + + s.last_cluster_size = nv; +} + diff --git a/lib/include/wolff/finite_states.hpp b/lib/include/wolff/finite_states.hpp new file mode 100644 index 0000000..426edad --- /dev/null +++ b/lib/include/wolff/finite_states.hpp @@ -0,0 +1,40 @@ + +#pragma once + +#include +#include +#include + +#define FINITE_STATES + +// must have N_STATES, states[N_STATES], and state_to_ind defined before +// invoking header + +std::array, N_STATES>, N_STATES> J_probs; +#ifndef NOFIELD +std::array, N_STATES> H_probs; +#endif + +template +#ifndef NOFIELD +void initialize_probs(std::function J, std::function H, double T) { + for (q_t i = 0; i < N_STATES; i++) { + for (q_t j = 0; j < N_STATES; j++) { + H_probs[i][j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T); + } + } +#else +void initialize_probs(std::function J, double T) { +#endif + for (q_t i = 0; i < N_STATES; i++) { + for (q_t j = 0; j < N_STATES; j++) { + for (q_t k = 0; k < N_STATES; k++) { + J_probs[i][j][k] = 1.0 - exp(-(J(states[i], states[k]) - J(states[j], states[k])) / T); + } + } + } +} + + + + diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp new file mode 100644 index 0000000..165544a --- /dev/null +++ b/lib/include/wolff/graph.hpp @@ -0,0 +1,21 @@ + +#pragma once + +#include +#include +#include +#include + +#include "types.h" + +class graph_t { + public: + v_t ne; + v_t nv; + std::vector> v_adj; + std::vector> coordinate; + + graph_t(D_t D, L_t L); + void add_ext(); +}; + diff --git a/lib/include/wolff/state.hpp b/lib/include/wolff/state.hpp new file mode 100644 index 0000000..1f5e359 --- /dev/null +++ b/lib/include/wolff/state.hpp @@ -0,0 +1,89 @@ + +#pragma once + +#include +#include + +#include "types.h" +#include "graph.hpp" + +template +class state_t { + private: + // updating fourier terms F requires many cos and sin calls, faster to do it beforehand. + std::vector> precomputed_cos; + std::vector> precomputed_sin; + public: + D_t D; + L_t L; + v_t nv; + v_t ne; + graph_t g; + double T; + std::vector spins; + R_t R; + double E; + typename X_t::M_t M; // the "sum" of the spins, like the total magnetization + v_t last_cluster_size; + std::vector ReF; + std::vector ImF; + + std::function J; +#ifndef NOFIELD + std::function H; + + state_t(D_t D, L_t L, double T, std::function J, std::function H) : D(D), L(L), g(D, L), T(T), R(), J(J), H(H) { +#else + state_t(D_t D, L_t L, double T, std::function J) : D(D), L(L), g(D, L), T(T), R(), J(J) { +#endif + nv = g.nv; + ne = g.ne; + spins.resize(nv); +#ifndef NOFIELD + g.add_ext(); + E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]); +#else + E = - (double)ne * J(spins[0], spins[0]); +#endif + M = spins[0] * nv; + last_cluster_size = 0; + ReF.resize(D); + ImF.resize(D); + for (D_t i = 0; i < D; i++) { + ReF[i] = spins[0] * 0.0; + ImF[i] = spins[0] * 0.0; + } + precomputed_cos.resize(nv); + precomputed_sin.resize(nv); + for (v_t i = 0; i < nv; i++) { + precomputed_cos[i].resize(D); + precomputed_sin[i].resize(D); + for (D_t j = 0; j < D; j++) { + precomputed_cos[i][j] = cos(2 * M_PI * g.coordinate[i][j] / (double)L); + precomputed_sin[i][j] = sin(2 * M_PI * g.coordinate[i][j] / (double)L); + } + } + } + + void update_magnetization(const X_t& s_old, const X_t& s_new) { + M += s_new - s_old; + } + + void update_energy(const double& dE) { + E += dE; + } + + void update_fourierZero(v_t v, const X_t& s_old, const X_t& s_new) { +#ifdef DIMENSION + for (D_t i = 0; i < DIMENSION; i++) { +#else + for (D_t i = 0; i < D; i++) { +#endif + ReF[i] += (s_new - s_old) * precomputed_cos[v][i]; + ImF[i] += (s_new - s_old) * precomputed_sin[v][i]; + } + } +}; + + + diff --git a/lib/include/wolff/types.h b/lib/include/wolff/types.h new file mode 100644 index 0000000..ec9efaf --- /dev/null +++ b/lib/include/wolff/types.h @@ -0,0 +1,36 @@ + +#include + +typedef uint_fast32_t v_t; +typedef uint_fast8_t q_t; +typedef uint_fast16_t R_t; +typedef uint_fast8_t D_t; +typedef uint_fast16_t L_t; +typedef uint_fast64_t count_t; +typedef int_fast64_t h_t; + +#define MAX_v UINT_FAST32_MAX +#define MAX_Q UINT_FAST8_MAX +#define MAX_R UINT_FAST16_MAX +#define MAX_D UINT_FAST8_MAX +#define MAX_L UINT_FAST16_MAX +#define MAX_COUNT UINT_FAST64_MAX +#define MAX_H INT_FAST64_MAX +#define MIN_H INT_FAST64_MIN + +#define PRIv PRIuFAST32 +#define PRIq PRIuFAST8 +#define PRIR PRIuFAST16 +#define PRID PRIuFAST8 +#define PRIL PRIuFAST16 +#define PRIcount PRIuFAST64 +#define PRIh PRIdFAST64 + +#define SCNv SCNuFAST32 +#define SCNq SCNuFAST8 +#define SCNR SCNuFAST16 +#define SCND SCNuFAST8 +#define SCNL SCNuFAST16 +#define SCNcount SCNuFAST64 +#define SCNh SCNdFAST64 + -- cgit v1.2.3-70-g09d2