fluid.c (9563B)
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 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 /* Determines the largest time step for the current simulation state. Note 19 * that the time step should be recalculated if cell sizes or diffusivities 20 * (i.e., permeabilities, porosities, viscosities, or compressibilities) 21 * change. The safety factor should be in ]0;1] */ 22 void 23 set_largest_fluid_timestep(struct simulation* sim, const double safety) 24 { 25 int i; 26 double dx_min, diff, diff_max; 27 double *dx; 28 29 dx = spacing(sim->z, sim->nz); 30 dx_min = INFINITY; 31 for (i=0; i<sim->nz-1; ++i) { 32 if (dx[i] < 0.0) { 33 fprintf(stderr, "error: cell spacing negative (%g) in cell %d\n", 34 dx[i], i); 35 free(dx); 36 exit(10); 37 } 38 if (dx[i] < dx_min) dx_min = dx[i]; 39 } 40 free(dx); 41 42 diff_max = -INFINITY; 43 for (i=0; i<sim->nz; ++i) { 44 diff = sim->k[i]/(sim->phi[i]*sim->beta_f*sim->mu_f); 45 if (diff > diff_max) diff_max = diff; 46 } 47 48 sim->dt = safety*0.5*dx_min*dx_min/diff_max; 49 if (sim->file_dt*0.5 < sim->dt) 50 sim->dt = sim->file_dt; 51 } 52 53 static double 54 sine_wave(const double time, 55 const double amplitude, 56 const double frequency, 57 const double phase) 58 { 59 return amplitude*sin(2.0*PI*frequency*time + phase); 60 } 61 62 static double 63 triangular_pulse(const double time, 64 const double peak_amplitude, 65 const double frequency, 66 const double peak_time) 67 { 68 if (peak_time - 1.0/(2.0*frequency) < time && time <= peak_time) 69 return peak_amplitude*2.0*frequency*(time - peak_time) + peak_amplitude; 70 else if (peak_time < time && time < peak_time + 1.0/(2.0*frequency)) 71 return peak_amplitude*2.0*frequency*(peak_time - time) + peak_amplitude; 72 else 73 return 0.0; 74 } 75 76 static double 77 square_pulse(const double time, 78 const double peak_amplitude, 79 const double frequency, 80 const double peak_time) 81 { 82 if (peak_time - 1.0/(2.0*frequency) < time && 83 time < peak_time + 1.0/(2.0*frequency)) 84 return peak_amplitude; 85 else 86 return 0.0; 87 } 88 89 static void 90 set_fluid_bcs(struct simulation* sim, const double p_f_top) 91 { 92 set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); 93 sim->p_f_ghost[sim->nz] = p_f_top; /* Include top node in BC */ 94 set_bc_neumann(sim->p_f_ghost, 95 sim->nz, 96 -1, 97 sim->phi[0]*sim->rho_f*sim->G, 98 sim->dz); 99 } 100 101 static double 102 darcy_pressure_change_1d(const int i, 103 const int nz, 104 const double* p_f_ghost_in, 105 const double* phi, 106 const double* k, 107 const double dz, 108 const double beta_f, 109 const double mu_f) 110 { 111 double k_, div_k_grad_p, k_zn, k_zp; 112 113 k_ = k[i]; 114 if (i==0) k_zn = k_; else k_zn = k[i-1]; 115 if (i==nz-1) k_zp = k_; else k_zp = k[i+1]; 116 #ifdef DEBUG 117 printf("%d->%d: p=[%g, %g, %g]\tk=[%g, %g, %g]\n", 118 i, i+1, 119 p_zn, p, p_zp, 120 k_zn, k_, k_zp); 121 #endif 122 123 div_k_grad_p = (2.0*k_zp*k_/(k_zp + k_) 124 * (p_f_ghost_in[i+2] 125 - p_f_ghost_in[i+1])/dz 126 - 2.0*k_zn*k_/(k_zn + k_) 127 * (p_f_ghost_in[i+1] 128 - p_f_ghost_in[i])/dz 129 )/dz; 130 131 #ifdef DEBUG 132 printf("phi[%d]=%g\tdiv_k_grad_p[%d]=%g\n", 133 i, phi[i], i, div_k_grad_p); 134 #endif 135 136 return 1.0/(beta_f*phi[i]*mu_f)*div_k_grad_p; 137 } 138 139 /*void 140 darcy_pressure_change_1d(double* dp_f_dt, 141 const int nz, 142 const double* p_f_ghost_in, 143 const double* phi, 144 const double* k, 145 const double dz, 146 const double beta_f, 147 const double mu_f) 148 { 149 int i; 150 double k_, div_k_grad_p, k_zn, k_zp; 151 152 #pragma omp parallel for private(i, k_, div_k_grad_p, k_zn, k_zp) if(nz>12) 153 for (i=0; i<nz; ++i) { 154 155 k_ = k[i]; 156 if (i==0) k_zn = k_; else k_zn = k[i-1]; 157 if (i==nz-1) k_zp = k_; else k_zp = k[i+1]; 158 159 div_k_grad_p = (2.0*k_zp*k_/(k_zp + k_) 160 * (p_f_ghost_in[i+2] 161 - p_f_ghost_in[i+1])/dz 162 - 2.0*k_zn*k_/(k_zn + k_) 163 * (p_f_ghost_in[i+1] 164 - p_f_ghost_in[i])/dz 165 )/dz; 166 167 dp_f_dt[i] = 1.0/(beta_f*phi[i]*mu_f)*div_k_grad_p; 168 } 169 }*/ 170 171 int 172 darcy_solver_1d(struct simulation* sim, 173 const int max_iter, 174 const double rel_tol) 175 { 176 int i, iter, solved; 177 double epsilon, theta, p_f_top, r_norm_max; 178 double *dp_f_dt_expl; 179 double *p_f_ghost_old, *dp_f_dt_impl, *p_f_ghost_new, *r_norm; 180 181 r_norm_max = NAN; 182 183 /* choose integration method, parameter in [0.0; 1.0] 184 * epsilon = 0.0: explicit 185 * epsilon = 0.5: Crank-Nicolson 186 * epsilon = 1.0: implicit */ 187 epsilon = 0.5; 188 189 /* choose relaxation factor, parameter in ]0.0; 1.0] 190 * theta in ]0.0; 1.0]: underrelaxation 191 * theta = 1.0: Gauss-Seidel 192 * theta > 1.0: overrelaxation */ 193 /* TODO: values other than 1.0 do not work! */ 194 theta = 1.0; 195 196 if (isnan(sim->p_f_mod_pulse_time)) 197 p_f_top = sim->p_f_top + sine_wave(sim->t, 198 sim->p_f_mod_ampl, 199 sim->p_f_mod_freq, 200 sim->p_f_mod_phase); 201 else 202 if (sim->p_f_mod_pulse_shape == 1) 203 p_f_top = sim->p_f_top + square_pulse(sim->t, 204 sim->p_f_mod_ampl, 205 sim->p_f_mod_freq, 206 sim->p_f_mod_pulse_time); 207 else 208 p_f_top = sim->p_f_top + triangular_pulse(sim->t, 209 sim->p_f_mod_ampl, 210 sim->p_f_mod_freq, 211 sim->p_f_mod_pulse_time); 212 213 /* set fluid BCs (1 of 2) */ 214 set_fluid_bcs(sim, p_f_top); 215 216 solved = 0; 217 218 if (epsilon < 1.0) { 219 /* compute explicit solution to pressure change */ 220 dp_f_dt_expl = zeros(sim->nz); 221 for (i=0; i<sim->nz; ++i) 222 dp_f_dt_expl[i] = darcy_pressure_change_1d(i, 223 sim->nz, 224 sim->p_f_ghost, 225 sim->phi, 226 sim->k, 227 sim->dz, 228 sim->beta_f, 229 sim->mu_f); 230 } 231 232 if (epsilon > 0.0) { 233 /* compute implicit solution with Jacobian iterations */ 234 dp_f_dt_impl = zeros(sim->nz); 235 p_f_ghost_new = zeros(sim->nz+2); 236 r_norm = zeros(sim->nz); 237 p_f_ghost_old = empty(sim->nz+2); 238 copy_values(sim->p_f_ghost, p_f_ghost_old, sim->nz+2); 239 240 for (iter=0; iter<max_iter; ++iter) { 241 242 /* set fluid BCs (2 of 2) */ 243 set_fluid_bcs(sim, p_f_top); 244 #ifdef DEBUG 245 puts(".. p_f_ghost after BC:"); 246 print_array(sim->p_f_ghost, sim->nz+2); 247 #endif 248 249 for (i=0; i<sim->nz-1; ++i) 250 dp_f_dt_impl[i] = darcy_pressure_change_1d(i, 251 sim->nz, 252 sim->p_f_ghost, 253 sim->phi, 254 sim->k, 255 sim->dz, 256 sim->beta_f, 257 sim->mu_f); 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 265 p_f_ghost_new[i+1] = p_f_ghost_old[i+1] 266 + epsilon*dp_f_dt_impl[i]*sim->dt; 267 268 if (epsilon < 1.0) 269 p_f_ghost_new[i+1] += (1.0 - epsilon) 270 *dp_f_dt_expl[i]*sim->dt; 271 272 p_f_ghost_new[i+1] = p_f_ghost_old[i+1]*(1.0 - theta) 273 + p_f_ghost_new[i+1]*theta; 274 275 r_norm[i] = (p_f_ghost_new[i+1] 276 - sim->p_f_ghost[i+1]) 277 /(sim->p_f_ghost[i+1] + 1e-16); 278 } 279 280 r_norm_max = max(r_norm, sim->nz); 281 #ifdef DEBUG 282 puts(".. p_f_ghost_new:"); 283 print_array(p_f_ghost_new, sim->nz+2); 284 #endif 285 286 copy_values(p_f_ghost_new, sim->p_f_ghost, sim->nz+2); 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 return solved - 1; 326 }