summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/cluster.h12
-rw-r--r--lib/correlation.h14
-rw-r--r--lib/measure.h4
-rw-r--r--lib/state.h16
-rw-r--r--lib/vector.h24
-rw-r--r--src/wolff_On.cpp48
6 files changed, 78 insertions, 40 deletions
diff --git a/lib/cluster.h b/lib/cluster.h
index 8061d34..2225e2f 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -94,12 +94,15 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
if (is_ext) {
X_t rs_old, rs_new;
+ v_t non_ghost;
if (vn == state->nv) {
rs_old = act_inverse (R_old, si_old);
rs_new = act_inverse (R_old, si_new);
+ non_ghost = v;
} else {
rs_old = act_inverse (R_old, sj);
rs_new = act_inverse (R_new, sj);
+ non_ghost = vn;
}
double dE = state->H(rs_old) - state->H(rs_new);
prob = 1.0 - exp(-dE / state->T);
@@ -108,6 +111,15 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
add (&(state->M), rs_new);
state->E += dE;
+ for (D_t i = 0; i < state->D; i++) {
+ double x = (double)((non_ghost / (v_t)pow(state->L, state->D - i - 1)) % state->L) / (double)state->L;
+ scalar_subtract(&(state->ReF[i]), cos(2 * M_PI * x), rs_old);
+ scalar_add(&(state->ReF[i]), cos(2 * M_PI * x), rs_new);
+
+ scalar_subtract(&(state->ImF[i]), sin(2 * M_PI * x), rs_old);
+ scalar_add(&(state->ImF[i]), sin(2 * M_PI * x), rs_new);
+ }
+
free_spin (rs_old);
free_spin (rs_new);
} else {
diff --git a/lib/correlation.h b/lib/correlation.h
index b3b7d86..26c3a99 100644
--- a/lib/correlation.h
+++ b/lib/correlation.h
@@ -7,15 +7,13 @@
#include <fftw3.h>
template <class R_t, class X_t>
-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++) {
- in[i] = correlation_component(s->spins[i]);
- }
-
- fftw_execute(plan);
+double correlation_length(const state_t <R_t, X_t> *s) {
+ double total = 0;
- double length = pow(out[0], 2);
+ for (D_t j = 0; j < s->D; j++) {
+ total += norm_squared(s->ReF[j]) + norm_squared(s->ImF[j]);
+ }
- return length;
+ return total / s->D;
}
diff --git a/lib/measure.h b/lib/measure.h
index 3474684..d20c081 100644
--- a/lib/measure.h
+++ b/lib/measure.h
@@ -36,9 +36,9 @@ std::function <void(const state_t <R_t, X_t> *)> measurement_magnetization_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) {
+std::function <void(const state_t <R_t, X_t> *)> measurement_fourier_file(FILE *file) {
return [=](const state_t <R_t, X_t> *s) {
- float smaller_X = (float)correlation_length(s, plan, fftw_in, fftw_out);
+ float smaller_X = (float)correlation_length(s);
fwrite(&smaller_X, sizeof(float), 1, file);
};
}
diff --git a/lib/state.h b/lib/state.h
index 1491938..c70459f 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -20,6 +20,8 @@ class state_t {
double E;
X_t M; // the "sum" of the spins, like the total magnetization
v_t last_cluster_size;
+ X_t *ReF;
+ X_t *ImF;
std::function <double(X_t, X_t)> J;
std::function <double(X_t)> H;
@@ -36,8 +38,14 @@ class state_t {
}
init (&R);
E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]);
- M = scalar_multiple (nv, spins[0]);
+ M = scalar_multiple(nv, spins[0]);
last_cluster_size = 0;
+ ReF = (X_t *)malloc(D * sizeof(X_t));
+ ImF = (X_t *)malloc(D * sizeof(X_t));
+ for (D_t i = 0; i < D; i++) {
+ ReF[i] = scalar_multiple(0, spins[0]);
+ ImF[i] = scalar_multiple(0, spins[0]);
+ }
}
~state_t() {
@@ -48,6 +56,12 @@ class state_t {
free(spins);
free_spin(R);
free_spin(M);
+ for (D_t i = 0; i < D; i++) {
+ free_spin(ReF[i]);
+ free_spin(ImF[i]);
+ }
+ free(ReF);
+ free(ImF);
}
};
diff --git a/lib/vector.h b/lib/vector.h
index 2099d28..a1084e2 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -37,6 +37,20 @@ void add (vector_t <q, T> *v1, vector_t <q, T> v2) {
}
template <q_t q, class T>
+void scalar_add (vector_t <q, T> *v1, T a, vector_t <q, T> v2) {
+ for (q_t i = 0; i < q; i++) {
+ v1->x[i] += a * v2.x[i];
+ }
+}
+
+template <q_t q, class T>
+void scalar_subtract (vector_t <q, T> *v1, T a, vector_t <q, T> v2) {
+ for (q_t i = 0; i < q; i++) {
+ v1->x[i] -= a * v2.x[i];
+ }
+}
+
+template <q_t q, class T>
void subtract (vector_t <q, T> *v1, vector_t <q, T> v2) {
for (q_t i = 0; i < q; i++) {
v1->x[i] -= v2.x[i];
@@ -44,6 +58,16 @@ void subtract (vector_t <q, T> *v1, vector_t <q, T> v2) {
}
template <q_t q, class T>
+T norm_squared (vector_t <q, T> v) {
+ T tmp = 0;
+ for (q_t i = 0; i < q; i++) {
+ tmp += v.x[i] * v.x[i];
+ }
+
+ return tmp;
+}
+
+template <q_t q, class T>
vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) {
vector_t <q, T> multiple;
multiple.x = (T *)malloc(q * sizeof(T));
diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index 329f468..d718440 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -115,24 +115,33 @@ int main(int argc, char *argv[]) {
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, pow(L, D), D * pow(L, D), 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, "<| \"ID\" -> %lu, \"MODEL\" -> \"%s\", \"q\" -> %d, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"FIELD_TYPE\" -> ", timestamp, ON_strings[N_COMP], N_COMP, D, L, (v_t)pow(L, D), D * (v_t)pow(L, D), T);
+ if (modulated_field) {
+ fprintf(outfile_info, "\"MODULATED\", \"ORDER\" -> %d, \"H\" -> %.15f, ", order, H_vec[0]);
+ } else {
+ fprintf(outfile_info, "\"VECTOR\", \"H\" -> {");
+ 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, "}, ");
}
- fprintf(outfile_info, "}, \"GENERATOR\" -> \"%s\", \"EPS\" -> %g |>\n", pert_type, epsilon);
+ fprintf(outfile_info, "\"GENERATOR\" -> \"%s\"", pert_type);
+
+ if (use_pert) {
+ fprintf(outfile_info, ", \"EPS\" -> %g", epsilon);
+ }
+
+ fprintf(outfile_info, " |>\n");
fclose(outfile_info);
unsigned int n_measurements = 0;
std::function <void(const On_t *)> *measurements = (std::function <void(const On_t *)> *)calloc(POSSIBLE_MEASUREMENTS, sizeof(std::function <void(const On_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));
@@ -166,22 +175,7 @@ int main(int argc, char *argv[]) {
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_R_t, vector_R_t> (outfile_F, plan, fftw_in, fftw_out);
+ measurements[n_measurements] = measurement_fourier_file<orthogonal_R_t, vector_R_t> (outfile_F);
n_measurements++;
}
@@ -208,10 +202,6 @@ int main(int argc, char *argv[]) {
}
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);