From a0e09642be5b63704ec48778b4f87e9843949ecf Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 10 Dec 2024 14:28:30 +0100 Subject: Initial commit --- percolation.cpp | 234 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 234 insertions(+) create mode 100644 percolation.cpp (limited to 'percolation.cpp') diff --git a/percolation.cpp b/percolation.cpp new file mode 100644 index 0000000..08c6970 --- /dev/null +++ b/percolation.cpp @@ -0,0 +1,234 @@ +#include +#include +#include +#include + +#include "eigen/Eigen/Dense" +#include "randutils/randutils.hpp" +#include "pcg-cpp/include/pcg_random.hpp" + +using Rng = randutils::random_generator; + +using Real = double; +using Vector = Eigen::Matrix; +using Matrix = Eigen::Matrix; + +class Vertex; +class Edge; + +class HalfEdge { +public: + Edge& parent; + Vertex& neighbor; + HalfEdge(Edge& e, Vertex& v) : parent(e), neighbor(v) {}; +}; + +class Vertex { +public: + std::vector bonds; + unsigned index; + unsigned cluster; +}; + +class Edge { +public: + bool active; +}; + +class Cluster : public std::list> {}; +class Clusters { +public: + std::vector lists; + std::vector indices; +}; + +class Graph { +public: + unsigned D; + unsigned L; + std::vector vertices; + std::vector edges; + std::vector clusters; + + unsigned Nv() const { + return vertices.size(); + } + + unsigned Ne() const { + return edges.size(); + } + + unsigned squaredDistance(unsigned vi, unsigned vj) const { + unsigned sd = 0; + for (unsigned i = 0; i < D; i++) { + int x1 = (unsigned)fmod((unsigned)(vi / (unsigned)pow(L, i)), pow(L, i + 1)); + int x2 = (unsigned)fmod((unsigned)(vj / (unsigned)pow(L, i)), pow(L, i + 1)); + unsigned Δx = abs(x1 - x2); + if (Δx > L / 2) { + Δx = L - Δx; + } + sd += pow(Δx, 2); + } + return sd; + } + + /* Initialize a square lattice */ + Graph(unsigned D, unsigned L) : D(D), L(L), vertices(pow(L, D)), edges(D * pow(L, D)) { + for (unsigned i = 0; i < Nv(); i++) { + vertices[i].index = i; + vertices[i].bonds.reserve(2 * D); + for (unsigned j = 0; j < D; j++) { + unsigned n1 = (unsigned)pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) + + (unsigned)fmod(i + pow(L, j), pow(L, j + 1)); + unsigned n2 = (unsigned)pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) + + (unsigned)fmod(pow(L, j + 1) + i - pow(L, j), pow(L, j + 1)); + + unsigned e1 = j * Nv() + i; + unsigned e2 = j * Nv() + n2; + + HalfEdge he1(edges[e1], vertices[n1]); + HalfEdge he2(edges[e2], vertices[n2]); + + vertices[i].bonds.push_back(he1); + vertices[i].bonds.push_back(he2); + } + } + } + + void dilute(Real p, Rng& r) { + for (Edge& e : edges) { + e.active = p > r.variate(0.0, 1.0); + } + } + + void markClusters() { + clusters = {}; + for (Vertex& v : vertices) { + v.cluster = 0; + }; + unsigned nClusters = 0; + for (Vertex& v : vertices) { + if (v.cluster == 0) { + nClusters++; + clusters.push_back({}); + std::stack> inCluster; + inCluster.push(v); + while (!inCluster.empty()) { + Vertex& vn = inCluster.top(); + inCluster.pop(); + if (vn.cluster == 0) { + vn.cluster = nClusters; + clusters[nClusters - 1].push_back(vn); + } + for (HalfEdge& e : vn.bonds) { + if (e.parent.active && e.neighbor.cluster == 0) { + inCluster.push(e.neighbor); + } + } + } + } + } + } + + Matrix laplacian() const { + Matrix L = Matrix::Zero(Nv(), Nv()); + + for (const Vertex& v : vertices) { + for (const HalfEdge& e : v.bonds) { + if (e.parent.active) { + L(v.index, v.index) += 1; + L(v.index, e.neighbor.index) -= 1; + } + } + } + + return L; + } +}; + +int main(int argc, char* argv[]) { + Rng r; + + unsigned D = 2; + unsigned L = 8; + double p = 0.5; + unsigned n = 10; + + int opt; + + while ((opt = getopt(argc, argv, "D:L:p:n:")) != -1) { + switch (opt) { + case 'D': + D = atoi(optarg); + break; + case 'L': + L = atoi(optarg); + break; + case 'p': + p = atof(optarg); + break; + case 'n': + n = (unsigned)atof(optarg); + break; + default: + exit(1); + } + } + + std::vector multiplicities(D * L * L / 4); + + Graph G(D, L); + + for (const Vertex& v1 : G.vertices) { + for (const Vertex& v2 : G.vertices) { + unsigned dist = G.squaredDistance(v1.index, v2.index); + if (dist > 0) { + multiplicities[dist - 1]++; + } + } + } + + for (unsigned trials = 0; trials < n; trials++) { + G.dilute(p, r); + G.markClusters(); + Matrix laplacian = G.laplacian(); + + std::vector conductivities(L * L * D / 4); + + for (const Cluster& c : G.clusters) { + std::vector inds; + inds.reserve(c.size()); + for (const Vertex& v : c){ + inds.push_back(v.index); + } + + Matrix subLaplacian = laplacian(inds, inds); + subLaplacian(0,0)++; /* Set voltage of first node to zero */ + Matrix subLaplacianInv = subLaplacian.inverse(); + + for (unsigned i = 0; i < c.size(); i++) { + for (unsigned j = 0; j < i; j++) { + Vector input = Vector::Zero(inds.size()); + input(i) = 1; + input(j) = -1; + + Vector output = subLaplacianInv * input; + double ΔV = std::abs(output(i) - output(j)); + + conductivities[G.squaredDistance(inds[i], inds[j]) - 1] += 1 / ΔV; + } + } + } + + for (unsigned i = 0; i < conductivities.size(); i++) { + if (multiplicities[i] != 0) { + std::cout << conductivities[i] / multiplicities[i] << " "; + } else { + std::cout << 0 << " "; + } + } + std::cout << std::endl; + } + + return 0; +} -- cgit v1.2.3-70-g09d2