From d6f6fab35801f4c8eb014736fabdc85b8b90a491 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 24 Jan 2019 14:32:47 -0500 Subject: removed bin counting measurement, added measurement of higher cumulants of the correlation functions --- src/fracture.cpp | 2 +- src/measurements.cpp | 87 ++++++++++++++++++++++++++++++++-------------------- src/measurements.hpp | 17 +++++----- 3 files changed, 63 insertions(+), 43 deletions(-) (limited to 'src') diff --git a/src/fracture.cpp b/src/fracture.cpp index 3505166..e260e14 100644 --- a/src/fracture.cpp +++ b/src/fracture.cpp @@ -57,7 +57,7 @@ int main(int argc, char* argv[]) { cholmod_common c; CHOL_F(start)(&c); - ma meas(Lx, Ly, ceil(Lx), ceil(Ly), beta); + ma meas(Lx, Ly, 2*ceil(Lx), 2*ceil(Ly), beta, 10); randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; diff --git a/src/measurements.cpp b/src/measurements.cpp index 72f074a..ffff06d 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -36,12 +36,17 @@ void update_field_file(std::string id, const std::vector& data, unsigned int std::ifstream file(filename); uint64_t N_old = 0; - std::vector data_old(data.size(), 0); + std::vector> data_old(data.size()); + for (unsigned j = 0; j < data.size(); j++) { + data_old[j].resize(data[j].size(), 0); + } if (file.is_open()) { file >> N_old; - for (unsigned int i = 0; i < data.size(); i++) { - file >> data_old[i]; + for (unsigned j = 0; j < data.size(); j++) { + for (unsigned int i = 0; i < data[j].size(); i++) { + file >> data_old[j][i]; + } } file.close(); } @@ -49,8 +54,11 @@ void update_field_file(std::string id, const std::vector& data, unsigned int std::ofstream file_out(filename); file_out <& data, const s } template -void autocorrelation(unsigned int Mx, unsigned int 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) { +void autocorrelation(unsigned int Mx, unsigned int 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); for (unsigned int i = 0; i < My * (Mx / 2 + 1); i++) { @@ -99,8 +107,10 @@ void autocorrelation(unsigned int Mx, unsigned int My, std::vector& out_data, fftw_execute(reverse_plan); - for (unsigned int i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) { - out_data[i] += (T)(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My)); + for (unsigned j = 0; j < out_data.size(); j++) { + for (unsigned int i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) { + out_data[j][i] += (T)pow(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My), j + 1); + } } } @@ -108,27 +118,35 @@ unsigned int edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned i return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly); } -ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta) : +ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta, unsigned Ncum) : Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(2 * (unsigned int)ceil(Lx * Ly / 2)), sc(2 * (unsigned int)ceil(Lx * Ly / 2), 0), sa(2 * (unsigned int)ceil(Lx * Ly / 2), 0), sC(2 * (unsigned int)ceil(Lx * Ly / 2), 0), sA(2 * (unsigned int)ceil(Lx * Ly / 2), 0), sd(3 * (unsigned int)ceil(Lx * Ly / 2), 0), - sb(log2(Mx < My ? Mx : My) + 1, 0), - Ccc((Mx / 2 + 1) * (My / 2 + 1), 0), - Css((Mx / 2 + 1) * (My / 2 + 1), 0), - Cmm((Mx / 2 + 1) * (My / 2 + 1), 0), - Caa((Mx / 2 + 1) * (My / 2 + 1), 0), - Cdd((Mx / 2 + 1) * (My / 2 + 1), 0), - Cll((Mx / 2 + 1) * (My / 2 + 1), 0), - Cee((Mx / 2 + 1) * (My / 2 + 1), 0) + Ccc(Ncum), + Css(Ncum), + Cmm(Ncum), + Caa(Ncum), + Cdd(Ncum), + Cll(Ncum), + Cee(Ncum) { N = 0; Nc = 0; Na = 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); + Cdd[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); + Cll[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); + Cee[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); + } // FFTW setup for correlation functions fftw_set_timelimit(1); @@ -140,6 +158,8 @@ ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta) : 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() { @@ -157,7 +177,6 @@ ma::~ma() { update_distribution_file("sA", sA, NA, Lx, Ly, beta); update_distribution_file("sC", sC, NC, Lx, Ly, beta); update_distribution_file("sd", sd, N, Lx, Ly, beta); - update_distribution_file("sb", sb, 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); @@ -166,6 +185,8 @@ ma::~ma() { 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("Cee", Cee, N, Lx, Ly, beta, Mx, My); + + //stress_file.close(); } void ma::pre_fracture(const network&) { @@ -243,22 +264,6 @@ void ma::post_fracture(network &n) { } } - unsigned int max_factor = log2(Mx < My ? Mx : My) + 1; - std::vector> bins(max_factor); - for (unsigned int i = 0; i < max_factor; i++) { - bins[i].resize(Mx * My / (unsigned int)pow(2, 2 * i)); - } - - for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) { - for (unsigned int i = 0; i < max_factor; i++) { - bins[i][edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx / pow(2, i), My / pow(2, i))] = 1; - } - } - - for (unsigned int i =0 ; i < max_factor; i++) { - sb[i] += bins[i].sum(); - } - // spanning cluster std::fill_n(fftw_forward_in, Mx * My, 0.0); @@ -311,6 +316,20 @@ void ma::post_fracture(network &n) { 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]; + + for (unsigned int i = 0; i < n.G.edges.size(); i++) { + double x; + if (!n.fuses[i]) { + stress_file << std::scientific << fabs(cur.currents[i]) << " "; + } + } + */ + sd[total_broken]++; N++; diff --git a/src/measurements.hpp b/src/measurements.hpp index 3d4ff3e..aa41fd1 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -32,6 +32,7 @@ class ma : public hooks { unsigned int My; double beta; Graph G; +// std::ofstream stress_file; // measurement storage std::vector sc; // cluster size distribution @@ -40,13 +41,13 @@ class ma : public hooks { std::vector sA; // avalanche size distribution std::vector sd; // avalanche size distribution std::vector sb; // bin size distribution - std::vector Ccc; // cluster-cluster correlations - std::vector Css; // surface-surface correlations - std::vector Cmm; // surface-surface correlations - std::vector Caa; // avalanche-avalanche correlations - std::vector Cdd; // damage-damage distribution - std::vector Cll; // damage-damage distribution - std::vector Cee; // damage-damage distribution + std::vector> Ccc; // cluster-cluster correlations + std::vector> Css; // surface-surface correlations + std::vector> Cmm; // surface-surface correlations + std::vector> Caa; // avalanche-avalanche correlations + std::vector> Cdd; // damage-damage distribution + std::vector> Cll; // damage-damage distribution + std::vector> Cee; // damage-damage distribution uint64_t Nc; uint64_t NC; uint64_t Na; @@ -58,7 +59,7 @@ class ma : public hooks { std::list> avalanches; - ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta); + ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta, unsigned Ncum); ~ma(); void pre_fracture(const network &) override; -- cgit v1.2.3-70-g09d2