diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/wolff.h | 11 | ||||
| -rw-r--r-- | lib/wolff_tools.c | 90 |
2 files changed, 90 insertions, 11 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); |