simulation.c - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
HTML git clone git://src.adamsgaard.dk/cngf-pf
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
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; }