From dd2c47db3512658858685c83dd772603203aaab1 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sun, 22 Jul 2018 01:20:29 -0400 Subject: potts now fully functional --- src/wolff_On.cpp | 33 ++------------------- src/wolff_potts.cpp | 82 +++++++++++++++++++++++++++++++++-------------------- 2 files changed, 54 insertions(+), 61 deletions(-) (limited to 'src') diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp index 32a6fc4..a59876f 100644 --- a/src/wolff_On.cpp +++ b/src/wolff_On.cpp @@ -1,5 +1,6 @@ #include + #ifdef HAVE_GLUT #include #endif @@ -7,6 +8,7 @@ #include #include #include +#include typedef orthogonal_t orthogonal_R_t; typedef vector_t vector_R_t; @@ -32,36 +34,6 @@ 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; @@ -81,7 +53,6 @@ int main(int argc, char *argv[]) { int order = 2; int opt; - q_t J_ind = 0; q_t H_ind = 0; double epsilon = 1; diff --git a/src/wolff_potts.cpp b/src/wolff_potts.cpp index 9d22ea4..3b55472 100644 --- a/src/wolff_potts.cpp +++ b/src/wolff_potts.cpp @@ -1,29 +1,37 @@ #include + +#ifdef HAVE_GLUT #include +#endif // include your group and spin space #include #include +#include // include wolff.h #include +typedef state_t , potts_t> sim_t; + 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; double T = 2.26918531421; - double H = 0.0; + double *H_vec = (double *)calloc(MAX_Q, sizeof(double)); bool silent = false; bool draw = false; + unsigned int window_size = 512; int opt; + q_t H_ind = 0; - 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); @@ -37,15 +45,24 @@ int main(int argc, char *argv[]) { case 'T': // temperature T = atof(optarg); break; - case 'H': // external field - H = atof(optarg); + case 'H': // external field. nth call couples to state n + H_vec[H_ind] = atof(optarg); + H_ind++; break; case 's': // don't print anything during simulation. speeds up slightly 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); } @@ -56,41 +73,49 @@ int main(int argc, char *argv[]) { gsl_rng_set(r, rand_seed()); // define spin-spin coupling - std::function Z = [] (ising_t s1, ising_t s2) -> double { + std::function , potts_t)> Z = [] (potts_t s1, potts_t s2) -> double { if (s1.x == s2.x) { return 1.0; } else { - return -1.0; + return 0.0; } }; // define spin-field coupling - std::function B = [=] (ising_t s) -> double { - if (s.x) { - return -H; - } else { - return H; - } + std::function )> B = [=] (potts_t s) -> double { + return H_vec[s.x]; }; // initialize state object - state_t s(D, L, T, Z, B); + state_t , potts_t> s(D, L, T, Z, B); // define function that generates self-inverse rotations - std::function *)> gen_R = [] (gsl_rng *, const state_t *) -> z2_t { - z2_t rot; - rot.x = true; + std::function (gsl_rng *, const sim_t *)> gen_R = [] (gsl_rng *r, const sim_t *s) -> symmetric_t { + symmetric_t rot; + init(&rot); + + for (int i = POTTSQ - 1; i >= 0; i--) { + if (rot.perm[i] == i) { + q_t j = gsl_rng_uniform_int(r, i + 1); + if (rot.perm[j] == j) { + q_t tmp = rot.perm[i]; + rot.perm[i] = rot.perm[j]; + rot.perm[j] = tmp; + } + } + } + return rot; }; // define function that updates any number of measurements - std::function *)> measurement; + std::function measurement; double average_M = 0; if (!draw) { // a very simple example: measure the average magnetization - measurement = [&] (const state_t *s) { - average_M += (double)s->M / (double)N / (double)s->nv; + measurement = [&] (const sim_t *s) { + average_M += (double)s->M[0] / (double)N / (double)s->nv; }; } else { // a more complex example: measure the average magnetization, and draw the spin configuration to the screen @@ -98,22 +123,19 @@ int main(int argc, char *argv[]) { // 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(); gluOrtho2D(0.0, L, 0.0, L); - measurement = [&] (const state_t *s) { - average_M += (double)s->M / (double)N / (double)s->nv; + measurement = [&] (const sim_t *s) { + average_M += (double)s->M[0] / (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); - } + potts_t tmp_s = act_inverse(s->R, s->spins[i]); + glColor3f(hue_to_R(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_G(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_B(tmp_s.x * 2 * M_PI / POTTSQ)); glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1); } glFlush(); @@ -124,7 +146,7 @@ int main(int argc, char *argv[]) { 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, = %g\n", N, D, L, T, H, average_M); + printf("%" PRIcount " %d-Potts runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, = %g\n", N, POTTSQ, D, L, T, H_vec[0], average_M); // free the random number generator gsl_rng_free(r); -- cgit v1.2.3-70-g09d2