summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitmodules9
-rw-r--r--compile_flags.txt3
m---------eigen0
m---------pcg-cpp0
-rw-r--r--percolation.cpp234
m---------randutils0
6 files changed, 246 insertions, 0 deletions
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..c1474be
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,9 @@
+[submodule "eigen"]
+ path = eigen
+ url = https://gitlab.com/libeigen/eigen
+[submodule "pcg-cpp"]
+ path = pcg-cpp
+ url = https://github.com/imneme/pcg-cpp
+[submodule "randutils"]
+ path = randutils
+ url = https://gist.github.com/imneme/540829265469e673d045
diff --git a/compile_flags.txt b/compile_flags.txt
new file mode 100644
index 0000000..4cb2489
--- /dev/null
+++ b/compile_flags.txt
@@ -0,0 +1,3 @@
+-std=c++17
+-Wno-mathematical-notation-identifier-extension
+-Wall
diff --git a/eigen b/eigen
new file mode 160000
+Subproject 3147391d946bb4b6c68edd901f2add6ac1f31f8
diff --git a/pcg-cpp b/pcg-cpp
new file mode 160000
+Subproject 428802d1a5634f96bcd0705fab379ff0113bcf1
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 <stack>
+#include <list>
+#include <cstdlib>
+#include <getopt.h>
+
+#include "eigen/Eigen/Dense"
+#include "randutils/randutils.hpp"
+#include "pcg-cpp/include/pcg_random.hpp"
+
+using Rng = randutils::random_generator<pcg32>;
+
+using Real = double;
+using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
+using Matrix = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
+
+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<HalfEdge> bonds;
+ unsigned index;
+ unsigned cluster;
+};
+
+class Edge {
+public:
+ bool active;
+};
+
+class Cluster : public std::list<std::reference_wrapper<Vertex>> {};
+class Clusters {
+public:
+ std::vector<Cluster> lists;
+ std::vector<unsigned> indices;
+};
+
+class Graph {
+public:
+ unsigned D;
+ unsigned L;
+ std::vector<Vertex> vertices;
+ std::vector<Edge> edges;
+ std::vector<Cluster> 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<Real, std::uniform_real_distribution>(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<std::reference_wrapper<Vertex>> 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<unsigned> 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<Real> conductivities(L * L * D / 4);
+
+ for (const Cluster& c : G.clusters) {
+ std::vector<unsigned> 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;
+}
diff --git a/randutils b/randutils
new file mode 160000
+Subproject 5afb4439c23a6c060eda72121bff8bf9da59591