diff options
| -rw-r--r-- | CMakeLists.txt | 5 | ||||
| -rw-r--r-- | lib/cluster.h | 68 | ||||
| -rw-r--r-- | lib/correlation.h | 4 | ||||
| -rw-r--r-- | lib/graph.cpp | 7 | ||||
| -rw-r--r-- | lib/graph.h | 1 | ||||
| -rw-r--r-- | lib/ising.h | 12 | ||||
| -rw-r--r-- | lib/state.h | 41 | ||||
| -rw-r--r-- | lib/vector.h | 6 | ||||
| -rw-r--r-- | src/wolff_On.cpp | 10 |
9 files changed, 105 insertions, 49 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index fe564f7..568585f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -23,6 +23,7 @@ add_executable(wolff_7potts src/wolff_potts.cpp ${CPPSOURCES} ${CSOURCES}) add_executable(wolff_3clock src/wolff_clock.cpp ${CPPSOURCES} ${CSOURCES}) add_executable(wolff_5clock src/wolff_clock.cpp ${CPPSOURCES} ${CSOURCES}) add_executable(wolff_planar src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES}) +add_executable(wolff_planar_2D src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES}) add_executable(wolff_heisenberg src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES}) add_executable(analyze_correlations src/analyze_correlations.cpp ${CPPSOURCES} ${CSOURCES}) @@ -32,6 +33,7 @@ SET_TARGET_PROPERTIES(wolff_7potts PROPERTIES COMPILE_FLAGS "-DPOTTSQ=7") SET_TARGET_PROPERTIES(wolff_3clock PROPERTIES COMPILE_FLAGS "-DPOTTSQ=3") SET_TARGET_PROPERTIES(wolff_5clock PROPERTIES COMPILE_FLAGS "-DPOTTSQ=5") SET_TARGET_PROPERTIES(wolff_planar PROPERTIES COMPILE_FLAGS "-DN_COMP=2") +SET_TARGET_PROPERTIES(wolff_planar_2D PROPERTIES COMPILE_FLAGS "-DN_COMP=2 -DDIMENSION=2") SET_TARGET_PROPERTIES(wolff_heisenberg PROPERTIES COMPILE_FLAGS "-DN_COMP=3") find_library(GSL REQUIRED NAMES gsl) @@ -54,6 +56,7 @@ if (${GLUT} MATCHES "GLUT-NOTFOUND") target_link_libraries(wolff_5clock cblas gsl m) target_link_libraries(wolff_heisenberg cblas gsl m) target_link_libraries(wolff_planar cblas gsl m) + target_link_libraries(wolff_planar_2D cblas gsl m) else() target_link_libraries(wolff_ising cblas gsl m glut GL GLU) target_link_libraries(wolff_dgm cblas gsl m glut GL GLU) @@ -65,6 +68,7 @@ else() target_link_libraries(wolff_5clock cblas gsl m glut GL GLU) target_link_libraries(wolff_heisenberg cblas gsl m glut GL GLU) target_link_libraries(wolff_planar cblas gsl m glut GL GLU) + target_link_libraries(wolff_planar_2D cblas gsl m glut GL GLU) target_compile_definitions(wolff_ising PUBLIC HAVE_GLUT) target_compile_definitions(wolff_dgm PUBLIC HAVE_GLUT) target_compile_definitions(wolff_cgm PUBLIC HAVE_GLUT) @@ -74,6 +78,7 @@ else() target_compile_definitions(wolff_3clock PUBLIC HAVE_GLUT) target_compile_definitions(wolff_5clock PUBLIC HAVE_GLUT) target_compile_definitions(wolff_planar PUBLIC HAVE_GLUT) + target_compile_definitions(wolff_planar_2D PUBLIC HAVE_GLUT) target_compile_definitions(wolff_heisenberg PUBLIC HAVE_GLUT) endif() diff --git a/lib/cluster.h b/lib/cluster.h index f3271b0..9dcf50d 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -25,78 +25,64 @@ void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, gsl_rng *rand) { stack.pop(); if (!marks[v]) { // don't reprocess anyone we've already visited! - X_t si_old, si_new; - R_t R_old, R_new; + X_t si_new; + R_t R_new; - R_old = s.R; marks[v] = true; bool v_is_ghost = (v == s.nv); // ghost site has the last index if (v_is_ghost) { - R_new = r.act(R_old); // if we are, then we're moving the transformation + R_new = r.act(s.R); // if we are, then we're moving the transformation } else { - si_old = s.spins[v]; - si_new = r.act(si_old); // otherwise, we're moving the spin at our site + 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]) { - X_t sj; bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost - if (!vn_is_ghost) { - sj = s.spins[vn]; - } - - double prob; + double dE, prob; - if (v_is_ghost || vn_is_ghost) { // if this is a ghost-involved bond... + 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) { - rs_old = R_old.act_inverse(si_old); - rs_new = R_old.act_inverse(si_new); + // 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 { - rs_old = R_old.act_inverse(sj); - rs_new = R_new.act_inverse(sj); + /* 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; } - double dE = s.H(rs_old) - s.H(rs_new); + 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)]; -#else - prob = 1.0 - exp(-dE / s.T); #endif - s.M -= rs_old; - s.M += rs_new; - - s.E += dE; - - for (D_t i = 0; i < s.D; i++) { - L_t x = (non_ghost / (v_t)pow(s.L, s.D - i - 1)) % s.L; - - s.ReF[i] -= rs_old * s.precomputed_cos[x]; - s.ReF[i] += rs_new * s.precomputed_cos[x]; - - s.ImF[i] -= rs_old * s.precomputed_sin[x]; - s.ImF[i] += rs_new * s.precomputed_sin[x]; - } - } else { // otherwise, we're at a perfectly normal bond! - double dE = s.J(si_old, sj) - s.J(si_new, sj); + s.update_magnetization(rs_old, rs_new); + s.update_fourierZero(non_ghost, rs_old, rs_new); + } else { // this is a perfectly normal bond! + 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(si_old)][state_to_ind(si_new)][state_to_ind(sj)]; -#else - prob = 1.0 - exp(-dE / s.T); + prob = J_probs[state_to_ind(s.spins[v])][state_to_ind(si_new)][state_to_ind(s.spins[vn])]; #endif - - s.E += dE; } + s.update_energy(dE); + +#ifndef FINITE_STATES + prob = 1.0 - exp(-dE / s.T); +#endif + if (gsl_rng_uniform(rand) < prob) { stack.push(vn); // push the neighboring vertex to the stack } diff --git a/lib/correlation.h b/lib/correlation.h index d08f713..29357a5 100644 --- a/lib/correlation.h +++ b/lib/correlation.h @@ -10,7 +10,11 @@ template <class R_t, class X_t> double correlation_length(const state_t <R_t, X_t>& s) { double total = 0; +#ifdef DIMENSION + for (D_t j = 0; j < DIMENSION; j++) { +#else for (D_t j = 0; j < s.D; j++) { +#endif total += norm_squared(s.ReF[j]) + norm_squared(s.ImF[j]); } diff --git a/lib/graph.cpp b/lib/graph.cpp index 8c97274..ca251f3 100644 --- a/lib/graph.cpp +++ b/lib/graph.cpp @@ -6,13 +6,17 @@ graph_t::graph_t(D_t D, L_t L) { ne = D * nv; v_adj.resize(nv); + coordinate.resize(nv); for (std::vector<v_t> v_adj_i : v_adj) { v_adj_i.reserve(2 * D); } for (v_t i = 0; i < nv; i++) { + coordinate[i].resize(D); for (D_t j = 0; j < D; j++) { + coordinate[i][j] = (i / (v_t)pow(L, D - j - 1)) % L; + v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1))); v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1))); } @@ -25,12 +29,15 @@ void graph_t::add_ext() { } v_adj.resize(nv + 1); + coordinate.resize(nv + 1); v_adj[nv].reserve(nv); for (v_t i = 0; i < nv; i++) { v_adj[nv].push_back(i); } + coordinate[nv].resize(coordinate[0].size()); + ne += nv; nv += 1; } diff --git a/lib/graph.h b/lib/graph.h index a4732fb..de06924 100644 --- a/lib/graph.h +++ b/lib/graph.h @@ -13,6 +13,7 @@ class graph_t { v_t ne; v_t nv; std::vector<std::vector<v_t>> v_adj; + std::vector<std::vector<double>> coordinate; graph_t(D_t D, L_t L); void add_ext(); diff --git a/lib/ising.h b/lib/ising.h index 473a18a..aaa0c52 100644 --- a/lib/ising.h +++ b/lib/ising.h @@ -53,6 +53,18 @@ class ising_t { return a; } } + + inline int operator-(const ising_t &s) const { + if (x == s.x) { + return 0; + } else { + if (x) { + return -2; + } else { + return 2; + } + } + } }; inline int& operator+=(int& M, const ising_t &s) { diff --git a/lib/state.h b/lib/state.h index 5abf65b..550d100 100644 --- a/lib/state.h +++ b/lib/state.h @@ -9,6 +9,10 @@ template <class R_t, class X_t> class state_t { + private: + // updating fourier terms F requires many cos and sin calls, faster to do it beforehand. + std::vector<std::vector<double>> precomputed_cos; + std::vector<std::vector<double>> precomputed_sin; public: D_t D; L_t L; @@ -23,9 +27,6 @@ class state_t { v_t last_cluster_size; std::vector<typename X_t::F_t> ReF; std::vector<typename X_t::F_t> ImF; - // updating fourier terms F requires many cos and sin calls, faster to do it beforehand. - std::vector<double> precomputed_cos; - std::vector<double> precomputed_sin; std::function <double(const X_t&, const X_t&)> J; std::function <double(const X_t&)> H; @@ -44,11 +45,35 @@ class state_t { ReF[i] = spins[0] * 0.0; ImF[i] = spins[0] * 0.0; } - precomputed_cos.resize(L); - precomputed_sin.resize(L); - for (L_t i = 0; i < L; i++) { - precomputed_cos[i] = cos(2 * M_PI * (double)i / (double)L); - precomputed_sin[i] = sin(2 * M_PI * (double)i / (double)L); + 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_old; + M += s_new; + } + + 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/vector.h b/lib/vector.h index 72633a8..7d0ee36 100644 --- a/lib/vector.h +++ b/lib/vector.h @@ -61,6 +61,12 @@ class vector_t : public std::array<T, q> { return result; } + + inline vector_t<q, T> operator-(const vector_t<q, T>& v) const { + vector_t<q, T> diff = *this; + diff -= v; + return diff; + } }; diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp index 3fa5840..6cd8777 100644 --- a/src/wolff_On.cpp +++ b/src/wolff_On.cpp @@ -42,7 +42,11 @@ int main(int argc, char *argv[]) { count_t N = (count_t)1e7; +#ifdef DIMENSION + D_t D = DIMENSION; +#else D_t D = 2; +#endif L_t L = 128; double T = 2.26918531421; double *H_vec = (double *)calloc(MAX_Q, sizeof(double)); @@ -69,9 +73,15 @@ int main(int argc, char *argv[]) { case 'N': // number of steps N = (count_t)atof(optarg); break; +#ifdef DIMENSION + case 'D': // dimension + printf("Dimension was specified at compile time, you can't change it now!\n"); + exit(EXIT_FAILURE); +#else case 'D': // dimension D = atoi(optarg); break; +#endif case 'L': // linear size L = atoi(optarg); break; |