From 5fe3e194a2d4690b541ecee8ea342681c4eeaa6c Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 21 Feb 2020 11:27:59 -0500 Subject: Fixed the periodic sphere code to work with the new backend and implemented saving cluster flip information --- euclidean.hpp | 10 +++ spheres.cpp | 173 ++++++++++++++++++++++++++++----------------------- spheres_infinite.cpp | 70 +++++++++++++-------- 3 files changed, 148 insertions(+), 105 deletions(-) diff --git a/euclidean.hpp b/euclidean.hpp index 49ec2d9..b1ce2c6 100644 --- a/euclidean.hpp +++ b/euclidean.hpp @@ -95,6 +95,16 @@ public: return s_new; } + Vector act(const Vector& s) const { + Vector s_new = t + r * s; + + for (unsigned i = 0; i < D; i++) { + s_new(i) = fmod(L + s_new(i), L); + } + + return s_new; + } + TorusGroup act(const TorusGroup& x) const { Vector tnew = r * x.t + t; Matrix rnew = r * x.r; diff --git a/spheres.cpp b/spheres.cpp index b7a9d94..4235d66 100644 --- a/spheres.cpp +++ b/spheres.cpp @@ -15,9 +15,13 @@ private: uint64_t t2; unsigned n; unsigned tmp; + unsigned wait; + double L; + TorusGroup s0_tmp; + std::ofstream& file; public: - animation(double L, unsigned w, int argc, char* argv[]) { + animation(double L, unsigned w, int argc, char* argv[], unsigned wait, std::ofstream& file) : s0_tmp(0), wait(wait), L(50 * L), file(file) { t1 = 0; t2 = 0; n = 0; @@ -25,37 +29,80 @@ public: glutInit(&argc, argv); glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); glutInitWindowSize(w, w); - glutCreateWindow("wolff"); + glutCreateWindow("wolffWindow"); glClearColor(0.0, 0.0, 0.0, 0.0); glMatrixMode(GL_PROJECTION); glLoadIdentity(); - gluOrtho2D(-1, L + 1, -1, L + 1); + gluOrtho2D(-1, L+1, -1, L+1); } - void pre_cluster(const model&, unsigned, const TorusGroup&) override { tmp = 0; } - - void plain_site_transformed(const model&, const Spin*, - const Spin&) override { - tmp++; - } - - void post_cluster(const model& m) override { + void pre_cluster(const model& m, unsigned, const Transformation, double>* t) override { + s0_tmp = m.s0.inverse(); + tmp = 0; glClearColor(1.0f, 1.0f, 1.0f, 1.0f); glClear(GL_COLOR_BUFFER_BIT); - for (const Spin& s : m.s) { + glColor3d(1.0, 0.0, 0.0); + glBegin(GL_LINES); + Vector r = t->r.t / 2; + double θ = atan2(r(1), r(0)); + Vector v1 = s0_tmp.act({r(0) + L * sin(θ), r(1) - L * cos(θ)}); + Vector v2 = s0_tmp.act({r(0) - L * sin(θ), r(1) + L * cos(θ)}); + glVertex2d(v1(0), v1(1)); + glVertex2d(v2(0), v2(1)); + if (n > wait) + file << v1.transpose() << " " << v2.transpose() << std::endl; + glEnd(); + for (const Spin* s : m.s) { glBegin(GL_POLYGON); unsigned n_points = 50; glColor3f(0.0f, 0.0f, 0.0f); + if (n > wait) + file << s << " " << s->s << " " << s0_tmp.act(*s).x.transpose() << " "; for (unsigned i = 0; i < n_points; i++) { - glVertex2d(m.s0.inverse().act(s).x(0) + s.s * cos(2 * i * M_PI / n_points), - m.s0.inverse().act(s).x(1) + s.s * sin(2 * i * M_PI / n_points)); + glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), + s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); } glEnd(); } + if (n > wait) + file << std::endl; + glLineWidth(3); glFlush(); + } + + void plain_site_transformed(const model& m, const Transformation, double>& t) override { + if (n > wait) { + if (t.current().size() > 1) { + glColor3f(0.0f, 0.0f, 1.0f); + } else { + glColor3f(0.0f, 1.0f, 0.0f); + } + for (const Spin* s : t.current()) { + if (s != NULL) { + file << s << " " << s0_tmp.act(t.r.act(*s)).x.transpose() << " "; + glBegin(GL_POLYGON); + unsigned n_points = 50; + for (unsigned i = 0; i < n_points; i++) { + glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), + s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); + } + glEnd(); + } else { + file << s << " " << s0_tmp.act({0, 0}).transpose() << " "; + } + } + } + tmp++; + } + void post_cluster(const model& m) override { t1 += tmp; t2 += tmp * tmp; + if (n > wait) { + file << std::endl; + glFlush(); + sleep(2); + } n++; } @@ -68,40 +115,29 @@ public: double var() { return (t2 - t1 * t1 / (double)n) / (double)n; } }; -std::function(const model&, randutils::mt19937_rng&)> -eGen(const std::vector>& mats, const std::vector>& vecs, - double ε, double L) { - return - [&mats, &vecs, L, ε](const model& M, randutils::mt19937_rng& rng) -> TorusGroup { - Vector t; - Matrix m; - unsigned flip = rng.uniform((unsigned)0, (unsigned)(mats.size() + vecs.size() - 1)); - if (flip < mats.size()) { - unsigned f_ind = rng.uniform((unsigned)0, (unsigned)M.s.size()); - t = M.s[f_ind].x; - for (unsigned j = 0; j < D; j++) { - t(j) = fmod(10 * L + t(j) + rng.variate(0.0, ε), L); - } - 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; - } - } - } +Gen, double> eGen(double L) { + std::vector> torusVectors = torus_vecs(L); + std::vector> torusMatrices = torus_mats(); + return [L, torusVectors, + torusMatrices](Model, double>& M, + Rng& r) -> Transformation, double>* { + Matrix m; + Vector t; - t = vecs[flip - mats.size()]; - } + m = r.pick(torusMatrices); + t(0) = r.uniform(0, L); + t(1) = r.uniform(0, L); + t = t - m * t; + + TorusGroup g = TorusGroup({(double)L, t, m}); - TorusGroup g(M.L, t, m); - return g; - }; + Spin* ss = r.pick(M.s); + + return new SpinFlip, double>(M, g, ss); + }; } + int main(int argc, char* argv[]) { const unsigned D = 2; @@ -135,8 +171,8 @@ int main(int argc, char* argv[]) { } } - double k = 1e2; - double a = 0.05; + double k = 1e8; + double a = 0.0; std::function&, const Spin&)> Z = [L, a, k](const Spin& s1, const Spin& s2) -> double { @@ -158,47 +194,26 @@ int main(int argc, char* argv[]) { return H * s.x(1); }; - std::vector> mats = torus_mats(); - std::vector> vecs = torus_vecs(L); - auto g = eGen(mats, vecs, L, L); - animation A(L, 750, argc, argv); + auto g = eGen(L); + std::ofstream ofile("test.dat"); + animation A(L, 750, argc, argv, 1000, ofile); model sphere(L, Z, B); randutils::mt19937_rng rng; - sphere.s.reserve(n); + sphere.s.resize(n); unsigned nx = floor(sqrt(n)); - for (unsigned i = 0; i < n; i++) { - Vector pos = {(i / nx) * L / nx, (i % nx) * L / nx}; - sphere.s.push_back({pos, 0.5}); - sphere.dict.insert(&sphere.s.back()); - } - - sphere.wolff(T, g, A, N); - std::ofstream outfile; - outfile.open("test.dat"); - - for (signed i = -10; i <= 10; i++) { - A.clear(); - double ε = pow(2, -4 + i / 2.0); - auto gn = eGen(mats, vecs, ε, L); - sphere.wolff(T, gn, A, N); - outfile << ε << " " << A.var() / sphere.s.size() << std::endl; - std::cout << ε << " " << A.var() / sphere.s.size() << std::endl; - } - - outfile.close(); - - std::ofstream snapfile; - snapfile.open("sphere_snap.dat"); - - for (Spin s : sphere.s) { - Spin rs = sphere.s0.inverse().act(s); - snapfile << rs.x.transpose() << "\n"; + for (unsigned i = 0; i < sphere.s.size(); i++) { + Spin* ss = new Spin(); + ss->x = {(i / nx) * L / nx, (i % nx) * L / nx}; + ss->s = rng.pick({0.45, 0.45}); + sphere.s[i] = ss; + sphere.dict.insert(ss); } - snapfile.close(); + sphere.wolff(T, {g}, A, N); + ofile.close(); return 0; } diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp index 3f8747b..56abe77 100644 --- a/spheres_infinite.cpp +++ b/spheres_infinite.cpp @@ -36,28 +36,6 @@ public: gluOrtho2D(-L, L, -L, L); } - void plain_site_transformed(const model& m, const Transformation, double>& t) override { - if (n > wait) { - if (t.current().size() > 1) { - glColor3f(0.0f, 0.0f, 1.0f); - } else { - glColor3f(0.0f, 1.0f, 0.0f); - } - for (const Spin* s : t.current()) { - if (s != NULL) { - glBegin(GL_POLYGON); - unsigned n_points = 50; - for (unsigned i = 0; i < n_points; i++) { - glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), - s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); - } - glEnd(); - } - } - } - tmp++; - } - void pre_cluster(const model& m, unsigned, const Transformation, double>* t) override { s0_tmp = m.s0.inverse(); tmp = 0; @@ -71,26 +49,58 @@ public: Vector v2 = s0_tmp.act({r(0) - L * sin(θ), r(1) + L * cos(θ)}); glVertex2d(v1(0), v1(1)); glVertex2d(v2(0), v2(1)); + if (n > wait) + file << v1.transpose() << " " << v2.transpose() << std::endl; glEnd(); for (const Spin* s : m.s) { glBegin(GL_POLYGON); unsigned n_points = 50; glColor3f(0.0f, 0.0f, 0.0f); + if (n > wait) + file << s << " " << s->s << " " << s0_tmp.act(*s).x.transpose() << " "; for (unsigned i = 0; i < n_points; i++) { glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); } glEnd(); } + if (n > wait) + file << std::endl; glLineWidth(3); glFlush(); } + void plain_site_transformed(const model& m, const Transformation, double>& t) override { + if (n > wait) { + if (t.current().size() > 1) { + glColor3f(0.0f, 0.0f, 1.0f); + } else { + glColor3f(0.0f, 1.0f, 0.0f); + } + for (const Spin* s : t.current()) { + if (s != NULL) { + file << s << " " << s0_tmp.act(t.r.act(*s)).x.transpose() << " "; + glBegin(GL_POLYGON); + unsigned n_points = 50; + for (unsigned i = 0; i < n_points; i++) { + glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points), + s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points)); + } + glEnd(); + } else { + file << s << " " << s0_tmp.act({0, 0}).transpose() << " "; + } + } + } + tmp++; + } + void post_cluster(const model& m) override { - glFlush(); t1 += tmp; t2 += tmp * tmp; if (n > wait) { + file << std::endl; + glFlush(); sleep(2); } n++; @@ -280,10 +290,18 @@ int main(int argc, char* argv[]) { return H * s.x.norm(); }; + std::function)> B_hard = [L, H](Spin s) -> double { + if (fabs(s.x(0)) < 3 * L / 4 && fabs(s.x(1)) < 3 * L / 4) { + return 0; + } else { + return std::numeric_limits::infinity(); + } + }; + auto g1 = eGen(1); auto g2 = mGen(0.1); auto g3 = tGen(0.1); - auto g4 = rGen(1); + auto g4 = rGen(0); auto tag = std::chrono::high_resolution_clock::now(); @@ -303,8 +321,8 @@ int main(int argc, char* argv[]) { unsigned nx = floor(sqrt(n)); for (unsigned i = 0; i < sphere.s.size(); i++) { Spin* ss = new Spin(); - ss->x = {(i / nx) * L / nx, (i % nx) * L / nx}; - ss->s = rng.pick({0.5, 0.5}); + ss->x = {(i / nx) * L / nx - L / 2, (i % nx) * L / nx - L / 2}; + ss->s = rng.pick({0.45, 0.45}); sphere.s[i] = ss; sphere.dict.insert(ss); } -- cgit v1.2.3