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