diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/fracture.cpp | 74 | ||||
| -rw-r--r-- | src/measurements.cpp | 100 | ||||
| -rw-r--r-- | src/measurements.hpp | 6 |
3 files changed, 130 insertions, 50 deletions
diff --git a/src/fracture.cpp b/src/fracture.cpp index 53a11a5..b1bde45 100644 --- a/src/fracture.cpp +++ b/src/fracture.cpp @@ -35,7 +35,11 @@ int main(int argc, char* argv[]) { double Ly = 16; double beta = 0.5; - while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) { + 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); @@ -49,6 +53,14 @@ int main(int argc, char* argv[]) { 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); } @@ -57,30 +69,54 @@ int main(int argc, char* argv[]) { cholmod_common c; CHOL_F(start)(&c); - // fourier transforms of powers of two are faster - unsigned Mx = pow(2, ceil(log2(4*Lx))); - unsigned My = pow(2, ceil(log2(4*Ly))); - - ma meas(Lx, Ly, Mx, My, beta); - randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; - for (unsigned trial = 0; trial < N; trial++) { - while (true) { - try { - graph G(Lx, Ly, rng); - network network(G, &c); - network.set_thresholds(beta, rng); - network.fracture(meas); - break; - } catch (std::exception &e) { - std::cout << e.what() << '\n'; + if (use_aN) { + unsigned Mx = pow(2, ceil(log2(4*2*n*a))); + unsigned My = pow(2, ceil(log2(4*2*n/a))); + + ma meas(n, a, Mx, My, beta); + + for (unsigned trial = 0; trial < N; trial++) { + while (true) { + try { + graph G(n, a, rng); + network network(G, &c); + network.set_thresholds(beta, rng); + network.fracture(meas); + break; + } catch (std::exception &e) { + std::cout << e.what() << '\n'; + } } + + if (quit.load()) + break; } + } else { + // fourier transforms of powers of two are faster + unsigned Mx = pow(2, ceil(log2(4*Lx))); + unsigned My = pow(2, ceil(log2(4*Ly))); + + ma meas(Lx, Ly, Mx, My, beta); + + for (unsigned trial = 0; trial < N; trial++) { + while (true) { + try { + graph G(Lx, Ly, rng); + network network(G, &c); + network.set_thresholds(beta, rng); + network.fracture(meas); + break; + } catch (std::exception &e) { + std::cout << e.what() << '\n'; + } + } - if (quit.load()) - break; + if (quit.load()) + break; + } } CHOL_F(finish)(&c); diff --git a/src/measurements.cpp b/src/measurements.cpp index ed4b8eb..0b570cd 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -3,8 +3,8 @@ #include <iostream> #include <cstdio> -void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned N, double Lx, double Ly, double beta) { - std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat"; +void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned N, std::string model_string) { + std::string filename = model_string + id + ".dat"; std::ifstream file(filename); uint64_t N_old = 0; @@ -32,8 +32,8 @@ void update_distribution_file(std::string id, const std::vector<uint64_t>& data, } template <class T> -void update_field_file(std::string id, const std::vector<T>& data, unsigned N, double Lx, double Ly, double beta, unsigned Mx, unsigned My) { - std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat"; +void update_field_file(std::string id, const std::vector<T>& data, unsigned N, std::string model_string, unsigned Mx, unsigned My) { + std::string filename = model_string + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat"; std::ifstream file(filename); uint64_t N_old = 0; @@ -95,7 +95,7 @@ void autocorrelation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vecto 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)]++; + data[(unsigned)(Mx * (Δx / Lx)) + (Mx / 2) * (unsigned)(My * (Δy / Ly))]++; } } } @@ -136,8 +136,48 @@ 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(unsigned n, double a, unsigned Mx, unsigned My, double beta) : + Mx(Mx), My(My), G(2 * n), + sc(3 * n, 0), + ss(3 * n, 0), + sm(3 * n, 0), + sa(3 * n, 0), + sl(3 * n, 0), + sD(3 * n, 0), + 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; + + if (beta != 0.0) { + model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_"; + } else { + model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_INF_"; + } + + // FFTW setup for correlation functions + /* + fftw_set_timelimit(1); + + fftw_forward_in = (double *)fftw_malloc(Mx * My * sizeof(double)); + fftw_forward_out = (fftw_complex *)fftw_malloc(Mx * My * sizeof(fftw_complex)); + fftw_reverse_in = (fftw_complex *)fftw_malloc(Mx * My * sizeof(fftw_complex)); + fftw_reverse_out = (double *)fftw_malloc(Mx * My * sizeof(double)); + + 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); + */ +} + 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)), + Mx(Mx), My(My), 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), @@ -156,6 +196,12 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta) : Nc = 0; Na = 0; + if (beta != 0.0) { + model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_"; + } else { + model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_INF_"; + } + // FFTW setup for correlation functions /* fftw_set_timelimit(1); @@ -182,20 +228,20 @@ ma::~ma() { 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_field_file("ccc", ccc, Nc, Lx, Ly, beta, Mx, My); - update_field_file("css", css, N, Lx, Ly, beta, Mx, My); - 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("csD", csD, N, Lx, Ly, beta, Mx, My); + update_distribution_file("sc", sc, Nc, model_string); + update_distribution_file("ss", ss, N, model_string); + update_distribution_file("sm", sm, N, model_string); + update_distribution_file("sa", sa, Na, model_string); + update_distribution_file("sl", sl, N, model_string); + update_distribution_file("sD", sD, N, model_string); + + update_field_file("ccc", ccc, Nc, model_string, Mx, My); + update_field_file("css", css, N, model_string, Mx, My); + update_field_file("cmm", cmm, N, model_string, Mx, My); + update_field_file("caa", caa, Na, model_string, Mx, My); + update_field_file("cll", cll, N, model_string, Mx, My); + update_field_file("cDD", cDD, N, model_string, Mx, My); + update_field_file("csD", csD, N, model_string, Mx, My); //stress_file.close(); } @@ -212,7 +258,7 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) { sa[avalanches.back().size() - 1]++; Na++; - autocorrelation2(Lx, Ly, Mx, My, caa, avalanches.back()); + autocorrelation2(net.G.L.x, net.G.L.y, Mx, My, caa, avalanches.back()); lv = c; avalanches.push_back({net.G.edges[i].r}); @@ -236,7 +282,7 @@ void ma::post_fracture(network &n) { sr.push_back(n.G.dual_edges[edge].r); } - autocorrelation2(Lx, Ly, Mx, My, css, sr); + autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, css, sr); unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]]; @@ -251,7 +297,7 @@ void ma::post_fracture(network &n) { sc[components[i].size() - 1]++; Nc++; - autocorrelation2(Lx, Ly, Mx, My, ccc, components[i]); + autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ccc, components[i]); } } @@ -259,7 +305,7 @@ void ma::post_fracture(network &n) { sm[components[crack_component].size() - 1]++; - autocorrelation2(Lx, Ly, Mx, My, cmm, components[crack_component]); + autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cmm, components[crack_component]); /// damage at end std::list<graph::coordinate> Dr; @@ -272,15 +318,15 @@ void ma::post_fracture(network &n) { sD[Dr.size() - 1]++; - autocorrelation2(Lx, Ly, Mx, My, cDD, Dr); - correlation2(Lx, Ly, Mx, My, csD, sr, Dr); + autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cDD, Dr); + correlation2(n.G.L.x, n.G.L.y, Mx, My, csD, sr, Dr); // ********************** LAST AVALANCHE ********************* // rewind the last avalanche sl[avalanches.back().size() - 1]++; - autocorrelation2(Lx, Ly, Mx, My, cll, avalanches.back()); + autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cll, avalanches.back()); N++; } diff --git a/src/measurements.hpp b/src/measurements.hpp index 0fd4644..600c6ac 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -28,11 +28,8 @@ class ma : public hooks { fftw_plan reverse_plan; */ unsigned N; - double Lx; - double Ly; unsigned Mx; unsigned My; - double beta; Graph G; // std::ofstream stress_file; @@ -56,10 +53,11 @@ class ma : public hooks { public: long double lv; - std::list<std::list<graph::coordinate>> avalanches; + std::string model_string; ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta); + ma(unsigned n, double a, unsigned Mx, unsigned My, double beta); ~ma(); void pre_fracture(const network &) override; |