From 9437ed889390e04a1e693a5dfa5290cd90ae150e Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 14 Feb 2020 17:04:17 -0500 Subject: Dynamic pairing method appears to be working with spheres, needs to be cleaned up and ising code needs to be updated with new transformation interfaces --- space_wolff.hpp | 8 ++++++-- spheres_infinite.cpp | 27 +++++++++++++++++++++++++-- transformation.hpp | 46 ++++++++++++++++++---------------------------- 3 files changed, 49 insertions(+), 32 deletions(-) diff --git a/space_wolff.hpp b/space_wolff.hpp index 141396a..63bf118 100644 --- a/space_wolff.hpp +++ b/space_wolff.hpp @@ -50,19 +50,23 @@ public: queue.push(t0); std::set*> visited; + std::set*> fused; while (!queue.empty()) { Transformation* t = queue.front(); queue.pop(); std::set*> c = t->current(); - if (!std::any_of(c.begin(), c.end(), [&visited](Spin* s) { return visited.contains(s); })) { + if (!std::any_of(c.begin(), c.end(), [&visited](Spin* s) { return visited.contains(s); }) && !(c.size() == 1 && std::any_of(c.begin(), c.end(), [&fused](Spin* s) { return fused.contains(s); }))) { visited.insert(c.begin(), c.end()); + fused.insert(c.begin(), c.end()); for (Spin* s : t->toConsider()) { + if (!fused.contains(s)) { double ΔE = t->ΔE(s); if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform(0, 1)) { - queue.push(t->createNew(T, visited, s, rng)); + queue.push(t->createNew(T, fused, s, rng)); + } } } diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp index 5287a84..b39f164 100644 --- a/spheres_infinite.cpp +++ b/spheres_infinite.cpp @@ -117,6 +117,28 @@ Gen, double> mGen(double ε) { }; } +Gen, double> rGen(double ε) { + return [ε](model& M, Rng& rng) -> Transformation, double>* { + Vector t; + Matrix m; + + double θ = rng.uniform((double)0.0, 2 * M_PI); + m(0, 0) = -cos(2 * θ); + m(1, 1) = cos(2 * θ); + m(0, 1) = -2 * cos(θ) * sin(θ); + m(1, 0) = -2 * cos(θ) * sin(θ); + + Spin* s = rng.pick(M.s); + t = M.s0.t; + for (unsigned j = 0; j < D; j++) { + t(j) += rng.variate(0.0, ε); + } + + Euclidean g(t - m * t, m); + return new SpinFlip, double>(M, g, s); + }; +} + int main(int argc, char* argv[]) { const unsigned D = 2; @@ -179,8 +201,9 @@ int main(int argc, char* argv[]) { return H * s.x.norm(); }; - auto g1 = eGen(1); + auto g1 = eGen(2); auto g2 = mGen(0.5); + auto g3 = rGen(5); animation A(L, 750, argc, argv); model sphere(1.0, Z, B); @@ -197,7 +220,7 @@ int main(int argc, char* argv[]) { sphere.dict.insert(ss); } - sphere.wolff(T, {g1, g2}, A, N); + sphere.wolff(T, {g1, g2, g3}, A, N); std::ofstream outfile; outfile.open("test.dat"); diff --git a/transformation.hpp b/transformation.hpp index fc1833e..ad3b49f 100644 --- a/transformation.hpp +++ b/transformation.hpp @@ -12,7 +12,7 @@ public: virtual std::set*> current() const { return {NULL}; } virtual std::set*> toConsider() const { return {}; } virtual double ΔE(const Spin*) const { return 0; } - virtual Transformation* createNew(double T, const std::set*>&, Spin*, Rng& rng) const { return new Transformation(); } + virtual Transformation* createNew(double T, std::set*>&, Spin*, Rng& rng) const { return new Transformation(); } virtual void apply(){}; }; @@ -43,6 +43,8 @@ public: neighbors.insert(new_neighbors.begin(), new_neighbors.end()); neighbors.insert(NULL); + neighbors.erase(sOld); + return neighbors; } @@ -56,7 +58,7 @@ public: } } - Transformation* createNew(double, const std::set*>&, Spin* s, Rng&) const override; + Transformation* createNew(double, std::set*>&, Spin* s, Rng&) const override; void apply() override { M.dict.remove(sOld); @@ -97,6 +99,9 @@ public: neighbors.insert(new_neighbors_2.begin(), new_neighbors_2.end()); neighbors.insert(NULL); + neighbors.erase(s1Old); + neighbors.erase(s2Old); + return neighbors; } @@ -112,8 +117,8 @@ public: } } - Transformation* createNew(double T, const std::set*>& v, Spin* s, Rng&) const; - Transformation* createNew(double T, const std::set*>& v, Spin* s, Rng&) const; + Transformation* createNew(double T, std::set*>& v, Spin* s, Rng&) const; + Transformation* createNew(double T, std::set*>& v, Spin* s, Rng&) const; void apply() override { M.dict.remove(s1Old); @@ -150,11 +155,11 @@ public: return M.B(s0s_new) - M.B(s0s_old); } - Transformation* createNew(double T, const std::set*>& v, Spin* s, Rng&) const { + Transformation* createNew(double T, std::set*>& v, Spin* s, Rng&) const { return new SpinFlip(M, r, s); } - Transformation* createNew(double, const std::set*>&, Spin* s, Rng&) const { + Transformation* createNew(double, std::set*>&, Spin* s, Rng&) const { Vector v = r.act(*s).x; std::set*> on_site = M.dict.at(v); if (on_site.empty()) { @@ -168,33 +173,16 @@ public: }; template -Transformation* SpinFlip::createNew(double T, const std::set*>& v, Spin* s, Rng& rng) const { +Transformation* SpinFlip::createNew(double T, std::set*>& v, Spin* s, Rng& rng) const { if (s == NULL) { return new FieldFlip(M, r); } else { - SpinFlip* t = new SpinFlip(M, r, s); - Spin* sMax = NULL; - double ΔEMax = 0.0; - for (Spin* ss : t->toConsider()) { - if (ss != NULL && ss != s && !v.contains(ss)) { - double ΔE = t->ΔE(ss); - if (ΔE > ΔEMax) { - sMax = ss; - ΔEMax = ΔE; - } - } - } - if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform(0, 1)) { - return t; - } else { - delete t; - return new PairFlip(M, r, s, sMax); - } + return new SpinFlip(M, r, s); } } template -Transformation* PairFlip::createNew(double T, const std::set*>& v, Spin* s, Rng& rng) const { +Transformation* PairFlip::createNew(double T, std::set*>& v, Spin* s, Rng& rng) const { if (s == NULL) { return new FieldFlip(M, r); } else { @@ -210,17 +198,19 @@ Transformation* PairFlip::createNew(double T, const std: } } } - if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform(0, 1)) { + if (sMax == NULL || 1.0 - exp(-ΔEMax / T) < rng.uniform(0, 1) || v.contains(sMax)) { return t; } else { delete t; + v.insert(s); + v.insert(sMax); return new PairFlip(M, r, s, sMax); } } } template -Transformation* PairFlip::createNew(double, const std::set*>&, Spin* s, Rng&) const { +Transformation* PairFlip::createNew(double, std::set*>&, Spin* s, Rng&) const { if (s == NULL) { return new FieldFlip(M, r); } else { -- cgit v1.2.3-70-g09d2