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