From 5fe3e194a2d4690b541ecee8ea342681c4eeaa6c Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
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<T, D> act(const Vector<T, D>& s) const {
+ Vector<T, D> 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<T, D> tnew = r * x.t + t;
Matrix<T, D> 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<double, D> 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<double, D>&) override { tmp = 0; }
-
- void plain_site_transformed(const model&, const Spin<double, D, double>*,
- const Spin<double, D, double>&) override {
- tmp++;
- }
-
- void post_cluster(const model& m) override {
+ void pre_cluster(const model& m, unsigned, const Transformation<double, D, TorusGroup<double, D>, 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<double, 2, double>& s : m.s) {
+ glColor3d(1.0, 0.0, 0.0);
+ glBegin(GL_LINES);
+ Vector<double, 2> r = t->r.t / 2;
+ double θ = atan2(r(1), r(0));
+ Vector<double, 2> v1 = s0_tmp.act({r(0) + L * sin(θ), r(1) - L * cos(θ)});
+ Vector<double, 2> 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<double, 2, double>* 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, D, TorusGroup<double, D>, 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<double, 2, double>* 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<TorusGroup<double, D>(const model&, randutils::mt19937_rng&)>
-eGen(const std::vector<Matrix<double, 2>>& mats, const std::vector<Vector<double, 2>>& vecs,
- double ε, double L) {
- return
- [&mats, &vecs, L, ε](const model& M, randutils::mt19937_rng& rng) -> TorusGroup<double, 2> {
- Vector<double, D> t;
- Matrix<double, D> 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<double, std::normal_distribution>(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, D, TorusGroup<double, D>, double> eGen(double L) {
+ std::vector<Vector<double, 2>> torusVectors = torus_vecs<double, 2>(L);
+ std::vector<Matrix<double, 2>> torusMatrices = torus_mats<double, 2>();
+ return [L, torusVectors,
+ torusMatrices](Model<double, D, TorusGroup<double, D>, double>& M,
+ Rng& r) -> Transformation<double, D, TorusGroup<double, D>, double>* {
+ Matrix<double, 2> m;
+ Vector<double, 2> t;
- t = vecs[flip - mats.size()];
- }
+ m = r.pick(torusMatrices);
+ t(0) = r.uniform<double>(0, L);
+ t(1) = r.uniform<double>(0, L);
+ t = t - m * t;
+
+ TorusGroup<double, 2> g = TorusGroup<double, 2>({(double)L, t, m});
- TorusGroup<double, D> g(M.L, t, m);
- return g;
- };
+ Spin<double, 2, double>* ss = r.pick(M.s);
+
+ return new SpinFlip<double, 2, TorusGroup<double, 2>, 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<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
[L, a, k](const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
@@ -158,47 +194,26 @@ int main(int argc, char* argv[]) {
return H * s.x(1);
};
- std::vector<Matrix<double, D>> mats = torus_mats<double, D>();
- std::vector<Vector<double, D>> vecs = torus_vecs<double, D>(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<double, D> 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<double, D, double> s : sphere.s) {
- Spin<double, D, double> rs = sphere.s0.inverse().act(s);
- snapfile << rs.x.transpose() << "\n";
+ for (unsigned i = 0; i < sphere.s.size(); i++) {
+ Spin<double, 2, double>* ss = new Spin<double, 2, double>();
+ 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, D, Euclidean<double, D>, 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<double, 2, double>* 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, D, Euclidean<double, D>, double>* t) override {
s0_tmp = m.s0.inverse();
tmp = 0;
@@ -71,26 +49,58 @@ public:
Vector<double, 2> 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<double, 2, double>* 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, D, Euclidean<double, D>, 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<double, 2, double>* 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<double(Spin<double, D, double>)> B_hard = [L, H](Spin<double, D, double> s) -> double {
+ if (fabs(s.x(0)) < 3 * L / 4 && fabs(s.x(1)) < 3 * L / 4) {
+ return 0;
+ } else {
+ return std::numeric_limits<double>::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<double, 2, double>* ss = new Spin<double, 2, double>();
- 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-70-g09d2