diff options
| -rw-r--r-- | makefile | 2 | ||||
| -rw-r--r-- | makefile_ept | 2 | ||||
| -rw-r--r-- | makefile_hal | 2 | ||||
| -rw-r--r-- | src/anal_process.c | 9 | ||||
| -rw-r--r-- | src/big_anal_process.c | 11 | ||||
| -rw-r--r-- | src/break_edge.c | 16 | ||||
| -rw-r--r-- | src/cen_anal_process.c | 154 | ||||
| -rw-r--r-- | src/corr_test.c | 2 | ||||
| -rw-r--r-- | src/fracture.c | 15 | ||||
| -rw-r--r-- | src/fracture.h | 8 | ||||
| -rw-r--r-- | src/net.c | 10 | ||||
| -rw-r--r-- | src/net_fracture.c | 15 | ||||
| -rw-r--r-- | src/rand.c | 9 |
13 files changed, 222 insertions, 33 deletions
@@ -4,7 +4,7 @@ CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-fi LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler -ltcmalloc OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand -BIN = corr_test fracture fitting_functions anal_process big_anal_process +BIN = corr_test fracture anal_process big_anal_process cen_anal_process all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} diff --git a/makefile_ept b/makefile_ept index 81ebb62..7325a88 100644 --- a/makefile_ept +++ b/makefile_ept @@ -4,7 +4,7 @@ CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-fi LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler -ltcmalloc OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand -BIN = corr_test fracture anal_process big_anal_process +BIN = corr_test fracture anal_process big_anal_process cen_anal_process all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} diff --git a/makefile_hal b/makefile_hal index 7ebc211..6fb1803 100644 --- a/makefile_hal +++ b/makefile_hal @@ -4,7 +4,7 @@ CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-fi LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler #-ltcmalloc OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand -BIN = corr_test fracture fitting_functions anal_process big_anal_process +BIN = corr_test fracture anal_process big_anal_process cen_anal_process all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} diff --git a/src/anal_process.c b/src/anal_process.c index c66aea0..0ebec84 100644 --- a/src/anal_process.c +++ b/src/anal_process.c @@ -37,6 +37,8 @@ int main(int argc, char *argv[]) { double *errors_c2 = (double *)malloc(nc * sizeof(double)); double *vals_c3 = (double *)malloc(nc * sizeof(double)); double *errors_c3 = (double *)malloc(nc * sizeof(double)); + double *vals_c4 = (double *)malloc(nc * sizeof(double)); + double *errors_c4 = (double *)malloc(nc * sizeof(double)); char *out_filename = (char *)malloc(100 * sizeof(char)); uint_t out_filename_len = 0; @@ -94,6 +96,11 @@ int main(int argc, char *argv[]) { mom(dist_len, 3, dist, mom3); vals_c3[i] = mom3[0]; errors_c3[i] = mom3[1]; + + double mom4[2]; + mom(dist_len, 4, dist, mom4); + vals_c4[i] = mom4[0]; + errors_c4[i] = mom4[1]; } free(dist); } @@ -107,7 +114,7 @@ int main(int argc, char *argv[]) { FILE *out_file = fopen(out_filename, "w"); for (uint_t i = 0; i < nc; i++) { - fprintf(out_file, "%u %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i]); + fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]); } fclose(out_file); diff --git a/src/big_anal_process.c b/src/big_anal_process.c index 8345785..0f7106f 100644 --- a/src/big_anal_process.c +++ b/src/big_anal_process.c @@ -53,6 +53,8 @@ int main(int argc, char *argv[]) { double *errors_c2 = (double *)malloc(nc * sizeof(double)); double *vals_c3 = (double *)malloc(nc * sizeof(double)); double *errors_c3 = (double *)malloc(nc * sizeof(double)); + double *vals_c4 = (double *)malloc(nc * sizeof(double)); + double *errors_c4 = (double *)malloc(nc * sizeof(double)); char *out_filename = (char *)malloc(100 * sizeof(char)); uint_t out_filename_len = 0; @@ -114,6 +116,11 @@ int main(int argc, char *argv[]) { get_mom(dist_len, 3, dist, mom1, mom3); vals_c3[i] = mom3[0]; errors_c3[i] = mom3[1]; + + double mom4[2]; + get_mom(dist_len, 4, dist, mom1, mom4); + vals_c4[i] = mom4[0]; + errors_c4[i] = mom4[1]; } free(dist); } @@ -127,7 +134,7 @@ int main(int argc, char *argv[]) { FILE *out_file = fopen(out_filename, "w"); for (uint_t i = 0; i < nc; i++) { - fprintf(out_file, "%u %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i]); + fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]); } fclose(out_file); @@ -141,6 +148,8 @@ int main(int argc, char *argv[]) { free(errors_c2); free(vals_c3); free(errors_c3); + free(vals_c4); + free(errors_c4); return 0; } diff --git a/src/break_edge.c b/src/break_edge.c index 01570d8..538b88d 100644 --- a/src/break_edge.c +++ b/src/break_edge.c @@ -1,13 +1,21 @@ #include "fracture.h" -bool break_edge(net_t *instance, uint_t edge, cholmod_common *c) { +bool break_edge(net_t *instance, uint_t edge, cholmod_common *c, bool refactor) { instance->fuses[edge] = true; if (instance->factor != NULL) { - uint_t w1 = instance->graph->ev_break[2 * edge]; - uint_t w2 = instance->graph->ev_break[2 * edge + 1]; - factor_update(instance->factor, w1, w2, c); + if (refactor) { + cholmod_sparse *laplacian = gen_laplacian(instance, c, true); + CHOL_F(free_factor)(&instance->factor, c); + instance->factor = CHOL_F(analyze)(laplacian, c); + CHOL_F(factorize)(laplacian, instance->factor, c); + CHOL_F(free_sparse)(&laplacian, c); + } else { + uint_t w1 = instance->graph->ev_break[2 * edge]; + uint_t w2 = instance->graph->ev_break[2 * edge + 1]; + factor_update(instance->factor, w1, w2, c); + } } uint_t v1, v2, s1, s2, dv1, dv2, ds1, ds2; diff --git a/src/cen_anal_process.c b/src/cen_anal_process.c new file mode 100644 index 0000000..98f68ac --- /dev/null +++ b/src/cen_anal_process.c @@ -0,0 +1,154 @@ + +#include "fracture.h" +#include <gsl/gsl_sf_erf.h> +#include <gsl/gsl_sf_laguerre.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_blas.h> + +void get_mean(uint_t len, uint32_t *data, double result[2]) { + uint_t total = 0; + double mean = 0; + double mean_err = 0; + + for (uint_t i = 0; i < len; i++) { + uint32_t datai = data[i]; + + total += datai; + mean += i * datai; + mean_err += gsl_pow_2(i) * datai; + } + + double meanf = mean / total; + double meanf_err = meanf * sqrt(mean_err / gsl_pow_2(mean) + 1 / total); + + result[0] = meanf; + result[1] = meanf_err; +} + +void get_mom(uint_t len, uint_t n, uint32_t *data, double mean[2], double result[2]) { + uint_t total = 0; + double mom = 0; + double mom_err = 0; + + for (uint_t i = 0; i < len; i++) { + uint32_t datai = data[i]; + double in = pow(fabs(((double)i) - mean[0]), n); + + total += datai; + mom += in * datai; + mom_err += gsl_pow_2(in) * datai; + } + + double momf = mom / total; + double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total); + + result[0] = momf; + result[1] = momf_err; +} + +int main(int argc, char *argv[]) { + uint_t nc = argc - 1; + uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t)); + double *betas_c = (double *)malloc(nc * sizeof(double)); + double *vals_c1 = (double *)malloc(nc * sizeof(double)); + double *errors_c1 = (double *)malloc(nc * sizeof(double)); + double *vals_c2 = (double *)malloc(nc * sizeof(double)); + double *errors_c2 = (double *)malloc(nc * sizeof(double)); + double *vals_c3 = (double *)malloc(nc * sizeof(double)); + double *errors_c3 = (double *)malloc(nc * sizeof(double)); + double *vals_c4 = (double *)malloc(nc * sizeof(double)); + double *errors_c4 = (double *)malloc(nc * sizeof(double)); + + char *out_filename = (char *)malloc(100 * sizeof(char)); + uint_t out_filename_len = 0; + + for (uint_t i = 0; i < nc; i++) { + char *fn = argv[1 + i]; + char *buff = (char *)malloc(20 * sizeof(char)); + uint_t pos = 0; uint_t c = 0; + while (c < 5) { + if (fn[pos] == '_') { + c++; + } + if (i == 0) { + out_filename[pos] = fn[pos]; + out_filename_len++; + } + pos++; + } + c = 0; + while (fn[pos] != '_') { + buff[c] = fn[pos]; + pos++; + c++; + } + buff[c] = '\0'; + Ls_c[i] = atoi(buff); + c = 0; + pos++; + while (fn[pos] != '_') { + buff[c] = fn[pos]; + pos++; + c++; + } + buff[c] = '\0'; + betas_c[i] = atof(buff); + free(buff); + + uint_t dist_len = 4 * pow(Ls_c[i], 2); + uint32_t *dist = malloc(dist_len * sizeof(uint32_t)); + FILE *dist_file = fopen(fn, "rb"); + fread(dist, sizeof(uint32_t), dist_len, dist_file); + fclose(dist_file); + { + double mom1[2]; + get_mean(dist_len, dist, mom1); + vals_c1[i] = mom1[0]; + errors_c1[i] = mom1[1]; + + double mom2[2]; + get_mom(dist_len, 2, dist, mom1, mom2); + vals_c2[i] = mom2[0]; + errors_c2[i] = mom2[1]; + + double mom3[2]; + get_mom(dist_len, 3, dist, mom1, mom3); + vals_c3[i] = mom3[0]; + errors_c3[i] = mom3[1]; + + double mom4[2]; + get_mom(dist_len, 4, dist, mom1, mom4); + vals_c4[i] = mom4[0]; + errors_c4[i] = mom4[1]; + } + free(dist); + } + + out_filename[out_filename_len-1] = '.'; + out_filename[out_filename_len] = 't'; + out_filename[out_filename_len+1] = 'x'; + out_filename[out_filename_len+2] = 't'; + out_filename[out_filename_len+3] = '\0'; + + FILE *out_file = fopen(out_filename, "w"); + + for (uint_t i = 0; i < nc; i++) { + fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]); + } + + fclose(out_file); + + + free(Ls_c); + free(betas_c); + free(vals_c1); + free(errors_c1); + free(vals_c2); + free(errors_c2); + free(vals_c3); + free(errors_c3); + + return 0; +} + diff --git a/src/corr_test.c b/src/corr_test.c index b6f8a18..4dc7efc 100644 --- a/src/corr_test.c +++ b/src/corr_test.c @@ -9,7 +9,7 @@ int main() { graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 128, false, &c); net_t *instance = net_create(network, 1e-14, 3, 0, true, &c); - net_fracture(instance, &c, 1e-10); + net_fracture(instance, &c, 1e-10, 40); double *corr = get_corr(instance, NULL, &c); diff --git a/src/fracture.c b/src/fracture.c index 11e6efd..7ee82fb 100644 --- a/src/fracture.c +++ b/src/fracture.c @@ -8,7 +8,7 @@ int main(int argc, char *argv[]) { // defining variables to be (potentially) set by command line flags uint8_t filename_len; uint32_t N; - uint_t L; + uint_t L, refactor_every; double beta, inf, cutoff, crack_len; bool crack_growth_crit, save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network, save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_energy, save_conductivity, @@ -25,6 +25,7 @@ int main(int argc, char *argv[]) { // set default values N = 100; + refactor_every = UINT_MAX; L = 16; crack_len = 0.; beta = .3; @@ -45,7 +46,7 @@ int main(int argc, char *argv[]) { save_damage_field = false; save_conductivity = false; save_toughness = false; - save_energy = false; + save_energy = true; save_threshold = false; @@ -58,11 +59,14 @@ int main(int argc, char *argv[]) { // get commandline options - while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:gTE")) != -1) { + while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:gTER:")) != -1) { switch (opt) { case 'n': N = atoi(optarg); break; + case 'R': + refactor_every = atoi(optarg); + break; case 'L': L = atoi(optarg); break; @@ -303,7 +307,7 @@ int main(int argc, char *argv[]) { graph_t *g = graph_create(lattice, boundary, L, use_dual, &c); net_t *net = 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); + data_t *data = net_fracture(tmp_net, &c, cutoff, refactor_every); net_free(tmp_net, &c); uint_t max_pos = 0; @@ -357,7 +361,8 @@ int main(int argc, char *argv[]) { break; } - break_edge(net, next_broken, &c); + bool refactor = ((j + 1) % refactor_every) == 0; + break_edge(net, next_broken, &c, refactor); double val = data->extern_field[j]; if (save_cluster_dist) { diff --git a/src/fracture.h b/src/fracture.h index 0b078ff..f122869 100644 --- a/src/fracture.h +++ b/src/fracture.h @@ -72,7 +72,7 @@ typedef struct { typedef struct { const graph_t *graph; bool *fuses; - double *thres; + long double *thres; double inf; cholmod_dense *boundary_cond; cholmod_factor *factor; @@ -117,7 +117,7 @@ double dual_vert_to_coord(uint_t width, bool periodic, uint_t vert, void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c); void net_notch(net_t *net, double notch_len, cholmod_common *c); -data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff); +data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refactor_every); double *net_voltages(const net_t *net, cholmod_common *c); double *net_currents(const net_t *net, const double *voltages, cholmod_common *c); double net_conductivity(const net_t *net, const double *voltages); @@ -148,7 +148,7 @@ graph_t *ini_voro_graph(uint_t L, bound_t boundary, bool use_dual, double *(*genfunc)(uint_t, bound_t, gsl_rng *, uint_t *), cholmod_common *c); -bool break_edge(net_t *instance, uint_t edge, cholmod_common *c); +bool break_edge(net_t *instance, uint_t edge, cholmod_common *c, bool refactor); void finish_instance(net_t *instance, cholmod_common *c); @@ -183,4 +183,4 @@ bool set_connected(const cholmod_sparse *laplacian, uint_t *marks, int vertex, i unsigned long int rand_seed(); -double rand_dist_pow(const gsl_rng *r, double beta); +long double rand_dist_pow(const gsl_rng *r, double beta); @@ -1,8 +1,8 @@ #include "fracture.h" -double *get_thres(uint_t ne, double beta) { - double *thres = (double *)malloc(ne * sizeof(double)); +long double *get_thres(uint_t ne, double beta) { + long double *thres = (long double *)malloc(ne * sizeof(long double)); assert(thres != NULL); gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN); @@ -39,7 +39,7 @@ void net_notch(net_t *net, double notch_len, cholmod_common *c) { correct_length = v1x < notch_len && v2x < notch_len; if (crosses_center && not_wrapping && correct_length) { - break_edge(net, i, c); + break_edge(net, i, c, false); } } } @@ -96,8 +96,8 @@ net_t *net_copy(const net_t *net, cholmod_common *c) { assert(net_copy->fuses != NULL); memcpy(net_copy->fuses, net->fuses, fuses_size); - size_t thres_size = (net->graph)->ne * sizeof(double); - net_copy->thres = (double *)malloc(thres_size); + size_t thres_size = (net->graph)->ne * sizeof(long double); + net_copy->thres = (long double *)malloc(thres_size); assert(net_copy->thres != NULL); memcpy(net_copy->thres, net->thres, thres_size); diff --git a/src/net_fracture.c b/src/net_fracture.c index dcf67c9..3bfcd83 100644 --- a/src/net_fracture.c +++ b/src/net_fracture.c @@ -3,14 +3,14 @@ uint_t get_next_broken(net_t *net, double *currents, double cutoff) { uint_t max_pos = UINT_MAX; - double max_val = 0; + long double max_val = 0; for (uint_t i = 0; i < net->graph->ne; i++) { - double current = fabs(currents[i]); + long double current = fabs(currents[i]); bool broken = net->fuses[i]; if (!broken && current > cutoff) { - double val = current / net->thres[i]; + long double val = current / net->thres[i]; if (val > max_val) { max_val = val; @@ -27,10 +27,12 @@ uint_t get_next_broken(net_t *net, double *currents, double cutoff) { return max_pos; } -data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) { +data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refactor_every) { data_t *data = data_create(net->graph->ne); + uint_t n = 0; while (true) { + n++; double *voltages = net_voltages(net, c); double *currents = net_currents(net, voltages, c); @@ -52,12 +54,13 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) { sim_current = 1; } - data_update(data, last_broke, fabs(sim_current * (net->thres)[last_broke] / currents[last_broke]), conductivity); + data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / currents[last_broke]), conductivity); free(voltages); free(currents); - break_edge(net, last_broke, c); + bool refactor = (n % refactor_every) == 0; + break_edge(net, last_broke, c, refactor); if (net->num_components > 1 && net->graph->boundary == TORUS_BOUND) { break; @@ -9,12 +9,15 @@ unsigned long int rand_seed() { return seed; } -double rand_dist_pow(const gsl_rng *r, double beta) { - double x = 0; +long double rand_dist_pow(const gsl_rng *r, double beta) { + long double x = 0; // underflow means that for very small beta x is sometimes identically zero, // which causes problems - while ((x = exp(log(gsl_rng_uniform_pos(r)) / beta)) == 0.0); + while (x == 0.0) { + long double y = logl(gsl_rng_uniform_pos(r)) / beta; + x = expl(y); + } return x; } |