tfluid.c - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
(HTM) git clone git://src.adamsgaard.dk/cngf-pf
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tfluid.c (13556B)
---
1 #include <stdlib.h>
2 #include <math.h>
3 #include <err.h>
4 #include "simulation.h"
5 #include "arrays.h"
6
7 void
8 hydrostatic_fluid_pressure_distribution(struct simulation *sim)
9 {
10 int i;
11
12 for (i = 0; i < sim->nz; ++i)
13 sim->p_f_ghost[i + 1] = sim->p_f_top
14 + sim->phi[i] * sim->rho_f * sim->G
15 * (sim->L_z - sim->z[i]);
16 }
17
18 static double
19 diffusivity(struct simulation *sim, int i) {
20 if (sim->D > 0.0)
21 return sim->D;
22 else
23 return sim->k[i]
24 / ((sim->alpha + sim->phi[i] * sim->beta_f) * sim->mu_f);
25 }
26
27 /* Determines the largest time step for the current simulation state. Note
28 * that the time step should be recalculated if cell sizes or
29 * diffusivities (i.e., permeabilities, porosities, viscosities, or
30 * compressibilities) change. The safety factor should be in ]0;1] */
31 int
32 set_largest_fluid_timestep(struct simulation *sim, const double safety)
33 {
34 int i;
35 double dx_min, diff, diff_max, *dx;
36
37 dx = spacing(sim->z, sim->nz);
38 dx_min = INFINITY;
39 for (i = 0; i < sim->nz - 1; ++i) {
40 if (dx[i] < 0.0) {
41 fprintf(stderr, "error: cell spacing negative (%g) in cell %d\n",
42 dx[i], i);
43 free(dx);
44 return 1;
45 }
46 if (dx[i] < dx_min)
47 dx_min = dx[i];
48 }
49 free(dx);
50
51 diff_max = -INFINITY;
52 for (i = 0; i < sim->nz; ++i) {
53 diff = diffusivity(sim, i);
54 if (diff > diff_max)
55 diff_max = diff;
56 }
57
58 sim->dt = safety * 0.5 * dx_min * dx_min / diff_max;
59 if (sim->file_dt * 0.5 < sim->dt)
60 sim->dt = sim->file_dt;
61
62 return 0;
63 }
64
65 static double
66 sine_wave(const double time,
67 const double ampl,
68 const double freq,
69 const double phase)
70 {
71 return ampl * sin(2.0 * PI * freq * time + phase);
72 }
73
74 static double
75 triangular_pulse(const double time,
76 const double peak_ampl,
77 const double freq,
78 const double peak_time)
79 {
80 if (peak_time - 1.0 / (2.0 * freq) < time && time <= peak_time)
81 return peak_ampl * 2.0 * freq * (time - peak_time) + peak_ampl;
82 else if (peak_time < time && time < peak_time + 1.0 / (2.0 * freq))
83 return peak_ampl * 2.0 * freq * (peak_time - time) + peak_ampl;
84 else
85 return 0.0;
86 }
87
88 static double
89 square_pulse(const double time,
90 const double peak_ampl,
91 const double freq,
92 const double peak_time)
93 {
94 if (peak_time - 1.0 / (2.0 * freq) < time &&
95 time < peak_time + 1.0 / (2.0 * freq))
96 return peak_ampl;
97 else
98 return 0.0;
99 }
100
101 static void
102 set_fluid_bcs(double *p_f_ghost, struct simulation *sim, const double p_f_top)
103 {
104 /* correct ghost node at top BC for hydrostatic pressure gradient */
105 set_bc_dirichlet(p_f_ghost, sim->nz, +1,
106 p_f_top - sim->phi[0] * sim->rho_f * sim->G * sim->dz);
107 p_f_ghost[sim->nz] = p_f_top; /* include top node in BC */
108 set_bc_neumann(p_f_ghost,
109 sim->nz,
110 -1,
111 sim->phi[0] * sim->rho_f * sim->G,
112 sim->dz);
113 }
114
115 static double
116 darcy_pressure_change_1d(const int i,
117 const int nz,
118 const double *p_f_ghost_in,
119 const double *phi,
120 const double *phi_dot,
121 const double *k,
122 const double dz,
123 const double beta_f,
124 const double alpha,
125 const double mu_f,
126 const double D)
127 {
128 double k_, div_k_grad_p, k_zn, k_zp;
129
130 if (D > 0.0)
131 return D * (p_f_ghost_in[i + 2]
132 - 2.0 * p_f_ghost_in[i + 1]
133 + p_f_ghost_in[i]) / (dz * dz);
134 else {
135 k_ = k[i];
136 if (i == 0)
137 k_zn = k_;
138 else
139 k_zn = k[i - 1];
140 if (i == nz - 1)
141 k_zp = k_;
142 else
143 k_zp = k[i + 1];
144
145 div_k_grad_p = (2.0 * k_zp * k_ / (k_zp + k_)
146 * (p_f_ghost_in[i + 2] - p_f_ghost_in[i + 1]) / dz
147 - 2.0 * k_zn * k_ / (k_zn + k_)
148 * (p_f_ghost_in[i + 1] - p_f_ghost_in[i]) / dz
149 ) / dz;
150
151 #ifdef DEBUG
152 printf("%s [%d]: phi=%g\tdiv_k_grad_p=%g\tphi_dot=%g\n",
153 __func__, i, phi[i], div_k_grad_p, phi_dot[i]);
154
155 printf(" p[%d, %d, %d] = [%g, %g, %g]\tk=[%g, %g, %g]\n",
156 i, i + 1, i + 2,
157 p_f_ghost_in[i], p_f_ghost_in[i + 1], p_f_ghost_in[i + 2],
158 k_zn, k_, k_zp);
159 #endif
160
161 return 1.0 / ((alpha + beta_f * phi[i]) * mu_f) * div_k_grad_p
162 - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[i])) * phi_dot[i];
163 }
164 }
165
166 static double
167 darcy_pressure_change_1d_impl(const int i,
168 const int nz,
169 const double dt,
170 const double *p_f_old_val,
171 const double *p_f_ghost_in,
172 double *p_f_ghost_out,
173 const double *phi,
174 const double *phi_dot,
175 const double *k,
176 const double dz,
177 const double beta_f,
178 const double alpha,
179 const double mu_f,
180 const double D)
181 {
182 double k_, div_k_grad_p, k_zn, k_zp, rhs_term, omega = 1.0;
183
184 if (D > 0.0)
185 return D * (p_f_ghost_in[i + 2]
186 - 2.0 * p_f_ghost_in[i + 1]
187 + p_f_ghost_in[i]) / (dz * dz);
188 else {
189 k_ = k[i];
190 if (i == 0)
191 k_zn = k_;
192 else
193 k_zn = k[i - 1];
194 if (i == nz - 1)
195 k_zp = k_;
196 else
197 k_zp = k[i + 1];
198
199 rhs_term = dt / ((alpha + beta_f * phi[i]) * mu_f)
200 * ((2.0 * k_zp * k_ / (k_zp + k_) / (dz * dz))
201 + (2.0 * k_zn * k_ / (k_zn + k_) / (dz * dz)));
202
203 p_f_ghost_out[i + 1] = 1.0 / (1.0 + rhs_term)
204 * (p_f_old_val[i + 1] + dt
205 * (1.0 / ((alpha + beta_f * phi[i]) * mu_f)
206 * (2.0 * k_zp * k_ / (k_zp + k_)
207 * (p_f_ghost_in[i + 2]) / dz
208 - 2.0 * k_zn * k_ / (k_zn + k_)
209 * -p_f_ghost_in[i] / dz) / dz
210 - 1.0 / ((alpha + beta_f * phi[i])
211 * (1.0 - phi[i])) * phi_dot[i]));
212 p_f_ghost_out[i + 1] = omega * p_f_ghost_out[i + 1]
213 + (1.0 - omega) * p_f_ghost_in[i + 1];
214
215 div_k_grad_p = (2.0 * k_zp * k_ / (k_zp + k_)
216 * (p_f_ghost_out[i + 2] - p_f_ghost_out[i + 1]) / dz
217 - 2.0 * k_zn * k_ / (k_zn + k_)
218 * (p_f_ghost_out[i + 1] - p_f_ghost_out[i]) / dz
219 ) / dz;
220 #ifdef DEBUG
221 printf("%s [%d]: phi=%g\tdiv_k_grad_p=%g\tphi_dot=%g\n",
222 __func__, i, phi[i], div_k_grad_p, phi_dot[i]);
223
224 printf(" p[%d, %d, %d] = [%g, %g, %g]\tk=[%g, %g, %g]\n",
225 i, i + 1, i + 2,
226 p_f_ghost_in[i], p_f_ghost_in[i + 1], p_f_ghost_in[i + 2],
227 k_zn, k_, k_zp);
228 #endif
229 /* use the values from the next time step as the time derivative for this iteration */
230 return 1.0 / ((alpha + beta_f * phi[i]) * mu_f) * div_k_grad_p
231 - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[i])) * phi_dot[i];
232 }
233 }
234
235 int
236 darcy_solver_1d(struct simulation *sim,
237 const int max_iter,
238 const double rel_tol)
239 {
240 int i, iter, solved = 0;
241 double epsilon, p_f_top, omega, r_norm_max = NAN, *k_n, *phi_n;
242
243 copy_values(sim->p_f_dot, sim->p_f_dot_old, sim->nz);
244
245 /* choose integration method, parameter in [0.0; 1.0]
246 * epsilon = 0.0: explicit
247 * epsilon = 0.5: Crank-Nicolson
248 * epsilon = 1.0: implicit */
249 epsilon = 0.5;
250
251 /* underrelaxation parameter in ]0.0; 1.0],
252 * where omega = 1.0: no underrelaxation */
253 omega = 1.0;
254
255 for (i = 0; i < sim->nz; ++i)
256 sim->p_f_dot_impl[i] = sim->p_f_dot_expl[i] = 0.0;
257
258 if (isnan(sim->p_f_mod_pulse_time))
259 p_f_top = sim->p_f_top
260 + sine_wave(sim->t,
261 sim->p_f_mod_ampl,
262 sim->p_f_mod_freq,
263 sim->p_f_mod_phase);
264 else if (sim->p_f_mod_pulse_shape == 1)
265 p_f_top = sim->p_f_top
266 + square_pulse(sim->t,
267 sim->p_f_mod_ampl,
268 sim->p_f_mod_freq,
269 sim->p_f_mod_pulse_time);
270 else
271 p_f_top = sim->p_f_top
272 + triangular_pulse(sim->t,
273 sim->p_f_mod_ampl,
274 sim->p_f_mod_freq,
275 sim->p_f_mod_pulse_time);
276
277 /* set fluid BCs (1 of 2) */
278 set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
279 set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top);
280
281 /* explicit solution to pressure change */
282 if (epsilon < 1.0) {
283 #ifdef DEBUG
284 printf("\nEXPLICIT SOLVER IN %s\n", __func__);
285 #endif
286 for (i = 0; i < sim->nz - 1; ++i)
287 sim->p_f_dot_expl[i] = darcy_pressure_change_1d(i,
288 sim->nz,
289 sim->p_f_ghost,
290 sim->phi,
291 sim->phi_dot,
292 sim->k,
293 sim->dz,
294 sim->beta_f,
295 sim->alpha,
296 sim->mu_f,
297 sim->D);
298 }
299
300 /* implicit solution with Jacobian iterations */
301 if (epsilon > 0.0) {
302
303 #ifdef DEBUG
304 printf("\nIMPLICIT SOLVER IN %s\n", __func__);
305 #endif
306
307 k_n = zeros(sim->nz);
308 phi_n = zeros(sim->nz);
309 if (sim->transient)
310 for (i = 0; i < sim->nz; ++i) {
311 /* grabbing the n + 1 iteration values for k and phi */
312 phi_n[i] = sim->phi[i] + sim->dt * sim->phi_dot[i];
313 k_n[i] = kozeny_carman(sim->d, phi_n[i]);
314 }
315 else
316 for (i = 0; i < sim->nz; ++i) {
317 phi_n[i] = sim->phi[i];
318 k_n[i] = sim->k[i];
319 }
320 copy_values(sim->p_f_next_ghost, sim->tmp_ghost, sim->nz + 2);
321
322 for (iter = 0; iter < max_iter; ++iter) {
323 copy_values(sim->p_f_dot_impl, sim->fluid_old_val, sim->nz);
324
325 #ifdef DEBUG
326 puts(".. p_f_ghost bfore BC:");
327 print_array(sim->tmp_ghost, sim->nz + 2);
328 #endif
329
330 /* set fluid BCs (2 of 2) */
331 set_fluid_bcs(sim->tmp_ghost, sim, p_f_top);
332
333 #ifdef DEBUG
334 puts(".. p_f_ghost after BC:");
335 print_array(sim->tmp_ghost, sim->nz + 2);
336 #endif
337
338 for (i = 0; i < sim->nz - 1; ++i)
339 sim->p_f_dot_impl[i] = darcy_pressure_change_1d_impl(i,
340 sim->nz,
341 sim->dt,
342 sim->p_f_ghost,
343 sim->tmp_ghost,
344 sim->p_f_next_ghost,
345 phi_n,
346 sim->phi_dot,
347 k_n,
348 sim->dz,
349 sim->beta_f,
350 sim->alpha,
351 sim->mu_f,
352 sim->D);
353
354 for (i = 0; i < sim->nz; ++i)
355 if (isnan(sim->p_f_dot_impl[i]))
356 errx(1, "NaN at sim->p_f_dot_impl[%d] (t = %g s, iter = %d)",
357 i, sim->t, iter);
358
359 set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top);
360 for (i = 0; i < sim->nz-1; ++i)
361 sim->p_f_dot_impl_r_norm[i] = fabs(residual(sim->p_f_next_ghost[i],
362 sim->tmp_ghost[i]));
363 r_norm_max = max(sim->p_f_dot_impl_r_norm, sim->nz - 1);
364 copy_values(sim->p_f_next_ghost, sim->tmp_ghost, sim->nz + 2);
365
366 #ifdef DEBUG
367 puts(".. p_f_ghost_new:");
368 print_array(sim->tmp_ghost, sim->nz + 2);
369 #endif
370
371 if (r_norm_max <= rel_tol) {
372 #ifdef DEBUG
373 printf(".. Iterative solution converged after %d iterations\n", iter);
374 #endif
375 solved = 1;
376 break;
377 }
378 }
379 free(k_n);
380 free(phi_n);
381 if (!solved) {
382 fprintf(stderr, "darcy_solver_1d: ");
383 fprintf(stderr, "Solution did not converge after %d iterations\n",
384 iter);
385 fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max);
386 }
387 } else
388 solved = 1;
389
390 for (i = 0; i < sim->nz; ++i)
391 sim->p_f_dot[i] = epsilon * sim->p_f_dot_impl[i]
392 + (1.0 - epsilon) * sim->p_f_dot_expl[i];
393
394 for (i = 0; i < sim->nz; ++i)
395 sim->p_f_dot[i] = omega * sim->p_f_dot[i]
396 + (1.0 - omega) * sim->p_f_dot_old[i];
397
398 for (i = 0; i < sim->nz-1; ++i)
399 sim->p_f_next_ghost[i + 1] = sim->p_f_dot[i] * sim->dt + sim->p_f_ghost[i + 1];
400
401 set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
402 set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top);
403 #ifdef DEBUG
404 printf(".. epsilon = %g\n", epsilon);
405 puts(".. p_f_dot_expl:");
406 print_array(sim->p_f_dot_expl, sim->nz);
407 puts(".. p_f_dot_impl:");
408 print_array(sim->p_f_dot_impl, sim->nz);
409 #endif
410
411 for (i = 0; i < sim->nz; ++i)
412 if (isnan(sim->p_f_dot_expl[i]) || isinf(sim->p_f_dot_expl[i]))
413 errx(1, "invalid: sim->p_f_dot_expl[%d] = %g (t = %g s)",
414 i, sim->p_f_dot_expl[i], sim->t);
415
416 for (i = 0; i < sim->nz; ++i)
417 if (isnan(sim->p_f_dot_impl[i]) || isinf(sim->p_f_dot_impl[i]))
418 errx(1, "invalid: sim->p_f_dot_impl[%d] = %g (t = %g s)",
419 i, sim->p_f_dot_impl[i], sim->t);
420
421 return solved - 1;
422 }