fluid.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 }