Introduction
Introduction Statistics Contact Development Disclaimer Help
tsimulation.c - cngf-pf - continuum model for granular flows with pore-pressure…
git clone git://src.adamsgaard.dk/cngf-pf
Log
Files
Refs
README
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 solv…
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_…
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 i…
38
39 sim->A = 0.40; /* Loose fit to Damsga…
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 201…
47 sim->transient = 0;
48 sim->phi_min = 0.20;
49 sim->phi_max = 0.55;
50 sim->dilatancy_constant = 4.09; /* Pailha & Pouliquen 200…
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 …
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", na…
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…
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, &retu…
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_w…
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…
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, &retu…
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, &retu…
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", …
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 numbe…
266 &return_status);
267
268 if (sim->nz <= 0)
269 warn_parameter_value("sim->nz is not a positive number",…
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 si…
276 sim->L_z, &return_status);
277
278 if (sim->nz <= 0)
279 warn_parameter_value("sim->nz is not a positive number",…
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",…
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, …
323 if (sim->dilatancy_constant < 0.0 || sim->dilatancy_constant > 1…
324 warn_parameter_value("sim->dilatancy_constant is not wit…
325 sim->dilatancy_constant, &return_st…
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 valu…
333 (double) sim->transient, &return_st…
334
335 if (sim->fluid) {
336 check_float("sim->p_f_mod_ampl", sim->p_f_mod_ampl, &ret…
337 if (sim->p_f_mod_ampl < 0.0)
338 warn_parameter_value("sim->p_f_mod_ampl is not a…
339 sim->p_f_mod_ampl, &return_…
340
341 check_float("sim->p_f_mod_freq", sim->p_f_mod_freq, &ret…
342 if (sim->p_f_mod_freq < 0.0)
343 warn_parameter_value("sim->p_f_mod_freq is not a…
344 sim->p_f_mod_freq, &return_…
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 positiv…
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 …
374 "transient simulations",
375 sim->D, &return_status);
376 }
377
378 if (return_status != 0) {
379 fprintf(stderr, "error: aborting due to invalid paramete…
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 * s…
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 …
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], SIG…
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…
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], SIGM…
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]) * si…
437 (sim->L_z - sim->z[i]) / sim->…
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[…
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] …
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_g…
501 sim->mu[…
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 fie…
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] * si…
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 val…
525 sim->sigma_n_eff[i] = sim->sigma_n[i]
526 - ((sim->p_f_ghost[i + 1]
527 + sim->p_f_next_ghost[…
528 / 2.0);
529 /* sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->…
530 if (sim->sigma_n_eff[i] < 0)
531 warnx("%s: sigma_n_eff[%d] is negative w…
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[…
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)…
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->d…
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 iteration…
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…
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", outfil…
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_E…
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…
770 i, sim->t);
771 } else {
772 sim->p_f_ghost[i + 1] += sim->p_f_dot[i]…
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_wa…
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, …
798
799 /* step 1 */
800 compute_inertia_number(sim); /* E…
801 /* step 2 */
802 compute_critical_state_porosity(sim); …
803 /* step 3 */
804 compute_tan_dilatancy_angle(sim); …
805 }
806 compute_critical_state_friction(sim); /* …
807
808 /* step 4 */
809 if (sim->transient) {
810 compute_porosity_change(sim); /* …
811 compute_permeability(sim); /* Eq.…
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,…
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.…
826
827 /* step 8, Eq. 11 */
828 if (implicit_1d_jacobian_poisson_solver(sim,
829 MAX_ITER…
830 rel_tol))
831 exit(12);
832
833 /* step 9 */
834 compute_shear_strain_rate_plastic(sim); /…
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->…
844 printf("sim->tan_psi[%d] = %g\n", i, sim…
845 printf("sim->mu_c[%d] = %g\n", i, sim->m…
846 printf("sim->phi[%d] = %g\n", i, sim->ph…
847 printf("sim->phi_dot[%d] = %g\n", i, sim…
848 printf("sim->k[%d] = %g\n", i, sim->k[i]…
849 printf("sim->mu[%d] = %g\n", i, sim->mu[…
850 }
851 #endif
852
853 /* stable porosity change field == coupled solut…
854 if (sim->transient) {
855 for (i = 0; i < sim->nz; ++i)
856 sim->g_r_norm[i] = fabs(residual…
857 r_norm_max = max(sim->g_r_norm, sim->nz);
858 if (r_norm_max <= rel_tol && coupled_ite…
859 break;
860 if (coupled_iter++ >= max_iter) {
861 fprintf(stderr, "coupled_shear_s…
862 fprintf(stderr, "Transient solut…
863 "after %d iterations\n",…
864 fprintf(stderr, ".. Residual nor…
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_…
874 / (sim->v_x[sim->nz - 1] …
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[…
879 / (sim->v_x[sim->nz - 1] …
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 …
885 fprintf(stderr, "v_x=%g, v_x_fix=%g, v_x_limit=%…
886 "vel_res_norm=%g, mu_wall=%g\n",
887 sim->v_x[sim->nz - 1], sim->v_x_fix, sim…
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 doubl…
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_am…
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_…
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 …
943 sim->rho_s);
944
945 dt = xi * (sim->alpha + sim->phi[sim->nz - 1] * sim->beta_f) * s…
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 }
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.