diff options
| -rw-r--r-- | CMakeLists.txt | 2 | ||||
| -rw-r--r-- | src/fracture.c | 316 |
2 files changed, 314 insertions, 4 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index ade950f..e924e98 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,7 +14,7 @@ file(GLOB EXE_SOURCES src/*.c) foreach( src_file ${EXE_SOURCES} ) string( REGEX REPLACE ".*/src/(.*)\.c" "\\1" exe_name ${src_file} ) add_executable( ${exe_name} ${src_file} ) - target_link_libraries(${exe_name} fracture_stuff gsl c cblas lapack dl pthread cholmod amd colamd suitesparseconfig camd ccolamd rt metis m jst tcmalloc profiler) + target_link_libraries(${exe_name} fracture_stuff gsl c cblas lapack dl pthread cholmod amd colamd suitesparseconfig camd ccolamd rt metis fftw3 m jst) install(TARGETS ${exe_name} DESTINATION bin) endforeach( src_file ${EXE_SOURCES} ) 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; |