diff options
| -rw-r--r-- | src/CMakeLists.txt | 9 | ||||
| -rw-r--r-- | src/fracture.cpp | 6 | ||||
| -rw-r--r-- | src/fracture_square.cpp | 2 | ||||
| -rw-r--r-- | src/measurements.cpp | 173 | ||||
| -rw-r--r-- | src/measurements.hpp | 22 | ||||
| -rw-r--r-- | src/sample.cpp | 69 | ||||
| -rw-r--r-- | src/sample.hpp | 30 | ||||
| -rw-r--r-- | src/sample_fracture.cpp | 65 | ||||
| -rw-r--r-- | src/sample_fracture_square.cpp | 67 |
9 files changed, 316 insertions, 127 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4d4099f..eecbf72 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,6 +8,9 @@ target_link_libraries(measurements frac analysis_tools) add_library(animate animate.cpp) target_link_libraries(animate frac analysis_tools GL GLU glut) +add_library(sample sample.cpp) +target_link_libraries(sample frac) + add_executable(fracture fracture.cpp) target_link_libraries(fracture frac measurements fftw3 cholmod profiler GL GLU glut) @@ -20,3 +23,9 @@ target_link_libraries(animate_fracture frac animate fftw3 cholmod profiler GL GL add_executable(animate_fracture_square animate_fracture_square.cpp) target_link_libraries(animate_fracture_square frac animate fftw3 cholmod profiler GL GLU glut) +add_executable(sample_fracture sample_fracture.cpp) +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) + diff --git a/src/fracture.cpp b/src/fracture.cpp index 595cf1d..ebacd28 100644 --- a/src/fracture.cpp +++ b/src/fracture.cpp @@ -61,10 +61,10 @@ int main(int argc, char* argv[]) { unsigned Mx = pow(2, ceil(log2(4*Lx))); unsigned My = pow(2, ceil(log2(4*Ly))); - ma meas(Lx, Ly, Mx, My, beta, 4); + ma meas(Lx, Ly, Mx, My, beta); - randutils::auto_seed_128 seeds; - std::mt19937 rng{seeds}; + //randutils::auto_seed_128 seeds; + std::mt19937 rng{683129629}; for (unsigned trial = 0; trial < N; trial++) { while (true) { diff --git a/src/fracture_square.cpp b/src/fracture_square.cpp index 47d7d67..dec4a2c 100644 --- a/src/fracture_square.cpp +++ b/src/fracture_square.cpp @@ -57,7 +57,7 @@ int main(int argc, char* argv[]) { cholmod_common c; CHOL_F(start)(&c); - ma meas(Lx, Ly, 2*Lx, 2*Ly, beta, 10); + ma meas(Lx, Ly, 2*Lx, 2*Ly, beta); graph G(Lx, Ly); network perm_network(G, &c); diff --git a/src/measurements.cpp b/src/measurements.cpp index ff217b2..ed4b8eb 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -37,17 +37,12 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned N, d std::ifstream file(filename); uint64_t N_old = 0; - std::vector<std::vector<uint64_t>> data_old(data.size()); - for (unsigned j = 0; j < data.size(); j++) { - data_old[j].resize(data[j].size(), 0); - } + std::vector<T> data_old(data.size(), 0); if (file.is_open()) { file >> N_old; for (unsigned j = 0; j < data.size(); j++) { - for (unsigned i = 0; i < data[j].size(); i++) { - file >> data_old[j][i]; - } + file >> data_old[j]; } file.close(); } @@ -56,10 +51,7 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned N, d file_out <<std::fixed<< N_old + N << "\n"; for (unsigned j = 0; j < data.size(); j++) { - for (unsigned i = 0; i < data[j].size(); i++) { - file_out << data_old[j][i] + data[j][i] << " "; - } - file_out << "\n"; + file_out << data_old[j] + data[j] << " "; } file_out.close(); @@ -94,6 +86,34 @@ void correlation(unsigned Mx, unsigned My, std::vector<std::vector<T>>& data, co } } +void autocorrelation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data, const std::list<graph::coordinate>& pos) { + for (std::list<graph::coordinate>::const_iterator it1 = pos.begin(); it1 != pos.end(); it1++) { + for (std::list<graph::coordinate>::const_iterator it2 = it1; it2 != pos.end(); it2++) { + double Δx_tmp = fabs(it1->x - it2->x); + double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp; + + double Δy_tmp = fabs(it1->y - it2->y); + double Δy = Δy_tmp < Ly / 2 ? Δy_tmp : Ly - Δy_tmp; + + data[floor((Mx * Δx) / Lx) + (Mx / 2) * floor((My * Δy) / Ly)]++; + } + } +} + +void correlation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data, const std::list<graph::coordinate>& pos1, const std::list<graph::coordinate>& pos2) { + for (std::list<graph::coordinate>::const_iterator it1 = pos1.begin(); it1 != pos1.end(); it1++) { + for (std::list<graph::coordinate>::const_iterator it2 = pos2.begin(); it2 != pos2.end(); it2++) { + double Δx_tmp = fabs(it1->x - it2->x); + double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp; + + double Δy_tmp = fabs(it1->y - it2->y); + double Δy = Δy_tmp < Ly / 2 ? Δy_tmp : Ly - Δy_tmp; + + data[floor((Mx * Δx) / Lx) + (Mx / 2) * floor((My * Δy) / Ly)]++; + } + } +} + template <class T> void autocorrelation(unsigned Mx, unsigned My, std::vector<std::vector<T>>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) { fftw_execute(forward_plan); @@ -116,39 +136,28 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly); } -ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncum) : +ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta) : Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(2 * (unsigned)ceil(Lx * Ly / 2)), sc(3 * (unsigned)ceil(Lx * Ly / 2), 0), ss(3 * (unsigned)ceil(Lx * Ly / 2), 0), sm(3 * (unsigned)ceil(Lx * Ly / 2), 0), sa(3 * (unsigned)ceil(Lx * Ly / 2), 0), sl(3 * (unsigned)ceil(Lx * Ly / 2), 0), - sd(3 * (unsigned)ceil(Lx * Ly / 2), 0), sD(3 * (unsigned)ceil(Lx * Ly / 2), 0), - ccc(Ncum), - css(Ncum), - cmm(Ncum), - caa(Ncum), - cll(Ncum), - cdd(Ncum), - cDD(Ncum), - csD(Ncum) + ccc((Mx / 2) * (My / 2), 0), + css((Mx / 2) * (My / 2), 0), + cmm((Mx / 2) * (My / 2), 0), + caa((Mx / 2) * (My / 2), 0), + cll((Mx / 2) * (My / 2), 0), + cDD((Mx / 2) * (My / 2), 0), + csD((Mx / 2) * (My / 2), 0) { N = 0; Nc = 0; Na = 0; - for (unsigned i = 0; i < Ncum; i++) { - ccc[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); - css[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); - cmm[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); - caa[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); - cll[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); - cdd[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); - cDD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); - csD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); - } // FFTW setup for correlation functions + /* fftw_set_timelimit(1); fftw_forward_in = (double *)fftw_malloc(Mx * My * sizeof(double)); @@ -158,12 +167,12 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncu forward_plan = fftw_plan_dft_r2c_2d(My, Mx, fftw_forward_in, fftw_forward_out, 0); reverse_plan = fftw_plan_dft_c2r_2d(My, Mx, fftw_reverse_in, fftw_reverse_out, 0); - std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_current.dat"; - //stress_file.open(filename, std::ifstream::app); + */ } ma::~ma() { // clean up FFTW objects + /* fftw_free(fftw_forward_in); fftw_free(fftw_forward_out); fftw_free(fftw_reverse_in); @@ -171,13 +180,13 @@ ma::~ma() { fftw_destroy_plan(forward_plan); fftw_destroy_plan(reverse_plan); fftw_cleanup(); + */ update_distribution_file("sc", sc, Nc, Lx, Ly, beta); update_distribution_file("ss", ss, N, Lx, Ly, beta); update_distribution_file("sm", sm, N, Lx, Ly, beta); update_distribution_file("sa", sa, Na, Lx, Ly, beta); update_distribution_file("sl", sl, N, Lx, Ly, beta); - update_distribution_file("sd", sd, N, Lx, Ly, beta); update_distribution_file("sD", sD, N, Lx, Ly, beta); update_field_file("ccc", ccc, Nc, Lx, Ly, beta, Mx, My); @@ -185,7 +194,6 @@ ma::~ma() { update_field_file("cmm", cmm, N, Lx, Ly, beta, Mx, My); update_field_file("caa", caa, Na, Lx, Ly, beta, Mx, My); update_field_file("cll", cll, N, Lx, Ly, beta, Mx, My); - update_field_file("cdd", cdd, N, Lx, Ly, beta, Mx, My); update_field_file("cDD", cDD, N, Lx, Ly, beta, Mx, My); update_field_file("csD", csD, N, Lx, Ly, beta, Mx, My); @@ -204,19 +212,12 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) { sa[avalanches.back().size() - 1]++; Na++; - std::fill_n(fftw_forward_in, Mx * My, 0.0); - - for (auto e : avalanches.back()) { - unsigned ind = edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My); - fftw_forward_in[ind] += 1.0; - } - - autocorrelation(Mx, My, caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation2(Lx, Ly, Mx, My, caa, avalanches.back()); lv = c; - avalanches.push_back({i}); + avalanches.push_back({net.G.edges[i].r}); } else { - avalanches.back().push_back(i); + avalanches.back().push_back(net.G.edges[i].r); } boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G); @@ -228,110 +229,58 @@ void ma::post_fracture(network &n) { std::list<unsigned> crack = find_minimal_crack(G, n); - // crack surface correlations - std::fill_n(fftw_forward_in, Mx * My, 0.0); - - ss[crack.size()]++; + ss[crack.size() - 1]++; + std::list<graph::coordinate> sr; for (auto edge : crack) { - fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My)] += 0.5; - fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My)] += 0.5; + sr.push_back(n.G.dual_edges[edge].r); } - fftw_complex *tss = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out); - - correlation(Mx, My, css, tss, tss, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation2(Lx, Ly, Mx, My, css, sr); unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]]; - std::vector<std::list<unsigned>> components(num); + std::vector<std::list<graph::coordinate>> components(num); for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) { - components[component[i]].push_back(i); + components[component[i]].push_back(n.G.dual_vertices[i].r); } - // non-spanning clusters - std::fill_n(fftw_forward_in, Mx * My, 0.0); - for (unsigned i = 0; i < num; i++) { if (i != crack_component && components[i].size() > 0) { sc[components[i].size() - 1]++; Nc++; - for (auto it = components[i].begin(); it != components[i].end(); it++) { - unsigned ind = edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My); - fftw_forward_in[ind] += 1.0; - } - - autocorrelation(Mx, My, ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); - - for (auto it = components[i].begin(); it != components[i].end(); it++) { - fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] = 0.0; - } + autocorrelation2(Lx, Ly, Mx, My, ccc, components[i]); } } // spanning cluster - // std::fill_n(fftw_forward_in, Mx * My, 0.0); we already reset in the last loop sm[components[crack_component].size() - 1]++; - for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) { - fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] += 1.0; - } - - autocorrelation(Mx, My, cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation2(Lx, Ly, Mx, My, cmm, components[crack_component]); /// damage at end - std::fill_n(fftw_forward_in, Mx * My, 0.0); - - unsigned final_broken = 0; + std::list<graph::coordinate> Dr; for (unsigned i = 0; i < n.G.edges.size(); i++) { if (n.fuses[i]) { - final_broken++; - fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0; + Dr.push_back(n.G.edges[i].r); } } - sD[final_broken]++; - - fftw_complex *tDD = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out); - - correlation(Mx, My, cDD, tDD, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out); - correlation(Mx, My, csD, tss, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out); + sD[Dr.size() - 1]++; - free(tss); - free(tDD); + autocorrelation2(Lx, Ly, Mx, My, cDD, Dr); + correlation2(Lx, Ly, Mx, My, csD, sr, Dr); - std::fill_n(fftw_forward_in, Mx * My, 0.0); + // ********************** LAST AVALANCHE ********************* // rewind the last avalanche sl[avalanches.back().size() - 1]++; - // rewind the last avalanche - for (auto e : avalanches.back()) { - boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G); - n.break_edge(e, true); - fftw_forward_in[edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0; - } - - autocorrelation(Mx, My, cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); - - // damage size distribution - unsigned total_broken = 0; - - std::fill_n(fftw_forward_in, Mx * My, 0.0); - - for (unsigned i = 0; i < n.G.edges.size(); i++) { - if (n.fuses[i]) { - fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0; - } - } - - autocorrelation(Mx, My, cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); - - sd[total_broken]++; + autocorrelation2(Lx, Ly, Mx, My, cll, avalanches.back()); N++; } diff --git a/src/measurements.hpp b/src/measurements.hpp index e267064..0fd4644 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -19,12 +19,14 @@ class ma : public hooks { // - interface for turning on and off specific measurements // private: + /* double *fftw_forward_in; fftw_complex *fftw_forward_out; fftw_complex *fftw_reverse_in; double *fftw_reverse_out; fftw_plan forward_plan; fftw_plan reverse_plan; + */ unsigned N; double Lx; double Ly; @@ -40,16 +42,14 @@ class ma : public hooks { std::vector<uint64_t> sm; // spanning cluster size distribution std::vector<uint64_t> sa; // non-final avalanche size distribution std::vector<uint64_t> sl; // final avalanche size distribution - std::vector<uint64_t> sd; // pre-fracture damage distribution std::vector<uint64_t> sD; // post-fracture damage distribution - std::vector<std::vector<uint64_t>> ccc; // cluster-cluster correlations - std::vector<std::vector<uint64_t>> css; // surface-surface correlations - std::vector<std::vector<uint64_t>> cmm; // surface-surface correlations - std::vector<std::vector<uint64_t>> caa; // avalanche-avalanche correlations - std::vector<std::vector<uint64_t>> cll; // damage-damage distribution - std::vector<std::vector<uint64_t>> cdd; // damage-damage distribution - std::vector<std::vector<uint64_t>> cDD; // damage-damage distribution - std::vector<std::vector<uint64_t>> csD; // damage-damage distribution + std::vector<uint64_t> ccc; // cluster-cluster correlations + std::vector<uint64_t> css; // surface-surface correlations + std::vector<uint64_t> cmm; // surface-surface correlations + std::vector<uint64_t> caa; // avalanche-avalanche correlations + std::vector<uint64_t> cll; // damage-damage distribution + std::vector<uint64_t> cDD; // damage-damage distribution + std::vector<uint64_t> csD; // damage-damage distribution uint64_t Nc; uint64_t Na; @@ -57,9 +57,9 @@ class ma : public hooks { long double lv; - std::list<std::list<unsigned>> avalanches; + std::list<std::list<graph::coordinate>> avalanches; - ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncum); + ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta); ~ma(); void pre_fracture(const network &) override; diff --git a/src/sample.cpp b/src/sample.cpp new file mode 100644 index 0000000..b8af937 --- /dev/null +++ b/src/sample.cpp @@ -0,0 +1,69 @@ + +#include "sample.hpp" +#include <iostream> +#include <cstdio> + +sample::sample(double Lx, double Ly, double beta) { + std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_sample.dat"; + sample_file.open(filename, std::ifstream::app); +} + +sample::~sample() { + sample_file.close(); +} + +void sample::pre_fracture(const network& n) { + sample_file << "<|"; + sample_file << "\"vr\"->{"; + for (const graph::vertex &v : n.G.vertices) { + sample_file << "{" << v.r.x << "," << v.r.y << "},"; + } + sample_file << "},\"vp\"->{"; + for (const graph::vertex &v : n.G.vertices) { + sample_file << "{"; + for (const graph::coordinate &r : v.polygon) { + sample_file << "{" << r.x << "," << r.y << "},"; + } + sample_file << "},"; + } + sample_file << "},\"ur\"->{"; + for (const graph::vertex &v : n.G.dual_vertices) { + sample_file << "{" << v.r.x << "," << v.r.y << "},"; + } + sample_file << "},\"up\"->{"; + for (const graph::vertex &v : n.G.dual_vertices) { + sample_file << "{"; + for (const graph::coordinate &r : v.polygon) { + sample_file << "{" << r.x << "," << r.y << "},"; + } + sample_file << "},"; + } + sample_file << "},\"e\"->{"; + for (const graph::edge &e : n.G.edges) { + sample_file << "{" << e.v[0] << "," << e.v[1] << "},"; + } + sample_file << "},\"ec\"->{"; + for (const graph::edge &e : n.G.edges) { + sample_file << "{" << e.crossings[0] << "," << e.crossings[1] << "},"; + } + sample_file << "},\"d\"->{"; + for (const graph::edge &e : n.G.dual_edges) { + sample_file << "{" << e.v[0] << "," << e.v[1] << "},"; + } + sample_file << "},\"dc\"->{"; + for (const graph::edge &e : n.G.dual_edges) { + sample_file << "{" << e.crossings[0] << "," << e.crossings[1] << "},"; + } + sample_file << "},\"data\"->{"; +} + +void sample::bond_broken(const network& net, const current_info& cur, unsigned i) { + long double c = logl(cur.conductivity / fabs(cur.currents[i])) + net.thresholds[i]; + sample_file << "{" << i << "," << c << "," << cur.conductivity << "},"; +} + +void sample::post_fracture(network &n) { + sample_file << "}"; + sample_file << "|>\n"; +} + diff --git a/src/sample.hpp b/src/sample.hpp new file mode 100644 index 0000000..d5d6fe5 --- /dev/null +++ b/src/sample.hpp @@ -0,0 +1,30 @@ + +#include <cstring> +#include <fstream> +#include <string> +#include <cinttypes> +#include <sstream> +#include <functional> +#include <iostream> +#include <array> + +#include <hooks.hpp> +#include <graph.hpp> +#include <network.hpp> + +class sample : public hooks { + // need: + // - interface for turning on and off specific measurements + // + private: + std::ofstream sample_file; + + public: + sample(double, double, double); + ~sample(); + + 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/sample_fracture.cpp b/src/sample_fracture.cpp new file mode 100644 index 0000000..42a6a5a --- /dev/null +++ b/src/sample_fracture.cpp @@ -0,0 +1,65 @@ + +#include <random> +#include <iostream> + +#include <cholmod.h> + +#include "randutils/randutils.hpp" + +#include <graph.hpp> +#include <network.hpp> +#include <hooks.hpp> +#include "sample.hpp" + +#include <csignal> +#include <cstring> +#include <atomic> + +int main(int argc, char* argv[]) { + + int opt; + + unsigned N = 1; + double Lx = 16; + double Ly = 16; + double beta = 0.5; + + while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -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; + default: + exit(1); + } + } + + cholmod_common c; + CHOL_F(start)(&c); + + sample meas(Lx, Ly, beta); + + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + for (unsigned trial = 0; trial < N; trial++) { + graph G(Lx, Ly, rng); + network network(G, &c); + network.set_thresholds(beta, rng); + network.fracture(meas); + } + + CHOL_F(finish)(&c); + + return 0; +} + diff --git a/src/sample_fracture_square.cpp b/src/sample_fracture_square.cpp new file mode 100644 index 0000000..fb2eb33 --- /dev/null +++ b/src/sample_fracture_square.cpp @@ -0,0 +1,67 @@ + +#include <random> +#include <iostream> + +#include <cholmod.h> + +#include "randutils/randutils.hpp" + +#include <graph.hpp> +#include <network.hpp> +#include <hooks.hpp> +#include "sample.hpp" + +#include <csignal> +#include <cstring> +#include <atomic> + +int main(int argc, char* argv[]) { + + int opt; + + unsigned N = 1; + unsigned Lx = 16; + unsigned Ly = 16; + double beta = 0.5; + + while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) { + switch (opt) { + case 'N': + N = (unsigned)atof(optarg); + break; + case 'X': + Lx = atoi(optarg); + break; + case 'Y': + Ly = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + default: + exit(1); + } + } + + cholmod_common c; + CHOL_F(start)(&c); + + sample meas(Lx, Ly, beta); + + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + graph G(Lx, Ly); + network perm_network(G, &c); + + for (unsigned trial = 0; trial < N; trial++) { + network tmp_network(perm_network); + tmp_network.set_thresholds(beta, rng); + tmp_network.fracture(meas); + } + + CHOL_F(finish)(&c); + + return 0; +} + |