From 14b533f9d8b58ada158dac15469da035fa2d8159 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 17 Oct 2017 17:56:09 -0400 Subject: many changes, including ability to compute approximate correlation time --- lib/wolff_tools.c | 30 ++++++++++--------- src/wolff.c | 86 +++++++++++++++++++++++++++---------------------------- 2 files changed, 59 insertions(+), 57 deletions(-) diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c index bf78fda..5e88a57 100644 --- a/lib/wolff_tools.c +++ b/lib/wolff_tools.c @@ -58,6 +58,16 @@ graph_t *graph_add_ext(const graph_t *g) { return h; } +uint32_t get_neighbor(const graph_t *g, uint32_t v, uint32_t i) { + uint32_t e, v1, v2; + + e = g->ve[g->vei[v] + i]; // select the ith bond connected to site + v1 = g->ev[2 * e]; + v2 = g->ev[2 * e + 1]; + + return v == v1 ? v2 : v1; // distinguish neighboring site from site itself +} + cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_on_ghost, gsl_rng *r) { uint32_t v0; @@ -86,17 +96,13 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_o for (uint32_t i = 0; i < nn; i++) { bool is_ext; - uint32_t e, v1, v2, vn; + uint32_t vn; int32_t *bond_counter; double prob; - e = g->ve[g->vei[v] + i]; // select the ith bond connected to site - v1 = g->ev[2 * e]; - v2 = g->ev[2 * e + 1]; - - vn = v == v1 ? v2 : v1; // distinguish neighboring site from site itself + vn = get_neighbor(g, v, i); - 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 + is_ext = (v == g->nv - 1 || vn == 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]; @@ -153,17 +159,13 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r, double dE = 0; for (uint32_t i = 0; i < nn; i++) { bool is_ext; - uint32_t e, v1, v2, vn; + uint32_t 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]; - - vn = v0 == v1 ? v2 : v1; // distinguish neighboring site from site itself + vn = get_neighbor(s->g, v0, i); - 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 + is_ext = (v0 == s->g->nv - 1 || vn == s->g->nv - 1); // our edge contained the "ghost" spin if either of its vertices was the last in the graph if (is_ext) { dE += 2 * H * spin_to_sign(s->spins[vn]) * spin_to_sign(s->spins[v0]); diff --git a/src/wolff.c b/src/wolff.c index 5091c0c..9082c81 100644 --- a/src/wolff.c +++ b/src/wolff.c @@ -8,27 +8,31 @@ int main(int argc, char *argv[]) { lattice_t lat; uint16_t L; uint32_t min_runs, lattice_i, sim_i; - uint64_t N; + uint64_t N, n; double T, H, eps; sim = WOLFF; L = 128; - N = 1000; + N = 100000000000000; + n = 3; lat = SQUARE_LATTICE; use_dual = false; M_stop = false; record_autocorrelation = false; T = 2.3; H = 0; - eps = 1e30; + eps = 1e-30; output_state = false; min_runs = 10; - while ((opt = getopt(argc, argv, "N:L:T:H:m:e:oq:DMas:")) != -1) { + while ((opt = getopt(argc, argv, "n:N:L:T:H:m:e:oq:DMas:")) != -1) { switch (opt) { case 'N': N = (uint64_t)atof(optarg); break; + case 'n': + n = (uint64_t)atof(optarg); + break; case 'L': L = atoi(optarg); break; @@ -124,17 +128,7 @@ int main(int argc, char *argv[]) { 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; + meas_t *M, *aM, *eM, *mM, *E, *eE, *mE, *corrE, *corrmE; M = calloc(1, sizeof(meas_t)); aM = calloc(1, sizeof(meas_t)); @@ -144,8 +138,15 @@ int main(int argc, char *argv[]) { eE = calloc(1, sizeof(meas_t)); mE = calloc(1, sizeof(meas_t)); + corrE = (meas_t *)calloc(1, sizeof(meas_t)); + corrmE = (meas_t *)calloc(1, sizeof(meas_t)); + double batch_mean_energy; + uint64_t batch_size = exp(0.667 * log(N)); + uint64_t n_batch = 0; + uint64_t batch_flips = 0; + printf("\n"); - while (diff > eps || diff == 0. || n_runs < min_runs) { + while ((diff > eps && n_steps < N) || n_runs < min_runs) { printf("\033[F\033[JWOLFF: sweep %" PRIu64 ", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n", n_runs, fabs(E->dx / E->x), M->dx / M->x, E->dc / E->c, M->dc / M->c, clust_per_sweep); @@ -153,9 +154,25 @@ int main(int argc, char *argv[]) { uint32_t n_flips = 0; uint32_t n_clust = 0; - while (n_flips / h->nv < N) { - n_flips += wolff_step(T, H, s, sim, r, bond_probs); + while (n_flips / h->nv < n) { + uint32_t tmp_flips = wolff_step(T, H, s, sim, r, bond_probs); + n_flips += tmp_flips; n_clust++; + + n_steps += tmp_flips; + + if (n_runs > 0) { + update_meas(corrE, s->H); + if (batch_flips <= batch_size && batch_flips + tmp_flips > batch_size) { + update_meas(corrmE, batch_mean_energy / n_batch); + batch_mean_energy = 0; + n_batch = 0; + batch_flips = (int64_t)(batch_flips + tmp_flips) - (int64_t)batch_size; + } + batch_mean_energy += s->H; + n_batch++; + batch_flips += tmp_flips; + } } double HH = 1; @@ -181,20 +198,9 @@ int main(int argc, char *argv[]) { diff = fabs(M->dc / M->c); } - clust_per_sweep = add_to_avg(clust_per_sweep, n_clust * 1. / N, n_runs); + 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 @@ -203,13 +209,16 @@ int main(int argc, char *argv[]) { FILE *outfile = fopen("out.dat", "a"); + double tau = batch_size * corrmE->c / corrE->c; + double dtau = tau * sqrt(pow(corrmE->dc / corrmE->c, 2) + pow(corrE->dc / corrE->c, 2)); + fprintf(outfile, - "%u %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", L, + "%u %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", L, T, H, n_runs, E->x / h->nv, E->dx / h->nv, M->x / h->nv, M->dx / h->nv, E->c / h->nv, E->dc / h->nv, M->c / h->nv, M->dc / h->nv, eE->n, eE->x / h->nv, eE->dx / h->nv, eM->x / h->nv, eM->dx / h->nv, eE->c / h->nv, eE->dc / h->nv, eM->c / h->nv, eM->dc / h->nv, mE->n, mE->x / h->nv, mE->dx / h->nv, mM->x / h->nv, mM->dx / h->nv, mE->c / h->nv, mE->dc / h->nv, mM->c / h->nv, mM->dc / h->nv, - aM->x / h->nv, aM->dx / h->nv, aM->c / h->nv, aM->dc / h->nv); + aM->x / h->nv, aM->dx / h->nv, aM->c / h->nv, aM->dc / h->nv, tau, dtau); fclose(outfile); if (output_state) { @@ -221,17 +230,6 @@ 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); @@ -243,6 +241,8 @@ int main(int argc, char *argv[]) { free(E); free(eE); free(mE); + free(corrE); + free(corrmE); free(bond_probs); graph_free(h); -- cgit v1.2.3-70-g09d2