simulation.c (23251B)
1 #include <stdio.h> 2 #include <stdlib.h> 3 #include <math.h> 4 #include "arrays.h" 5 #include "simulation.h" 6 #include "fluid.h" 7 8 9 /* #define SHOW_PARAMETERS */ 10 11 /* relative tolerance criteria for granular fluidity solver */ 12 #define RTOL_GRANULAR 1e-5 13 #define MAX_ITER_GRANULAR 10000 14 15 /* relative tolerance criteria for fluid-pressure solver */ 16 #define RTOL_DARCY 1e-5 17 #define MAX_ITER_DARCY 10000 18 19 /* relative tolerance criteria when shear velocity is restricted */ 20 #define RTOL_STRESS 1e-3 21 #define MAX_ITER_STRESS 20000 22 23 24 /* Simulation settings */ 25 void 26 init_sim(struct simulation *sim) 27 { 28 snprintf(sim->name, sizeof(sim->name), DEFAULT_SIMULATION_NAME); 29 30 sim->G = 9.81; 31 32 sim->P_wall = 120e3; 33 sim->mu_wall = 0.45; 34 sim->v_x_bot = 0.0; 35 sim->v_x_fix = NAN; 36 sim->v_x_limit = NAN; 37 38 sim->nz = -1; /* cell size equals grain size if negative */ 39 40 /* lower values of A mean that the velocity curve can have sharper 41 * curves, e.g. at the transition from μ ≈ μ_s */ 42 sim->A = 0.40; /* Loose fit to Damsgaard et al 2013 */ 43 44 /* lower values of b mean larger shear velocity for a given stress 45 * ratio above yield */ 46 sim->b = 0.9377; /* Henann and Kamrin 2016 */ 47 48 /* Henann and Kamrin 2016 */ 49 /* sim->mu_s = 0.3819; */ 50 /* sim->C = 0.0; */ 51 52 /* Damsgaard et al 2013 */ 53 sim->mu_s = tan(DEG2RAD(22.0)); 54 sim->C = 0.0; 55 sim->phi = initval(0.25, 1); 56 sim->d = 0.04; /* Damsgaard et al 2013 */ 57 58 sim->transient = 0; 59 sim->phi_min = 0.20; 60 sim->phi_max = 0.55; 61 sim->dilatancy_angle = 1.0; 62 63 /* Iverson et al 1997, 1998: Storglaciaren till */ 64 /* sim->mu_s = tan(DEG2RAD(26.3)); */ 65 /* sim->C = 5.0e3; */ 66 /* sim->phi = initval(0.22, 1); */ 67 /* sim->d = ??; */ 68 69 /* Iverson et al 1997, 1998: Two Rivers till */ 70 /* sim->mu_s = tan(DEG2RAD(17.8)); */ 71 /* sim->C = 14.0e3; */ 72 /* sim->phi = initval(0.37, 1); */ 73 /* sim->d = ??; */ 74 75 /* Tulaczyk et al 2000a: Upstream B till */ 76 /* sim->mu_s = tan(DEG2RAD(23.9)); */ 77 /* sim->C = 3.0e3; */ 78 /* sim->phi = initval(0.35, 1); */ 79 /* sim->d = ??; */ 80 81 /* grain material density [kg/m^3] */ 82 sim->rho_s = 2.6e3; /* Damsgaard et al 2013 */ 83 84 /* spatial settings */ 85 sim->origo_z = 0.0; 86 sim->L_z = 1.0; 87 88 /* temporal settings */ 89 sim->t = 0.0; 90 sim->dt = 2.0; 91 sim->t_end = 1.0; 92 sim->file_dt = 1.0; 93 sim->n_file = 0; 94 95 /* fluid settings */ 96 sim->fluid = 0; 97 98 sim->rho_f = 1e3; 99 100 /* Water at 20 deg C, Goren et al 2011 */ 101 /* sim->beta_f = 4.5e-10; */ 102 /* sim->mu_f = 1.0-3; */ 103 104 /* Water at 0 deg C */ 105 sim->beta_f = 3.9e-10; /* doi:10.1063/1.1679903 */ 106 sim->mu_f = 1.787e-3; /* Cuffey and Paterson 2010 */ 107 108 sim->alpha = 1e-8; 109 110 sim->D = -1.0; /* disabled when negative */ 111 112 /* Damsgaard et al 2015 */ 113 sim->k = initval(1.9e-15, 1); 114 115 /* Iverson et al 1994: Storglaciaren */ 116 /* sim->k = initval(1.3e-14, 1); */ 117 118 /* Engelhardt et al 1990: Upstream B */ 119 /* sim->k = initval(2.0e-16, 1); */ 120 121 /* Leeman et al 2016: Upstream B till */ 122 /* sim->k = initval(4.9e-17, 1); */ 123 124 /* no fluid-pressure variations */ 125 sim->p_f_top = 0.0; 126 sim->p_f_mod_ampl = 0.0; 127 sim->p_f_mod_freq = 1.0; 128 sim->p_f_mod_phase = 0.0; 129 sim->p_f_mod_pulse_time = NAN; 130 sim->p_f_mod_pulse_shape = 0; 131 } 132 133 void 134 prepare_arrays(struct simulation *sim) 135 { 136 if (sim->nz < 2) { 137 fprintf(stderr, "error: grid size (nz) must be at least 2 but is %d\n", 138 sim->nz); 139 exit(1); 140 } 141 free(sim->phi); 142 free(sim->k); 143 144 sim->z = linspace(sim->origo_z, /* spatial coordinates */ 145 sim->origo_z + sim->L_z, 146 sim->nz); 147 sim->dz = sim->z[1] - sim->z[0]; /* cell spacing */ 148 sim->mu = zeros(sim->nz); /* stress ratio */ 149 sim->mu_c = zeros(sim->nz); /* critical-state stress ratio */ 150 sim->sigma_n_eff = zeros(sim->nz); /* effective normal stress */ 151 sim->sigma_n = zeros(sim->nz); /* normal stess */ 152 sim->p_f_ghost = zeros(sim->nz + 2); /* fluid pressure with ghost nodes */ 153 sim->phi = zeros(sim->nz); /* porosity */ 154 sim->phi_c = zeros(sim->nz); /* critical-state porosity */ 155 sim->phi_dot = zeros(sim->nz); /* rate of porosity change */ 156 sim->k = zeros(sim->nz);/* permeability */ 157 sim->xi = zeros(sim->nz); /* cooperativity length */ 158 sim->gamma_dot_p = zeros(sim->nz); /* shear velocity */ 159 sim->v_x = zeros(sim->nz); /* shear velocity */ 160 sim->g_local = zeros(sim->nz); /* local fluidity */ 161 sim->g_ghost = zeros(sim->nz + 2); /* fluidity with ghost nodes */ 162 sim->I = zeros(sim->nz);/* inertia number */ 163 sim->tan_psi = zeros(sim->nz); /* tan(dilatancy_angle) */ 164 } 165 166 void 167 free_arrays(struct simulation *sim) 168 { 169 free(sim->z); 170 free(sim->mu); 171 free(sim->mu_c); 172 free(sim->sigma_n_eff); 173 free(sim->sigma_n); 174 free(sim->p_f_ghost); 175 free(sim->k); 176 free(sim->phi); 177 free(sim->phi_c); 178 free(sim->phi_dot); 179 free(sim->xi); 180 free(sim->gamma_dot_p); 181 free(sim->v_x); 182 free(sim->g_local); 183 free(sim->g_ghost); 184 free(sim->I); 185 free(sim->tan_psi); 186 } 187 188 static void 189 warn_parameter_value(const char message[], 190 const double value, 191 int *return_status) 192 { 193 fprintf(stderr, 194 "check_simulation_parameters: %s (%.17g)\n", 195 message, value); 196 *return_status = 1; 197 } 198 199 static void 200 check_float(const char name[], const double value, int *return_status) 201 { 202 #ifdef SHOW_PARAMETERS 203 printf("%30s: %.17g\n", name, value); 204 #endif 205 if (isnan(value)) { 206 char message[100]; 207 208 snprintf(message, sizeof(message), "%s is NaN", name); 209 warn_parameter_value(message, value, return_status); 210 *return_status = 1; 211 } else if (isinf(value)) { 212 char message[100]; 213 214 snprintf(message, sizeof(message), "%s is infinite", name); 215 warn_parameter_value(message, value, return_status); 216 *return_status = 1; 217 } 218 } 219 220 void 221 check_simulation_parameters(struct simulation *sim) 222 { 223 int return_status = 0; 224 225 check_float("sim->G", sim->G, &return_status); 226 if (sim->G < 0.0) 227 warn_parameter_value("sim->G is negative", sim->G, &return_status); 228 229 check_float("sim->P_wall", sim->P_wall, &return_status); 230 if (sim->P_wall < 0.0) 231 warn_parameter_value("sim->P_wall is negative", sim->P_wall, 232 &return_status); 233 234 check_float("sim->v_x_bot", sim->v_x_bot, &return_status); 235 236 check_float("sim->mu_wall", sim->mu_wall, &return_status); 237 if (sim->mu_wall < 0.0) 238 warn_parameter_value("sim->mu_wall is negative", sim->mu_wall, 239 &return_status); 240 241 check_float("sim->A", sim->A, &return_status); 242 if (sim->A < 0.0) 243 warn_parameter_value("sim->A is negative", sim->A, &return_status); 244 245 check_float("sim->b", sim->b, &return_status); 246 if (sim->b < 0.0) 247 warn_parameter_value("sim->b is negative", sim->b, &return_status); 248 249 check_float("sim->mu_s", sim->mu_s, &return_status); 250 if (sim->mu_s < 0.0) 251 warn_parameter_value("sim->mu_s is negative", sim->mu_s, 252 &return_status); 253 254 check_float("sim->C", sim->C, &return_status); 255 256 check_float("sim->d", sim->d, &return_status); 257 if (sim->d <= 0.0) 258 warn_parameter_value("sim->d is not a positive number", sim->d, 259 &return_status); 260 261 check_float("sim->rho_s", sim->rho_s, &return_status); 262 if (sim->rho_s <= 0.0) 263 warn_parameter_value("sim->rho_s is not a positive number", sim->rho_s, 264 &return_status); 265 266 if (sim->nz <= 0) 267 warn_parameter_value("sim->nz is not a positive number", sim->nz, 268 &return_status); 269 270 check_float("sim->origo_z", sim->origo_z, &return_status); 271 check_float("sim->L_z", sim->L_z, &return_status); 272 if (sim->L_z <= sim->origo_z) 273 warn_parameter_value("sim->L_z is smaller or equal to sim->origo_z", 274 sim->L_z, &return_status); 275 276 if (sim->nz <= 0) 277 warn_parameter_value("sim->nz is not a positive number", sim->nz, 278 &return_status); 279 280 check_float("sim->dz", sim->dz, &return_status); 281 if (sim->dz <= 0.0) 282 warn_parameter_value("sim->dz is not a positive number", sim->dz, 283 &return_status); 284 285 check_float("sim->t", sim->t, &return_status); 286 if (sim->t < 0.0) 287 warn_parameter_value("sim->t is a negative number", 288 sim->t, &return_status); 289 290 check_float("sim->t_end", sim->t_end, &return_status); 291 if (sim->t > sim->t_end) 292 warn_parameter_value("sim->t_end is smaller than sim->t", 293 sim->t, &return_status); 294 295 check_float("sim->dt", sim->dt, &return_status); 296 if (sim->dt <= 0.0) 297 warn_parameter_value("sim->dt is not a positive number", 298 sim->dt, &return_status); 299 300 check_float("sim->file_dt", sim->file_dt, &return_status); 301 if (sim->file_dt < 0.0) 302 warn_parameter_value("sim->file_dt is a negative number", 303 sim->file_dt, &return_status); 304 305 check_float("sim->phi[0]", sim->phi[0], &return_status); 306 if (sim->phi[0] < 0.0 || sim->phi[0] > 1.0) 307 warn_parameter_value("sim->phi[0] is not within [0;1]", 308 sim->phi[0], &return_status); 309 310 check_float("sim->phi_min", sim->phi_min, &return_status); 311 if (sim->phi_min < 0.0 || sim->phi_min > 1.0) 312 warn_parameter_value("sim->phi_min is not within [0;1]", 313 sim->phi_min, &return_status); 314 315 check_float("sim->phi_max", sim->phi_max, &return_status); 316 if (sim->phi_max < 0.0 || sim->phi_max > 1.0) 317 warn_parameter_value("sim->phi_max is not within [0;1]", 318 sim->phi_max, &return_status); 319 320 check_float("sim->dilatancy_angle", sim->dilatancy_angle, &return_status); 321 if (sim->dilatancy_angle < 0.0 || sim->dilatancy_angle > 100.0) 322 warn_parameter_value("sim->dilatancy_angle is not within [0;100]", 323 sim->dilatancy_angle, &return_status); 324 325 if (sim->fluid != 0 && sim->fluid != 1) 326 warn_parameter_value("sim->fluid has an invalid value", 327 (double) sim->fluid, &return_status); 328 329 if (sim->transient != 0 && sim->transient != 1) 330 warn_parameter_value("sim->transient has an invalid value", 331 (double) sim->transient, &return_status); 332 333 if (sim->fluid) { 334 335 check_float("sim->p_f_mod_ampl", sim->p_f_mod_ampl, &return_status); 336 if (sim->p_f_mod_ampl < 0.0) 337 warn_parameter_value("sim->p_f_mod_ampl is not a zero or positive", 338 sim->p_f_mod_ampl, &return_status); 339 340 if (sim->P_wall - sim->p_f_mod_ampl < 0.0) 341 warn_parameter_value("sim->P_wall - sim->p_f_mod_ampl is negative", 342 sim->P_wall - sim->p_f_mod_ampl, 343 &return_status); 344 345 check_float("sim->p_f_mod_freq", sim->p_f_mod_freq, &return_status); 346 if (sim->p_f_mod_freq < 0.0) 347 warn_parameter_value("sim->p_f_mod_freq is not a zero or positive", 348 sim->p_f_mod_freq, &return_status); 349 350 check_float("sim->beta_f", sim->beta_f, &return_status); 351 if (sim->beta_f <= 0.0) 352 warn_parameter_value("sim->beta_f is not positive", 353 sim->beta_f, &return_status); 354 355 check_float("sim->alpha", sim->alpha, &return_status); 356 if (sim->alpha <= 0.0) 357 warn_parameter_value("sim->alpha is not positive", 358 sim->alpha, &return_status); 359 360 check_float("sim->mu_f", sim->mu_f, &return_status); 361 if (sim->mu_f <= 0.0) 362 warn_parameter_value("sim->mu_f is not positive", 363 sim->mu_f, &return_status); 364 365 check_float("sim->rho_f", sim->rho_f, &return_status); 366 if (sim->rho_f <= 0.0) 367 warn_parameter_value("sim->rho_f is not positive", 368 sim->rho_f, &return_status); 369 370 check_float("sim->k[0]", sim->k[0], &return_status); 371 if (sim->k[0] <= 0.0) 372 warn_parameter_value("sim->k[0] is not positive", 373 sim->k[0], &return_status); 374 } 375 if (return_status != 0) { 376 fprintf(stderr, "error: aborting due to invalid parameter choices\n"); 377 exit(return_status); 378 } 379 } 380 381 void 382 lithostatic_pressure_distribution(struct simulation *sim) 383 { 384 int i; 385 386 for (i = 0; i < sim->nz; ++i) 387 sim->sigma_n[i] = sim->P_wall + 388 (1.0 - sim->phi[i]) * sim->rho_s * sim->G * 389 (sim->L_z - sim->z[i]); 390 } 391 392 /* NEW FUNCTIONS START */ 393 394 void 395 compute_inertia_number(struct simulation *sim) 396 { 397 int i; 398 399 for (i = 0; i < sim->nz; ++i) 400 sim->I[i] = sim->gamma_dot_p[i] * sim->d 401 / sqrt(sim->sigma_n_eff[i] / sim->rho_s); 402 } 403 404 void 405 compute_critical_state_porosity(struct simulation *sim) 406 { 407 int i; 408 409 for (i = 0; i < sim->nz; ++i) 410 sim->phi_c[i] = sim->phi_min + (sim->phi_max - sim->phi_min) *sim->I[i]; 411 } 412 413 void 414 compute_critical_state_friction(struct simulation *sim) 415 { 416 int i; 417 418 if (sim->fluid) 419 for (i = 0; i < sim->nz; ++i) 420 sim->mu_c[i] = sim->mu_wall / 421 (sim->sigma_n_eff[i] / (sim->P_wall - sim->p_f_top)); 422 else 423 for (i = 0; i < sim->nz; ++i) 424 sim->mu_c[i] = sim->mu_wall / 425 (1.0 + (1.0 - sim->phi[i]) * sim->rho_s * sim->G * 426 (sim->L_z - sim->z[i]) / sim->P_wall); 427 } 428 429 void 430 compute_friction(struct simulation *sim) 431 { 432 int i; 433 434 if (sim->transient) 435 for (i = 0; i < sim->nz; ++i) 436 sim->mu[i] = sim->tan_psi[i] + sim->mu_c[i]; 437 else 438 for (i = 0; i < sim->nz; ++i) 439 sim->mu[i] = sim->mu_c[i]; 440 } 441 442 void 443 compute_tan_dilatancy_angle(struct simulation *sim) 444 { 445 int i; 446 447 for (i = 0; i < sim->nz; ++i) 448 sim->tan_psi[i] = sim->dilatancy_angle * (sim->phi_c[i] - sim->phi[i]); 449 } 450 451 void 452 compute_porosity_change(struct simulation *sim) 453 { 454 int i; 455 456 for (i = 0; i < sim->nz; ++i) { 457 sim->phi_dot[i] = sim->tan_psi[i] * sim->gamma_dot_p[i] * sim->phi[i]; 458 sim->phi[i] += sim->phi_dot[i] * sim->dt; 459 } 460 } 461 462 void 463 compute_permeability(struct simulation *sim) 464 { 465 int i; 466 467 for (i = 0; i < sim->nz; ++i) 468 sim->k[i] = pow(sim->d, 2.0) / 180.0 469 * pow(sim->phi[i], 3.0) 470 / pow(1.0 - sim->phi[i], 2.0); 471 } 472 473 /* NEW FUNCTIONS END */ 474 475 static double 476 shear_strain_rate_plastic(const double fluidity, const double friction) 477 { 478 return fluidity * friction; 479 } 480 481 void 482 compute_shear_strain_rate_plastic(struct simulation *sim) 483 { 484 int i; 485 486 for (i = 0; i < sim->nz; ++i) 487 sim->gamma_dot_p[i] = shear_strain_rate_plastic(sim->g_ghost[i + 1], 488 sim->mu[i]); 489 } 490 491 void 492 compute_shear_velocity(struct simulation *sim) 493 { 494 int i; 495 496 /* TODO: implement iterative solver for v_x from gamma_dot_p field */ 497 /* Dirichlet BC at bottom */ 498 sim->v_x[0] = sim->v_x_bot; 499 500 for (i = 1; i < sim->nz; ++i) 501 sim->v_x[i] = sim->v_x[i - 1] + sim->gamma_dot_p[i] * sim->dz; 502 } 503 504 void 505 compute_effective_stress(struct simulation *sim) 506 { 507 int i; 508 509 if (sim->fluid) 510 for (i = 0; i < sim->nz; ++i) 511 sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[i + 1]; 512 else 513 for (i = 0; i < sim->nz; ++i) 514 sim->sigma_n_eff[i] = sim->sigma_n[i]; 515 } 516 517 static double 518 cooperativity_length(const double A, 519 const double d, 520 const double mu, 521 const double p, 522 const double mu_s, 523 const double C) 524 { 525 return A * d / sqrt(fabs((mu - C / p) - mu_s)); 526 } 527 528 void 529 compute_cooperativity_length(struct simulation *sim) 530 { 531 int i; 532 533 for (i = 0; i < sim->nz; ++i) 534 sim->xi[i] = cooperativity_length(sim->A, 535 sim->d, 536 sim->mu[i], 537 sim->sigma_n_eff[i], 538 sim->mu_s, 539 sim->C); 540 } 541 542 static double 543 local_fluidity(const double p, 544 const double mu, 545 const double mu_s, 546 const double C, 547 const double b, 548 const double rho_s, 549 const double d) 550 { 551 if (mu - C / p <= mu_s) 552 return 0.0; 553 else 554 return sqrt(p / rho_s * d * d) * ((mu - C / p) - mu_s) / (b * mu); 555 } 556 557 void 558 compute_local_fluidity(struct simulation *sim) 559 { 560 int i; 561 562 for (i = 0; i < sim->nz; ++i) 563 sim->g_local[i] = local_fluidity(sim->sigma_n_eff[i], 564 sim->mu[i], 565 sim->mu_s, 566 sim->C, 567 sim->b, 568 sim->rho_s, 569 sim->d); 570 } 571 572 void 573 set_bc_neumann(double *a, 574 const int nz, 575 const int boundary, 576 const double df, 577 const double dx) 578 { 579 if (boundary == -1) 580 a[0] = a[1] + df * dx; 581 else if (boundary == +1) 582 a[nz + 1] = a[nz] - df * dx; 583 else { 584 fprintf(stderr, "set_bc_neumann: Unknown boundary %d\n", boundary); 585 exit(1); 586 } 587 } 588 589 void 590 set_bc_dirichlet(double *a, 591 const int nz, 592 const int boundary, 593 const double value) 594 { 595 if (boundary == -1) 596 a[0] = value; 597 else if (boundary == +1) 598 a[nz + 1] = value; 599 else { 600 fprintf(stderr, "set_bc_dirichlet: Unknown boundary %d\n", boundary); 601 exit(1); 602 } 603 } 604 605 double 606 residual(double new, double old) 607 { 608 return (new - old) / (old + 1e-16); 609 } 610 611 static void 612 poisson_solver_1d_cell_update(int i, 613 const double *g_in, 614 const double *g_local, 615 double *g_out, 616 double *r_norm, 617 const double dz, 618 const double *xi) 619 { 620 double coorp_term = dz * dz / (2.0 * pow(xi[i], 2.0)); 621 622 g_out[i + 1] = 1.0 / (1.0 + coorp_term) 623 * (coorp_term * g_local[i] + g_in[i + 2] / 2.0 + g_in[i] / 2.0); 624 625 r_norm[i] = residual(g_out[i + 1], g_in[i + 1]); 626 627 #ifdef DEBUG 628 printf("-- %d --------------\n", i); 629 printf("coorp_term: %g\n", coorp_term); 630 printf(" g_local: %g\n", g_local[i]); 631 printf(" g_in: %g\n", g_in[i + 1]); 632 printf(" g_out: %g\n", g_out[i + 1]); 633 printf(" r_norm: %g\n", r_norm[i]); 634 #endif 635 } 636 637 static int 638 implicit_1d_jacobian_poisson_solver(struct simulation *sim, 639 const int max_iter, 640 const double rel_tol) 641 { 642 int iter, i; 643 double r_norm_max = NAN; 644 double *tmp, *g_ghost_out = empty(sim->nz + 2), *r_norm = empty(sim->nz); 645 646 for (iter = 0; iter < max_iter; ++iter) { 647 #ifdef DEBUG 648 printf("\n@@@ ITERATION %d @@@\n", iter); 649 #endif 650 /* Dirichlet BCs resemble fixed particle velocities */ 651 set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0); 652 set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0); 653 654 /* Neumann BCs resemble free surfaces */ 655 /* set_bc_neumann(sim->g_ghost, sim->nz, +1, 0.0, sim->dz); */ 656 657 for (i = 0; i < sim->nz; ++i) 658 poisson_solver_1d_cell_update(i, 659 sim->g_ghost, 660 sim->g_local, 661 g_ghost_out, 662 r_norm, 663 sim->dz, 664 sim->xi); 665 r_norm_max = max(r_norm, sim->nz); 666 667 tmp = g_ghost_out; 668 g_ghost_out = sim->g_ghost; 669 sim->g_ghost = tmp; 670 671 if (r_norm_max <= rel_tol) { 672 set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0); 673 set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0); 674 free(g_ghost_out); 675 free(r_norm); 676 /* printf(".. Solution converged after %d 677 * iterations\n", iter); */ 678 return 0; 679 } 680 } 681 682 free(g_ghost_out); 683 free(r_norm); 684 fprintf(stderr, "implicit_1d_jacobian_poisson_solver: "); 685 fprintf(stderr, "Solution did not converge after %d iterations\n", iter); 686 fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max); 687 return 1; 688 } 689 690 void 691 write_output_file(struct simulation *sim, const int normalize) 692 { 693 char outfile[200]; 694 FILE *fp; 695 696 snprintf(outfile, sizeof(outfile), "%s.output%05d.txt", 697 sim->name, sim->n_file++); 698 699 if ((fp = fopen(outfile, "w")) != NULL) { 700 print_output(sim, fp, normalize); 701 fclose(fp); 702 } else { 703 fprintf(stderr, "could not open output file: %s", outfile); 704 exit(1); 705 } 706 } 707 708 void 709 print_output(struct simulation *sim, FILE *fp, const int norm) 710 { 711 int i; 712 double *v_x_out; 713 714 if (norm) 715 v_x_out = normalize(sim->v_x, sim->nz); 716 else 717 v_x_out = copy(sim->v_x, sim->nz); 718 719 for (i = 0; i < sim->nz; ++i) 720 fprintf(fp, 721 "%.17g\t%.17g\t%.17g\t" 722 "%.17g\t%.17g\t%.17g\t" 723 "%.17g\t%.17g\t%.17g\n", 724 sim->z[i], 725 v_x_out[i], 726 sim->sigma_n_eff[i], 727 sim->p_f_ghost[i + 1], 728 sim->mu[i], 729 sim->gamma_dot_p[i], 730 sim->phi[i], 731 sim->I[i], 732 sim->mu[i] * sim->sigma_n_eff[i]); 733 734 free(v_x_out); 735 } 736 737 int 738 coupled_shear_solver(struct simulation *sim, 739 const int max_iter, 740 const double rel_tol) 741 { 742 int i, coupled_iter, stress_iter = 0; 743 double *r_norm, *oldval; 744 double r_norm_max, vel_res_norm = NAN, mu_wall_orig = sim->mu_wall; 745 746 if (sim->transient) { 747 r_norm = empty(sim->nz); 748 oldval = empty(sim->nz); 749 } 750 do { /* stress iterations */ 751 752 coupled_iter = 0; 753 do { /* coupled iterations */ 754 755 if (sim->transient) { 756 copy_values(sim->gamma_dot_p, oldval, sim->nz); 757 758 /* step 1 */ 759 compute_inertia_number(sim); /* Eq. 1 */ 760 /* step 2 */ 761 compute_critical_state_porosity(sim); /* Eq. 2 */ 762 /* step 3 */ 763 compute_tan_dilatancy_angle(sim); /* Eq. 5 */ 764 } 765 compute_critical_state_friction(sim); /* Eq. 7 */ 766 767 /* step 4 */ 768 if (sim->transient) { 769 compute_porosity_change(sim); /* Eq. 3 */ 770 compute_permeability(sim); /* Eq. 6 */ 771 } 772 compute_friction(sim); /* Eq. 4 */ 773 774 /* step 5, Eq. 13 */ 775 if (sim->fluid) 776 if (darcy_solver_1d(sim, MAX_ITER_DARCY, RTOL_DARCY)) 777 exit(11); 778 779 /* step 6 */ 780 compute_effective_stress(sim); /* Eq. 9 */ 781 782 /* step 7 */ 783 compute_local_fluidity(sim); /* Eq. 10 */ 784 compute_cooperativity_length(sim); /* Eq. 12 */ 785 786 /* step 8, Eq. 11 */ 787 if (implicit_1d_jacobian_poisson_solver(sim, 788 MAX_ITER_GRANULAR, 789 RTOL_GRANULAR)) 790 exit(12); 791 792 /* step 9 */ 793 compute_shear_strain_rate_plastic(sim); /* Eq. 8 */ 794 compute_inertia_number(sim); 795 compute_shear_velocity(sim); 796 797 #ifdef DEBUG 798 for (i = 0; i < sim->nz; ++i) { 799 printf("\nsim->t = %g s\n", sim->t); 800 printf("sim->I[%d] = %g\n", i, sim->I[i]); 801 printf("sim->phi_c[%d] = %g\n", i, sim->phi_c[i]); 802 printf("sim->tan_psi[%d] = %g\n", i, sim->tan_psi[i]); 803 printf("sim->mu_c[%d] = %g\n", i, sim->mu_c[i]); 804 printf("sim->phi_dot[%d] = %g\n", i, sim->phi_dot[i]); 805 printf("sim->k[%d] = %g\n", i, sim->k[i]); 806 printf("sim->mu[%d] = %g\n", i, sim->mu[i]); 807 } 808 #endif 809 /* stable velocity field == coupled solution 810 * converged */ 811 if (sim->transient) { 812 for (i = 0; i < sim->nz; ++i) 813 r_norm[i] = fabs(residual(sim->gamma_dot_p[i], oldval[i])); 814 r_norm_max = max(r_norm, sim->nz); 815 if (r_norm_max <= rel_tol) 816 break; 817 818 if (coupled_iter++ >= max_iter) { 819 free(r_norm); 820 free(oldval); 821 fprintf(stderr, "coupled_shear_solver: "); 822 fprintf(stderr, "Transient solution did not converge " 823 "after %d iterations\n", coupled_iter); 824 fprintf(stderr, ".. Residual normalized error: %g\n", 825 r_norm_max); 826 return 1; 827 } 828 } 829 } while (sim->transient); 830 831 if (!isnan(sim->v_x_limit) || !isnan(sim->v_x_fix)) { 832 if (!isnan(sim->v_x_limit)) { 833 vel_res_norm = (sim->v_x_limit - sim->v_x[sim->nz - 1]) 834 / (sim->v_x[sim->nz - 1] + 1e-12); 835 if (vel_res_norm > 0.0) 836 vel_res_norm = 0.0; 837 } else { 838 vel_res_norm = (sim->v_x_fix - sim->v_x[sim->nz - 1]) 839 / (sim->v_x[sim->nz - 1] + 1e-12); 840 } 841 sim->mu_wall *= 1.0 + (vel_res_norm * 1e-2); 842 } 843 if (++stress_iter > MAX_ITER_STRESS) { 844 fprintf(stderr, "error: stress solution did not converge:\n"); 845 fprintf(stderr, "v_x=%g, v_x_fix=%g, v_x_limit=%g, " 846 "vel_res_norm=%g, mu_wall=%g\n", 847 sim->v_x[sim->nz - 1], sim->v_x_fix, sim->v_x_limit, 848 vel_res_norm, sim->mu_wall); 849 return 10; 850 } 851 } while ((!isnan(sim->v_x_fix) || !isnan(sim->v_x_limit)) 852 && fabs(vel_res_norm) > RTOL_STRESS); 853 854 if (!isnan(sim->v_x_limit)) 855 sim->mu_wall = mu_wall_orig; 856 857 if (sim->transient) { 858 free(r_norm); 859 free(oldval); 860 } 861 return 0; 862 } 863 864 double 865 find_flux(const struct simulation *sim) 866 { 867 int i; 868 double flux = 0.0; 869 870 for (i = 1; i < sim->nz; ++i) 871 flux += (sim->v_x[i] + sim->v_x[i - 1]) / 2.0 * sim->dz; 872 873 return flux; 874 }