summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/anal_process.c9
-rw-r--r--src/big_anal_process.c11
-rw-r--r--src/break_edge.c16
-rw-r--r--src/cen_anal_process.c154
-rw-r--r--src/corr_test.c2
-rw-r--r--src/fracture.c15
-rw-r--r--src/fracture.h8
-rw-r--r--src/net.c10
-rw-r--r--src/net_fracture.c15
-rw-r--r--src/rand.c9
10 files changed, 219 insertions, 30 deletions
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);
diff --git a/src/net.c b/src/net.c
index f3aeda9..444e76c 100644
--- a/src/net.c
+++ b/src/net.c
@@ -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;
diff --git a/src/rand.c b/src/rand.c
index 1a6a3d4..75722ac 100644
--- a/src/rand.c
+++ b/src/rand.c
@@ -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;
}