diff options
| -rw-r--r-- | CMakeLists.txt | 8 | ||||
| -rw-r--r-- | lib/correlation.h | 21 | ||||
| -rw-r--r-- | lib/measure.h | 48 | ||||
| -rw-r--r-- | lib/vector.h | 3 | ||||
| -rw-r--r-- | src/wolff_On.cpp | 219 | ||||
| -rw-r--r-- | src/wolff_heisenberg.cpp | 145 | ||||
| -rw-r--r-- | src/wolff_planar.cpp | 183 |
7 files changed, 279 insertions, 348 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 4f29504..98ef947 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,6 +5,7 @@ set(CMAKE_CXX_FLAGS_DEBUG "-g") set(CMAKE_CXX_FLAGS_RELEASE "-O3") set (CMAKE_CXX_STANDARD 17) +set (CMAKE_C_STANDARD 11) include_directories(lib ~/.local/include) link_directories(~/.local/lib) @@ -13,10 +14,13 @@ file(GLOB CSOURCES lib/*.c) file(GLOB CPPSOURCES lib/*.cpp) add_executable(wolff_finite src/wolff_finite.c ${CSOURCES}) -add_executable(wolff_heisenberg src/wolff_heisenberg.cpp ${CPPSOURCES} ${CSOURCES}) -add_executable(wolff_planar src/wolff_planar.cpp ${CPPSOURCES} ${CSOURCES}) +add_executable(wolff_planar src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES}) +add_executable(wolff_heisenberg src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES}) add_executable(analyze_correlations src/analyze_correlations.cpp ${CPPSOURCES} ${CSOURCES}) +SET_TARGET_PROPERTIES(wolff_planar PROPERTIES COMPILE_FLAGS "-DN_COMP=2") +SET_TARGET_PROPERTIES(wolff_heisenberg PROPERTIES COMPILE_FLAGS "-DN_COMP=3") + find_package(OpenMP) if (OPENMP_FOUND) set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") diff --git a/lib/correlation.h b/lib/correlation.h index 5722aab..b3b7d86 100644 --- a/lib/correlation.h +++ b/lib/correlation.h @@ -7,29 +7,14 @@ #include <fftw3.h> template <class R_t, class X_t> -double correlation_length(const state_t <R_t, X_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); - +double correlation_length(const state_t <R_t, X_t> *s, fftw_plan plan, double *in, double *out) { for (v_t i = 0; i < s->nv; i++) { - data[i] = correlation_component(s->spins[i]); + in[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); + double length = pow(out[0], 2); return length; } diff --git a/lib/measure.h b/lib/measure.h new file mode 100644 index 0000000..3474684 --- /dev/null +++ b/lib/measure.h @@ -0,0 +1,48 @@ + +#pragma once + +#define POSSIBLE_MEASUREMENTS 4 +const unsigned char measurement_energy = 1 << 0; +const unsigned char measurement_clusterSize = 1 << 1; +const unsigned char measurement_magnetization = 1 << 2; +const unsigned char measurement_fourierZero = 1 << 3; + +#ifdef __cplusplus + +#include "state.h" +#include "correlation.h" +#include <functional> + +template <class R_t, class X_t> +std::function <void(const state_t <R_t, X_t> *)> measurement_energy_file(FILE *file) { + return [=](const state_t <R_t, X_t> *s) { + float smaller_E = (float)s->E; + fwrite(&smaller_E, sizeof(float), 1, file); + }; +} + +template <class R_t, class X_t> +std::function <void(const state_t <R_t, X_t> *)> measurement_cluster_file(FILE *file) { + return [=](const state_t <R_t, X_t> *s) { + fwrite(&(s->last_cluster_size), sizeof(uint32_t), 1, file); + }; +} + +template <class R_t, class X_t> +std::function <void(const state_t <R_t, X_t> *)> measurement_magnetization_file(FILE *file) { + return [=](const state_t <R_t, X_t> *s) { + write_magnetization(s->M, file); + }; +} + +template <class R_t, class X_t> +std::function <void(const state_t <R_t, X_t> *)> measurement_fourier_file(FILE *file, fftw_plan plan, double *fftw_in, double *fftw_out) { + return [=](const state_t <R_t, X_t> *s) { + float smaller_X = (float)correlation_length(s, plan, fftw_in, fftw_out); + fwrite(&smaller_X, sizeof(float), 1, file); + }; +} + +#endif + + diff --git a/lib/vector.h b/lib/vector.h index 6c72fd1..2099d28 100644 --- a/lib/vector.h +++ b/lib/vector.h @@ -94,3 +94,6 @@ double H_vector(vector_t <q, T> v1, T *H) { H_vec.x = H; return (double)(dot <q, T> (v1, H_vec)); } + +char const *ON_strings[] = {"TRIVIAL", "ISING", "PLANAR", "HEISENBERG"}; + diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp new file mode 100644 index 0000000..336869f --- /dev/null +++ b/src/wolff_On.cpp @@ -0,0 +1,219 @@ + +#include <getopt.h> + +#include <wolff.h> +#include <correlation.h> +#include <measure.h> + +typedef state_t <orthogonal_t <N_COMP, double>, vector_t <N_COMP, double>> planar_t; + +// angle from the x-axis of a two-vector +double theta(vector_t <N_COMP, double> v) { + double x = v.x[0]; + double y = v.x[1]; + + double val = atan(y / x); + + if (x < 0.0 && y > 0.0) { + return M_PI + val; + } else if ( x < 0.0 && y < 0.0 ) { + return - M_PI + val; + } else { + return val; + } +} + +double H_modulated(vector_t <N_COMP, double> v, int order, double mag) { + return mag * cos(order * theta(v)); +} + +int main(int argc, char *argv[]) { + + count_t N = (count_t)1e7; + + D_t D = 2; + L_t L = 128; + double T = 2.26918531421; + double *H_vec = (double *)calloc(MAX_Q, sizeof(double)); + + bool silent = false; + bool use_pert = false; + + bool modulated_field = false; + int order = 2; + + int opt; + q_t J_ind = 0; + q_t H_ind = 0; + double epsilon = 1; + + unsigned char measurement_flags = 0; + + while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:M:")) != -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. 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 'p': + use_pert = true; + break; + case 'e': + epsilon = atof(optarg); + break; + case 'm': + modulated_field = true; + break; + case 'M': + measurement_flags |= 1 << atoi(optarg); + break; + case 'o': + order = atoi(optarg); + break; + default: + exit(EXIT_FAILURE); + } + } + + unsigned long timestamp; + + { + struct timespec spec; + clock_gettime(CLOCK_REALTIME, &spec); + timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; + } + + const char *pert_type; + + std::function <orthogonal_t <N_COMP, double>(gsl_rng *, const planar_t *)> gen_R; + + if (use_pert) { + gen_R = std::bind(generate_rotation_perturbation <N_COMP>, std::placeholders::_1, std::placeholders::_2, epsilon); + pert_type = "PERTURB"; + } else { + gen_R = generate_rotation_uniform <N_COMP>; + pert_type = "UNIFORM"; + } + + + FILE *outfile_info = fopen("wolff_metadata.txt", "a"); + + fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"%s\", \"q\" -> %d, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, ON_strings[N_COMP], N_COMP, D, L, L * L, D * L * L, T); + + for (q_t i = 0; i < N_COMP; i++) { + fprintf(outfile_info, "%.15f", H_vec[i]); + if (i < N_COMP - 1) { + fprintf(outfile_info, ", "); + } + } + + fprintf(outfile_info, "}, \"GENERATOR\" -> \"%s\", \"EPS\" -> %g |>\n", pert_type, epsilon); + + fclose(outfile_info); + + unsigned int n_measurements = 0; + std::function <void(const planar_t *)> *measurements = (std::function <void(const planar_t *)> *)calloc(POSSIBLE_MEASUREMENTS, sizeof(std::function <void(const planar_t *)>)); + FILE *outfile_M, *outfile_E, *outfile_S, *outfile_F; + double *fftw_in, *fftw_out; + fftw_plan plan; + + if (measurement_flags & measurement_energy) { + char *filename_E = (char *)malloc(255 * sizeof(char)); + sprintf(filename_E, "wolff_%lu_E.dat", timestamp); + outfile_E = fopen(filename_E, "wb"); + free(filename_E); + measurements[n_measurements] = measurement_energy_file<orthogonal_t <N_COMP, double>, vector_t <N_COMP, double>> (outfile_E); + n_measurements++; + } + + if (measurement_flags & measurement_clusterSize) { + char *filename_S = (char *)malloc(255 * sizeof(char)); + sprintf(filename_S, "wolff_%lu_S.dat", timestamp); + outfile_S = fopen(filename_S, "wb"); + free(filename_S); + measurements[n_measurements] = measurement_cluster_file<orthogonal_t <N_COMP, double>, vector_t <N_COMP, double>> (outfile_S); + n_measurements++; + } + + if (measurement_flags & measurement_magnetization) { + char *filename_M = (char *)malloc(255 * sizeof(char)); + sprintf(filename_M, "wolff_%lu_M.dat", timestamp); + outfile_M = fopen(filename_M, "wb"); + free(filename_M); + measurements[n_measurements] = measurement_magnetization_file<orthogonal_t <N_COMP, double>, vector_t <N_COMP, double>> (outfile_M); + n_measurements++; + } + + if (measurement_flags & measurement_fourierZero) { + char *filename_F = (char *)malloc(255 * sizeof(char)); + sprintf(filename_F, "wolff_%lu_F.dat", timestamp); + outfile_F = fopen(filename_F, "wb"); + free(filename_F); + + fftw_in = (double *)fftw_malloc(pow(L, D) * sizeof(double)); + fftw_out = (double *)fftw_malloc(pow(L, D) * sizeof(double)); + int rank = 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] = L; + kind[i] = FFTW_R2HC; + } + plan = fftw_plan_r2r(rank, n, fftw_in, fftw_out, kind, 0); + + free(n); + free(kind); + + measurements[n_measurements] = measurement_fourier_file<orthogonal_t <N_COMP, double>, vector_t <N_COMP, double>> (outfile_F, plan, fftw_in, fftw_out); + n_measurements++; + } + + std::function <double(vector_t <N_COMP, double>)> H; + + if (modulated_field) { + H = std::bind(H_modulated, std::placeholders::_1, order, H_vec[0]); + } else { + H = std::bind(H_vector <N_COMP, double>, std::placeholders::_1, H_vec); + } + + wolff <orthogonal_t <N_COMP, double>, vector_t <N_COMP, double>> (N, D, L, T, dot <N_COMP, double>, H, gen_R, n_measurements, measurements, silent); + + free(measurements); + + if (measurement_flags & measurement_energy) { + fclose(outfile_E); + } + if (measurement_flags & measurement_clusterSize) { + fclose(outfile_S); + } + if (measurement_flags & measurement_magnetization) { + fclose(outfile_M); + } + if (measurement_flags & measurement_fourierZero) { + fclose(outfile_F); + fftw_destroy_plan(plan); + fftw_free(fftw_in); + fftw_free(fftw_out); + fftw_cleanup(); // fftw is only used if fourier modes are measured! + } + + free(H_vec); + + return 0; +} + diff --git a/src/wolff_heisenberg.cpp b/src/wolff_heisenberg.cpp deleted file mode 100644 index e05453f..0000000 --- a/src/wolff_heisenberg.cpp +++ /dev/null @@ -1,145 +0,0 @@ - -#include <getopt.h> - -#include <correlation.h> -#include <wolff.h> - -typedef state_t <orthogonal_t <3, double>, vector_t <3, double>> heisenberg_t; - -int main(int argc, char *argv[]) { - - count_t N = (count_t)1e7; - - D_t D = 2; - L_t L = 128; - double T = 2.26918531421; - double *H = (double *)calloc(MAX_Q, sizeof(double)); - - bool silent = false; - bool use_pert = false; - - int opt; - q_t J_ind = 0; - q_t H_ind = 0; - double epsilon = 1; - - while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:")) != -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. nth call couples to state n - H[H_ind] = atof(optarg); - H_ind++; - break; - case 's': // don't print anything during simulation. speeds up slightly - silent = true; - break; - case 'p': - use_pert = true; - break; - case 'e': - epsilon = atof(optarg); - break; - default: - exit(EXIT_FAILURE); - } - } - - unsigned long timestamp; - - { - struct timespec spec; - clock_gettime(CLOCK_REALTIME, &spec); - timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; - } - - std::function <orthogonal_t <3, double>(gsl_rng *, const heisenberg_t *)> gen_R; - - const char *pert_type; - - if (use_pert) { - gen_R = std::bind(generate_rotation_perturbation <3>, std::placeholders::_1, std::placeholders::_2, epsilon); - pert_type = "PERTURB"; - } else { - gen_R = generate_rotation_uniform <3>; - pert_type = "UNIFORM"; - } - - FILE *outfile_info = fopen("wolff_metadata.txt", "a"); - - fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"HEISENBERG\", \"q\" -> 3, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, D, L, L * L, D * L * L, T); - - for (q_t i = 0; i < 3; i++) { - fprintf(outfile_info, "%.15f", H[i]); - if (i < 3 - 1) { - fprintf(outfile_info, ", "); - } - } - - fprintf(outfile_info, "}, \"GENERATOR\" -> \"%s\", \"EPS\" -> %g |>\n", pert_type, epsilon); - - fclose(outfile_info); - - char *filename_M = (char *)malloc(255 * sizeof(char)); - char *filename_E = (char *)malloc(255 * sizeof(char)); - char *filename_S = (char *)malloc(255 * sizeof(char)); - char *filename_X = (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); - sprintf(filename_X, "wolff_%lu_X.dat", timestamp); - - FILE *outfile_M = fopen(filename_M, "wb"); - FILE *outfile_E = fopen(filename_E, "wb"); - FILE *outfile_S = fopen(filename_S, "wb"); - FILE *outfile_X = fopen(filename_X, "wb"); - - free(filename_M); - free(filename_E); - free(filename_S); - free(filename_X); - - std::function <void(const heisenberg_t *)> *measurements = (std::function <void(const heisenberg_t *)> *)malloc(4 * sizeof(std::function <void(const heisenberg_t *)>)); - - measurements[0] = [&](const heisenberg_t *s) { - float smaller_E = (float)s->E; - fwrite(&smaller_E, sizeof(float), 1, outfile_E); - }; - measurements[1] = [&](const heisenberg_t *s) { - float smaller_X = (float)correlation_length(s); - fwrite(&smaller_X, sizeof(float), 1, outfile_X); - }; - measurements[2] = [&](const heisenberg_t *s) { - write_magnetization(s->M, outfile_M); - }; - measurements[3] = [&](const heisenberg_t *s) { - fwrite(&(s->last_cluster_size), sizeof(uint32_t), 1, outfile_S); - }; - - wolff <orthogonal_t <3, double>, vector_t <3, double>> (N, D, L, T, dot <3, double>, std::bind(H_vector <3, double>, std::placeholders::_1, H), gen_R, 4, measurements, silent); - - free(measurements); - - fclose(outfile_M); - fclose(outfile_E); - fclose(outfile_S); - fclose(outfile_X); - - free(H); - - fftw_cleanup(); - - return 0; -} - diff --git a/src/wolff_planar.cpp b/src/wolff_planar.cpp deleted file mode 100644 index 4b9b5f0..0000000 --- a/src/wolff_planar.cpp +++ /dev/null @@ -1,183 +0,0 @@ - -#include <getopt.h> - -#include <wolff.h> -#include <correlation.h> - -typedef state_t <orthogonal_t <2, double>, vector_t <2, double>> planar_t; - -// angle from the x-axis of a two-vector -double theta(vector_t <2, double> v) { - double x = v.x[0]; - double y = v.x[1]; - - double val = atan(y / x); - - if (x < 0.0 && y > 0.0) { - return M_PI + val; - } else if ( x < 0.0 && y < 0.0 ) { - return - M_PI + val; - } else { - return val; - } -} - -double H_modulated(vector_t <2, double> v, int order, double mag) { - return mag * cos(order * theta(v)); -} - -int main(int argc, char *argv[]) { - - count_t N = (count_t)1e7; - - D_t D = 2; - L_t L = 128; - double T = 2.26918531421; - double *H_vec = (double *)calloc(MAX_Q, sizeof(double)); - - bool silent = false; - bool use_pert = false; - - bool modulated_field = false; - int order = 2; - - int opt; - q_t J_ind = 0; - q_t H_ind = 0; - double epsilon = 1; - - while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:")) != -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. 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 'p': - use_pert = true; - break; - case 'e': - epsilon = atof(optarg); - break; - case 'm': - modulated_field = true; - break; - case 'o': - order = atoi(optarg); - break; - default: - exit(EXIT_FAILURE); - } - } - - unsigned long timestamp; - - { - struct timespec spec; - clock_gettime(CLOCK_REALTIME, &spec); - timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; - } - - const char *pert_type; - - std::function <orthogonal_t <2, double>(gsl_rng *, const planar_t *)> gen_R; - - if (use_pert) { - gen_R = std::bind(generate_rotation_perturbation <2>, std::placeholders::_1, std::placeholders::_2, epsilon); - pert_type = "PERTURB"; - } else { - gen_R = generate_rotation_uniform <2>; - pert_type = "UNIFORM"; - } - - - FILE *outfile_info = fopen("wolff_metadata.txt", "a"); - - fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"PLANAR\", \"q\" -> 2, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, D, L, L * L, D * L * L, T); - - for (q_t i = 0; i < 2; i++) { - fprintf(outfile_info, "%.15f", H_vec[i]); - if (i < 2 - 1) { - fprintf(outfile_info, ", "); - } - } - - fprintf(outfile_info, "}, \"GENERATOR\" -> \"%s\", \"EPS\" -> %g |>\n", pert_type, epsilon); - - fclose(outfile_info); - - char *filename_M = (char *)malloc(255 * sizeof(char)); - char *filename_E = (char *)malloc(255 * sizeof(char)); - char *filename_S = (char *)malloc(255 * sizeof(char)); - char *filename_X = (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); - sprintf(filename_X, "wolff_%lu_X.dat", timestamp); - - FILE *outfile_M = fopen(filename_M, "wb"); - FILE *outfile_E = fopen(filename_E, "wb"); - FILE *outfile_S = fopen(filename_S, "wb"); - FILE *outfile_X = fopen(filename_X, "wb"); - - free(filename_M); - free(filename_E); - free(filename_S); - free(filename_X); - - std::function <void(const planar_t *)> *measurements = (std::function <void(const planar_t *)> *)calloc(4, sizeof(std::function <void(const planar_t *)>)); - - measurements[0] = (std::function <void(const planar_t *)>)[&](const planar_t *s) { - float smaller_E = (float)s->E; - fwrite(&smaller_E, sizeof(float), 1, outfile_E); - }; - measurements[1] = [&](const planar_t *s) { - float smaller_X = (float)correlation_length(s); - fwrite(&smaller_X, sizeof(float), 1, outfile_X); - }; - measurements[2] = [&](const planar_t *s) { - write_magnetization(s->M, outfile_M); - }; - measurements[3] = [&](const planar_t *s) { - fwrite(&(s->last_cluster_size), sizeof(uint32_t), 1, outfile_S); - }; - - std::function <double(vector_t <2, double>)> H; - - if (modulated_field) { - H = std::bind(H_modulated, std::placeholders::_1, order, H_vec[0]); - } else { - H = std::bind(H_vector <2, double>, std::placeholders::_1, H_vec); - } - - wolff <orthogonal_t <2, double>, vector_t <2, double>> (N, D, L, T, dot <2, double>, H, gen_R, 4, measurements, silent); - - free(measurements); - - fclose(outfile_M); - fclose(outfile_E); - fclose(outfile_S); - fclose(outfile_X); - - free(H_vec); - - fftw_cleanup(); - - return 0; -} - |