From 1e939e597964fa081b347e40af2be1069867b906 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 17 May 2019 16:38:24 -0400 Subject: added support for running percolation using the fuse network rules --- lib/include/network.hpp | 15 +++++ lib/src/network.cpp | 80 +++++++++++++++++++++++ src/CMakeLists.txt | 5 ++ src/perc_meas.cpp | 167 ++++++++++++++++++++++++++++++++++++++++++++++++ src/perc_meas.hpp | 52 +++++++++++++++ src/percolation.cpp | 100 +++++++++++++++++++++++++++++ 6 files changed, 419 insertions(+) create mode 100644 src/perc_meas.cpp create mode 100644 src/perc_meas.hpp create mode 100644 src/percolation.cpp diff --git a/lib/include/network.hpp b/lib/include/network.hpp index dd11342..15cb8ea 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -82,3 +82,18 @@ class elastic_network : public network { current_info get_current_info() override; }; +class percolation_network : public network { + private: + double weight; + public: + problem hook_x; + problem hook_y; + + percolation_network(const graph&, cholmod_common*); + percolation_network(const percolation_network&); + + void fracture(hooks&, double weight = 0.5, double cutoff = 1e-10); + void break_edge(unsigned, bool unbreak = false) override; + current_info get_current_info() override; +}; + diff --git a/lib/src/network.cpp b/lib/src/network.cpp index ed51c8e..1d9b3f9 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -363,3 +363,83 @@ void elastic_network::break_edge(unsigned e, bool unbreak) { hook_y.break_edge(e, unbreak); } +percolation_network::percolation_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) { +} + +percolation_network::percolation_network(const percolation_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) { +} + +void percolation_network::fracture(hooks& m, double weight, double cutoff) { + this->weight = weight; + m.pre_fracture(*this); + + while (true) { + current_info ctot = this->get_current_info(); + + if (ctot.conductivity == 0) { + break; + } + + unsigned max_pos = UINT_MAX; + long double max_val = std::numeric_limits::lowest(); + + for (unsigned i = 0; i < G.edges.size(); i++) { + if (!fuses[i] && ctot.currents[i] > 1.0 / G.edges.size()) { + long double val = - thresholds[i]; + if (val > max_val) { + max_val = val; + max_pos = i; + } + } + } + + if (max_pos == UINT_MAX) { + throw nofuseex; + } + + this->break_edge(max_pos); + + m.bond_broken(*this, ctot, max_pos); + } + + m.post_fracture(*this); +} + +current_info percolation_network::get_current_info() { + current_info cx = hook_x.solve(fuses); + current_info cy = hook_y.solve(fuses); + + bool done_x = cx.conductivity < 1.0 / G.edges.size(); + bool done_y = cy.conductivity < 1.0 / G.edges.size(); + + current_info ctot; + ctot.currents.resize(G.edges.size()); + + if (done_x && done_y) { + ctot.conductivity = 0; + } else if (done_x) { + ctot.conductivity = cy.conductivity; + for (unsigned i = 0; i < G.edges.size(); i++) { + ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity; + } + } else if (done_y) { + ctot.conductivity = cx.conductivity; + for (unsigned i = 0; i < G.edges.size(); i++) { + ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity; + } + } else { + ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); + for (unsigned i = 0; i < G.edges.size(); i++) { + ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2)); + } + } + + return ctot; +} + +void percolation_network::break_edge(unsigned e, bool unbreak) { + fuses[e] = !unbreak; + hook_x.break_edge(e, unbreak); + hook_y.break_edge(e, unbreak); +} + 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 +#include + +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& data, std::string model_string) { + std::string filename = model_string + id + ".dat"; + std::ifstream file(filename); + + std::vector 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 <>& data, std::string model_string) { + std::string filename = model_string + id + ".dat"; + std::ifstream file(filename); + + std::vector> 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 &x : sn) { + x.resize(2 * n); + } + for (std::vector &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 components(parent.begin(), parent.end()); + std::vector 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 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 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> 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 +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include "analysis_tools.hpp" +typedef boost::graph_traits::vertex_descriptor Vertex; +typedef boost::graph_traits::vertices_size_type VertexIndex; + + +class pm : public hooks { + private: + unsigned p; + Graph G; + std::vector rank; + std::vector parent; + boost::disjoint_sets ds; + + // measurement storage + std::vector sc; // non-spanning cluster size distribution + std::vector> sn; // non-spanning cluster size distribution + std::vector ss; // minimal spanning cluster size distribution + std::vector sm; // spanning cluster size distribution + std::vector sl; // final avalanche size distribution + std::vector> sb; // final avalanche size distribution + std::vector 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 +#include + +#include + +#include "randutils/randutils.hpp" + +#include +#include +#include +#include "perc_meas.hpp" + +#include +#include +#include + +std::atomic 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; +} + -- cgit v1.2.3-54-g00ecf