From 5ffaf0a1bb0f0b47d57d0f24ee1134659775dacb Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 20 Jul 2018 22:57:39 -0400 Subject: added ising example to cpp collection --- src/wolff_On.cpp | 39 ++++++++++++- src/wolff_ising.cpp | 160 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 196 insertions(+), 3 deletions(-) create mode 100644 src/wolff_ising.cpp (limited to 'src') diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp index d718440..76831ba 100644 --- a/src/wolff_On.cpp +++ b/src/wolff_On.cpp @@ -40,6 +40,7 @@ int main(int argc, char *argv[]) { bool silent = false; bool use_pert = false; + bool N_is_sweeps = false; bool modulated_field = false; int order = 2; @@ -51,7 +52,7 @@ int main(int argc, char *argv[]) { unsigned char measurement_flags = measurement_energy | measurement_clusterSize; - while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:M:")) != -1) { + while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:M:S")) != -1) { switch (opt) { case 'N': // number of steps N = (count_t)atof(optarg); @@ -82,11 +83,14 @@ int main(int argc, char *argv[]) { modulated_field = true; break; case 'M': - measurement_flags |= 1 << atoi(optarg); + measurement_flags ^= 1 << atoi(optarg); break; case 'o': order = atoi(optarg); break; + case 'S': + N_is_sweeps = true; + break; default: exit(EXIT_FAILURE); } @@ -179,6 +183,13 @@ int main(int argc, char *argv[]) { n_measurements++; } + meas_t *meas_sweeps; + if (N_is_sweeps) { + meas_sweeps = (meas_t *)calloc(1, sizeof(meas_t)); + measurements[n_measurements] = measurement_average_cluster (meas_sweeps); + n_measurements++; + } + std::function H; if (modulated_field) { @@ -187,7 +198,24 @@ int main(int argc, char *argv[]) { H = std::bind(H_vector , std::placeholders::_1, H_vec); } - wolff (N, D, L, T, dot , H, gen_R, n_measurements, measurements, silent); + // initialize random number generator + gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); + gsl_rng_set(r, rand_seed()); + + state_t s(D, L, T, dot , H); + + if (N_is_sweeps) { + count_t N_rounds = 0; + printf("\n"); + while (N_rounds * N * meas_sweeps->x < N * s.nv) { + printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size); + wolff (N, &s, gen_R, n_measurements, measurements, r, silent); + N_rounds++; + } + printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size); + } else { + wolff (N, &s, gen_R, n_measurements, measurements, r, silent); + } free(measurements); @@ -204,7 +232,12 @@ int main(int argc, char *argv[]) { fclose(outfile_F); } + if (N_is_sweeps) { + free(meas_sweeps); + } + free(H_vec); + gsl_rng_free(r); return 0; } diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp new file mode 100644 index 0000000..4d4c791 --- /dev/null +++ b/src/wolff_ising.cpp @@ -0,0 +1,160 @@ + +#include + +#include +#include +#include + +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 = 0.0; + + bool silent = false; + bool N_is_sweeps = false; + + int opt; + + unsigned char measurement_flags = measurement_energy | measurement_clusterSize; + + while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:sM: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. nth call couples to state n + H = atof(optarg); + break; + case 's': // don't print anything during simulation. speeds up slightly + silent = true; + break; + case 'M': + measurement_flags ^= 1 << atoi(optarg); + break; + case 'S': + N_is_sweeps = true; + break; + default: + exit(EXIT_FAILURE); + } + } + + unsigned long timestamp; + + { + struct timespec spec; + clock_gettime(CLOCK_REALTIME, &spec); + timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec; + } + + 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, D, L, (v_t)pow(L, D), D * (v_t)pow(L, D), T, H); + + fclose(outfile_info); + + unsigned int n_measurements = 0; + std::function *)> *measurements = (std::function *)> *)calloc(POSSIBLE_MEASUREMENTS, sizeof(std::function *)>)); + FILE *outfile_M, *outfile_E, *outfile_S, *outfile_F; + + 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 (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 (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 (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); + measurements[n_measurements] = measurement_fourier_file (outfile_F); + n_measurements++; + } + + meas_t *meas_sweeps; + if (N_is_sweeps) { + meas_sweeps = (meas_t *)calloc(1, sizeof(meas_t)); + measurements[n_measurements] = measurement_average_cluster (meas_sweeps); + n_measurements++; + } + + // initialize random number generator + gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); + gsl_rng_set(r, rand_seed()); + + state_t s(D, L, T, ising_dot, std::bind(scalar_field, std::placeholders::_1, H)); + + std::function *)> gen = generate_ising_rotation; + + if (N_is_sweeps) { + count_t N_rounds = 0; + printf("\n"); + while (N_rounds * N * meas_sweeps->x < N * s.nv) { + printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size); + wolff(N, &s, gen, n_measurements, measurements, r, silent); + N_rounds++; + } + printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size); + } else { + wolff(N, &s, gen, n_measurements, measurements, r, 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); + } + + if (N_is_sweeps) { + free(meas_sweeps); + } + + gsl_rng_free(r); + + return 0; +} + -- cgit v1.2.3-70-g09d2