summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2019-05-17 16:38:24 -0400
committerJaron Kent-Dobias <jaron@kent-dobias.com>2019-05-17 16:38:24 -0400
commit1e939e597964fa081b347e40af2be1069867b906 (patch)
tree395de78337d7ed85971fa3e84860ee4f4e5a5956 /src
parent456d748b8cbf141e1fe6e82a7c0106f469e28ba6 (diff)
downloadfuse_networks-1e939e597964fa081b347e40af2be1069867b906.tar.gz
fuse_networks-1e939e597964fa081b347e40af2be1069867b906.tar.bz2
fuse_networks-1e939e597964fa081b347e40af2be1069867b906.zip
added support for running percolation using the fuse network rules
Diffstat (limited to 'src')
-rw-r--r--src/CMakeLists.txt5
-rw-r--r--src/perc_meas.cpp167
-rw-r--r--src/perc_meas.hpp52
-rw-r--r--src/percolation.cpp100
4 files changed, 324 insertions, 0 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index eecbf72..9a7ce20 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -5,6 +5,9 @@ target_link_libraries(analysis_tools frac)
add_library(measurements measurements.cpp)
target_link_libraries(measurements frac analysis_tools)
+add_library(perc_meas perc_meas.cpp)
+target_link_libraries(perc_meas frac analysis_tools)
+
add_library(animate animate.cpp)
target_link_libraries(animate frac analysis_tools GL GLU glut)
@@ -29,3 +32,5 @@ target_link_libraries(sample_fracture frac sample cholmod)
add_executable(sample_fracture_square sample_fracture_square.cpp)
target_link_libraries(sample_fracture_square frac sample cholmod)
+add_executable(percolation percolation.cpp)
+target_link_libraries(percolation frac perc_meas fftw3 cholmod profiler GL GLU glut)
diff --git a/src/perc_meas.cpp b/src/perc_meas.cpp
new file mode 100644
index 0000000..4220a93
--- /dev/null
+++ b/src/perc_meas.cpp
@@ -0,0 +1,167 @@
+
+#include "perc_meas.hpp"
+#include <iostream>
+#include <cstdio>
+
+class badcycleException: public std::exception
+{
+ virtual const char* what() const throw()
+ {
+ return "Could not find a valid cycle on the broken system.";
+ }
+} badcycleex;
+
+void update_distribution_file(std::string id, const std::vector<uint64_t>& data, std::string model_string) {
+ std::string filename = model_string + id + ".dat";
+ std::ifstream file(filename);
+
+ std::vector<uint64_t> data_old(data.size(), 0);
+
+ if (file.is_open()) {
+ for (unsigned i = 0; i < data.size(); i++) {
+ uint64_t num;
+ file >> num;
+ data_old[i] = num;
+ }
+
+ file.close();
+ }
+
+ std::ofstream file_out(filename);
+
+ for (unsigned i = 0; i < data.size(); i++) {
+ file_out <<std::fixed<< data_old[i] + data[i] << " ";
+ }
+
+ file_out.close();
+}
+
+void update_distribution_file(std::string id, const std::vector<std::vector<uint64_t>>& data, std::string model_string) {
+ std::string filename = model_string + id + ".dat";
+ std::ifstream file(filename);
+
+ std::vector<std::vector<uint64_t>> data_old(data.size());
+
+ for (unsigned i = 0; i < data.size(); i++) {
+ data_old[i].resize(data[0].size());
+ }
+
+ if (file.is_open()) {
+ for (unsigned i = 0; i < data.size(); i++) {
+ for (unsigned j = 0; j < data[0].size(); j++) {
+ uint64_t num;
+ file >> num;
+ data_old[i][j] = num;
+ }
+ }
+
+ file.close();
+ }
+
+ std::ofstream file_out(filename);
+
+ for (unsigned i = 0; i < data.size(); i++) {
+ for (unsigned j = 0; j < data[0].size(); j++) {
+ file_out << std::fixed << data_old[i][j] + data[i][j] << " ";
+ }
+ file_out << "\n";
+ }
+
+ file_out.close();
+}
+
+pm::pm(unsigned n, double a) :
+ G(2 * n),
+ rank(2 * n),
+ parent(2 * n),
+ ds(&rank[0], &parent[0]),
+ sm(2 * n),
+ sl(2 * n),
+ sn(3 * n),
+ ss(2 * n),
+ sb(3 * n),
+ sp(3 * n)
+{
+ model_string = "percolation_" + std::to_string(n) + "_" + std::to_string(a) + "_";
+ for (std::vector<uint64_t> &x : sn) {
+ x.resize(2 * n);
+ }
+ for (std::vector<uint64_t> &x : sb) {
+ x.resize(3 * n);
+ }
+}
+
+pm::~pm() {
+ update_distribution_file("sm", sm, model_string);
+ update_distribution_file("sl", sl, model_string);
+ update_distribution_file("sn", sn, model_string);
+ update_distribution_file("ss", ss, model_string);
+ update_distribution_file("sb", sb, model_string);
+ update_distribution_file("sp", sp, model_string);
+}
+
+void pm::pre_fracture(const network&) {
+ boost::remove_edge_if(trivial, G);
+ initialize_incremental_components(G, ds);
+ incremental_components(G, ds);
+ p = 0;
+}
+
+void pm::bond_broken(const network& net, const current_info& cur, unsigned i) {
+ boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
+ ds.union_set(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1]);
+
+ boost::component_index<VertexIndex> components(parent.begin(), parent.end());
+ std::vector<unsigned> counts(components.size());
+
+ for (unsigned i = 0; i < components.size(); i++) {
+ unsigned comp_size = 0;
+ BOOST_FOREACH(VertexIndex child_index, components[i]) {
+ comp_size++;
+ }
+
+ sn[p][comp_size - 1]++;
+ }
+
+ unsigned bb_size = 0;
+
+ for (unsigned i = 0; i <net.G.edges.size(); i++) {
+ if (cur.currents[i] > 1.0 / net.G.edges.size()) {
+ bb_size++;
+ }
+ }
+
+ sb[p][bb_size - 1]++;
+ p++;
+}
+
+void pm::post_fracture(network &n) {
+ auto post_cracks = find_minimal_crack(G, n);
+ std::vector<unsigned> component(boost::num_vertices(G));
+ unsigned num = boost::connected_components(G, &component[0]);
+ if (post_cracks.size() > 2 || post_cracks.size() == 0) {
+ throw badcycleex;
+ }
+ for (auto c : post_cracks) {
+ sl[c.second.size() - 1]++;
+ }
+ unsigned crack_component = component[n.G.dual_edges[post_cracks.front().second.front()].v[0]];
+
+ std::vector<std::list<graph::coordinate>> components(num);
+
+ for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
+ components[component[i]].push_back(n.G.dual_vertices[i].r);
+ }
+
+ for (unsigned i = 0; i < num; i++) {
+ if (i != crack_component) {
+ sm[components[i].size() - 1]++;
+ } else {
+ ss[components[i].size() - 1]++;
+ }
+ sn[p][components[i].size() - 1]++;
+ }
+
+ sp[p - 1]++;
+}
+
diff --git a/src/perc_meas.hpp b/src/perc_meas.hpp
new file mode 100644
index 0000000..b85e27e
--- /dev/null
+++ b/src/perc_meas.hpp
@@ -0,0 +1,52 @@
+
+#include <cstring>
+#include <fstream>
+#include <string>
+#include <cinttypes>
+#include <sstream>
+#include <functional>
+#include <iostream>
+#include <array>
+
+#include <fftw3.h>
+
+#include <hooks.hpp>
+#include <boost/graph/graph_utility.hpp>
+#include <boost/graph/incremental_components.hpp>
+#include <boost/pending/disjoint_sets.hpp>
+#include <boost/foreach.hpp>
+
+#include "analysis_tools.hpp"
+typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
+typedef boost::graph_traits<Graph>::vertices_size_type VertexIndex;
+
+
+class pm : public hooks {
+ private:
+ unsigned p;
+ Graph G;
+ std::vector<VertexIndex> rank;
+ std::vector<Vertex> parent;
+ boost::disjoint_sets<VertexIndex*, Vertex*> ds;
+
+ // measurement storage
+ std::vector<uint64_t> sc; // non-spanning cluster size distribution
+ std::vector<std::vector<uint64_t>> sn; // non-spanning cluster size distribution
+ std::vector<uint64_t> ss; // minimal spanning cluster size distribution
+ std::vector<uint64_t> sm; // spanning cluster size distribution
+ std::vector<uint64_t> sl; // final avalanche size distribution
+ std::vector<std::vector<uint64_t>> sb; // final avalanche size distribution
+ std::vector<uint64_t> sp;
+
+ public:
+ std::string model_string;
+
+ pm(unsigned Lx, unsigned Ly);
+ pm(unsigned n, double a);
+ ~pm();
+
+ void pre_fracture(const network &) override;
+ void bond_broken(const network& net, const current_info& cur, unsigned i) override;
+ void post_fracture(network &n) override;
+};
+
diff --git a/src/percolation.cpp b/src/percolation.cpp
new file mode 100644
index 0000000..03da117
--- /dev/null
+++ b/src/percolation.cpp
@@ -0,0 +1,100 @@
+
+#include <random>
+#include <iostream>
+
+#include <cholmod.h>
+
+#include "randutils/randutils.hpp"
+
+#include <graph.hpp>
+#include <network.hpp>
+#include <hooks.hpp>
+#include "perc_meas.hpp"
+
+#include <csignal>
+#include <cstring>
+#include <atomic>
+
+std::atomic<bool> quit(false); // signal flag
+
+void got_signal(int) {
+ quit.store(true);
+}
+
+int main(int argc, char* argv[]) {
+ struct sigaction sa;
+ memset( &sa, 0, sizeof(sa) );
+ sa.sa_handler = got_signal;
+ sigfillset(&sa.sa_mask);
+ sigaction(SIGINT, &sa, NULL);
+
+ int opt;
+
+ unsigned N = 1;
+ double Lx = 16;
+ double Ly = 16;
+ double beta = 0.0001;
+
+ unsigned n = 128;
+ double a = 1.0;
+ bool use_aN = false;
+
+ while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:")) != -1) {
+ switch (opt) {
+ case 'N':
+ N = (unsigned)atof(optarg);
+ break;
+ case 'X':
+ Lx = atof(optarg);
+ break;
+ case 'Y':
+ Ly = atof(optarg);
+ break;
+ case 'b':
+ beta = atof(optarg);
+ break;
+ case 'n':
+ n = (unsigned)atof(optarg);
+ use_aN = true;
+ break;
+ case 'a':
+ a = atof(optarg);
+ use_aN = true;
+ break;
+ default:
+ exit(1);
+ }
+ }
+
+ cholmod_common c;
+ CHOL_F(start)(&c);
+
+ randutils::auto_seed_128 seeds;
+ std::mt19937 rng{seeds};
+
+ if (use_aN) {
+ pm meas(n, a);
+
+ for (unsigned trial = 0; trial < N; trial++) {
+ while (true) {
+ try {
+ graph G(n, a, rng);
+ percolation_network fuse_network(G, &c);
+ fuse_network.set_thresholds(beta, rng);
+ fuse_network.fracture(meas);
+ break;
+ } catch (std::exception &e) {
+ std::cout << e.what() << '\n';
+ }
+ }
+
+ if (quit.load())
+ break;
+ }
+ }
+
+ CHOL_F(finish)(&c);
+
+ return 0;
+}
+