From 6e8b19e1f1a244ef09e1b63d7593250d6ce01692 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sun, 14 Oct 2018 20:54:01 -0400 Subject: added site and bond dependence, activated by compiler flags --- examples/src/models/On/wolff_On.cpp | 4 + examples/src/models/ising/CMakeLists.txt | 6 +- examples/src/models/ising/wolff_ising.cpp | 4 + .../src/models/ising/wolff_random-field_ising.cpp | 207 +++++++++++++++++++++ lib/include/wolff/cluster.hpp | 9 + lib/include/wolff/state.hpp | 76 ++++++-- 6 files changed, 291 insertions(+), 15 deletions(-) create mode 100644 examples/src/models/ising/wolff_random-field_ising.cpp diff --git a/examples/src/models/On/wolff_On.cpp b/examples/src/models/On/wolff_On.cpp index e3568c7..ad2ac77 100644 --- a/examples/src/models/On/wolff_On.cpp +++ b/examples/src/models/On/wolff_On.cpp @@ -212,7 +212,11 @@ int main(int argc, char *argv[]) { other_f = [&] (const On_t& s) { glClear(GL_COLOR_BUFFER_BIT); for (v_t i = 0; i < pow(L, 2); i++) { +#ifdef NOFIELD + vector_R_t v_tmp = s.spins[i]; +#else vector_R_t v_tmp = s.R.act_inverse(s.spins[i]); +#endif double thetai = fmod(2 * M_PI + theta(v_tmp), 2 * M_PI); double saturation = 0.7; double value = 0.9; diff --git a/examples/src/models/ising/CMakeLists.txt b/examples/src/models/ising/CMakeLists.txt index 495c7e5..9f4acd4 100644 --- a/examples/src/models/ising/CMakeLists.txt +++ b/examples/src/models/ising/CMakeLists.txt @@ -2,6 +2,7 @@ add_executable(wolff_ising wolff_ising.cpp) add_executable(wolff_ising_2d wolff_ising.cpp) add_executable(wolff_ising_2d_no-field wolff_ising.cpp) +add_executable(wolff_random-field_ising wolff_random-field_ising.cpp) set_target_properties(wolff_ising_2d PROPERTIES COMPILE_FLAGS "-DDIMENSION=2") set_target_properties(wolff_ising_2d_no-field PROPERTIES COMPILE_FLAGS "-DDIMENSION=2 -DNOFIELD") @@ -14,15 +15,18 @@ if (${GLUT} MATCHES "GLUT-NOTFOUND") target_link_libraries(wolff_ising wolff wolff_examples) target_link_libraries(wolff_ising_2d wolff wolff_examples) target_link_libraries(wolff_ising_2d_no-field wolff wolff_examples) + target_link_libraries(wolff_random-field_ising wolff wolff_examples) else() target_compile_definitions(wolff_ising PUBLIC HAVE_GLUT) target_compile_definitions(wolff_ising_2d PUBLIC HAVE_GLUT) target_compile_definitions(wolff_ising_2d_no-field PUBLIC HAVE_GLUT) + target_compile_definitions(wolff_random-field_ising PUBLIC HAVE_GLUT) target_link_libraries(wolff_ising wolff wolff_examples glut GL GLU) target_link_libraries(wolff_ising_2d wolff wolff_examples glut GL GLU) target_link_libraries(wolff_ising_2d_no-field wolff wolff_examples glut GL GLU) + target_link_libraries(wolff_random-field_ising wolff wolff_examples glut GL GLU) endif() -install(TARGETS wolff_ising wolff_ising_2d wolff_ising_2d_no-field DESTINATION ${CMAKE_INSTALL_BINDIR} OPTIONAL) +install(TARGETS wolff_ising wolff_ising_2d wolff_ising_2d_no-field wolff_random-field_ising DESTINATION ${CMAKE_INSTALL_BINDIR} OPTIONAL) diff --git a/examples/src/models/ising/wolff_ising.cpp b/examples/src/models/ising/wolff_ising.cpp index 5bdaa82..f7a3ada 100644 --- a/examples/src/models/ising/wolff_ising.cpp +++ b/examples/src/models/ising/wolff_ising.cpp @@ -151,7 +151,11 @@ int main(int argc, char *argv[]) { other_f = [] (const state_t & s) { glClear(GL_COLOR_BUFFER_BIT); for (v_t i = 0; i < pow(s.L, 2); i++) { +#ifdef NOFIELD + if (s.spins[i].x == false) { +#else if (s.spins[i].x == s.R.x) { +#endif glColor3f(0.0, 0.0, 0.0); } else { glColor3f(1.0, 1.0, 1.0); diff --git a/examples/src/models/ising/wolff_random-field_ising.cpp b/examples/src/models/ising/wolff_random-field_ising.cpp new file mode 100644 index 0000000..ce26b88 --- /dev/null +++ b/examples/src/models/ising/wolff_random-field_ising.cpp @@ -0,0 +1,207 @@ + +#define SITE_DEPENDENCE + +#include +#include + +// if you have GLUT installed, you can see graphics! +#ifdef HAVE_GLUT +#include +#endif + +// include your group and spin space +#include "z2.hpp" +#include "ising.hpp" + +#include + +// measure.hpp contains useful functions for saving timeseries to files +#include + +// include wolff.hpp +#include + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e4; + + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double H = 0.0; + + bool silent = false; + bool draw = false; + bool N_is_sweeps = false; + unsigned int window_size = 512; + + // don't measure anything by default + unsigned char measurement_flags = 0; + + int opt; + + while ((opt = getopt(argc, argv, "N:D:L:T:H:sdw:M:S")) != -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 '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; + case 'M': + measurement_flags ^= 1 << atoi(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + // get nanosecond timestamp for unique run id + unsigned long timestamp; + + { + struct timespec spec; + clock_gettime(CLOCK_REALTIME, &spec); + timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; + } + + // initialize random number generator + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + // define spin-spin coupling + std::function Z = [] (const ising_t& s1, const ising_t& s2) -> double { + if (s1.x == s2.x) { + return 1.0; + } else { + return -1.0; + } + }; + + // create random field + std::vector random_field_values(pow(L, D)); + std::normal_distribution distribution(0.0, H); + for (v_t i = 0; i < pow(L, D); i++) { + random_field_values[i] = distribution(rng); + } + + // define spin-field coupling + std::function B = [&] (v_t v, const ising_t& s) -> double { + if (s.x) { + return -random_field_values[v]; + } else { + return random_field_values[v]; + } + }; + + // initialize state object +#ifndef NOFIELD + state_t s(D, L, T, Z, B); +#else + state_t s(D, L, T, Z); +#endif + + // define function that generates self-inverse rotations + std::function gen_R = [] (std::mt19937&, const ising_t& s) -> z2_t { + return z2_t(true); + }; + + FILE **outfiles = measure_setup_files(measurement_flags, timestamp); + + std::function &)> other_f; + uint64_t sum_of_clusterSize = 0; + + if (N_is_sweeps) { + other_f = [&] (const state_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 state_t & s) { + glClear(GL_COLOR_BUFFER_BIT); + for (v_t i = 0; i < pow(s.L, 2); i++) { +#ifdef NOFIELD + if (s.spins[i].x == false) { +#else + if (s.spins[i].x == s.R.x) { +#endif + glColor3f(0.0, 0.0, 0.0); + } else { + glColor3f(1.0, 1.0, 1.0); + } + glRecti(i / s.L, i % s.L, (i / s.L) + 1, (i % s.L) + 1); + } + glFlush(); + }; +#endif + } else { + other_f = [] (const state_t& s) {}; + } + + std::function &)> measurements = measure_function_write_files(measurement_flags, outfiles, other_f); + + // add line to metadata file with run info + { + FILE *outfile_info = fopen("wolff_metadata.txt", "a"); + + fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"ISING\", \"q\" -> 2, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> %.15f |>\n", timestamp, s.D, s.L, s.nv, s.ne, T, H); + + fclose(outfile_info); + } + + // run wolff for N cluster flips + if (N_is_sweeps) { + count_t N_rounds = 0; + printf("\n"); + while (sum_of_clusterSize < N * s.nv) { + printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size); + wolff(N, s, gen_R, measurements, rng, silent); + N_rounds++; + } + printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)((double)sum_of_clusterSize / (double)s.nv), N, s.E, s.last_cluster_size); + } else { + wolff(N, s, gen_R, measurements, rng, silent); + } + + measure_free_files(measurement_flags, outfiles); + + return 0; + +} + diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp index 104f3c2..e9dff7b 100644 --- a/lib/include/wolff/cluster.hpp +++ b/lib/include/wolff/cluster.hpp @@ -65,7 +65,11 @@ void flip_cluster(state_t& s, v_t v0, const R_t& r, std::mt19937& rand non_ghost = vn; } +#ifdef SITE_DEPENDENCE + dE = s.H(non_ghost, rs_old) - s.H(non_ghost, rs_new); +#else dE = s.H(rs_old) - s.H(rs_new); +#endif #ifdef FINITE_STATES prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)]; @@ -76,7 +80,12 @@ void flip_cluster(state_t& s, v_t v0, const R_t& r, std::mt19937& rand } else // this is a perfectly normal bond! #endif { +#ifdef BOND_DEPENDENCE + dE = s.J(v, s.spins[v], vn, s.spins[vn]) - s.J(v, si_new, vn, s.spins[vn]); +#else dE = s.J(s.spins[v], s.spins[vn]) - s.J(si_new, s.spins[vn]); +#endif + #ifdef FINITE_STATES prob = J_probs[state_to_ind(s.spins[v])][state_to_ind(si_new)][state_to_ind(s.spins[vn])]; diff --git a/lib/include/wolff/state.hpp b/lib/include/wolff/state.hpp index 1f5e359..e7c0ac3 100644 --- a/lib/include/wolff/state.hpp +++ b/lib/include/wolff/state.hpp @@ -16,35 +16,83 @@ class state_t { 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; + v_t nv; // the number of vertices in the lattice + v_t ne; // the number of edges in the lattice + graph_t g; // the graph defining the lattice without ghost + double T; // the temperature + std::vector spins; // the state of the ordinary spins +#ifndef NOFIELD + R_t R; // the current state of the ghost site +#endif + double E; // the system's total energy typename X_t::M_t M; // the "sum" of the spins, like the total magnetization - v_t last_cluster_size; + v_t last_cluster_size; // the size of the last cluster std::vector ReF; std::vector ImF; - std::function J; +#ifdef BOND_DEPENDENCE + std::function J; // coupling between sites +#else + std::function J; // coupling between sites +#endif + #ifndef NOFIELD - std::function H; +#ifdef SITE_DEPENDENCE + std::function H; // coupling with the external field +#else + std::function H; // coupling with the external field +#endif +#endif - 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) { + state_t(D_t D, L_t L, double T, +#ifdef BOND_DEPENDENCE + std::function J +#else + std::function J +#endif +#ifndef NOFIELD +#ifdef SITE_DEPENDENCE + , std::function 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) { + , std::function H +#endif +#endif + ) : D(D), L(L), g(D, L), T(T), +#ifndef NOFIELD + R(), #endif + J(J) +#ifndef NOFIELD + , H(H) +#endif + { nv = g.nv; ne = g.ne; spins.resize(nv); +#ifdef BOND_DEPENDENCE + E = 0; + for (v_t v = 0; v < nv; v++) { + for (const v_t &vn : g.v_adj[v]) { + if (v < vn) { + E -= J(v, spins[v], vn, spins[vn]); + } + } + } +#else + E = - (double)ne * J(spins[0], spins[0]); +#endif + #ifndef NOFIELD g.add_ext(); - E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]); +#ifdef SITE_DEPENDENCE + for (v_t i = 0; i < nv; i++) { + E -= H(i, spins[i]); + } #else - E = - (double)ne * J(spins[0], spins[0]); + E -= (double)nv * H(spins[0]); +#endif #endif + M = spins[0] * nv; last_cluster_size = 0; ReF.resize(D); -- cgit v1.2.3-70-g09d2