From 614575bb88a2cadc9e35b684d0f1712de822ef0d Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 14 Jan 2020 18:00:07 -0500 Subject: generalized code somewhat to accomodate simulation of infinite space, sample central field sim --- space_wolff.hpp | 157 ++++++++++++++++++++++++++++++++++++++---- spheres.cpp | 18 ++--- spheres_infinite.cpp | 189 +++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 340 insertions(+), 24 deletions(-) create mode 100644 spheres_infinite.cpp diff --git a/space_wolff.hpp b/space_wolff.hpp index 615fb93..7345726 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -9,6 +9,25 @@ #include #include #include +#include +#include +#include + +namespace std { +template struct hash> { + typedef array argument_type; + typedef size_t result_type; + + result_type operator()(const argument_type& a) const { + hash hasher; + result_type h = 0; + for (result_type i = 0; i < N; ++i) { + h = h * 31 + hasher(a[i]); + } + return h; + } +}; +} // namespace std #include "randutils/randutils.hpp" @@ -56,6 +75,54 @@ public: }; template class Euclidean { +public: + Vector t; + Matrix r; + + Euclidean(T L) { + for (unsigned i = 0; i < D; i++) { + t(i) = 0; + r(i, i) = 1; + for (unsigned j = 1; j < D; j++) { + r(i, (i + j) % D) = 0; + } + } + } + + Euclidean(Vector t0, Matrix r0) { + t = t0; + r = r0; + } + + template Spin act(const Spin& s) const { + Spin s_new; + + s_new.x = t + r * s.x; + s_new.s = s.s; + + return s_new; + } + + Euclidean act(const Euclidean& x) const { + Vector tnew = r * x.t + t; + Matrix rnew = r * x.r; + + Euclidean pnew(tnew, rnew); + + return pnew; + } + + Euclidean inverse() const { + Vector tnew = -r.transpose() * t; + Matrix rnew = r.transpose(); + + Euclidean pnew(tnew, rnew); + + return pnew; + } +}; + +template class TorusGroup { private: T L; @@ -63,9 +130,9 @@ public: Vector t; Matrix r; - /** brief Euclidean - default constructor, constructs the identity + /** brief TorusGroup - default constructor, constructs the identity */ - Euclidean(T L) : L(L) { + TorusGroup(T L) : L(L) { for (unsigned i = 0; i < D; i++) { t(i) = 0; r(i, i) = 1; @@ -75,7 +142,7 @@ public: } } - Euclidean(T L, Vector t0, Matrix r0) : L(L) { + TorusGroup(T L, Vector t0, Matrix r0) : L(L) { t = t0; r = r0; } @@ -93,7 +160,7 @@ public: return s_new; } - Euclidean act(const Euclidean& x) const { + TorusGroup act(const TorusGroup& x) const { Vector tnew = r * x.t + t; Matrix rnew = r * x.r; @@ -101,21 +168,82 @@ public: tnew(i) = fmod(L + tnew(i), L); } - Euclidean pnew(this->L, tnew, rnew); + TorusGroup pnew(this->L, tnew, rnew); return pnew; } - Euclidean inverse() const { + TorusGroup inverse() const { Vector tnew = -r.transpose() * t; Matrix rnew = r.transpose(); - Euclidean pnew(this->L, tnew, rnew); + TorusGroup pnew(this->L, tnew, rnew); return pnew; } }; +template class Octree { +private: + unsigned N; + T L; + std::unordered_map, std::set*>> data; + +public: + Octree(T L, unsigned N) { + L = L; + N = N; + } + + std::array ind(Vector x) const { + std::array ind; + + for (unsigned i = 0; i < D; i++) { + ind[i] = (signed)std::floor(N * x(i) / L); + } + + return ind; + } + + void insert(Spin* s) { data[ind(s->x)].insert(s); }; + + void remove(Spin* s) { + data[ind(s->x)].erase(s); + if (data[ind(s->x)].empty()) { + data.erase(ind(s->x)); + } + }; + + std::set*> neighbors(const Vector& x) const { + std::array i0 = ind(x); + std::set*> ns; + + nearest_neighbors_of(i0, D + 1, ns); + + return ns; + }; + + void nearest_neighbors_of(std::array i0, unsigned depth, std::set*>& ns) const { + if (depth == 0) { + auto it = data.find(i0); + if (it != data.end()) { + ns.insert(it->second.begin(), it->second.end()); + } + } else { + for (signed j : {-1, 0, 1}) { + std::array i1 = i0; + if (N < 2) { + i1[depth - 1] += j; + } else { + i1[depth - 1] = (N + i1[depth - 1] + j) % N; + } + nearest_neighbors_of(i1, depth - 1, ns); + } + } + }; +}; + +/* template class Octree { private: T L; @@ -168,6 +296,7 @@ public: } } }; +*/ class quantity { private: @@ -360,18 +489,16 @@ public: R s0; std::vector> s; Octree dict; - unsigned dDepth; - unsigned nDepth; std::function&, const Spin&)> Z; std::function&)> B; randutils::mt19937_rng rng; Model(U L, std::function&, const Spin&)> Z, - std::function&)> B, unsigned dDepth, unsigned nDepth) - : L(L), s0(L), dict(L, dDepth), Z(Z), B(B), dDepth(dDepth), nDepth(nDepth), rng() {} + std::function&)> B) + : L(L), s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {} - void step(double T, unsigned ind, Euclidean r, measurement& A) { + void step(double T, unsigned ind, R& r, measurement& A) { unsigned cluster_size = 0; std::queue*> queue; @@ -387,7 +514,7 @@ public: visited.insert(si); if (si == NULL) { - Euclidean s0_new = r.act(s0); + R s0_new = r.act(s0); for (Spin& ss : s) { Spin s0s_old = s0.inverse().act(ss); Spin s0s_new = s0_new.inverse().act(ss); @@ -402,8 +529,8 @@ public: } else { Spin si_new = r.act(*si); std::set*> all_neighbors; - std::set*> current_neighbors = dict.neighbors(si->x, nDepth); - std::set*> new_neighbors = dict.neighbors(si_new.x, nDepth); + std::set*> current_neighbors = dict.neighbors(si->x); + std::set*> new_neighbors = dict.neighbors(si_new.x); all_neighbors.insert(current_neighbors.begin(), current_neighbors.end()); all_neighbors.insert(new_neighbors.begin(), new_neighbors.end()); diff --git a/spheres.cpp b/spheres.cpp index f3dde42..97c1e18 100644 --- a/spheres.cpp +++ b/spheres.cpp @@ -3,9 +3,9 @@ #include const unsigned D = 2; -typedef Model, double> model; +typedef Model, double> model; -class animation : public measurement, double> { +class animation : public measurement, double> { private: uint64_t t1; uint64_t t2; @@ -27,7 +27,7 @@ class animation : public measurement, double> { gluOrtho2D(-1, L + 1, -1 , L + 1); } - void pre_cluster(const model&, unsigned, const Euclidean&) override { + void pre_cluster(const model&, unsigned, const TorusGroup&) override { tmp = 0; } @@ -65,8 +65,8 @@ class animation : public measurement, double> { } }; -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) -> Euclidean { +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)); @@ -91,7 +91,7 @@ std::function(const model&, randutils::mt19937_rng&)> eGen( t = vecs[flip - mats.size()]; } - Euclidean g(M.L, t, m); + TorusGroup g(M.L, t, m); return g; }; @@ -130,8 +130,8 @@ int main(int argc, char* argv[]) { } } - double k = 100.0; - double a = 0.2; + double k = 1e2; + double a = 0.05; std::function&, const Spin&)> Z = [L, a, k] (const Spin& s1, const Spin& s2) -> double { @@ -158,7 +158,7 @@ int main(int argc, char* argv[]) { std::vector> vecs = torus_vecs(L); auto g = eGen(mats, vecs, L, L); animation A(L, 750, argc, argv); - model sphere(L, Z, B, std::floor(log2(L)), 2); + model sphere(L, Z, B); randutils::mt19937_rng rng; diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp new file mode 100644 index 0000000..f2f998a --- /dev/null +++ b/spheres_infinite.cpp @@ -0,0 +1,189 @@ + +#include "space_wolff.hpp" +#include + +const unsigned D = 2; +typedef Model, double> model; + +class animation : public measurement, double> { + private: + uint64_t t1; + uint64_t t2; + unsigned n; + unsigned tmp; + public: + animation(double L, unsigned w, int argc, char *argv[]) { + t1 = 0; + t2 = 0; + n = 0; + + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(w, w); + glutCreateWindow("wolff"); + glClearColor(0.0,0.0,0.0,0.0); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(-L, L, -L , L); + } + + void pre_cluster(const model&, unsigned, const Euclidean&) override { + tmp = 0; + } + + void plain_site_transformed(const model&, const Spin*, const Spin&) override { + tmp++; + } + + void post_cluster(const model& m) override { + glClearColor(1.0f, 1.0f, 1.0f, 1.0f ); + glClear(GL_COLOR_BUFFER_BIT); + for (const Spin& s : m.s) { + glBegin(GL_POLYGON); + unsigned n_points = 50; + glColor3f(0.0f, 0.0f, 0.0f); + 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)); + } + glEnd(); + } + glFlush(); + + t1 += tmp; + t2 += tmp * tmp; + n++; + } + + void clear() { + t1 = 0; + t2 = 0; + n = 0; + } + + double var() { + return (t2 - t1 * t1 / (double)n) / (double)n; + } +}; + +std::function(const model&, randutils::mt19937_rng&)> eGen(double ε) { + return [ε] (const model& M, randutils::mt19937_rng& rng) -> Euclidean { + Vector t; + Matrix m; + + double θ = rng.uniform((double)0.0, 2 * M_PI); + m(0, 0) = cos(θ); + m(1, 1) = cos(θ); + m(0, 1) = sin(θ); + m(1, 0) = -sin(θ); + + 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) = rng.variate(0.0, ε); + } + + Euclidean g(t, m); + return g; + }; + +} + +int main(int argc, char* argv[]) { + const unsigned D = 2; + + double L = 32; + unsigned N = 1000; + double T = 2.0 / log(1.0 + sqrt(2.0)); + double H = 1.0; + unsigned n = 25; + + int opt; + + while ((opt = getopt(argc, argv, "n:N:L:T:H:")) != -1) { + switch (opt) { + case 'n': + n = (unsigned)atof(optarg); + break; + case 'N': + N = (unsigned)atof(optarg); + break; + case 'L': + L = atof(optarg); + break; + case 'T': + T = atof(optarg); + break; + case 'H': + H = atof(optarg); + break; + default: + exit(1); + } + } + + double k = 1e2; + double a = 0.05; + + std::function&, const Spin&)> Z = + [L, a, k] (const Spin& s1, const Spin& s2) -> double { + Vector d = s1.x - s2.x; + + double σ = s1.s + s2.s; + double δ = σ - sqrt(d.transpose() * d); + + if (δ > - a * σ) { + return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2)); + } else if (δ > - 2 * a * σ) { + return 0.5 * k * pow(δ + 2 * a * σ, 2); + } else { + return 0; + } + }; + + std::function)> B = + [L, H] (Spin s) -> double { + return H * s.x.norm(); + }; + + auto g = eGen(1); + animation A(L, 750, argc, argv); + model sphere(1, Z, B); + + randutils::mt19937_rng rng; + + sphere.s.reserve(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(ε); + 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"; + } + + snapfile.close(); + + return 0; +} -- cgit v1.2.3-54-g00ecf