fluid.c (8928B)
1 #include <stdlib.h> 2 #include <math.h> 3 #include "simulation.h" 4 #include "arrays.h" 5 6 /* #define DEBUG true */ 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]/((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; 36 double *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(struct simulation *sim, const double p_f_top) 104 { 105 set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); 106 sim->p_f_ghost[sim->nz] = p_f_top; /* Include top node in BC */ 107 set_bc_neumann(sim->p_f_ghost, 108 sim->nz, 109 -1, 110 sim->phi[0] * sim->rho_f * sim->G, 111 sim->dz); 112 } 113 114 static double 115 darcy_pressure_change_1d(const int i, 116 const int nz, 117 const double *p_f_ghost_in, 118 const double *phi, 119 const double *phi_dot, 120 const double *k, 121 const double dz, 122 const double beta_f, 123 const double alpha, 124 const double mu_f, 125 const double D) 126 { 127 double k_, div_k_grad_p, k_zn, k_zp; 128 129 if (D > 0.0) 130 return D*(p_f_ghost_in[i + 2] - 2.0 * p_f_ghost_in[i + 1] + p_f_ghost_in[i])/(dz*dz); 131 132 else { 133 134 k_ = k[i]; 135 if (i == 0) 136 k_zn = k_; 137 else 138 k_zn = k[i - 1]; 139 if (i == nz - 1) 140 k_zp = k_; 141 else 142 k_zp = k[i + 1]; 143 #ifdef DEBUG 144 printf("%d->%d: p=[%g, %g, %g]\tk=[%g, %g, %g]\n", 145 i, i + 1, 146 p_f_ghost_in[i], p_f_ghost_in[i + 1], p_f_ghost_in[i + 2], 147 k_zn, k_, k_zp); 148 #endif 149 150 div_k_grad_p = (2.0 * k_zp * k_ / (k_zp + k_) 151 * (p_f_ghost_in[i + 2] 152 - p_f_ghost_in[i + 1]) / dz 153 - 2.0 * k_zn * k_ / (k_zn + k_) 154 * (p_f_ghost_in[i + 1] 155 - p_f_ghost_in[i]) / dz 156 ) / dz; 157 158 #ifdef DEBUG 159 printf("darcy_pressure_change [%d]: phi=%g\tdiv_k_grad_p=%g\tphi_dot=%g\n", 160 i, phi[i], div_k_grad_p, phi_dot[i]); 161 #endif 162 163 /* TODO: add advective term */ 164 return 1.0 / ((alpha + beta_f * phi[i]) * mu_f) * div_k_grad_p 165 - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[i])) * phi_dot[i]; 166 } 167 } 168 169 int 170 darcy_solver_1d(struct simulation *sim, 171 const int max_iter, 172 const double rel_tol) 173 { 174 int i, iter, solved = 0; 175 double epsilon, theta, p_f_top, r_norm_max = NAN; 176 double *dp_f_dt_expl, *p_f_ghost_old, *dp_f_dt_impl, *p_f_ghost_new; 177 double *tmp, *r_norm; 178 179 /* choose integration method, parameter in [0.0; 1.0] 180 * epsilon = 0.0: explicit 181 * epsilon = 0.5: Crank-Nicolson 182 * epsilon = 1.0: * implicit */ 183 epsilon = 0.5; 184 185 /* choose relaxation factor, parameter in ]0.0; 1.0] theta in 186 * ]0.0; 1.0]: underrelaxation theta = 1.0: Gauss-Seidel theta > 187 * 1.0: overrelaxation */ 188 /* TODO: values other than 1.0 do not work! */ 189 theta = 1.0; 190 191 if (isnan(sim->p_f_mod_pulse_time)) 192 p_f_top = sim->p_f_top 193 + sine_wave(sim->t, 194 sim->p_f_mod_ampl, 195 sim->p_f_mod_freq, 196 sim->p_f_mod_phase); 197 else if (sim->p_f_mod_pulse_shape == 1) 198 p_f_top = sim->p_f_top 199 + square_pulse(sim->t, 200 sim->p_f_mod_ampl, 201 sim->p_f_mod_freq, 202 sim->p_f_mod_pulse_time); 203 else 204 p_f_top = sim->p_f_top 205 + triangular_pulse(sim->t, 206 sim->p_f_mod_ampl, 207 sim->p_f_mod_freq, 208 sim->p_f_mod_pulse_time); 209 210 /* set fluid BCs (1 of 2) */ 211 set_fluid_bcs(sim, p_f_top); 212 213 if (epsilon < 1.0) { 214 /* compute explicit solution to pressure change */ 215 dp_f_dt_expl = zeros(sim->nz); 216 for (i = 0; i < sim->nz; ++i) 217 dp_f_dt_expl[i] = darcy_pressure_change_1d(i, 218 sim->nz, 219 sim->p_f_ghost, 220 sim->phi, 221 sim->phi_dot, 222 sim->k, 223 sim->dz, 224 sim->beta_f, 225 sim->alpha, 226 sim->mu_f, 227 sim->D); 228 } 229 if (epsilon > 0.0) { 230 /* compute implicit solution with Jacobian iterations */ 231 dp_f_dt_impl = zeros(sim->nz); 232 p_f_ghost_new = zeros(sim->nz + 2); 233 r_norm = zeros(sim->nz); 234 p_f_ghost_old = empty(sim->nz + 2); 235 copy_values(sim->p_f_ghost, p_f_ghost_old, sim->nz + 2); 236 237 for (iter = 0; iter < max_iter; ++iter) { 238 239 /* set fluid BCs (2 of 2) */ 240 set_fluid_bcs(sim, p_f_top); 241 #ifdef DEBUG 242 puts(".. p_f_ghost after BC:"); 243 print_array(sim->p_f_ghost, sim->nz + 2); 244 #endif 245 246 for (i = 0; i < sim->nz - 1; ++i) 247 dp_f_dt_impl[i] = darcy_pressure_change_1d(i, 248 sim->nz, 249 sim->p_f_ghost, 250 sim->phi, 251 sim->phi_dot, 252 sim->k, 253 sim->dz, 254 sim->beta_f, 255 sim->alpha, 256 sim->mu_f, 257 sim->D); 258 259 for (i = 0; i < sim->nz - 1; ++i) { 260 #ifdef DEBUG 261 printf("dp_f_expl[%d] = %g\ndp_f_impl[%d] = %g\n", 262 i, dp_f_dt_expl[i], i, dp_f_dt_impl[i]); 263 #endif 264 p_f_ghost_new[i + 1] = p_f_ghost_old[i + 1] 265 + epsilon * dp_f_dt_impl[i] * sim->dt; 266 267 if (epsilon < 1.0) 268 p_f_ghost_new[i + 1] += (1.0 - epsilon) 269 * dp_f_dt_expl[i] * sim->dt; 270 271 p_f_ghost_new[i + 1] = p_f_ghost_old[i + 1] * (1.0 - theta) 272 + p_f_ghost_new[i + 1] * theta; 273 274 r_norm[i] = fabs(residual(p_f_ghost_new[i + 1], 275 sim->p_f_ghost[i + 1])); 276 } 277 278 r_norm_max = max(r_norm, sim->nz); 279 #ifdef DEBUG 280 puts(".. p_f_ghost_new:"); 281 print_array(p_f_ghost_new, sim->nz + 2); 282 #endif 283 284 tmp = p_f_ghost_new; 285 p_f_ghost_new = sim->p_f_ghost; 286 sim->p_f_ghost = tmp; 287 #ifdef DEBUG 288 puts(".. p_f_ghost after update:"); 289 print_array(sim->p_f_ghost, sim->nz + 2); 290 #endif 291 292 if (r_norm_max <= rel_tol) { 293 #ifdef DEBUG 294 printf(".. Solution converged after %d iterations\n", iter); 295 #endif 296 solved = 1; 297 break; 298 } 299 } 300 free(dp_f_dt_impl); 301 free(p_f_ghost_new); 302 free(r_norm); 303 free(p_f_ghost_old); 304 if (!solved) { 305 fprintf(stderr, "darcy_solver_1d: "); 306 fprintf(stderr, "Solution did not converge after %d iterations\n", 307 iter); 308 fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max); 309 } 310 } else { 311 for (i = 0; i < sim->nz; ++i) 312 sim->p_f_ghost[i + 1] += dp_f_dt_expl[i] * sim->dt; 313 solved = 1; 314 #ifdef DEBUG 315 puts(".. dp_f_dt_expl:"); 316 print_array(dp_f_dt_expl, sim->nz); 317 puts(".. p_f_ghost after explicit solution:"); 318 print_array(sim->p_f_ghost, sim->nz + 2); 319 #endif 320 } 321 set_fluid_bcs(sim, p_f_top); 322 323 if (epsilon < 1.0) 324 free(dp_f_dt_expl); 325 326 return solved - 1; 327 }