summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--ising.cpp37
-rw-r--r--space_wolff.hpp107
2 files changed, 113 insertions, 31 deletions
diff --git a/ising.cpp b/ising.cpp
index 645c759..e4ab65c 100644
--- a/ising.cpp
+++ b/ising.cpp
@@ -1,18 +1,25 @@
#include "space_wolff.hpp"
+std::function<double(spin<signed, 2, signed>)> B_sin(unsigned L, unsigned n, double H) {
+ return [n, H, L] (spin<signed, 2, signed> s) -> double {
+ return H * cos(2 * M_PI * n * s.x[0] / ((double)L));
+ };
+}
+
int main(int argc, char* argv[]) {
const unsigned D = 2;
unsigned L = 32;
unsigned N = 1000;
+ unsigned mod = 0;
double T = 2.0 / log(1.0 + sqrt(2.0));
double H = 1.0;
double ε = 0.1;
int opt;
- while ((opt = getopt(argc, argv, "N:L:T:H:e:")) != -1) {
+ while ((opt = getopt(argc, argv, "N:L:T:H:e:m:")) != -1) {
switch (opt) {
case 'N':
N = (unsigned)atof(optarg);
@@ -29,6 +36,9 @@ int main(int argc, char* argv[]) {
case 'e':
ε = atof(optarg);
break;
+ case 'm':
+ mod = atoi(optarg);
+ break;
default:
exit(1);
}
@@ -62,11 +72,19 @@ int main(int argc, char* argv[]) {
}
};
- std::function<double(spin<signed, D, signed>)> B =
+ std::function<double(spin<signed, D, signed>)> B_face =
[L, H] (spin<signed, D, signed> s) -> double {
return H * s.s * smiley[s.x(1) * 16 / L][s.x(0) * 16 / L];
};
+ std::function<double(spin<signed, D, signed>)> B;
+
+ if (mod > 0) {
+ B = B_sin(L, mod, H);
+ } else {
+ B = B_face;
+ }
+
std::function<std::set<unsigned>(model<signed, D, signed>&, unsigned, spin<signed, D, signed>)> neighbors =
[] (model<signed, D, signed>& m, unsigned i0, spin<signed, D, signed> s1) -> std::set<unsigned> {
std::set<unsigned> nn;
@@ -127,7 +145,7 @@ int main(int argc, char* argv[]) {
ising.wolff(T, N, rng);
std::array<double, 2> τ = ising.Eq.τ();
std::cout << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n";
- if (τ[1] / τ[0] < ε) {
+ if (τ[1] / τ[0] < ε && τ[0] * 1e4 < ising.Eq.num_added()) {
break;
}
}
@@ -140,20 +158,11 @@ int main(int argc, char* argv[]) {
}
std::ofstream outfile;
- outfile.open("snap.dat");
-
- for (unsigned i = 0; i < L; i++) {
- for (unsigned j = 0; j < L; j++) {
- unsigned out = output[L * i + j] == 1 ? 1 : 0;
- outfile << out << " ";
- }
- outfile << "\n";
- }
- outfile.close();
+ outfile.open("out.dat", std::ios::app);
std::array<double, 2> act = ising.Eq.τ();
- std::cout << act[0] << " " << act[1] << "\n";
+ outfile << L << " " << T << " " << mod << " " << H << " " << ising.Cq.avg() << " " << ising.Cq.serr() << " " << act[0] << " " << act[1] << "\n";
return 0;
}
diff --git a/space_wolff.hpp b/space_wolff.hpp
index b60e274..9a4b4b4 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -237,6 +237,10 @@ class quantity {
double serr() const {
return sqrt(this->σ());
}
+
+ unsigned num_added() const {
+ return n - wait;
+ }
};
template <class U, int D, class state>
@@ -249,16 +253,87 @@ class model {
std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> neighbors;
std::function<double(spin<U, D, state>, spin<U, D, state>)> Z;
std::function<double(spin<U, D, state>)> B;
+ std::vector<matrix<U, D>> mats;
+ std::vector<vector<U, D>> steps;
long double E;
quantity Eq;
quantity Cq;
+ void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
+ if (level > 0) {
+ unsigned new_level = level - 1;
+ unsigned old_length = sequences.size();
+ for (std::array<double, D>& sequence : sequences) {
+ std::array<double, D> new_sequence = sequence;
+ new_sequence[new_level] = -1;
+ sequences.push_front(new_sequence);
+ }
+ one_sequences(sequences, new_level);
+ }
+ }
+
model(unsigned L, std::function<double(spin<U, D, state>, spin<U, D, state>)> Z,
std::function<double(spin<U, D, state>)> B,
std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> ns) :
L(L), s0(L), dict(L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) {
+ std::array<double, D> ini_sequence;
+ ini_sequence.fill(1);
+ std::list<std::array<double, D>> sequences;
+ sequences.push_back(ini_sequence);
+
+ one_sequences(sequences, D);
+
+ sequences.pop_front(); // don't want the identity matrix!
+
+ for (std::array<double, D> sequence : sequences) {
+ matrix<U, D> m;
+ for (unsigned i = 0; i < D; i++) {
+ for (unsigned j = 0; j < D; j++) {
+ if (i == j) {
+ m(i, j) = sequence[i];
+ } else {
+ m(i, j) = 0;
+ }
+ }
+ }
+
+ mats.push_back(m);
+
+ vector<U, D> v;
+ for (unsigned i = 0; i < D; i++) {
+ if (sequence[i] == 1) {
+ v(i) = 0;
+ } else {
+ v(i) = L / 2;
+ }
+ }
+
+ steps.push_back(v);
}
+ for (unsigned i = 0; i < D; i++) {
+ for (unsigned j = 0; j < D; j++) {
+ if (i != j) {
+ matrix<U, D> m;
+ for (unsigned k = 0; k < D; k++) {
+ for (unsigned l = 0; l < D; l++) {
+ if ((k == i && l == j) || (k == j && l == i)) {
+ if (i < j) {
+ m(k, l) = 1;
+ } else {
+ m(k, l) = -1;
+ }
+ } else {
+ m(k, l) = 0;
+ }
+ }
+ }
+ mats.push_back(m);
+ }
+ }
+ }
+ }
+
void update_energy() {
E = 0;
for (unsigned i = 0; i < s.size(); i++) {
@@ -347,31 +422,29 @@ class model {
std::uniform_real_distribution<double> t_dist(0, L);
std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1);
+ std::uniform_int_distribution<unsigned> coin(0, mats.size() + steps.size() - 1);
for (unsigned i = 0; i < N; i++) {
vector<U, D> t;
matrix<U, D> m;
- for (unsigned j = 0; j < D; j++) {
- t(j) = (U)t_dist(rng);
- }
-
- unsigned flip_D1 = r_dist(rng);
- unsigned flip_D2 = r_dist(rng);
-
- for (unsigned j = 0; j < D; j++) {
- for (unsigned k = 0; k < D; k++) {
- if ((j == flip_D1 && k == flip_D2) || (j == flip_D2 && k == flip_D1)) {
- if (flip_D1 <= flip_D2) {
- m(j, k) = -1;
- } else {
+ unsigned flip = coin(rng);
+ if (flip < mats.size()) {
+ for (unsigned j = 0; j < D; j++) {
+ t(j) = (U)t_dist(rng);
+ }
+ m = mats[flip];
+ } else {
+ for (unsigned j = 0; j < D; j++) {
+ for (unsigned k = 0; k < D; k++) {
+ if (j == k) {
m(j, k) = 1;
+ } else {
+ m(j, k) = 0;
}
- } else if ((j == k && j != flip_D1) && j != flip_D2) {
- m(j, k) = 1;
- } else {
- m(j, k) = 0;
}
}
+
+ t = steps[flip - mats.size()];
}
euclidean<U, D> g(L, t, m);