From cea2351e7283099ebfd7d9d29688fe6c817bf4b8 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 16 Oct 2017 18:18:55 -0400 Subject: added new algorithm methods, specifically metropolis and wolff with a static ghost spin --- lib/wolff.h | 11 +++++-- lib/wolff_tools.c | 90 +++++++++++++++++++++++++++++++++++++++++++++++++------ src/wolff.c | 65 +++++++++++++++++++++++++++++++++++++--- 3 files changed, 151 insertions(+), 15 deletions(-) diff --git a/lib/wolff.h b/lib/wolff.h index 82ccb9e..bd49bf0 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -18,6 +18,12 @@ #include "queue.h" +typedef enum { + WOLFF, + WOLFF_GHOST, + METROPOLIS +} sim_t; + typedef struct { graph_t *g; bool *spins; @@ -29,6 +35,7 @@ typedef struct { uint32_t nv; int32_t dJb; int32_t dHb; + bool hit_ghost; } cluster_t; typedef struct { @@ -44,12 +51,12 @@ typedef struct { int32_t sign(double x); -cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, +cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_on_ghost, gsl_rng *r); graph_t *graph_add_ext(const graph_t *g); -uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, +uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r, double *ps); void update_meas(meas_t *m, double x); diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c index 59d04fc..bf78fda 100644 --- a/lib/wolff_tools.c +++ b/lib/wolff_tools.c @@ -2,6 +2,14 @@ #include "queue.h" #include "wolff.h" +int32_t spin_to_sign(bool spin) { + if (spin) { + return -1; + } else { + return 1; + } +} + int32_t sign(double x) { return x > 0 ? 1 : -1; } @@ -50,7 +58,7 @@ graph_t *graph_add_ext(const graph_t *g) { return h; } -cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, +cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_on_ghost, gsl_rng *r) { uint32_t v0; int32_t n_h_bonds, n_bonds; @@ -73,7 +81,7 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, v = stack_pop(&stack); nn = g->vei[v + 1] - g->vei[v]; - if (x[v] == x0) { // if the vertex hasn't already been flipped + if (x[v] == x0 && !(c->hit_ghost && stop_on_ghost)) { // if the vertex hasn't already been flipped x[v] = !x[v]; // flip the vertex for (uint32_t i = 0; i < nn; i++) { @@ -88,7 +96,7 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, vn = v == v1 ? v2 : v1; // distinguish neighboring site from site itself - is_ext = (v1 == g->nv - 1 || v2 == g->nv - 1); + is_ext = (v1 == g->nv - 1 || v2 == g->nv - 1); // our edge contained the "ghost" spin if either of its vertices was the last in the graph bond_counter = is_ext ? &(c->dHb) : &(c->dJb); prob = is_ext ? ps[1] : ps[0]; @@ -98,6 +106,9 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, (*bond_counter)++; if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]... + if (is_ext && stop_on_ghost) { + c->hit_ghost = true; + } stack_push(&stack, vn); // push the neighboring vertex to the stack } } else { @@ -114,8 +125,9 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, return c; } -uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, +uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r, double *ps) { + uint32_t n_flips; bool no_r, no_ps; no_r = false; no_ps = false; @@ -133,14 +145,74 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, ps[1] = 1 - exp(-2 * fabs(H) / T); } - cluster_t *c = flip_cluster(s->g, ps, s->spins, r); + switch (sim) { + case METROPOLIS: { + uint32_t v0 = gsl_rng_uniform_int(r, s->g->nv); // pick a random vertex + uint32_t nn = s->g->vei[v0 + 1] - s->g->vei[v0]; + + double dE = 0; + for (uint32_t i = 0; i < nn; i++) { + bool is_ext; + uint32_t e, v1, v2, vn; + int32_t *bond_counter; + double prob; + + e = s->g->ve[s->g->vei[v0] + i]; // select the ith bond connected to site + v1 = s->g->ev[2 * e]; + v2 = s->g->ev[2 * e + 1]; - s->M += - sign(H) * 2 * c->dHb; - s->H += 2 * (c->dJb + sign (H) * H * c->dHb); + vn = v0 == v1 ? v2 : v1; // distinguish neighboring site from site itself - uint32_t n_flips = c->nv; + is_ext = (v1 == s->g->nv - 1 || v2 == s->g->nv - 1); // our edge contained the "ghost" spin if either of its vertices was the last in the graph - free(c); + if (is_ext) { + dE += 2 * H * spin_to_sign(s->spins[vn]) * spin_to_sign(s->spins[v0]); + } else { + dE += 2 * spin_to_sign(s->spins[vn]) * spin_to_sign(s->spins[v0]); + } + } + + double p = exp(-dE / T); + if (gsl_rng_uniform(r) < p) { + s->M += - sign(H) * 2 * spin_to_sign(s->spins[v0]) * spin_to_sign(s->spins[s->g->nv - 1]); + s->H += dE; + s->spins[v0] = !s->spins[v0]; + } + + n_flips = 1; + } + break; + + case WOLFF: { + cluster_t *c = flip_cluster(s->g, ps, s->spins, false, r); + s->M += - sign(H) * 2 * c->dHb; + s->H += 2 * (c->dJb + sign (H) * H * c->dHb); + n_flips = c->nv; + + free(c); + } + break; + case WOLFF_GHOST: { + bool *spins_bak; + spins_bak = (bool *)malloc(s->g->nv * sizeof(bool)); + memcpy(spins_bak, s->spins, s->g->nv * sizeof(bool)); + + cluster_t *c = flip_cluster(s->g, ps, s->spins, true, r); + + if (c->hit_ghost) { + memcpy(s->spins, spins_bak, s->g->nv * sizeof(bool)); + } else { + s->M += - sign(H) * 2 * c->dHb; + s->H += 2 * (c->dJb + sign (H) * H * c->dHb); + } + free(spins_bak); + + n_flips = c->nv; + + free(c); + } + break; + } if (no_ps) { free(ps); diff --git a/src/wolff.c b/src/wolff.c index da57057..5091c0c 100644 --- a/src/wolff.c +++ b/src/wolff.c @@ -3,25 +3,28 @@ int main(int argc, char *argv[]) { int opt; - bool output_state, use_dual, M_stop; + sim_t sim; + bool output_state, use_dual, M_stop, record_autocorrelation; lattice_t lat; uint16_t L; - uint32_t min_runs, lattice_i; + uint32_t min_runs, lattice_i, sim_i; uint64_t N; double T, H, eps; + sim = WOLFF; L = 128; N = 1000; lat = SQUARE_LATTICE; use_dual = false; M_stop = false; + record_autocorrelation = false; T = 2.3; H = 0; eps = 1e30; output_state = false; min_runs = 10; - while ((opt = getopt(argc, argv, "N:L:T:H:m:e:oq:DM")) != -1) { + while ((opt = getopt(argc, argv, "N:L:T:H:m:e:oq:DMas:")) != -1) { switch (opt) { case 'N': N = (uint64_t)atof(optarg); @@ -50,6 +53,27 @@ int main(int argc, char *argv[]) { case 'D': use_dual = true; break; + case 'a': + record_autocorrelation = true; + break; + case 's': + sim_i = atoi(optarg); + switch (sim_i) { + case 0: + sim = WOLFF; + break; + case 1: + sim = WOLFF_GHOST; + break; + case 2: + sim = METROPOLIS; + break; + default: + printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 " + "(DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n"); + exit(EXIT_FAILURE); + } + break; case 'q': lattice_i = atoi(optarg); switch (lattice_i) { @@ -97,8 +121,19 @@ int main(int argc, char *argv[]) { double diff = 1e31; uint64_t n_runs = 0; + uint64_t n_steps = 0; double clust_per_sweep = 0; + uint64_t auto_n = 10 * h->nv; + uint64_t *auto_time; + double *auto_energy; + if (record_autocorrelation) { + auto_time = (uint64_t *)malloc(auto_n * sizeof(uint64_t)); + auto_energy = (double *)malloc(auto_n * sizeof(double)); + auto_time[0] = 0; + auto_energy[0] = s->H / h->nv; + } + meas_t *M, *aM, *eM, *mM, *E, *eE, *mE; M = calloc(1, sizeof(meas_t)); @@ -119,7 +154,7 @@ int main(int argc, char *argv[]) { uint32_t n_clust = 0; while (n_flips / h->nv < N) { - n_flips += wolff_step(T, H, s, r, bond_probs); + n_flips += wolff_step(T, H, s, sim, r, bond_probs); n_clust++; } @@ -149,6 +184,17 @@ int main(int argc, char *argv[]) { clust_per_sweep = add_to_avg(clust_per_sweep, n_clust * 1. / N, n_runs); n_runs++; + n_steps += n_flips; + + if (record_autocorrelation) { + if (n_runs == auto_n) { + auto_n *= 2; + auto_time = (uint64_t *)realloc(auto_time, auto_n * sizeof(uint64_t)); + auto_energy = (double *)realloc(auto_energy, auto_n * sizeof(double)); + } + auto_time[n_runs] = n_steps; + auto_energy[n_runs] = s->H / h->nv; + } } printf("\033[F\033[JWOLFF: sweep %" PRIu64 @@ -175,6 +221,17 @@ int main(int argc, char *argv[]) { fclose(state_file); } + if (record_autocorrelation) { + FILE *auto_file = fopen("autocorrelation.dat", "a"); + for (uint64_t i = 0; i < n_runs + 1; i++) { + fprintf(auto_file, "%" PRIu64 " %.15f ", auto_time[i], auto_energy[i]); + } + fprintf(auto_file, "\n"); + fclose(auto_file); + free(auto_time); + free(auto_energy); + } + gsl_rng_free(r); graph_free(s->g); free(s->spins); -- cgit v1.2.3-70-g09d2