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