diff options
| -rw-r--r-- | CMakeLists.txt | 31 | ||||
| -rw-r--r-- | lib/ising.h | 8 | ||||
| -rw-r--r-- | lib/potts.h | 100 | ||||
| -rw-r--r-- | lib/state.h | 4 | ||||
| -rw-r--r-- | lib/symmetric.h | 88 | ||||
| -rw-r--r-- | lib/vector.h | 11 | ||||
| -rw-r--r-- | src/wolff_On.cpp | 75 | ||||
| -rw-r--r-- | src/wolff_ising.cpp | 21 | ||||
| -rw-r--r-- | src/wolff_potts.cpp | 138 |
9 files changed, 459 insertions, 17 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index b225d39..1bec5c6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,7 @@ cmake_minimum_required(VERSION 3.0) project(wolff) + set(CMAKE_CXX_FLAGS_DEBUG "-g") set(CMAKE_CXX_FLAGS_RELEASE "-O3") @@ -22,17 +23,27 @@ add_executable(analyze_correlations src/analyze_correlations.cpp ${CPPSOURCES} $ SET_TARGET_PROPERTIES(wolff_planar PROPERTIES COMPILE_FLAGS "-DN_COMP=2") SET_TARGET_PROPERTIES(wolff_heisenberg PROPERTIES COMPILE_FLAGS "-DN_COMP=3") -find_package(OpenMP) -if (OPENMP_FOUND) - set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") - set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") +find_library(GSL gsl) +find_library(FFTW fftw3) +find_library(M m) +FIND_LIBRARY(GL GL) +FIND_LIBRARY(GLU GLU) +FIND_LIBRARY(GLUT glut) + +target_link_libraries(wolff_finite cblas gsl m) +target_link_libraries(analyze_correlations cblas gsl fftw3 m) +if (GL_FOUND AND GLU_FOUND AND GLUT_FOUND) + target_link_libraries(wolff_ising 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_compile_definitions(wolff_ising HAVE_GLUT) + target_compile_definitions(wolff_planar HAVE_GLUT) + target_compile_definitions(wolff_heisenberg HAVE_GLUT) +else () + target_link_libraries(wolff_ising cblas gsl m) + target_link_libraries(wolff_heisenberg cblas gsl m) + target_link_libraries(wolff_planar cblas gsl m) endif() -target_link_libraries(wolff_finite gsl cblas fftw3 m) -target_link_libraries(wolff_ising gsl cblas fftw3 m glut GL GLU) -target_link_libraries(wolff_heisenberg gsl cblas fftw3 m) -target_link_libraries(wolff_planar gsl cblas fftw3 m) -target_link_libraries(analyze_correlations gsl cblas fftw3 m) - install(TARGETS wolff_finite wolff_ising wolff_heisenberg wolff_planar analyze_correlations DESTINATION bin) diff --git a/lib/ising.h b/lib/ising.h index 4ad88f4..b4856c3 100644 --- a/lib/ising.h +++ b/lib/ising.h @@ -76,6 +76,14 @@ int scalar_multiple(int factor, ising_t s) { } } +double scalar_multiple(double factor, ising_t s) { + if (s.x) { + return -factor; + } else { + return factor; + } +} + double norm_squared(double s) { return pow(s, 2); } diff --git a/lib/potts.h b/lib/potts.h new file mode 100644 index 0000000..e7f0899 --- /dev/null +++ b/lib/potts.h @@ -0,0 +1,100 @@ +#pragma once + +#include <cmath> +#include <stdio.h> + +#include "types.h" + +/* The following is the minimum definition of a spin class. + * + * The class must contain an M_t and an F_t for holding the sum of an + * integer number of spins and a double-weighted number of spins, + * respectively. + * + * void init(X_t *p); + * void free_spin(X_t p); + * void free_spin(M_t p); + * void free_spin(F_t p); + * X_t copy(X_t x); + * void add(M_t *x1, int factor, X_t x2); + * void add(F_t *x1, double factor, X_t x2); + * M_t scalar_multiple(int factor, X_t x); + * F_t scalar_multiple(double factor, X_t x); + * double norm_squared(F_t x); + * void write_magnetization(M_t M, FILE *outfile); + * + */ + +template <q_t q> +class potts_t { + public: + q_t x; + + typedef int *M_t; + typedef double *F_t; +}; + +template <q_t q> +void init(potts_t <q> *p) { + p->x = 0; +} + +template <q_t q> +void free_spin(potts_t <q> s) { + // do nothing! +} + +template <q_t q> +void free_spin(typename potts_t<q>::M_t s) { + free(s); +} + +template <q_t q> +void free_spin(typename potts_t<q>::F_t s) { + free(s); +} + +template <q_t q> +potts_t <q> copy(potts_t <q> s) { + return s; +} + +template <q_t q> +void add(typename potts_t<q>::M_t s1, int a, potts_t <q> s2) { + s1[s2.x] += a; +} + +template <q_t q> +void add(typename potts_t<q>::F_t s1, double a, potts_t <q> s2) { + s1[s2.x] += a; +} + +template <q_t q> +typename potts_t<q>::M_t scalar_multiple(int factor, potts_t <q> s) { + int *M = (int *)calloc(q, sizeof(int)); + M[s.x] += factor; + return M; +} + +template <q_t q> +typename potts_t<q>::F_t scalar_multiple(double factor, potts_t <q> s) { + int *F = (double *)calloc(q, sizeof(double)); + M[s.x] += factor; + return M; +} + +template <q_t q> +double norm_squared(typename potts<q>::F_t s) { + double total = 0; + for (q_t i = 0; i < q; i++) { + total += pow(s[i], 2); + } + + return total * (double)q / ((double)q - 1.0); +} + +template <q_t q> +void write_magnetization(typename potts_t<q>::M_t M, FILE *outfile) { + fwrite(&M, sizeof(int), q, outfile); +} + diff --git a/lib/state.h b/lib/state.h index 76d3e5a..8630810 100644 --- a/lib/state.h +++ b/lib/state.h @@ -47,8 +47,8 @@ class state_t { ReF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t)); ImF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t)); for (D_t i = 0; i < D; i++) { - ReF[i] = scalar_multiple(0, spins[0]); - ImF[i] = scalar_multiple(0, spins[0]); + ReF[i] = scalar_multiple(0.0, spins[0]); + ImF[i] = scalar_multiple(0.0, spins[0]); } precomputed_cos = (double *)malloc(L * sizeof(double)); precomputed_sin = (double *)malloc(L * sizeof(double)); diff --git a/lib/symmetric.h b/lib/symmetric.h index c71521d..0b292a6 100644 --- a/lib/symmetric.h +++ b/lib/symmetric.h @@ -5,6 +5,10 @@ #include "types.h" +#ifdef __cplusplus +extern "C" { +#endif + q_t *symmetric_compose(q_t q, const q_t *g1, const q_t *g2); q_t symmetric_act(const q_t *g, q_t s); @@ -13,3 +17,87 @@ q_t *symmetric_invert(q_t q, const q_t *g); q_t *symmetric_gen_transformations(q_t q); +#ifdef __cplusplus +} +#endif + +#ifdef __cplusplus +template <q_t q> +class symmetric_t { + public: + q_t *perm; +}; + +template <q_t q> +void init(symmetric_t<q> *p) { + p->perm = (q_t *)malloc(q * sizeof(q_t)); + + for (q_t i = 0; i < q; i++) { + p->perm[i] = i; + } +} + +template <q_t q> +void free_spin(symmetric_t<q> p) { + free(p->perm); +} + +template <q_t q> +symmetric_t<q_t> copy(symmetric_t<q_t> x) { + symmetric_t<q> x2; + x2.perm = (q_t *)malloc(q * sizeof(q_t)); + + for (q_t i = 0; i < q; i++) { + x2.perm[i] = x.perm[i]; + } + + return x2; +} + +template <q_t q> +potts_t<q> act(symmetric_t<q> r, potts_t<q> s) { + potts_t<q> s2; + s2.x = r.perm[s.x]; + return s2; +} + +template <q_t q> +symmetric_t<q> act(symmetric_t<q> r1, symmetric_t<q> r2) { + symmetric_t<q> r3; + r3.perm = (q_t *)malloc(q * sizeof(q_t)); + for (q_t i = 0; i < q; i++) { + r3.perm[i] = r1.perm[r2.perm[i]]; + } + + return r3; +} + +template <q_t q> +potts_t<q> act_inverse(symmetric_t<q> r, potts_t<q> s) { + potts_t<q> s2; + + q_t i; + + for (i = 0; i < q; i++) { + if (r.perm[i] == s.x) { + break; + } + } + + s2.x = i; + + return s2; +} + +template <q_t q> +symmetric_t<q> act_inverse(symmetric_t<q> r1, symmetric_t<q> r2) { + symmetric_t<q> r3; + r3.perm = (q_t *)malloc(q * sizeof(q_t)); + for (q_t i = 0; i < q; i++) { + r3.perm[r1.perm[i]] = r2.perm[i]; + } + + return r3; +} +#endif + diff --git a/lib/vector.h b/lib/vector.h index c478618..2fe6ab8 100644 --- a/lib/vector.h +++ b/lib/vector.h @@ -80,6 +80,17 @@ vector_t <q, T> scalar_multiple(int a, vector_t <q, T> v) { } template <q_t q, class T> +vector_t <q, T> scalar_multiple(double a, vector_t <q, T> v) { + vector_t <q, T> multiple; + multiple.x = (T *)malloc(q * sizeof(T)); + for (q_t i = 0; i < q; i++) { + multiple.x[i] = a * v.x[i]; + } + + return multiple; +} + +template <q_t q, class T> double norm_squared (vector_t <q, T> v) { double tmp = 0; diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp index 491fb27..32a6fc4 100644 --- a/src/wolff_On.cpp +++ b/src/wolff_On.cpp @@ -1,5 +1,8 @@ #include <getopt.h> +#ifdef HAVE_GLUT +#include <GL/glut.h> +#endif #include <wolff.h> #include <correlation.h> @@ -29,6 +32,36 @@ double H_modulated(vector_R_t v, int order, double mag) { return mag * cos(order * theta(v)); } +double hue_to_R(double theta) { + if (((M_PI / 3 <= theta) && (theta < 2 * M_PI / 3)) || ((4 * M_PI / 3 <= theta) && (theta < 5 * M_PI / 3))) { + return 1.0 - fabs(fmod(theta / (2 * M_PI / 6), 2) - 1.0); + } else if (((0 <= theta) && (theta < M_PI / 3)) || ((5 * M_PI / 3 <= theta) && (theta <= 2 * M_PI))) { + return 1.0; + } else { + return 0.0; + } +} + +double hue_to_G(double theta) { + if (((0 <= theta) && (theta < M_PI / 3)) || ((M_PI <= theta) && (theta < 4 * M_PI / 3))) { + return 1.0 - fabs(fmod(theta / (2 * M_PI / 6), 2) - 1.0); + } else if (((M_PI / 3 <= theta) && (theta < 2 * M_PI / 3)) || ((2 * M_PI / 3 <= theta) && (theta < M_PI))) { + return 1.0; + } else { + return 0.0; + } +} + +double hue_to_B(double theta) { + if (((2 * M_PI / 3 <= theta) && (theta < M_PI)) || ((5 * M_PI / 3 <= theta) && (theta <= 2 * M_PI))) { + return 1.0 - fabs(fmod(theta / (2 * M_PI / 6), 2) - 1.0); + } else if (((M_PI <= theta) && (theta < 4 * M_PI / 3)) || ((4 * M_PI / 3 <= theta) && (theta < 5 * M_PI / 3))) { + return 1.0; + } else { + return 0.0; + } +} + int main(int argc, char *argv[]) { count_t N = (count_t)1e7; @@ -41,6 +74,8 @@ int main(int argc, char *argv[]) { bool silent = false; bool use_pert = false; bool N_is_sweeps = false; + bool draw = false; + unsigned int window_size = 512; bool modulated_field = false; int order = 2; @@ -54,7 +89,7 @@ int main(int argc, char *argv[]) { unsigned char measurement_flags = 0; - while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:M:S")) != -1) { + while ((opt = getopt(argc, argv, "N:D:L:T:H:spe:mo:M:Sdw:")) != -1) { switch (opt) { case 'N': // number of steps N = (count_t)atof(optarg); @@ -93,6 +128,17 @@ int main(int argc, char *argv[]) { case 'S': N_is_sweeps = true; break; + case 'd': +#ifdef HAVE_GLUT + draw = true; + break; +#else + printf("You didn't compile this with the glut library installed!\n"); + exit(EXIT_FAILURE); +#endif + case 'w': + window_size = atoi(optarg); + break; default: exit(EXIT_FAILURE); } @@ -153,6 +199,33 @@ int main(int argc, char *argv[]) { other_f = [&] (const On_t *s) { sum_of_clusterSize += s->last_cluster_size; }; + } else if (draw) { +#ifdef HAVE_GLUT + // initialize glut + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(window_size, window_size); + glutCreateWindow("wolff"); + glClearColor(0.0,0.0,0.0,0.0); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(0.0, L, 0.0, L); + + other_f = [&] (const On_t *s) { + glClear(GL_COLOR_BUFFER_BIT); + for (v_t i = 0; i < pow(L, 2); i++) { + vector_R_t v_tmp = act_inverse(s->R, s->spins[i]); + double thetai = fmod(2 * M_PI + theta(v_tmp), 2 * M_PI); + free_spin(v_tmp); + double saturation = 0.7; + double value = 0.9; + double chroma = saturation * value; + glColor3f(chroma * hue_to_R(thetai) + (value - chroma), chroma * hue_to_G(thetai) + (value - chroma), chroma * hue_to_B(thetai) + (value - chroma)); + glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); + } + glFlush(); + }; +#endif } else { other_f = [] (const On_t *s) {}; } diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp index f4073da..83b6448 100644 --- a/src/wolff_ising.cpp +++ b/src/wolff_ising.cpp @@ -1,6 +1,8 @@ #include <getopt.h> +#ifdef HAVE_GLUT #include <GL/glut.h> +#endif // include your group and spin space #include <z2.h> @@ -11,7 +13,7 @@ int main(int argc, char *argv[]) { - count_t N = (count_t)1e7; + count_t N = (count_t)1e4; D_t D = 2; L_t L = 128; @@ -20,10 +22,11 @@ int main(int argc, char *argv[]) { bool silent = false; bool draw = false; + unsigned int window_size = 512; int opt; - while ((opt = getopt(argc, argv, "N:D:L:T:H:sd")) != -1) { + while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:")) != -1) { switch (opt) { case 'N': // number of steps N = (count_t)atof(optarg); @@ -44,8 +47,16 @@ int main(int argc, char *argv[]) { silent = true; break; case 'd': +#ifdef HAVE_GLUT draw = true; break; +#else + printf("You didn't compile this with the glut library installed!\n"); + exit(EXIT_FAILURE); +#endif + case 'w': + window_size = atoi(optarg); + break; default: exit(EXIT_FAILURE); } @@ -95,11 +106,12 @@ int main(int argc, char *argv[]) { } else { // a more complex example: measure the average magnetization, and draw the spin configuration to the screen +#ifdef HAVE_GLUT // initialize glut glutInit(&argc, argv); glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); - glutInitWindowSize(L,L); - glutCreateWindow("null"); + glutInitWindowSize(window_size, window_size); + glutCreateWindow("wolff"); glClearColor(0.0,0.0,0.0,0.0); glMatrixMode(GL_PROJECTION); glLoadIdentity(); @@ -118,6 +130,7 @@ int main(int argc, char *argv[]) { } glFlush(); }; +#endif } // run wolff for N cluster flips diff --git a/src/wolff_potts.cpp b/src/wolff_potts.cpp new file mode 100644 index 0000000..9d22ea4 --- /dev/null +++ b/src/wolff_potts.cpp @@ -0,0 +1,138 @@ + +#include <getopt.h> +#include <GL/glut.h> + +// include your group and spin space +#include <symmetric.h> +#include <potts.h> + +// include wolff.h +#include <wolff.h> + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e7; + + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double H = 0.0; + + bool silent = false; + bool draw = false; + + int opt; + + while ((opt = getopt(argc, argv, "N:D:L:T:H:sd")) != -1) { + switch (opt) { + case 'N': // number of steps + N = (count_t)atof(optarg); + break; + case 'D': // dimension + D = atoi(optarg); + break; + case 'L': // linear size + L = atoi(optarg); + break; + case 'T': // temperature + T = atof(optarg); + break; + case 'H': // external field + H = atof(optarg); + break; + case 's': // don't print anything during simulation. speeds up slightly + silent = true; + break; + case 'd': + draw = true; + break; + default: + exit(EXIT_FAILURE); + } + } + + // initialize random number generator + gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); + gsl_rng_set(r, rand_seed()); + + // define spin-spin coupling + std::function <double(potts_t, ising_t)> Z = [] (ising_t s1, ising_t s2) -> double { + if (s1.x == s2.x) { + return 1.0; + } else { + return -1.0; + } + }; + + // define spin-field coupling + std::function <double(ising_t)> B = [=] (ising_t s) -> double { + if (s.x) { + return -H; + } else { + return H; + } + }; + + // initialize state object + state_t <z2_t, ising_t> s(D, L, T, Z, B); + + // define function that generates self-inverse rotations + std::function <z2_t(gsl_rng *, const state_t <z2_t, ising_t> *)> gen_R = [] (gsl_rng *, const state_t <z2_t, ising_t> *) -> z2_t { + z2_t rot; + rot.x = true; + return rot; + }; + + // define function that updates any number of measurements + std::function <void(const state_t <z2_t, ising_t> *)> measurement; + + double average_M = 0; + if (!draw) { + // a very simple example: measure the average magnetization + measurement = [&] (const state_t <z2_t, ising_t> *s) { + average_M += (double)s->M / (double)N / (double)s->nv; + }; + } else { + // a more complex example: measure the average magnetization, and draw the spin configuration to the screen + + // initialize glut + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(L,L); + glutCreateWindow("null"); + glClearColor(0.0,0.0,0.0,0.0); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(0.0, L, 0.0, L); + + measurement = [&] (const state_t <z2_t, ising_t> *s) { + average_M += (double)s->M / (double)N / (double)s->nv; + glClear(GL_COLOR_BUFFER_BIT); + for (v_t i = 0; i < pow(L, 2); i++) { + if (s->spins[i].x == s->R.x) { + glColor3f(0.0, 0.0, 0.0); + } else { + glColor3f(1.0, 1.0, 1.0); + } + glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); + } + glFlush(); + }; + } + + // run wolff for N cluster flips + wolff(N, &s, gen_R, measurement, r, silent); + + // tell us what we found! + printf("%" PRIcount " Ising runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, <M> = %g\n", N, D, L, T, H, average_M); + + // free the random number generator + gsl_rng_free(r); + + if (draw) { + } + + return 0; + +} + |