added support for measuring correlation functions

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;