From bb519f0c9f5252bf6eb39aeb05add14f8af7be86 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 11 Dec 2024 14:05:01 +0100 Subject: Changed many variable names, and now save multiple moments to separate files. --- percolation.cpp | 169 +++++++++++++++++++++++++++++--------------------------- 1 file changed, 87 insertions(+), 82 deletions(-) (limited to 'percolation.cpp') diff --git a/percolation.cpp b/percolation.cpp index fac04ec..c1237e1 100644 --- a/percolation.cpp +++ b/percolation.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include "eigen/Eigen/Dense" #include "randutils/randutils.hpp" @@ -9,6 +10,7 @@ using Rng = randutils::random_generator; +using Bool = short unsigned; using Real = double; using Vector = Eigen::Matrix; using Matrix = Eigen::Matrix; @@ -24,8 +26,8 @@ public: class Vertex { public: - std::vector bonds; unsigned index; + std::vector outgoingEdges; std::vector coordinates; }; @@ -58,73 +60,72 @@ public: return D * Nv(); } - unsigned squaredDistance(unsigned vi, unsigned vj) const { - unsigned sd = 0; - for (unsigned i = 0; i < D; i++) { - int x1 = vertices[vi].coordinates[i]; - int x2 = vertices[vj].coordinates[i]; - unsigned Δx = abs(x1 - x2); + unsigned squaredDistance(unsigned i₁, unsigned i₂) const { + unsigned Δx² = 0; + for (unsigned d = 0; d < D; d++) { + int x₁ = vertices[i₁].coordinates[d]; + int x₂ = vertices[i₂].coordinates[d]; + unsigned Δx = abs(x₁ - x₂); if (Δx > L / 2) { Δx = L - Δx; } - sd += Δx * Δx; + Δx² += Δx * Δx; } - return sd; + return Δx²; } - /* Initialize a square lattice */ Graph(unsigned D, unsigned L) : D(D), L(L), vertices(uPow(L, D)), multiplicities((D * L * L) / 4) { for (unsigned i = 0; i < Nv(); i++) { - vertices[i].coordinates.resize(D); vertices[i].index = i; - vertices[i].bonds.reserve(2 * D); - for (unsigned j = 0; j < D; j++) { - vertices[i].coordinates[j] = (i / uPow(L, j)) % L; + vertices[i].coordinates.reserve(D); + vertices[i].outgoingEdges.reserve(2 * D); + for (unsigned d = 0; d < D; d++) { + vertices[i].coordinates.push_back((i / uPow(L, d)) % L); - unsigned n1 = uPow(L, j + 1) * (i / uPow(L, j + 1)) + (i + uPow(L, j)) % uPow(L, j + 1); - unsigned n2 = uPow(L, j + 1) * (i / uPow(L, j + 1)) + (uPow(L, j + 1) + i - uPow(L, j)) % uPow(L, j + 1); + unsigned j₁ = uPow(L, d + 1) * (i / uPow(L, d + 1)) + (i + uPow(L, d)) % uPow(L, d + 1); + unsigned j₂ = uPow(L, d + 1) * (i / uPow(L, d + 1)) + (uPow(L, d + 1) + i - uPow(L, d)) % uPow(L, d + 1); - unsigned e1 = j * Nv() + i; - unsigned e2 = j * Nv() + n2; + unsigned e₁ = d * Nv() + i; + unsigned e₂ = d * Nv() + j₂; - HalfEdge he1(e1, vertices[n1]); - HalfEdge he2(e2, vertices[n2]); + HalfEdge he₁(e₁, vertices[j₁]); + HalfEdge he₂(e₂, vertices[j₂]); - vertices[i].bonds.push_back(he1); - vertices[i].bonds.push_back(he2); + vertices[i].outgoingEdges.push_back(he₁); + vertices[i].outgoingEdges.push_back(he₂); } } - for (const Vertex& v1 : vertices) { - for (const Vertex& v2 : vertices) { - unsigned dist = squaredDistance(v1.index, v2.index); - if (dist > 0) { - multiplicities[dist - 1]++; + for (const Vertex& v₁ : vertices) { + for (const Vertex& v₂ : vertices) { + unsigned Δx² = squaredDistance(v₁.index, v₂.index); + if (Δx² > 0) { + multiplicities[Δx² - 1]++; } } } } - std::vector markClusters(const std::vector& activeEdges) const { - std::vector clusters; - std::vector indicies(Nv()); + std::list findClusters(const std::vector& activeEdges) const { + std::list clusters; + std::vector clusterIndicies(Nv()); unsigned nClusters = 0; - for (const Vertex& v : vertices) { - if (indicies[v.index] == 0) { + for (const Vertex& v₀ : vertices) { + if (clusterIndicies[v₀.index] == 0) { nClusters++; clusters.push_back({}); - std::stack> inCluster; - inCluster.push(v); - while (!inCluster.empty()) { - const Vertex& vn = inCluster.top(); - inCluster.pop(); - if (indicies[vn.index] == 0) { - indicies[vn.index] = nClusters; - clusters[nClusters - 1].push_back(vn); - } - for (const HalfEdge& e : vn.bonds) { - if (activeEdges[e.index] && indicies[e.neighbor.index] == 0) { - inCluster.push(e.neighbor); + std::stack> verticiesToProcess; + verticiesToProcess.push(v₀); + while (!verticiesToProcess.empty()) { + const Vertex& v₁ = verticiesToProcess.top(); + verticiesToProcess.pop(); + if (clusterIndicies[v₁.index] == 0) { + clusterIndicies[v₁.index] = nClusters; + clusters.back().push_back(v₁); + for (const HalfEdge& e : v₁.outgoingEdges) { + if (activeEdges[e.index] && clusterIndicies[e.neighbor.index] == 0) { + verticiesToProcess.push(e.neighbor); + } } } } @@ -134,11 +135,11 @@ public: return clusters; } - Matrix laplacian(const std::vector& activeEdges) const { + Matrix laplacian(const std::vector& activeEdges) const { Matrix L = Matrix::Zero(Nv(), Nv()); for (const Vertex& v : vertices) { - for (const HalfEdge& e : v.bonds) { + for (const HalfEdge& e : v.outgoingEdges) { if (activeEdges[e.index]) { L(v.index, v.index) += 1; L(v.index, e.neighbor.index) -= 1; @@ -153,14 +154,15 @@ public: int main(int argc, char* argv[]) { Rng r; - unsigned D = 2; - unsigned L = 8; - double p = 0.5; - unsigned n = 10; + unsigned D = 2; // Dimension of space + unsigned L = 8; // Size of lattice + Real p = 0.5; // Probability of bond + unsigned M = 1; // Moments of conductivity + unsigned n = 1; // Number of experiments int opt; - while ((opt = getopt(argc, argv, "D:L:p:n:")) != -1) { + while ((opt = getopt(argc, argv, "D:L:p:n:M:")) != -1) { switch (opt) { case 'D': D = atoi(optarg); @@ -171,6 +173,9 @@ int main(int argc, char* argv[]) { case 'p': p = atof(optarg); break; + case 'M': + M = atoi(optarg); + break; case 'n': n = (unsigned)atof(optarg); break; @@ -181,59 +186,59 @@ int main(int argc, char* argv[]) { Graph G(D, L); - for (unsigned trials = 0; trials < n; trials++) { - std::vector activeEdges(G.Ne()); - for (short unsigned& edge : activeEdges) { + for (unsigned trial = 0; trial < n; trial++) { + std::vector activeEdges(G.Ne()); + for (Bool& edge : activeEdges) { edge = r.uniform(0.0, 1.0) < p; } - std::vector clusters = G.markClusters(activeEdges); + + std::list clusters = G.findClusters(activeEdges); Matrix laplacian = G.laplacian(activeEdges); - std::vector conductivities((L * L * D) / 4); - std::vector probabilities((L * L * D) / 4); + std::vector> conductivities(M + 1); + + for (std::vector& dat : conductivities) { + dat.resize((L * L * D) / 4); + } - for (const Cluster& c : clusters) { - std::vector inds; - inds.reserve(c.size()); - for (const Vertex& v : c){ - inds.push_back(v.index); + for (const Cluster& cluster : clusters) { + std::vector clusterIndicies; + clusterIndicies.reserve(cluster.size()); + for (const Vertex& v : cluster){ + clusterIndicies.push_back(v.index); } - Matrix subLaplacian = laplacian(inds, inds); + Matrix subLaplacian = laplacian(clusterIndicies, clusterIndicies); subLaplacian(0,0)++; /* Set voltage of first node to zero */ Matrix subLaplacianInv = subLaplacian.inverse(); - for (unsigned i = 0; i < c.size(); i++) { + for (unsigned i = 0; i < cluster.size(); i++) { for (unsigned j = 0; j < i; j++) { - Vector input = Vector::Zero(inds.size()); + Vector input = Vector::Zero(clusterIndicies.size()); input(i) = 1; input(j) = -1; Vector output = subLaplacianInv * input; - double ΔV = std::abs(output(i) - output(j)); + Real ΔV = std::abs(output(i) - output(j)); - conductivities[G.squaredDistance(inds[i], inds[j]) - 1] += 1 / ΔV; - probabilities[G.squaredDistance(inds[i], inds[j]) - 1]++; + for (unsigned m = 0; m <= M; m++) { + conductivities[m][G.squaredDistance(clusterIndicies[i], clusterIndicies[j]) - 1] += pow(1 / ΔV, m); + } } } } - for (unsigned i = 0; i < conductivities.size(); i++) { - if (G.multiplicities[i] != 0) { - std::cout << conductivities[i] / G.multiplicities[i] << " "; - } else { - std::cout << 0 << " "; - } - } - std::cout << std::endl; - for (unsigned i = 0; i < probabilities.size(); i++) { - if (G.multiplicities[i] != 0) { - std::cout << ((Real)probabilities[i]) / G.multiplicities[i] << " "; - } else { - std::cout << 0 << " "; + for (unsigned m = 0; m <= M; m++) { + std::ofstream file("perc_" + std::to_string(D) + "_" + std::to_string(L) + "_" + std::to_string(p) + "_" + std::to_string(m) + ".dat", std::ios_base::app); + for (unsigned i = 0; i < conductivities[m].size(); i++) { + if (G.multiplicities[i] != 0) { + file << conductivities[m][i] / G.multiplicities[i] << " "; + } else { + file << 0 << " "; + } } + file << std::endl; } - std::cout << std::endl; } return 0; -- cgit v1.2.3-70-g09d2