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