From 578652bf9d57a94a362b4c8a6d36c1755a214863 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sat, 9 Feb 2019 22:13:14 -0500
Subject: now collect correlatino functions and size distributions for
truncated and untruncated measurements
---
src/measurements.cpp | 200 ++++++++++++++++++++++++++++++++++++++++++++-------
src/measurements.hpp | 15 ++++
2 files changed, 191 insertions(+), 24 deletions(-)
(limited to 'src')
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;
--
cgit v1.2.3-70-g09d2