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