Introduction
Introduction Statistics Contact Development Disclaimer Help
tfluid.c - cngf-pf - continuum model for granular flows with pore-pressure dyna…
git clone git://src.adamsgaard.dk/cngf-pf
Log
Files
Refs
README
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…
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…
25 }
26
27 /* Determines the largest time step for the current simulation state. No…
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 (%…
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) + pea…
82 else if (peak_time < time && time < peak_time + 1.0 / (2.0 * fre…
83 return peak_ampl * 2.0 * freq * (peak_time - time) + pea…
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_…
103 {
104 /* correct ghost node at top BC for hydrostatic pressure gradien…
105 set_bc_dirichlet(p_f_ghost, sim->nz, +1,
106 p_f_top - sim->phi[0] * sim->rho_f * sim->G * s…
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 …
147 - 2.0 * k_zn * k_ / (k_zn + k_)
148 * (p_f_ghost_in[i + 1] - p_f_ghost_in[i]…
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_i…
158 k_zn, k_, k_zp);
159 #endif
160
161 return 1.0 / ((alpha + beta_f * phi[i]) * mu_f) * div_k_…
162 - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[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 * d…
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…
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 …
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…
217 - 2.0 * k_zn * k_ / (k_zn + k_)
218 * (p_f_ghost_out[i + 1] - p_f_ghost_out…
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_i…
227 k_zn, k_, k_zp);
228 #endif
229 /* use the values from the next time step as the time de…
230 return 1.0 / ((alpha + beta_f * phi[i]) * mu_f) * div_k_…
231 - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[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(…
288 …
289 …
290 …
291 …
292 …
293 …
294 …
295 …
296 …
297 …
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 f…
312 phi_n[i] = sim->phi[i] + sim->dt * sim->…
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…
321
322 for (iter = 0; iter < max_iter; ++iter) {
323 copy_values(sim->p_f_dot_impl, sim->fluid_old_va…
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_ch…
340 …
341 …
342 …
343 …
344 …
345 …
346 …
347 …
348 …
349 …
350 …
351 …
352 …
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_imp…
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(resid…
362 …
363 r_norm_max = max(sim->p_f_dot_impl_r_norm, sim->…
364 copy_values(sim->p_f_next_ghost, sim->tmp_ghost,…
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 …
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…
384 iter);
385 fprintf(stderr, ".. Residual normalized error: %…
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[…
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…
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 +…
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_ex…
413 errx(1, "invalid: sim->p_f_dot_expl[%d] = %g (t …
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_im…
418 errx(1, "invalid: sim->p_f_dot_impl[%d] = %g (t …
419 i, sim->p_f_dot_impl[i], sim->t);
420
421 return solved - 1;
422 }
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.