summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--lib/ising.h55
-rw-r--r--lib/z2.h25
-rw-r--r--src/wolff_ising.cpp36
3 files changed, 74 insertions, 42 deletions
diff --git a/lib/ising.h b/lib/ising.h
index 4e5164b..e09b39e 100644
--- a/lib/ising.h
+++ b/lib/ising.h
@@ -1,10 +1,28 @@
#pragma once
#include <cmath>
-#include <stdlib.h>
#include "types.h"
+/* The following is the minimum definition of a spin class.
+ *
+ * The class must contain an M_t and an F_t for holding the sum of an
+ * integer number of spins and a double-weighted number of spins,
+ * respectively.
+ *
+ * void init(X_t *p);
+ * void free_spin(X_t p);
+ * void free_spin(M_t p);
+ * void free_spin(F_t p);
+ * X_t copy(X_t x);
+ * void add(M_t *x1, int factor, X_t x2);
+ * void add(F_t *x1, double factor, X_t x2);
+ * M_t scalar_multiple(int factor, X_t x);
+ * double norm_squared(F_t x);
+ * void write_magnetization(M_t M, FILE *outfile);
+ *
+ */
+
class ising_t {
public:
bool x;
@@ -33,8 +51,15 @@ ising_t copy(ising_t s) {
return s;
}
-template <class T>
-void add(T *s1, T a, ising_t s2) {
+void add(int *s1, int a, ising_t s2) {
+ if (s2.x) {
+ *s1 -= a;
+ } else {
+ *s1 += a;
+ }
+}
+
+void add(double *s1, double a, ising_t s2) {
if (s2.x) {
*s1 -= a;
} else {
@@ -50,31 +75,11 @@ int scalar_multiple(int factor, ising_t s) {
}
}
-
double norm_squared(double s) {
return pow(s, 2);
}
-template <class T>
-void write_magnetization(T M, FILE *outfile) {
- fwrite(&M, sizeof(T), 1, outfile);
-}
-
-// below this line is unnecessary, but convenient
-
-double ising_dot(ising_t s1, ising_t s2) {
- if (s1.x == s2.x) {
- return 1.0;
- } else {
- return -1.0;
- }
-}
-
-double scalar_field(ising_t s, double H) {
- if (s.x) {
- return -H;
- } else {
- return H;
- }
+void write_magnetization(int M, FILE *outfile) {
+ fwrite(&M, sizeof(int), 1, outfile);
}
diff --git a/lib/z2.h b/lib/z2.h
index 599a6a5..dae1eb7 100644
--- a/lib/z2.h
+++ b/lib/z2.h
@@ -1,10 +1,21 @@
#pragma once
-#include <gsl/gsl_rng.h>
#include "types.h"
#include "ising.h"
-#include "state.h"
+
+/* The minimum definition for a group type R_t to act on a spin type X_t is
+ * given by the following.
+ *
+ * void init(R_t *p);
+ * void free_spin(R_t r);
+ * R_t copy(R_t r);
+ * X_t act(R_t r, X_t x);
+ * R_t act(R_t r, R_t r);
+ * X_t act_inverse(R_t r, X_t x);
+ * R_t act_inverse(R_t r, R_t r);
+ *
+ */
struct z2_t { bool x; };
@@ -52,13 +63,3 @@ z2_t act_inverse(z2_t r1, z2_t r2) {
return act(r1, r2);
}
-// these are all functions necessary for wolff.h
-
-z2_t generate_ising_rotation(gsl_rng *r, const state_t <z2_t, ising_t> *s) {
- z2_t rot;
- rot.x = true;
- return rot;
-}
-
-
-
diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp
index fd2f0bb..e9d0185 100644
--- a/src/wolff_ising.cpp
+++ b/src/wolff_ising.cpp
@@ -45,24 +45,50 @@ int main(int argc, char *argv[]) {
gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
gsl_rng_set(r, rand_seed());
- state_t <z2_t, ising_t> s(D, L, T, ising_dot, std::bind(scalar_field, std::placeholders::_1, H));
+ // define spin-spin coupling
+ std::function <double(ising_t, ising_t)> Z = [] (ising_t s1, ising_t s2) -> double {
+ if (s1.x == s2.x) {
+ return 1.0;
+ } else {
+ return -1.0;
+ }
+ };
- std::function <z2_t(gsl_rng *, const state_t <z2_t, ising_t> *)> gen_R = generate_ising_rotation;
+ // define spin-field coupling
+ std::function <double(ising_t)> B = [=] (ising_t s) -> double {
+ if (s.x) {
+ return -H;
+ } else {
+ return H;
+ }
+ };
- double average_M = 0;
+ // initialize state object
+ state_t <z2_t, ising_t> s(D, L, T, Z, B);
- typedef std::function <void(const state_t <z2_t, ising_t> *)> meas_func;
+ // define function that generates self-inverse rotations
+ std::function <z2_t(gsl_rng *, const state_t <z2_t, ising_t> *)> gen_R = [] (gsl_rng *, const state_t <z2_t, ising_t> *) -> z2_t {
+ z2_t rot;
+ rot.x = true;
+ return rot;
+ };
- meas_func measurement = [&] (const state_t <z2_t, ising_t> *s) {
+ // define function that updates any number of measurements
+ double average_M = 0;
+ std::function <void(const state_t <z2_t, ising_t> *)> measurement = [&] (const state_t <z2_t, ising_t> *s) {
average_M += (double)s->M / (double)N / (double)s->nv;
};
+ // run wolff for N cluster flips
wolff(N, &s, gen_R, measurement, r, silent);
+ // tell us what we found!
printf("%" PRIcount " Ising runs completed. D = %" PRID ", L = %" PRIL ", T = %g, H = %g, <M> = %g\n", N, D, L, T, H, average_M);
+ // free the random number generator
gsl_rng_free(r);
return 0;
+
}