diff options
| -rw-r--r-- | ising.cpp | 37 | ||||
| -rw-r--r-- | space_wolff.hpp | 107 |
2 files changed, 113 insertions, 31 deletions
@@ -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); |