From d944bcc3df0a8d7a10b755b5858c85e61a835a35 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 21 Feb 2019 14:01:24 -0500 Subject: simplified and sped up measurementsn substantially, correlation functions and some distributions again incompatible --- src/measurements.cpp | 173 ++++++++++++++++++--------------------------------- 1 file changed, 61 insertions(+), 112 deletions(-) (limited to 'src/measurements.cpp') 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& data, unsigned N, d std::ifstream file(filename); uint64_t N_old = 0; - std::vector> data_old(data.size()); - for (unsigned j = 0; j < data.size(); j++) { - data_old[j].resize(data[j].size(), 0); - } + std::vector 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& data, unsigned N, d file_out <>& data, co } } +void autocorrelation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector& data, const std::list& pos) { + for (std::list::const_iterator it1 = pos.begin(); it1 != pos.end(); it1++) { + for (std::list::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& data, const std::list& pos1, const std::list& pos2) { + for (std::list::const_iterator it1 = pos1.begin(); it1 != pos1.end(); it1++) { + for (std::list::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 void autocorrelation(unsigned Mx, unsigned My, std::vector>& 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 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 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> components(num); + std::vector> 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 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++; } -- cgit v1.2.3-70-g09d2