diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/measurements.cpp | 200 | ||||
| -rw-r--r-- | src/measurements.hpp | 15 |
2 files changed, 191 insertions, 24 deletions
diff --git a/src/measurements.cpp b/src/measurements.cpp index 86c7dfd..324313c 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -125,6 +125,13 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncu sl(2 * (unsigned)ceil(Lx * Ly / 2), 0), sd(3 * (unsigned)ceil(Lx * Ly / 2), 0), sD(3 * (unsigned)ceil(Lx * Ly / 2), 0), + Sc(2 * (unsigned)ceil(Lx * Ly / 2), 0), + Ss(2 * (unsigned)ceil(Lx * Ly / 2), 0), + Sm(2 * (unsigned)ceil(Lx * Ly / 2), 0), + Sa(2 * (unsigned)ceil(Lx * Ly / 2), 0), + Sl(2 * (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), @@ -132,7 +139,15 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncu Cll(Ncum), Cdd(Ncum), CDD(Ncum), - CsD(Ncum) + CsD(Ncum), + ccc(Ncum), + css(Ncum), + cmm(Ncum), + caa(Ncum), + cll(Ncum), + cdd(Ncum), + cDD(Ncum), + csD(Ncum) { N = 0; Nc = 0; @@ -146,6 +161,14 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncu 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); + 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 @@ -180,6 +203,14 @@ ma::~ma() { update_distribution_file("sd", sd, N, Lx, Ly, beta); update_distribution_file("sD", sD, N, Lx, Ly, beta); + 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); update_field_file("Css", Css, N, Lx, Ly, beta, Mx, My); update_field_file("Cmm", Cmm, N, Lx, Ly, beta, Mx, My); @@ -189,6 +220,15 @@ ma::~ma() { update_field_file("CDD", CDD, N, Lx, Ly, beta, Mx, My); update_field_file("CsD", CsD, N, Lx, Ly, beta, Mx, My); + 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("cDD", cDD, N, Lx, Ly, beta, Mx, My); + update_field_file("csD", csD, N, Lx, Ly, beta, Mx, My); + //stress_file.close(); } @@ -206,12 +246,28 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) { std::fill_n(fftw_forward_in, Mx * My, 0.0); + unsigned avalanches_trunc = 0; + for (auto e : avalanches.back()) { - fftw_forward_in[edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My)] = 1.0; + unsigned ind = edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My); + if (fftw_forward_in[ind] == 0.0) { + fftw_forward_in[ind] = 1.0; + avalanches_trunc++; + } } + Sa[avalanches_trunc]++; autocorrelation(Mx, My, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + 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); + lv = c; avalanches.push_back({i}); } else { @@ -230,16 +286,39 @@ void ma::post_fracture(network &n) { // crack surface correlations std::fill_n(fftw_forward_in, Mx * My, 0.0); + unsigned surface_size_trunc = 0; + 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)] = 1.0; - fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My)] = 1.0; + unsigned ind1 = edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My); + unsigned ind2 = edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My); + if (fftw_forward_in[ind1] == 0.0) { + fftw_forward_in[ind1] = 1.0; + surface_size_trunc++; + } + if (fftw_forward_in[ind2] == 0.0) { + fftw_forward_in[ind2] = 1.0; + surface_size_trunc++; + } } ss[crack.size()]++; + 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); + + std::fill_n(fftw_forward_in, Mx * My, 0.0); + + 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; + } + + Ss[surface_size_trunc]++; + 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); + correlation(Mx, My, css, tss, tss, reverse_plan, fftw_reverse_in, fftw_reverse_out); unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]]; @@ -251,50 +330,105 @@ void ma::post_fracture(network &n) { // 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) { + unsigned cluster_size_trunc = 0; 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)] = 1.0; + unsigned ind = edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My); + if (fftw_forward_in[ind] == 0.0) { + cluster_size_trunc++; + fftw_forward_in[ind] = 1.0; + } } + Sc[cluster_size_trunc - 1]++; sc[components[i].size() - 1]++; + autocorrelation(Mx, My, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); Nc++; 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; } + + 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; + } } } // spanning cluster - std::fill_n(fftw_forward_in, Mx * My, 0.0); + // std::fill_n(fftw_forward_in, Mx * My, 0.0); we already reset in the last loop sm[components[crack_component].size() - 1]++; + unsigned spanning_cluster_trunc = 0; 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; + unsigned ind = edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My); + if (fftw_forward_in[ind] == 0.0) { + fftw_forward_in[ind] = 1.0; + spanning_cluster_trunc++; + } } autocorrelation(Mx, My, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + std::fill_n(fftw_forward_in, Mx * My, 0.0); + + Sm[spanning_cluster_trunc - 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); + /// damage at end std::fill_n(fftw_forward_in, Mx * My, 0.0); unsigned final_broken = 0; + unsigned final_trunc = 0; 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; + unsigned ind = edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My); + if (fftw_forward_in[ind] == 0.0) { + final_trunc++; + fftw_forward_in[ind] = 1.0; + } } } - fftw_complex *tDD = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out); + fftw_complex *TDD = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out); + SD[final_trunc]++; sD[final_broken]++; - 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); + 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); + + free(Tss); + free(TDD); + + 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; + } + } + + 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); free(tss); free(tDD); @@ -303,41 +437,59 @@ void ma::post_fracture(network &n) { // rewind the last avalanche sl[avalanches.back().size() - 1]++; + unsigned last_trunc = 0; 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; + unsigned ind = edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My); + if (fftw_forward_in[ind] == 0.0) { + fftw_forward_in[ind] = 1.0; + last_trunc++; + } } + Sl[last_trunc - 1]++; autocorrelation(Mx, My, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + std::fill_n(fftw_forward_in, Mx * My, 0.0); + + // rewind the last avalanche + for (auto e : avalanches.back()) { + 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); + unsigned damage_trunc = 0; + for (unsigned i = 0; i < n.G.edges.size(); i++) { if (n.fuses[i]) { total_broken++; - fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] = 1.0; + unsigned ind = edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My); + if (fftw_forward_in[ind] == 0.0) { + fftw_forward_in[ind] = 1.0; + damage_trunc++; + } } } + Sd[damage_trunc]++; autocorrelation(Mx, My, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); - /* - current_info cur = n.get_current_info(); - - unsigned ii = avalanches.back().front(); - long double c = logl(cur.conductivity / fabs(cur.currents[ii])) + n.thresholds[ii]; + std::fill_n(fftw_forward_in, Mx * My, 0.0); for (unsigned i = 0; i < n.G.edges.size(); i++) { - double x; - if (!n.fuses[i]) { - stress_file << std::scientific << fabs(cur.currents[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]++; diff --git a/src/measurements.hpp b/src/measurements.hpp index 1d263c3..c735163 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -42,6 +42,13 @@ class ma : public hooks { 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<uint64_t> Sc; // non-spanning cluster size distribution + std::vector<uint64_t> Ss; // minimal spanning cluster size distribution + 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 @@ -50,6 +57,14 @@ class ma : public hooks { 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<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 uint64_t Nc; uint64_t Na; |