summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2018-07-31 15:44:26 -0400
committerJaron Kent-Dobias <jaron@kent-dobias.com>2018-07-31 15:44:26 -0400
commit5ff02c5e6608afe0d5c5e58acb57311a4b03544e (patch)
treec5c8369e4cabc55ddbfae0bf04de40b7204637c7 /src
parent22d55ab7279affffce87cdbdf173aa1928a4ff75 (diff)
downloadfuse_networks-5ff02c5e6608afe0d5c5e58acb57311a4b03544e.tar.gz
fuse_networks-5ff02c5e6608afe0d5c5e58acb57311a4b03544e.tar.bz2
fuse_networks-5ff02c5e6608afe0d5c5e58acb57311a4b03544e.zip
added support for measuring correlation functions
Diffstat (limited to 'src')
-rw-r--r--src/fracture.c316
1 files changed, 313 insertions, 3 deletions
diff --git a/src/fracture.c b/src/fracture.c
index 41bd056..49fb039 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -1,10 +1,13 @@
+#include <fftw3.h>
#include "fracture.h"
int main(int argc, char *argv[]) {
int opt;
+ fftw_set_timelimit(1);
+
// defining variables to be (potentially) set by command line flags
uint8_t filename_len;
uint32_t N;
@@ -12,7 +15,7 @@ int main(int argc, char *argv[]) {
double beta, inf, cutoff, crack_len;
bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual,
save_network, save_crit_stress, save_energy, save_conductivity,
- save_damage, save_threshold, save_current_load;
+ save_damage, save_threshold, save_current_load, save_correlations;
bound_t boundary;
lattice_t lattice;
@@ -41,6 +44,7 @@ int main(int argc, char *argv[]) {
save_energy = false;
save_threshold = false;
save_current_load = false;
+ save_correlations = false;
uint8_t bound_i;
char boundc2 = 'f';
@@ -50,7 +54,7 @@ int main(int argc, char *argv[]) {
// get commandline options
- while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TE")) != -1) {
+ while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TEz")) != -1) {
switch (opt) {
case 'n':
N = atoi(optarg);
@@ -155,6 +159,9 @@ int main(int argc, char *argv[]) {
case 'C':
save_current_load = true;
break;
+ case 'z':
+ save_correlations = true;
+ break;
default: /* '?' */
exit(EXIT_FAILURE);
}
@@ -166,6 +173,72 @@ int main(int argc, char *argv[]) {
} else {
boundc = 'c';
}
+
+ double *dd_correlations; // damage-damage
+ double *dp_correlations; // damage-precursor
+ double *ds_correlations; // damage-stress
+ double *pp_correlations; // precursor-precursor
+ double *ps_correlations; // precursor-stress
+ double *ss_correlations; // stress-stress
+ double *DD_correlations; // after-crack damage-damage
+ double *DP_correlations; // after-crack damage-precursor
+ double *PP_correlations; // after-crack precursor-precursor
+ double *fftw_forward_in;
+ fftw_complex *fftw_forward_out;
+ fftw_complex *fftw_reverse_in;
+ double *fftw_reverse_out;
+ fftw_plan forward_plan;
+ fftw_plan reverse_plan;
+ uint64_t N_averaged = 0;
+ double mean_D = 0;
+ double mean_P = 0;
+ double mean_d = 0;
+ double mean_p = 0;
+ double mean_s = 0;
+ char *correlations_filename;
+ if (save_correlations) {
+ assert(lattice == DIAGONAL_LATTICE);
+ dd_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+ dp_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+ ds_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+ pp_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+ ps_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+ ss_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+ DD_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+ DP_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+ PP_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+ fftw_forward_in = (double *)fftw_malloc(pow(L, 2) * sizeof(double));
+ fftw_forward_out = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex));
+ fftw_reverse_in = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex));
+ fftw_reverse_out = (double *)fftw_malloc(pow(L, 2) * sizeof(double));
+ forward_plan = fftw_plan_dft_r2c_2d(L, L, fftw_forward_in, fftw_forward_out, 0);
+ reverse_plan = fftw_plan_dft_c2r_2d(L, L, fftw_reverse_in, fftw_reverse_out, 0);
+
+ correlations_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(correlations_filename, filename_len, "corr_%c_%c_%c_%c_%d_%g_%g.dat",
+ lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+
+ FILE *correlations_out = fopen(correlations_filename, "rb");
+
+ if (correlations_out != NULL) {
+ fread(&N_averaged, sizeof(uint64_t), 1, correlations_out);
+ fread(&mean_d, sizeof(double), 1, correlations_out);
+ fread(&mean_p, sizeof(double), 1, correlations_out);
+ fread(&mean_s, sizeof(double), 1, correlations_out);
+ fread(&mean_D, sizeof(double), 1, correlations_out);
+ fread(&mean_P, sizeof(double), 1, correlations_out);
+ fread(dd_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fread(dp_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fread(ds_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fread(pp_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fread(ps_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fread(ss_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fread(DD_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fread(DP_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fread(PP_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fclose(correlations_out);
+ }
+ }
FILE *data_out;
@@ -284,7 +357,6 @@ int main(int argc, char *argv[]) {
net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c);
net_t *tmp_net = net_copy(net, &c);
data_t *data = net_fracture(tmp_net, &c, cutoff);
- net_free(tmp_net, &c);
uint_t max_pos = 0;
long double max_val = 0;
@@ -327,6 +399,205 @@ int main(int argc, char *argv[]) {
}
}
+ if (save_correlations) {
+ uint32_t damage1 = 0;
+ for (uint32_t j = 0; j < g->ne; j++) {
+ if (tmp_net->fuses[j]) {
+ fftw_forward_in[j] = 1.0;
+ damage1++;
+ } else {
+ fftw_forward_in[j] = 0.0;
+ }
+ }
+
+ fftw_execute(forward_plan);
+ fftw_complex *D_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
+ memcpy(D_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
+
+ dll_t *cycle = find_cycles(g, tmp_net->fuses);
+ components_t *comp = get_clusters(tmp_net);
+
+ uint32_t in_crack1 = 0;
+ for (uint32_t j = 0; j < g->ne; j++) {
+ if (tmp_net->fuses[j] && comp->labels[g->dev[2 * cycle->e]] == comp->labels[g->dev[2 * j]]) {
+ fftw_forward_in[j] = 1.0;
+ in_crack1++;
+ } else {
+ fftw_forward_in[j] = 0.0;
+ }
+ }
+
+ fftw_execute(forward_plan);
+ fftw_complex *P_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
+ memcpy(P_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
+
+ uint32_t damage2 = 0;
+ for (uint32_t j = 0; j < g->ne; j++) {
+ if (net->fuses[j]) {
+ fftw_forward_in[j] = 1.0;
+ damage2++;
+ } else {
+ fftw_forward_in[j] = 0.0;
+ }
+ }
+
+ fftw_execute(forward_plan);
+ fftw_complex *d_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
+ memcpy(d_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
+
+ uint32_t in_crack2 = 0;
+ for (uint32_t j = 0; j < g->ne; j++) {
+ if (net->fuses[j] && comp->labels[2 * g->dev[cycle->e]] == comp->labels[g->dev[2 * j]]) {
+ fftw_forward_in[j] = 1.0;
+ in_crack2++;
+ } else {
+ fftw_forward_in[j] = 0.0;
+ }
+ }
+
+ fftw_execute(forward_plan);
+ fftw_complex *p_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
+ memcpy(p_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
+
+ graph_components_free(comp);
+ int test = 0;
+ while (cycle != NULL) {
+ dll_t *old = cycle;
+ cycle = cycle->right;
+ free(old);
+ }
+
+ double *tmp_voltage = net_voltages(net, &c);
+ double *tmp_current = net_currents(net, tmp_voltage, &c);
+ free(tmp_voltage);
+
+ double t_stress = 0;
+ for (uint32_t j = 0; j < g->ne; j++) {
+ fftw_forward_in[j] = tmp_current[j];
+ t_stress += tmp_current[j];
+ }
+
+ free(tmp_current);
+
+ fftw_execute(forward_plan);
+ fftw_complex *s_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
+ memcpy(s_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
+
+ mean_D = (double)damage1 / (1.0 + N_averaged) + (double)N_averaged * mean_D / (N_averaged + 1.0);
+ mean_P = (double)in_crack1 / (1.0 + N_averaged) + (double)N_averaged * mean_P / (N_averaged + 1.0);
+ mean_d = (double)damage2 / (1.0 + N_averaged) + (double)N_averaged * mean_d / (N_averaged + 1.0);
+ mean_p = (double)in_crack2 / (1.0 + N_averaged) + (double)N_averaged * mean_p / (N_averaged + 1.0);
+ mean_s = (double)t_stress / (1.0 + N_averaged) + (double)N_averaged * mean_s / (N_averaged + 1.0);
+
+ for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+ fftw_reverse_in[j][0] = D_transform[j][0] * D_transform[j][0] + D_transform[j][1] * D_transform[j][1];
+ fftw_reverse_in[j][1] = 0.0;
+ }
+
+ fftw_execute(reverse_plan);
+
+ for (uint32_t j = 0; j < g->ne; j++) {
+ DD_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DD_correlations[j] / (N_averaged + 1.0);
+ }
+
+ for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+ fftw_reverse_in[j][0] = D_transform[j][0] * P_transform[j][0] + D_transform[j][1] * P_transform[j][1];
+ fftw_reverse_in[j][1] = D_transform[j][0] * P_transform[j][1] - D_transform[j][1] * P_transform[j][0];
+ }
+
+ fftw_execute(reverse_plan);
+
+ for (uint32_t j = 0; j < g->ne; j++) {
+ DP_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DP_correlations[j] / (N_averaged + 1.0);
+ }
+
+ for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+ fftw_reverse_in[j][0] = P_transform[j][0] * P_transform[j][0] + P_transform[j][1] * P_transform[j][1];
+ fftw_reverse_in[j][1] = 0.0;
+ }
+
+ fftw_execute(reverse_plan);
+
+ for (uint32_t j = 0; j < g->ne; j++) {
+ PP_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * PP_correlations[j] / (N_averaged + 1.0);
+ }
+
+ for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+ fftw_reverse_in[j][0] = d_transform[j][0] * d_transform[j][0] + d_transform[j][1] * d_transform[j][1];
+ fftw_reverse_in[j][1] = 0.0;
+ }
+
+ fftw_execute(reverse_plan);
+
+ for (uint32_t j = 0; j < g->ne; j++) {
+ dd_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dd_correlations[j] / (N_averaged + 1.0);
+ }
+
+ for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+ fftw_reverse_in[j][0] = d_transform[j][0] * p_transform[j][0] + d_transform[j][1] * p_transform[j][1];
+ fftw_reverse_in[j][1] = d_transform[j][0] * p_transform[j][1] - d_transform[j][1] * p_transform[j][0];
+ }
+
+ fftw_execute(reverse_plan);
+
+ for (uint32_t j = 0; j < g->ne; j++) {
+ dp_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dp_correlations[j] / (N_averaged + 1.0);
+ }
+
+ for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+ fftw_reverse_in[j][0] = d_transform[j][0] * s_transform[j][0] + d_transform[j][1] * s_transform[j][1];
+ fftw_reverse_in[j][1] = d_transform[j][0] * s_transform[j][1] - d_transform[j][1] * s_transform[j][0];
+ }
+
+ fftw_execute(reverse_plan);
+
+ for (uint32_t j = 0; j < g->ne; j++) {
+ ds_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * ds_correlations[j] / (N_averaged + 1.0);
+ }
+
+ for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+ fftw_reverse_in[j][0] = p_transform[j][0] * p_transform[j][0] + p_transform[j][1] * p_transform[j][1];
+ fftw_reverse_in[j][1] = 0.0;
+ }
+
+ fftw_execute(reverse_plan);
+
+ for (uint32_t j = 0; j < g->ne; j++) {
+ pp_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * pp_correlations[j] / (N_averaged + 1.0);
+ }
+
+ for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+ fftw_reverse_in[j][0] = p_transform[j][0] * s_transform[j][0] + p_transform[j][1] * s_transform[j][1];
+ fftw_reverse_in[j][1] = p_transform[j][0] * s_transform[j][1] - p_transform[j][1] * s_transform[j][0];
+ }
+
+ fftw_execute(reverse_plan);
+
+ for (uint32_t j = 0; j < g->ne; j++) {
+ ps_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * ps_correlations[j] / (N_averaged + 1.0);
+ }
+
+ for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+ fftw_reverse_in[j][0] = s_transform[j][0] * s_transform[j][0] + s_transform[j][1] * s_transform[j][1];
+ fftw_reverse_in[j][1] = 0.0;
+ }
+
+ fftw_execute(reverse_plan);
+
+ for (uint32_t j = 0; j < g->ne; j++) {
+ ss_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * ss_correlations[j] / (N_averaged + 1.0);
+ }
+
+ free(D_transform);
+ free(P_transform);
+ free(d_transform);
+ free(p_transform);
+ free(s_transform);
+ N_averaged++;
+ }
+
+ net_free(tmp_net, &c);
+
if (save_crit_stress) {
crit_stress[i] = data->extern_field[max_pos];
}
@@ -514,6 +785,45 @@ int main(int argc, char *argv[]) {
free(crit_stress);
}
+ if (save_correlations) {
+ FILE *correlations_out = fopen(correlations_filename, "wb");
+ fwrite(&N_averaged, sizeof(uint64_t), 1, correlations_out);
+ fwrite(&mean_d, sizeof(double), 1, correlations_out);
+ fwrite(&mean_p, sizeof(double), 1, correlations_out);
+ fwrite(&mean_s, sizeof(double), 1, correlations_out);
+ fwrite(&mean_D, sizeof(double), 1, correlations_out);
+ fwrite(&mean_P, sizeof(double), 1, correlations_out);
+ fwrite(dd_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fwrite(dp_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fwrite(ds_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fwrite(pp_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fwrite(ps_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fwrite(ss_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fwrite(DD_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fwrite(DP_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fwrite(PP_correlations, sizeof(double), pow(L, 2), correlations_out);
+ fclose(correlations_out);
+
+ free(dd_correlations);
+ free(dp_correlations);
+ free(ds_correlations);
+ free(pp_correlations);
+ free(ps_correlations);
+ free(ss_correlations);
+ free(DD_correlations);
+ free(DP_correlations);
+ free(PP_correlations);
+ fftw_free(fftw_forward_in);
+ fftw_free(fftw_forward_out);
+ fftw_free(fftw_reverse_in);
+ fftw_free(fftw_reverse_out);
+ fftw_destroy_plan(forward_plan);
+ fftw_destroy_plan(reverse_plan);
+ free(correlations_filename);
+ }
+
+ fftw_cleanup();
+
CHOL_F(finish)(&c);
return 0;