diff options
-rw-r--r-- | CMakeLists.txt | 6 | ||||
-rw-r--r-- | lib/include/current_info.hpp | 10 | ||||
-rw-r--r-- | lib/include/hooks.hpp | 4 | ||||
-rw-r--r-- | lib/include/network.hpp | 8 | ||||
-rw-r--r-- | lib/src/network.cpp | 82 | ||||
-rw-r--r-- | src/fracture.cpp | 23 | ||||
-rw-r--r-- | src/measurements.cpp | 15 | ||||
-rw-r--r-- | src/measurements.hpp | 10 |
8 files changed, 69 insertions, 89 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 6746aa9..a804ac5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,8 +1,8 @@ -cmake_minimum_required(VERSION 3.9) +cmake_minimum_required(VERSION 3.8) -set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall") -set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto") +set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall -stdlib=libc++") +set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -stdlib=libc++") set (CMAKE_CXX_STANDARD 17) set (CMAKE_C_STANDARD 11) diff --git a/lib/include/current_info.hpp b/lib/include/current_info.hpp new file mode 100644 index 0000000..90dd3c7 --- /dev/null +++ b/lib/include/current_info.hpp @@ -0,0 +1,10 @@ + +#pragma once + +#include <vector> + +typedef struct current_info { + double conductivity; + std::vector<double> currents; +} current_info; + diff --git a/lib/include/hooks.hpp b/lib/include/hooks.hpp index 350f318..e617ac3 100644 --- a/lib/include/hooks.hpp +++ b/lib/include/hooks.hpp @@ -1,12 +1,14 @@ #pragma once +#include "current_info.hpp" + class network; class hooks { public: virtual void pre_fracture(const network&) {}; - virtual void bond_broken(const network&, const std::pair<double, std::vector<double>>&, unsigned int) {}; + virtual void bond_broken(const network&, const current_info&, unsigned int) {}; virtual void post_fracture(network&) {}; // post fracture hook can be destructive }; diff --git a/lib/include/network.hpp b/lib/include/network.hpp index 5cbba0d..e6d67b9 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -6,10 +6,11 @@ #include <utility> #include <random> -#include "cholmod.h" +#include <cholmod.h> #include "graph.hpp" #include "hooks.hpp" +#include "current_info.hpp" #ifdef FRACTURE_LONGINT @@ -40,9 +41,8 @@ class network { ~network(); void set_thresholds(double, std::mt19937&); - void break_edge(unsigned int); - void add_edge(unsigned int); - std::pair<double, std::vector<double>> currents(); + void break_edge(unsigned int, bool unbreak = false); + current_info get_current_info(); void fracture(hooks&, double cutoff = 1e-14); }; diff --git a/lib/src/network.cpp b/lib/src/network.cpp index a27ca39..9f6cb0d 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -56,7 +56,6 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges. t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c); - t->nnz = 2 * G.edges.size(); for (unsigned int i = 0; i < G.edges.size(); i++) { @@ -89,68 +88,17 @@ network::~network() { void network::set_thresholds(double beta, std::mt19937& rng) { std::uniform_real_distribution<long double> d(0.0, 1.0); - for (unsigned int i = 0; i < G.edges.size(); i++) { - long double x = 0.0; + for (long double& threshold : thresholds) { + threshold = 0.0; - while (x == 0.0) { - x = expl(logl(d(rng)) / (long double)beta); + while (threshold == 0.0) { + threshold = expl(logl(d(rng)) / (long double)beta); } - - thresholds[i] = x; - } -} - -void network::add_edge(unsigned int e) { - fuses[e] = false; - unsigned int v0 = G.edges[e][0]; - unsigned int v1 = G.edges[e][1]; - - unsigned int n = factor->n; - - cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); - - unsigned int s1, s2; - s1 = v0 < v1 ? v0 : v1; - s2 = v0 < v1 ? v1 : v0; - - CHOL_INT *pp = (CHOL_INT *)update_mat->p; - CHOL_INT *ii = (CHOL_INT *)update_mat->i; - double *xx = (double *)update_mat->x; - - for (unsigned int i = 0; i <= s1; i++) { - pp[i] = 0; - } - - for (unsigned int i = s1 + 1; i <= n; i++) { - pp[i] = 2; - } - - ii[0] = s1; - ii[1] = s2; - xx[0] = -1; - xx[1] = 1; - - cholmod_sparse *perm_update_mat = CHOL_F(submatrix)( - update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c); - - CHOL_F(updown)(false, perm_update_mat, factor, c); - - CHOL_F(free_sparse)(&perm_update_mat, c); - CHOL_F(free_sparse)(&update_mat, c); - - double v0y = G.vertices[v0].r.y; - double v1y = G.vertices[v1].r.y; - - if (fabs(v0y - v1y) > 0.5) { - bool ind = v0y < v1y ? 0 : 1; - - ((double *)b->x)[G.edges[e][ind]] += 1.0; - ((double *)b->x)[G.edges[e][!ind]] -= 1.0; } } -void network::break_edge(unsigned int e) { - fuses[e] = true; +void network::break_edge(unsigned int e, bool unbreak) { + fuses[e] = !unbreak; unsigned int v0 = G.edges[e][0]; unsigned int v1 = G.edges[e][1]; @@ -182,7 +130,7 @@ void network::break_edge(unsigned int e) { cholmod_sparse *perm_update_mat = CHOL_F(submatrix)( update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c); - CHOL_F(updown)(false, perm_update_mat, factor, c); + CHOL_F(updown)(unbreak, perm_update_mat, factor, c); CHOL_F(free_sparse)(&perm_update_mat, c); CHOL_F(free_sparse)(&update_mat, c); @@ -191,14 +139,14 @@ void network::break_edge(unsigned int e) { double v1y = G.vertices[v1].r.y; if (fabs(v0y - v1y) > 0.5) { - bool ind = v0y < v1y ? 0 : 1; + bool ind = v0y < v1y ? unbreak : !unbreak; ((double *)b->x)[G.edges[e][ind]] -= 1.0; ((double *)b->x)[G.edges[e][!ind]] += 1.0; } } -std::pair<double, std::vector<double>> network::currents() { +current_info network::get_current_info() { cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); if (((double *)x->x)[0] != ((double *)x->x)[0]) { @@ -240,16 +188,16 @@ std::pair<double, std::vector<double>> network::currents() { CHOL_F(free_dense)(&x, c); CHOL_F(free_dense)(&y, c); - return std::make_pair(total_current, currents); + return {total_current, currents}; } void network::fracture(hooks& m, double cutoff) { m.pre_fracture(*this); while (true) { - auto currents = this->currents(); + current_info ci = this->get_current_info(); - if (currents.first < cutoff * G.vertices.size()) { + if (ci.conductivity < cutoff * G.vertices.size()) { break; } @@ -257,8 +205,8 @@ void network::fracture(hooks& m, double cutoff) { long double max_val = 0; for (unsigned int i = 0; i < G.edges.size(); i++) { - if (!fuses[i] && fabs(currents.second[i]) > cutoff) { - long double val = (long double)fabs(currents.second[i]) / thresholds[i]; + if (!fuses[i] && fabs(ci.currents[i]) > cutoff) { + long double val = (long double)fabs(ci.currents[i]) / thresholds[i]; if (val > max_val) { max_val = val; max_pos = i; @@ -272,7 +220,7 @@ void network::fracture(hooks& m, double cutoff) { this->break_edge(max_pos); - m.bond_broken(*this, currents, max_pos); + m.bond_broken(*this, ci, max_pos); } m.post_fracture(*this); diff --git a/src/fracture.cpp b/src/fracture.cpp index 77af253..0aad56e 100644 --- a/src/fracture.cpp +++ b/src/fracture.cpp @@ -11,7 +11,23 @@ #include <hooks.hpp> #include "measurements.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 int N = 1; @@ -36,12 +52,10 @@ int main(int argc, char* argv[]) { cholmod_common c; CHOL_F(start)(&c); - c.supernodal = CHOLMOD_SUPERNODAL; - graph G(L); network base_network(G, &c); - ma meas(N, L, beta); + ma meas(L, beta); randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; @@ -50,6 +64,9 @@ int main(int argc, char* argv[]) { network tmp_network(base_network); tmp_network.set_thresholds(beta, rng); tmp_network.fracture(meas); + + if (quit.load()) + break; } CHOL_F(finish)(&c); diff --git a/src/measurements.cpp b/src/measurements.cpp index 0bd72c4..ed96855 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -38,8 +38,8 @@ std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) { class find_cycle : public boost::default_dfs_visitor { public: - unsigned int end; std::list<unsigned int>& E; + unsigned int end; struct done{}; find_cycle(std::list<unsigned int>& E, unsigned int end) : E(E), end(end) {} @@ -226,7 +226,7 @@ void autocorrelation(unsigned int L, std::vector<T>& out_data, fftw_plan forward } } -ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2)), N(N), beta(beta), +ma::ma(unsigned int L, double beta) : L(L), beta(beta), G(2 * pow(L / 2, 2)), sc(pow(L, 2), 0), sa(pow(L, 2), 0), sd(pow(L, 2), 0), @@ -239,6 +239,7 @@ ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2)) Cll(pow(L / 2 + 1, 2), 0), Cee(pow(L / 2 + 1, 2), 0) { + N = 0; Nc = 0; Na = 0; @@ -284,8 +285,9 @@ void ma::pre_fracture(const network &) { boost::remove_edge_if(trivial, G); } -void ma::bond_broken(const network& net, const std::pair<double, std::vector<double>>& cur, unsigned int i) { - if (cur.first / fabs(cur.second[i]) * net.thresholds[i] > lv) { +void ma::bond_broken(const network& net, const current_info& cur, unsigned int i) { + double c = cur.conductivity / fabs(cur.currents[i]) * net.thresholds[i]; + if (c > lv) { sa[avalanches.back().size()]++; Na++; @@ -297,7 +299,7 @@ void ma::bond_broken(const network& net, const std::pair<double, std::vector<dou autocorrelation(L, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); - lv = cur.first / fabs(cur.second[i]) * net.thresholds[i]; + lv = c; avalanches.push_back({i}); } else { avalanches.back().push_back(i); @@ -373,7 +375,7 @@ void ma::post_fracture(network &n) { // rewind the last avalanche for (auto e : avalanches.back()) { boost::remove_edge(n.G.dual_edges[e][0], n.G.dual_edges[e][1], G); - n.add_edge(e); + n.break_edge(e, true); fftw_forward_in[e] = 1; } @@ -417,5 +419,6 @@ void ma::post_fracture(network &n) { sd[total_broken]++; + N++; } diff --git a/src/measurements.hpp b/src/measurements.hpp index 621c425..ab31ccd 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -41,6 +41,9 @@ class ma : public hooks { double *fftw_reverse_out; fftw_plan forward_plan; fftw_plan reverse_plan; + unsigned int N; + unsigned int L; + double beta; Graph G; // measurement storage @@ -59,19 +62,16 @@ class ma : public hooks { uint64_t Na; public: - unsigned int N; - unsigned int L; double lv; - double beta; std::list<std::list<unsigned int>> avalanches; - ma(unsigned int N, unsigned int L, double beta); + ma(unsigned int L, double beta); ~ma(); void pre_fracture(const network &) override; - void bond_broken(const network& net, const std::pair<double, std::vector<double>>& cur, unsigned int i) override; + void bond_broken(const network& net, const current_info& cur, unsigned int i) override; void post_fracture(network &n) override; }; |