simulation.c (17572B)
1 #include <stdio.h> 2 #include <stdlib.h> 3 #include <math.h> 4 #include "arrays.h" 5 #include "simulation.h" 6 7 /* #define SHOW_PARAMETERS */ 8 9 void 10 prepare_arrays(struct simulation* sim) 11 { 12 sim->z = linspace(sim->origo_z, /* spatial coordinates */ 13 sim->origo_z + sim->L_z, 14 sim->nz); 15 sim->dz = sim->z[1] - sim->z[0]; /* cell spacing */ 16 sim->mu = zeros(sim->nz); /* stress ratio */ 17 sim->mu_c = zeros(sim->nz); /* critical-state stress ratio */ 18 sim->tau = zeros(sim->nz); /* shear stress */ 19 sim->tau_c = zeros(sim->nz); /* critical-state shear stress */ 20 sim->sigma_n_eff = zeros(sim->nz); /* effective normal stress */ 21 sim->sigma_n = zeros(sim->nz); /* normal stess */ 22 sim->p_f_ghost = zeros(sim->nz+2); /* fluid pressure with ghost nodes */ 23 free(sim->phi); 24 sim->phi = zeros(sim->nz); /* porosity */ 25 sim->phi_c = zeros(sim->nz); /* critical-state porosity */ 26 free(sim->k); 27 sim->k = zeros(sim->nz); /* permeability */ 28 sim->xi = zeros(sim->nz); /* cooperativity length */ 29 sim->gamma_dot_p = zeros(sim->nz); /* shear velocity */ 30 sim->v_x = zeros(sim->nz); /* shear velocity */ 31 sim->g_ghost = zeros(sim->nz+2); /* fluidity with ghost nodes */ 32 sim->I = zeros(sim->nz); /* inertia number */ 33 sim->tan_dilatancy_angle = zeros(sim->nz); /* tan(dilatancy_angle) */ 34 } 35 36 void 37 free_arrays(struct simulation* sim) 38 { 39 free(sim->z); 40 free(sim->mu); 41 free(sim->mu_c); 42 free(sim->tau); 43 free(sim->tau_c); 44 free(sim->sigma_n_eff); 45 free(sim->sigma_n); 46 free(sim->p_f_ghost); 47 free(sim->k); 48 free(sim->phi); 49 free(sim->phi_c); 50 free(sim->xi); 51 free(sim->gamma_dot_p); 52 free(sim->v_x); 53 free(sim->g_ghost); 54 free(sim->I); 55 free(sim->tan_dilatancy_angle); 56 } 57 58 static void 59 warn_parameter_value(const char message[], 60 const double value, 61 int* return_status) 62 { 63 fprintf(stderr, 64 "check_simulation_parameters: %s (%.17g)\n", 65 message, 66 value); 67 *return_status = 1; 68 } 69 70 static void 71 check_float(const char name[], const double value, int* return_status) 72 { 73 #ifdef SHOW_PARAMETERS 74 printf("%30s: %.17g\n", name, value); 75 #endif 76 if (isnan(value)) { 77 char message[100]; 78 snprintf(message, sizeof(message), "%s is NaN", name); 79 warn_parameter_value(message, value, return_status); 80 *return_status = 1; 81 } else if (isinf(value)) { 82 char message[100]; 83 snprintf(message, sizeof(message), "%s is infinite", name); 84 warn_parameter_value(message, value, return_status); 85 *return_status = 1; 86 } 87 } 88 89 void 90 check_simulation_parameters(const struct simulation* sim) 91 { 92 int return_status; 93 94 return_status = 0; 95 96 check_float("sim.G", sim->G, &return_status); 97 if (sim->G < 0.0) 98 warn_parameter_value("sim.G is negative", sim->G, &return_status); 99 100 check_float("sim.P_wall", sim->P_wall, &return_status); 101 if (sim->P_wall < 0.0) 102 warn_parameter_value("sim.P_wall is negative", sim->P_wall, 103 &return_status); 104 105 check_float("sim.v_x_bot", sim->v_x_bot, &return_status); 106 107 check_float("sim.mu_wall", sim->mu_wall, &return_status); 108 if (sim->mu_wall < 0.0) 109 warn_parameter_value("sim.mu_wall is negative", sim->mu_wall, 110 &return_status); 111 112 check_float("sim.A", sim->A, &return_status); 113 if (sim->A < 0.0) 114 warn_parameter_value("sim.A is negative", sim->A, &return_status); 115 116 check_float("sim.b", sim->b, &return_status); 117 if (sim->b < 0.0) 118 warn_parameter_value("sim.b is negative", sim->b, &return_status); 119 120 check_float("sim.mu_s", sim->mu_s, &return_status); 121 if (sim->mu_s < 0.0) 122 warn_parameter_value("sim.mu_s is negative", sim->mu_s, 123 &return_status); 124 125 check_float("sim.C", sim->C, &return_status); 126 if (sim->C < 0.0) 127 warn_parameter_value("sim.C is negative", sim->C, 128 &return_status); 129 130 check_float("sim.d", sim->d, &return_status); 131 if (sim->d <= 0.0) 132 warn_parameter_value("sim.d is not a positive number", sim->d, 133 &return_status); 134 135 check_float("sim.rho_s", sim->rho_s, &return_status); 136 if (sim->rho_s <= 0.0) 137 warn_parameter_value("sim.rho_s is not a positive number", sim->rho_s, 138 &return_status); 139 140 if (sim->nz <= 0) 141 warn_parameter_value("sim.nz is not a positive number", sim->nz, 142 &return_status); 143 144 check_float("sim.origo_z", sim->origo_z, &return_status); 145 check_float("sim.L_z", sim->L_z, &return_status); 146 if (sim->L_z <= sim->origo_z) 147 warn_parameter_value("sim.L_z is smaller or equal to sim.origo_z", 148 sim->L_z, &return_status); 149 150 if (sim->nz <= 0) 151 warn_parameter_value("sim.nz is not a positive number", sim->nz, 152 &return_status); 153 154 check_float("sim.dz", sim->dz, &return_status); 155 if (sim->dz <= 0.0) 156 warn_parameter_value("sim.dz is not a positive number", sim->dz, 157 &return_status); 158 159 check_float("sim.t", sim->t, &return_status); 160 if (sim->t < 0.0) 161 warn_parameter_value("sim.t is a negative number", 162 sim->t, &return_status); 163 164 check_float("sim.t_end", sim->t_end, &return_status); 165 if (sim->t > sim->t_end) 166 warn_parameter_value("sim.t_end is smaller than sim.t", 167 sim->t, &return_status); 168 169 check_float("sim.dt", sim->dt, &return_status); 170 if (sim->dt <= 0.0) 171 warn_parameter_value("sim.dt is not a positive number", 172 sim->dt, &return_status); 173 174 check_float("sim.file_dt", sim->file_dt, &return_status); 175 if (sim->file_dt < 0.0) 176 warn_parameter_value("sim.file_dt is a negative number", 177 sim->file_dt, &return_status); 178 179 check_float("sim.phi[0]", sim->phi[0], &return_status); 180 if (sim->phi[0] < 0.0 || sim->phi[0] > 1.0) 181 warn_parameter_value("sim.phi[0] is not within [0;1]", 182 sim->phi[0], &return_status); 183 184 check_float("sim.phi_min", sim->phi_min, &return_status); 185 if (sim->phi_min <= 0.0) 186 warn_parameter_value("sim.phi_min is not a positive number", 187 sim->phi_min, &return_status); 188 189 check_float("sim.phi_max", sim->phi_max, &return_status); 190 if (sim->phi_max <= 0.0) 191 warn_parameter_value("sim.phi_max is not a positive number", 192 sim->phi_max, &return_status); 193 194 check_float("sim.K", sim->K, &return_status); 195 if (sim->K <= 0.0) 196 warn_parameter_value("sim.K is not a positive number", sim->K, 197 &return_status); 198 199 if (sim->fluid != 0 && sim->fluid != 1) 200 warn_parameter_value("sim.fluid has an invalid value", 201 (double)sim->fluid, &return_status); 202 203 if (sim->fluid) { 204 205 check_float("sim.p_f_mod_ampl", sim->p_f_mod_ampl, &return_status); 206 if (sim->p_f_mod_ampl < 0.0) 207 warn_parameter_value("sim.p_f_mod_ampl is not a zero or positive", 208 sim->p_f_mod_ampl, &return_status); 209 210 if (sim->P_wall - sim->p_f_mod_ampl < 0.0) 211 warn_parameter_value("sim.P_wall - sim.p_f_mod_ampl is negative", 212 sim->P_wall - sim->p_f_mod_ampl, 213 &return_status); 214 215 check_float("sim.p_f_mod_freq", sim->p_f_mod_freq, &return_status); 216 if (sim->p_f_mod_freq < 0.0) 217 warn_parameter_value("sim.p_f_mod_freq is not a zero or positive", 218 sim->p_f_mod_freq, &return_status); 219 220 check_float("sim.beta_f", sim->beta_f, &return_status); 221 if (sim->beta_f <= 0.0) 222 warn_parameter_value("sim.beta_f is not positive", 223 sim->beta_f, &return_status); 224 225 check_float("sim.mu_f", sim->mu_f, &return_status); 226 if (sim->mu_f <= 0.0) 227 warn_parameter_value("sim.mu_f is not positive", 228 sim->mu_f, &return_status); 229 230 check_float("sim.rho_f", sim->rho_f, &return_status); 231 if (sim->rho_f <= 0.0) 232 warn_parameter_value("sim.rho_f is not positive", 233 sim->rho_f, &return_status); 234 235 check_float("sim.k[0]", sim->k[0], &return_status); 236 if (sim->k[0] <= 0.0) 237 warn_parameter_value("sim.k[0] is not positive", 238 sim->k[0], &return_status); 239 } 240 241 if (return_status != 0) { 242 fprintf(stderr, "error: aborting due to invalid parameter choices\n"); 243 exit(return_status); 244 } 245 } 246 247 void 248 lithostatic_pressure_distribution(struct simulation* sim) 249 { 250 int i; 251 for (i=0; i<sim->nz; ++i) 252 sim->sigma_n[i] = sim->P_wall + 253 (1.0 - sim->phi[i])*sim->rho_s*sim->G* 254 (sim->L_z - sim->z[i]); 255 } 256 257 /* NEW FUNCTIONS START */ 258 259 void 260 compute_inertia_number(struct simulation* sim) 261 { 262 int i; 263 for (i=0; i<sim->nz; ++i) 264 sim->I[i] = sim->gamma_dot_p[i]*sim->d 265 /sqrt(sim->sigma_n_eff[i]/sim->rho_s); 266 } 267 268 void 269 compute_critical_state_porosity(struct simulation* sim) 270 { 271 int i; 272 for (i=0; i<sim->nz; ++i) 273 sim->phi_c[i] = sim->phi_min + (sim->phi_max - sim->phi_min)*sim->I[i]; 274 } 275 276 void 277 compute_critical_state_friction(struct simulation* sim) 278 { 279 int i; 280 if (sim->fluid) 281 for (i=0; i<sim->nz; ++i) 282 sim->mu_c[i] = sim->mu_wall/ 283 (sim->sigma_n_eff[i]/(sim->P_wall - sim->p_f_top)); 284 else 285 for (i=0; i<sim->nz; ++i) 286 sim->mu_c[i] = sim->mu_wall/ 287 (1.0 + (1.0 - sim->phi[i])*sim->rho_s*sim->G* 288 (sim->L_z - sim->z[i])/sim->P_wall); 289 } 290 291 void 292 compute_friction(struct simulation* sim) 293 { 294 int i; 295 for (i=0; i<sim->nz; ++i) 296 sim->mu[i] = sim->tau[i]/sim->sigma_n_eff[i]; 297 } 298 299 void 300 compute_critical_state_shear_stress(struct simulation* sim) 301 { 302 int i; 303 for (i=0; i<sim->nz; ++i) 304 sim->tau_c[i] = sim->mu_c[i]*sim->sigma_n_eff[i]; 305 } 306 307 void 308 compute_shear_stress(struct simulation* sim) 309 { 310 int i; 311 for (i=0; i<sim->nz; ++i) 312 sim->tau[i] = sim->tan_dilatancy_angle[i]*sim->sigma_n_eff[i] + sim->tau_c[i]; 313 } 314 315 void 316 compute_tan_dilatancy_angle(struct simulation* sim) 317 { 318 int i; 319 for (i=0; i<sim->nz; ++i) 320 sim->tan_dilatancy_angle[i] = sim->K*(sim->phi_c[i] - sim->phi[i]); 321 } 322 323 void 324 compute_porosity_change(struct simulation* sim) 325 { 326 int i; 327 for (i=0; i<sim->nz; ++i) 328 sim->phi[i] += sim->tan_dilatancy_angle[i]*sim->gamma_dot_p[i]*sim->phi[i]*sim->dt; 329 } 330 331 void 332 compute_permeability(struct simulation* sim) 333 { 334 int i; 335 for (i=0; i<sim->nz; ++i) 336 sim->k[i] = pow(sim->d, 2.0)/180.0 * pow(sim->phi[i], 3.0)/pow(1.0 - sim->phi[i], 2.0); 337 } 338 339 /* NEW FUNCTIONS END */ 340 341 double 342 shear_strain_rate_plastic(const double fluidity, const double friction) 343 { 344 return fluidity*friction; 345 } 346 347 void 348 compute_shear_strain_rate_plastic(struct simulation* sim) 349 { 350 int i; 351 for (i=0; i<sim->nz; ++i) 352 sim->gamma_dot_p[i] = shear_strain_rate_plastic(sim->g_ghost[i+1], 353 sim->mu[i]); 354 } 355 356 void 357 compute_shear_velocity(struct simulation* sim) 358 { 359 int i; 360 361 // TODO: implement iterative solver 362 // Dirichlet BC at bottom 363 sim->v_x[0] = sim->v_x_bot; 364 365 for (i=1; i<sim->nz; ++i) 366 sim->v_x[i] = sim->v_x[i-1] + sim->gamma_dot_p[i]*sim->dz; 367 } 368 369 void 370 compute_effective_stress(struct simulation* sim) 371 { 372 int i; 373 if (sim->fluid) 374 for (i=0; i<sim->nz; ++i) 375 sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[i+1]; 376 else 377 for (i=0; i<sim->nz; ++i) 378 sim->sigma_n_eff[i] = sim->sigma_n[i]; 379 } 380 381 static double 382 cooperativity_length(const double A, 383 const double d, 384 const double mu, 385 const double p, 386 const double mu_s, 387 const double C) 388 { 389 return A*d/sqrt(fabs((mu - C/p) - mu_s)); 390 } 391 392 void 393 compute_cooperativity_length(struct simulation* sim) 394 { 395 int i; 396 for (i=0; i<sim->nz; ++i) 397 sim->xi[i] = cooperativity_length(sim->A, 398 sim->d, 399 sim->mu[i], 400 sim->sigma_n_eff[i], 401 sim->mu_s, 402 sim->C); 403 } 404 405 static double 406 local_fluidity(const double p, 407 const double mu, 408 const double mu_s, 409 const double C, 410 const double b, 411 const double rho_s, 412 const double d) 413 { 414 if (mu - C/p <= mu_s) 415 return 0.0; 416 else 417 return sqrt(p/rho_s*d*d) * ((mu - C/p) - mu_s)/(b*mu); 418 } 419 420 void 421 compute_local_fluidity(struct simulation* sim) 422 { 423 int i; 424 for (i=0; i<sim->nz; ++i) 425 sim->g_ghost[i+1] = local_fluidity(sim->sigma_n_eff[i], 426 sim->mu[i], 427 sim->mu_s, 428 sim->C, 429 sim->b, 430 sim->rho_s, 431 sim->d); 432 } 433 434 void 435 set_bc_neumann(double* g_ghost, 436 const int nz, 437 const int boundary, 438 const double df, 439 const double dx) 440 { 441 if (boundary == -1) 442 g_ghost[0] = g_ghost[1] + df*dx; 443 else if (boundary == +1) 444 g_ghost[nz+1] = g_ghost[nz] - df*dx; 445 else { 446 fprintf(stderr, "set_bc_neumann: Unknown boundary %d\n", boundary); 447 exit(1); 448 } 449 } 450 451 void 452 set_bc_dirichlet(double* g_ghost, 453 const int nz, 454 const int boundary, 455 const double value) 456 { 457 if (boundary == -1) 458 g_ghost[0] = value; 459 else if (boundary == +1) 460 g_ghost[nz+1] = value; 461 else { 462 fprintf(stderr, "set_bc_dirichlet: Unknown boundary %d\n", boundary); 463 exit(1); 464 } 465 } 466 467 static void 468 poisson_solver_1d_cell_update(int i, 469 const double* g_in, 470 double* g_out, 471 double* r_norm, 472 const double dz, 473 const double* mu, 474 const double* p, 475 const double* xi, 476 const double mu_s, 477 const double C, 478 const double b, 479 const double rho_s, 480 const double d) 481 { 482 double coorp_term; 483 484 coorp_term = dz*dz/(2.0*pow(xi[i], 2.0)); 485 g_out[i+1] = 1.0/(1.0 + coorp_term)*(coorp_term* 486 local_fluidity(p[i], mu[i], mu_s, C, b, rho_s, d) 487 + g_in[i+2]/2.0 488 + g_in[i]/2.0); 489 490 r_norm[i] = pow(g_out[i+1] - g_in[i+1], 2.0) / (pow(g_out[i+1], 2.0) + 1e-16); 491 492 #ifdef DEBUG 493 printf("-- %d --------------\n", i); 494 printf("coorp_term: %g\n", coorp_term); 495 printf(" g_local: %g\n", local_fluidity(p[i], mu[i], mu_s, b, rho_s, d)); 496 printf(" g_in: %g\n", g_in[i+1]); 497 printf(" g_out: %g\n", g_out[i+1]); 498 printf(" r_norm: %g\n", r_norm[i]); 499 #endif 500 } 501 502 int 503 implicit_1d_jacobian_poisson_solver(struct simulation* sim, 504 const int max_iter, 505 const double rel_tol) 506 { 507 double r_norm_max; 508 double *g_ghost_out, *r_norm; 509 int iter, i; 510 511 r_norm_max = NAN; 512 g_ghost_out = empty(sim->nz+2); 513 r_norm = empty(sim->nz); 514 515 for (iter=0; iter<max_iter; ++iter) { 516 #ifdef DEBUG 517 printf("\n@@@ ITERATION %d @@@\n", iter); 518 #endif 519 /* Dirichlet BCs resemble fixed particle velocities */ 520 set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0); 521 set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0); 522 523 /* Neumann BCs resemble free surfaces */ 524 /* set_bc_neumann(sim->g_ghost, sim->nz, +1, 0.0, sim->dz); */ 525 526 for (i=0; i<sim->nz; ++i) 527 poisson_solver_1d_cell_update(i, 528 sim->g_ghost, 529 g_ghost_out, 530 r_norm, 531 sim->dz, 532 sim->mu, 533 sim->sigma_n_eff, 534 sim->xi, 535 sim->mu_s, 536 sim->C, 537 sim->b, 538 sim->rho_s, 539 sim->d); 540 r_norm_max = max(r_norm, sim->nz); 541 542 copy_values(g_ghost_out, sim->g_ghost, sim->nz+2); 543 544 if (r_norm_max <= rel_tol) { 545 set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0); 546 set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0); 547 free(g_ghost_out); 548 free(r_norm); 549 /* printf(".. Solution converged after %d iterations\n", iter); */ 550 return 0; 551 } 552 } 553 554 free(g_ghost_out); 555 free(r_norm); 556 fprintf(stderr, "implicit_1d_jacobian_poisson_solver: "); 557 fprintf(stderr, "Solution did not converge after %d iterations\n", iter); 558 fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max); 559 return 1; 560 } 561 562 void 563 write_output_file(struct simulation* sim, const int normalize) 564 { 565 char outfile[200]; 566 FILE *fp; 567 568 snprintf(outfile, sizeof(outfile), "%s.output%05d.txt", 569 sim->name, sim->n_file++); 570 571 fp = fopen(outfile, "w"); 572 if (sim->fluid) 573 print_wet_output(fp, sim, normalize); 574 else 575 print_dry_output(fp, sim, normalize); 576 577 fclose(fp); 578 } 579 580 void 581 print_dry_output(FILE* fp, struct simulation* sim, const int norm) 582 { 583 int i; 584 double *v_x_out; 585 586 if (norm) 587 v_x_out = normalize(sim->v_x, sim->nz); 588 else 589 v_x_out = copy(sim->v_x, sim->nz); 590 591 for (i=0; i<sim->nz; ++i) 592 fprintf(fp, "%.17g\t%.17g\t%.17g\t%.17g\t%.17g\n", 593 sim->z[i], 594 v_x_out[i], 595 sim->sigma_n_eff[i], 596 sim->mu[i], 597 sim->gamma_dot_p[i]); 598 599 free(v_x_out); 600 } 601 602 void 603 print_wet_output(FILE* fp, struct simulation* sim, const int norm) 604 { 605 int i; 606 double *v_x_out; 607 608 if (norm) 609 v_x_out = normalize(sim->v_x, sim->nz); 610 else 611 v_x_out = copy(sim->v_x, sim->nz); 612 613 for (i=0; i<sim->nz; ++i) 614 fprintf(fp, "%.17g\t%.17g\t%.17g\t%.17g\t%.17g\t%.17g\n", 615 sim->z[i], 616 v_x_out[i], 617 sim->sigma_n_eff[i], 618 sim->p_f_ghost[i+1], 619 sim->mu[i], 620 sim->gamma_dot_p[i]); 621 622 free(v_x_out); 623 }