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