From cd85d3696be4a7ee65b0cf6669fc62dc4841aef3 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 18 Jul 2018 15:37:27 -0400 Subject: - added support for computing spatial fourier transforms - measurements now are functions passed to wolff in array --- lib/cluster.h | 8 +++---- lib/correlation.h | 36 ++++++++++++++++++++++++++++++++ lib/height.h | 62 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ lib/orthogonal.h | 3 ++- lib/state.h | 2 ++ lib/vector.h | 15 +++++++++----- lib/wolff.h | 43 +++++++------------------------------- 7 files changed, 123 insertions(+), 46 deletions(-) create mode 100644 lib/correlation.h create mode 100644 lib/height.h (limited to 'lib') diff --git a/lib/cluster.h b/lib/cluster.h index 5f770ad..8061d34 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -52,7 +52,7 @@ template T subtract(T, T); template -v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { +void flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { v_t nv = 0; ll_t *stack = NULL; // create a new stack @@ -104,8 +104,8 @@ v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { double dE = state->H(rs_old) - state->H(rs_new); prob = 1.0 - exp(-dE / state->T); - subtract (state->M, rs_old); - add (state->M, rs_new); + subtract (&(state->M), rs_old); + add (&(state->M), rs_new); state->E += dE; free_spin (rs_old); @@ -137,6 +137,6 @@ v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { free(marks); - return nv; + state->last_cluster_size = nv; } diff --git a/lib/correlation.h b/lib/correlation.h new file mode 100644 index 0000000..5722aab --- /dev/null +++ b/lib/correlation.h @@ -0,0 +1,36 @@ + +#pragma once + +#include "types.h" +#include "state.h" + +#include + +template +double correlation_length(const state_t *s) { + double *data = (double *)fftw_malloc(s->nv * sizeof(double)); + int rank = s->D; + int *n = (int *)malloc(rank * sizeof(int)); + fftw_r2r_kind *kind = (fftw_r2r_kind *)malloc(rank * sizeof(fftw_r2r_kind)); + for (D_t i = 0; i < rank; i++) { + n[i] = s->L; + kind[i] = FFTW_R2HC; + } + fftw_plan plan = fftw_plan_r2r(rank, n, data, data, kind, 0); + + for (v_t i = 0; i < s->nv; i++) { + data[i] = correlation_component(s->spins[i]); + } + + fftw_execute(plan); + + double length = pow(data[0], 2); + + fftw_destroy_plan(plan); + fftw_free(data); + free(n); + free(kind); + + return length; +} + diff --git a/lib/height.h b/lib/height.h new file mode 100644 index 0000000..16c4dc5 --- /dev/null +++ b/lib/height.h @@ -0,0 +1,62 @@ + +#pragma once + +#include + +#include "types.h" + +// object definition +template +struct height_t { T x; }; + +// init, copy, add, subtract, scalar_multiple, free_spin, and +// write_magnetization are necessary for the operation of wolff.h +template +void init(height_t *ptr) { + ptr->x = (T)0; +} + +template +height_t copy (height_t h) { + return h; +} + +template +void add (height_t *h1, height_t h2) { + h1->x += h2.x; +} + +template +void subtract (height_t *h1, height_T h2) { + h1->x -= h2.x; +} + +template +height_t scalar_multiple(v_t a, height_t h) { + height_t hm; + hm.x = a * h.x; + + return hm; +} + +template +void free_spin (height_t h) { +} + +template +void write_magnetization(height_t M, FILE *outfile) { + fwrite(&(M.x), sizeof(T), 1, outfile); +} + +template +double correlation_component(height_t h) { + return (double)h.x; +} + +// below here are not necessary for operation + +template +T dot(height_t h1, height_t h2) { + return (h1.x - h2.x) * (h1.x - h2.x); +} + diff --git a/lib/orthogonal.h b/lib/orthogonal.h index d8ad33d..340ee2c 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -184,9 +184,10 @@ orthogonal_t generate_rotation_perturbation (gsl_rng *r, const state } double v2 = 0; + double factor = gsl_ran_ugaussian(r); for (q_t i = 0; i < q; i++) { - tmp_v.x[i] = (tmp_v.x[i] - tmpM * s->M.x[i] / M2) + epsilon * gsl_ran_ugaussian(r); + tmp_v.x[i] = (tmp_v.x[i] - tmpM * s->M.x[i] / M2) + epsilon * factor * s->M.x[i] / sqrt(M2); v2 += tmp_v.x[i] * tmp_v.x[i]; } diff --git a/lib/state.h b/lib/state.h index 8d4265b..1491938 100644 --- a/lib/state.h +++ b/lib/state.h @@ -19,6 +19,7 @@ class state_t { R_t R; double E; X_t M; // the "sum" of the spins, like the total magnetization + v_t last_cluster_size; std::function J; std::function H; @@ -36,6 +37,7 @@ class state_t { init (&R); E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]); M = scalar_multiple (nv, spins[0]); + last_cluster_size = 0; } ~state_t() { diff --git a/lib/vector.h b/lib/vector.h index 62ce59e..88934c1 100644 --- a/lib/vector.h +++ b/lib/vector.h @@ -30,16 +30,16 @@ vector_t copy (vector_t v) { } template -void add (vector_t v1, vector_t v2) { +void add (vector_t *v1, vector_t v2) { for (q_t i = 0; i < q; i++) { - v1.x[i] += v2.x[i]; + v1->x[i] += v2.x[i]; } } template -void subtract (vector_t v1, vector_t v2) { +void subtract (vector_t *v1, vector_t v2) { for (q_t i = 0; i < q; i++) { - v1.x[i] -= v2.x[i]; + v1->x[i] -= v2.x[i]; } } @@ -72,7 +72,7 @@ void free_spin (vector_t v) { template void write_magnetization(vector_t M, FILE *outfile) { - fwrite(M.x, sizeof(double), q, outfile); + fwrite(M.x, sizeof(T), q, outfile); } template // save some space and don't write whole doubles @@ -83,3 +83,8 @@ void write_magnetization(vector_t M, FILE *outfile) { } } +template +double correlation_component(vector_t v) { + return (double)v.x[0]; +} + diff --git a/lib/wolff.h b/lib/wolff.h index 90470d5..2d7e8b2 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -1,7 +1,4 @@ -#include -#include - #include "cluster.h" #include "state.h" @@ -13,7 +10,7 @@ double H_vector(vector_t v1, T *H) { } template -void wolff(count_t N, D_t D, L_t L, double T, std::function J, std::function H, std::function *)> gen_R, unsigned long timestamp, bool silent) { +void wolff(count_t N, D_t D, L_t L, double T, std::function J, std::function H, std::function *)> gen_R, unsigned int n_measurements, std::function *)> *measurements, bool silent) { state_t s(D, L, T, J, H); @@ -21,52 +18,26 @@ void wolff(count_t N, D_t D, L_t L, double T, std::function J gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); gsl_rng_set(r, rand_seed()); - char *filename_M = (char *)malloc(255 * sizeof(char)); - char *filename_E = (char *)malloc(255 * sizeof(char)); - char *filename_S = (char *)malloc(255 * sizeof(char)); - - sprintf(filename_M, "wolff_%lu_M.dat", timestamp); - sprintf(filename_E, "wolff_%lu_E.dat", timestamp); - sprintf(filename_S, "wolff_%lu_S.dat", timestamp); - - FILE *outfile_M = fopen(filename_M, "wb"); - FILE *outfile_E = fopen(filename_E, "wb"); - FILE *outfile_S = fopen(filename_S, "wb"); - - free(filename_M); - free(filename_E); - free(filename_S); - - v_t cluster_size = 0; - if (!silent) printf("\n"); for (count_t steps = 0; steps < N; steps++) { - if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s.E, cluster_size); + if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s.E, s.last_cluster_size); v_t v0 = gsl_rng_uniform_int(r, s.nv); - R_t step = gen_R(r, &s); - - cluster_size = flip_cluster (&s, v0, step, r); - + flip_cluster (&s, v0, step, r); free_spin(step); - { - float smaller_E = (float)s.E; - fwrite(&smaller_E, sizeof(float), 1, outfile_E); + for (unsigned int i = 0; i < n_measurements; i++) { + measurements[i](&s); } - write_magnetization(s.M, outfile_M); - fwrite(&cluster_size, sizeof(uint32_t), 1, outfile_S); } if (!silent) { printf("\033[F\033[J"); } - printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, cluster_size); + printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, s.last_cluster_size); - fclose(outfile_M); - fclose(outfile_E); - fclose(outfile_S); + fftw_cleanup(); gsl_rng_free(r); } -- cgit v1.2.3-54-g00ecf