From 2bb0740b68fdb62d45adc00204b3990ca42ade77 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 22 Aug 2016 10:11:14 -0400 Subject: started repo again without all the data files gunking the works --- src/current_scaling.c | 300 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 300 insertions(+) create mode 100644 src/current_scaling.c (limited to 'src/current_scaling.c') diff --git a/src/current_scaling.c b/src/current_scaling.c new file mode 100644 index 0000000..6837f6b --- /dev/null +++ b/src/current_scaling.c @@ -0,0 +1,300 @@ + +#include "fracture.h" + +int main(int argc, char *argv[]) { + int opt; + + // defining variables to be (potentially) set by command line flags + int iter = 1; + int num = 100; + int width = 16; + double crack_len = 8; + double beta = .3; + double inf = 1e-8; + double cutoff = 1e-8; + bool beta_shift = false; + bool supplied_bound = false; + bool ash_beta = false; + char *bound_file; + bool voltage_bound = false; + bool use_first = false; + bool save_stress = false; + bool save_bound = false; + bool save_damage = false; + bool save_strength = false; + + while ((opt = getopt(argc, argv, "n:w:b:l:i:Bf:aVFsSed")) != -1) { + switch (opt) { + case 'n': + num = atoi(optarg); + break; + case 'w': + width = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + case 'l': + crack_len = atof(optarg); + break; + case 'i': + iter = atoi(optarg); + break; + case 'B': + beta_shift = true; + break; + case 'a': + ash_beta = true; + break; + case 'V': + voltage_bound = true; + break; + case 'F': + use_first = true; + break; + case 's': + save_stress = true; + break; + case 'd': + save_damage = true; + break; + case 'e': + save_bound = true; + break; + case 'S': + save_strength = true; + break; + case 'f': + supplied_bound = true; + bound_file = optarg; + break; + default: /* '?' */ + exit(EXIT_FAILURE); + } + } + + // start cholmod + cholmod_common c; + CHOL_F(start)(&c); + + /* if we use voltage boundary conditions, the laplacian matrix is positive + * definite and we can use a supernodal LL decomposition. otherwise we need + * to use the simplicial LDL decomposition + */ + if (voltage_bound) { + (&c)->supernodal = CHOLMOD_SUPERNODAL; + } else { + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + } + + fnet *network = ini_square_network(width, false, true, &c); + finst *perm_instance = + create_instance(network, inf, voltage_bound, false, &c); + gen_crack(perm_instance, crack_len, 1, &c); + finish_instance(perm_instance, &c); + + if (voltage_bound) { + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + finst *tmp_instance = create_instance(network, inf, false, false, &c); + gen_crack(tmp_instance, crack_len, 1, &c); + finish_instance(tmp_instance, &c); + double *voltage = get_voltage(tmp_instance, &c); + + for (int i = 0; i < network->num_bounds; i++) { + for (int j = 0; j < network->bound_inds[i + 1] - network->bound_inds[i]; + j++) { + ((double *)perm_instance->boundary_cond + ->x)[network->bound_verts[network->bound_inds[i] + j]] = + voltage[network->bound_verts[network->bound_inds[i] + j]]; + } + } + (&c)->supernodal = CHOLMOD_SUPERNODAL; + } + + if (supplied_bound) { + FILE *bound_f = fopen(bound_file, "r"); + for (int i = 0; i < network->num_verts; i++) { + double tmp; + fscanf(bound_f, "%lg ", &tmp); + ((double *)perm_instance->boundary_cond->x)[i] = tmp; + } + + ((double *)perm_instance->boundary_cond->x)[network->num_verts] = 0; + ((double *)perm_instance->boundary_cond->x)[network->num_verts + 1] = 0; + fclose(bound_f); + } + + printf("\n"); + for (int DUMB2 = 0; DUMB2 < iter; DUMB2++) { + + double *strength; + if (save_strength) { + strength = (double *)malloc(num * sizeof(double)); + } + + double *damage; + if (save_damage) { + damage = (double *)calloc(network->num_edges, sizeof(double)); + } + double *avg_current = (double *)calloc(network->num_edges, sizeof(double)); + unsigned int *num_current_skipped = + (unsigned int *)calloc(network->num_edges, sizeof(unsigned int)); + double *avg_voltage = (double *)calloc(network->num_verts, sizeof(double)); + unsigned int *num_voltage_skipped = + (unsigned int *)calloc(network->num_verts, sizeof(unsigned int)); + + for (int DUMB = 0; DUMB < num; DUMB++) { + printf("\033[F\033[JCURRENT_SCALING: ITERATION %0*d: %0*d / %d\n", + (int)log10(iter) + 1, DUMB2 + 1, (int)log10(num) + 1, DUMB + 1, + num); + + break_data *breaking_data = NULL; + while (breaking_data == NULL) { + double *fuse_thres = gen_fuse_thres( + network->num_edges, network->edge_coords, beta, beta_scaling_flat); + finst *instance = copy_instance(perm_instance, &c); + breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); + free_instance(instance, &c); + free(fuse_thres); + } + + unsigned int min_pos = 0; + double min_val = DBL_MAX; + + for (unsigned int j = 0; j < breaking_data->num_broken; j++) { + double val = fabs(breaking_data->extern_field[j]); + if (val < min_val) { + min_pos = j; + min_val = val; + } + } + + if (save_strength) { + strength[DUMB] = fabs(breaking_data->extern_field[min_pos]); + } + + finst *tmp_instance = copy_instance(perm_instance, &c); + + unsigned int until = min_pos; + if (use_first) + until = 1; + for (unsigned int i = 0; i < until; i++) { + break_edge(tmp_instance, breaking_data->break_list[i], &c); + if (save_damage) { + damage[breaking_data->break_list[i]] += 1. / num; + } + } + + double *voltage = get_voltage(tmp_instance, &c); + double *current = get_current(tmp_instance, &c); + + for (unsigned int i = 0; i < network->num_edges; i++) { + avg_current[i] += current[i]; + if (current[i] == 0) + num_current_skipped[i]++; + } + + for (unsigned int i = 0; i < network->num_verts; i++) { + if (tmp_instance->marks[i] == tmp_instance->marks[network->num_verts]) { + avg_voltage[i] += voltage[i]; + } else { + num_voltage_skipped[i]++; + } + } + + free(current); + free(voltage); + free_instance(tmp_instance, &c); + free(breaking_data->break_list); + free(breaking_data->extern_field); + free(breaking_data); + } + + for (int i = 0; i < network->num_edges; i++) { + if (num_current_skipped[i] < num) { + avg_current[i] /= num - num_current_skipped[i]; + } + } + + for (int i = 0; i < network->num_verts; i++) { + if (num_voltage_skipped[i] < num) { + avg_voltage[i] /= num - num_voltage_skipped[i]; + } + } + + double *avg_field; + if (voltage_bound) + avg_field = avg_voltage; + else + avg_field = avg_current; + + update_boundary(perm_instance, avg_field); + + if (save_stress) { + char *c_filename = (char *)malloc(100 * sizeof(char)); + snprintf(c_filename, 100, "current_%d_%g_%d_%g.txt", width, crack_len, + iter, beta); + FILE *outfile = fopen(c_filename, "w"); + for (int i = 0; i < network->num_edges; i++) { + fprintf(outfile, "%g ", avg_current[i]); + } + fclose(outfile); + free(c_filename); + } + + if (save_damage) { + char *c_filename = (char *)malloc(100 * sizeof(char)); + snprintf(c_filename, 100, "damage_%d_%g_%d_%g.txt", width, crack_len, + iter, beta); + FILE *outfile = fopen(c_filename, "w"); + for (int i = 0; i < network->num_edges; i++) { + fprintf(outfile, "%g ", damage[i]); + } + fclose(outfile); + free(c_filename); + } + + if (save_strength) { + char *s_filename = (char *)malloc(100 * sizeof(char)); + snprintf(s_filename, 100, "strength_%d_%g_%d_%g.txt", width, crack_len, iter, beta); + FILE *f = fopen(s_filename, "a"); + for (int i = 0; i < num; i++) { + fprintf(f, "%g ", strength[i]); + } + fclose(f); + free(s_filename); + } + + if (save_bound) { + char *b_filename = (char *)malloc(100 * sizeof(char)); + snprintf(b_filename, 100, "bounds_%d_%g_%d_%g.txt", width, crack_len, + iter, beta); + FILE *outfile = fopen(b_filename, "w"); + for (int i = 0; i < network->num_verts; i++) { + fprintf(outfile, "%g ", ((double *)perm_instance->boundary_cond->x)[i]); + } + fclose(outfile); + free(b_filename); + } + + free(avg_current); + free(avg_voltage); + if (save_damage) free(damage); + free(num_current_skipped); + free(num_voltage_skipped); + if (save_strength) { + free(strength); + } + + printf( + "\033[F\033[JCURRENT_SCALING: ITERATION %0*d COMPLETE, BETA = %.2f\n\n", + (int)log10(iter) + 1, DUMB2 + 1, beta); + } + + free_instance(perm_instance, &c); + free_fnet(network, &c); + + CHOL_F(finish)(&c); + + return 0; +} -- cgit v1.2.3-70-g09d2