fluid.c (8391B)
1 #include "fluid.h" 2 #include "arrays.h" 3 #include "simulation.h" 4 #include <err.h> 5 #include <math.h> 6 #include <stdlib.h> 7 8 void hydrostatic_fluid_pressure_distribution(struct simulation *sim) { 9 int i; 10 11 for (i = 0; i < sim->nz; ++i) 12 sim->p_f_ghost[i + 1] = sim->p_f_top + sim->phi[i] * sim->rho_f * sim->G * 13 (sim->L_z - sim->z[i]); 14 } 15 16 static double diffusivity(struct simulation *sim, int i) { 17 if (sim->D > 0.0) 18 return sim->D; 19 else 20 return sim->k[i] / ((sim->alpha + sim->phi[i] * sim->beta_f) * sim->mu_f); 21 } 22 23 /* Determines the largest time step for the current simulation state. Note 24 * that the time step should be recalculated if cell sizes or 25 * diffusivities (i.e., permeabilities, porosities, viscosities, or 26 * compressibilities) change. The safety factor should be in ]0;1] */ 27 int set_largest_fluid_timestep(struct simulation *sim, const double safety) { 28 int i; 29 double dx_min, diff, diff_max, *dx; 30 31 dx = spacing(sim->z, sim->nz); 32 dx_min = INFINITY; 33 for (i = 0; i < sim->nz - 1; ++i) { 34 if (dx[i] < 0.0) { 35 fprintf(stderr, "error: cell spacing negative (%g) in cell %d\n", dx[i], 36 i); 37 free(dx); 38 return 1; 39 } 40 if (dx[i] < dx_min) 41 dx_min = dx[i]; 42 } 43 free(dx); 44 45 diff_max = -INFINITY; 46 for (i = 0; i < sim->nz; ++i) { 47 diff = diffusivity(sim, i); 48 if (diff > diff_max) 49 diff_max = diff; 50 } 51 52 sim->dt = safety * 0.5 * dx_min * dx_min / diff_max; 53 if (sim->file_dt < sim->dt) 54 sim->dt = sim->file_dt; 55 56 return 0; 57 } 58 59 static double sine_wave(const double time, const double ampl, const double freq, 60 const double phase) { 61 return ampl * sin(2.0 * PI * freq * time + phase); 62 } 63 64 static double triangular_pulse(const double time, const double peak_ampl, 65 const double freq, const double peak_time) { 66 if (peak_time - 1.0 / (2.0 * freq) < time && time <= peak_time) 67 return peak_ampl * 2.0 * freq * (time - peak_time) + peak_ampl; 68 else if (peak_time < time && time < peak_time + 1.0 / (2.0 * freq)) 69 return peak_ampl * 2.0 * freq * (peak_time - time) + peak_ampl; 70 else 71 return 0.0; 72 } 73 74 static double square_pulse(const double time, const double peak_ampl, 75 const double freq, const double peak_time) { 76 if (peak_time - 1.0 / (2.0 * freq) < time && 77 time < peak_time + 1.0 / (2.0 * freq)) 78 return peak_ampl; 79 else 80 return 0.0; 81 } 82 83 static void set_fluid_bcs(double *p_f_ghost, struct simulation *sim, 84 const double p_f_top) { 85 /* correct ghost node at top BC for hydrostatic pressure gradient */ 86 set_bc_dirichlet(p_f_ghost, sim->nz, +1, 87 p_f_top - sim->phi[0] * sim->rho_f * sim->G * sim->dz); 88 /* p_f_ghost[sim->nz] = p_f_top; */ /* REMOVED: Do not pin top physical node, 89 let solver handle it */ 90 set_bc_neumann(p_f_ghost, sim->nz, -1, sim->phi[0] * sim->rho_f * sim->G, 91 sim->dz); 92 } 93 94 int darcy_solver_1d(struct simulation *sim) { 95 int i, solved = 0; 96 double p_f_top; 97 double *a = sim->tdma_a; 98 double *b = sim->tdma_b; 99 double *c = sim->tdma_c; 100 double *d = sim->tdma_d; 101 double *x = sim->tdma_x; 102 double *k_n = sim->darcy_k_n; 103 double *phi_n = sim->darcy_phi_n; 104 105 for (i = 0; i < sim->nz; ++i) 106 sim->p_f_dot_impl[i] = sim->p_f_dot_expl[i] = 0.0; 107 108 if (isnan(sim->p_f_mod_pulse_time)) 109 p_f_top = sim->p_f_top + sine_wave(sim->t, sim->p_f_mod_ampl, 110 sim->p_f_mod_freq, sim->p_f_mod_phase); 111 else if (sim->p_f_mod_pulse_shape == 1) 112 p_f_top = 113 sim->p_f_top + square_pulse(sim->t, sim->p_f_mod_ampl, 114 sim->p_f_mod_freq, sim->p_f_mod_pulse_time); 115 else 116 p_f_top = sim->p_f_top + triangular_pulse(sim->t, sim->p_f_mod_ampl, 117 sim->p_f_mod_freq, 118 sim->p_f_mod_pulse_time); 119 120 /* set fluid BCs (1 of 2) */ 121 set_fluid_bcs(sim->p_f_ghost, sim, p_f_top); 122 set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top); 123 124 /* implicit solution with TDMA */ 125 { 126 #ifdef DEBUG 127 printf("\nIMPLICIT SOLVER IN %s\n", __func__); 128 #endif 129 /* Predictor step for nonlinear coefficients */ 130 if (sim->transient) 131 for (i = 0; i < sim->nz; ++i) { 132 phi_n[i] = sim->phi[i] + sim->dt * sim->phi_dot[i]; 133 k_n[i] = kozeny_carman(sim->d, phi_n[i]); 134 } 135 else 136 for (i = 0; i < sim->nz; ++i) { 137 phi_n[i] = sim->phi[i]; 138 k_n[i] = sim->k[i]; 139 } 140 141 /* Build Tridiagonal System */ 142 for (i = 0; i < sim->nz; ++i) { 143 if (sim->D > 0.0) { 144 /* Constant diffusivity mode */ 145 double coeff = sim->D * sim->dt / (sim->dz * sim->dz); 146 a[i] = -coeff; 147 b[i] = 1.0 + 2.0 * coeff; 148 c[i] = -coeff; 149 /* RHS is just p_old (plus potential source terms if applicable, 150 but explicit function for D > 0 uses only diffusion term) */ 151 d[i] = sim->p_f_ghost[i + 1]; 152 } else { 153 /* Coefficients calculation matches the discretized equation */ 154 double k_i = k_n[i]; 155 double k_zn, k_zp; 156 157 if (i == 0) 158 k_zn = k_i; 159 else 160 k_zn = k_n[i - 1]; 161 162 if (i == sim->nz - 1) 163 k_zp = k_i; 164 else 165 k_zp = k_n[i + 1]; 166 167 /* Harmonic means for conductivity at faces */ 168 double k_harm_p = 2.0 * k_zp * k_i / fmax(k_zp + k_i, 1e-30); 169 double k_harm_n = 2.0 * k_zn * k_i / fmax(k_zn + k_i, 1e-30); 170 171 /* Diffusion terms */ 172 double coupling = 173 1.0 / ((sim->alpha + sim->beta_f * phi_n[i]) * sim->mu_f); 174 double porosity_term = 175 -1.0 / ((sim->alpha + sim->beta_f * phi_n[i]) * (1.0 - phi_n[i])); 176 177 /* Matrix coefficients (LHS) */ 178 /* term: - dt * coupling * (k_n * (p_i - p_{i-1}) / dz^2) */ 179 double alpha_i = -sim->dt * coupling * k_harm_n / 180 (sim->dz * sim->dz); // coeff for p_{i-1} 181 double gamma_i = -sim->dt * coupling * k_harm_p / 182 (sim->dz * sim->dz); // coeff for p_{i+1} 183 double beta_i = 184 1.0 - 185 (alpha_i + gamma_i); // coeff for p_i (sum of abs(off-diags) + 1) 186 187 a[i] = alpha_i; 188 b[i] = beta_i; 189 c[i] = gamma_i; 190 191 /* RHS: p_old + source term */ 192 d[i] = 193 sim->p_f_ghost[i + 1] + sim->dt * porosity_term * sim->phi_dot[i]; 194 } 195 } 196 197 /* Apply Boundary Conditions to Linear System */ 198 /* Bottom (i=0): Neumann. p_{-1} = p_0 + C. */ 199 /* Bottom (i=0): Neumann. p_{-1} = p_0 + C. */ 200 double bc_neumann_val = sim->phi[0] * sim->rho_f * sim->G * sim->dz; 201 202 b[0] += a[0]; 203 d[0] -= a[0] * bc_neumann_val; 204 a[0] = 0.0; 205 206 /* Top (i=nz-1): Dirichlet. p_{n} = p_ghost_top. 207 Term c[n-1] * p_{n} becomes c[n-1] * p_ghost_top. 208 Subtract from d[n-1]. 209 FIX: Use sim->p_f_ghost[sim->nz + 1] which has correct BC correction. 210 */ 211 d[sim->nz - 1] -= c[sim->nz - 1] * sim->p_f_ghost[sim->nz + 1]; 212 c[sim->nz - 1] = 0.0; 213 214 /* Solve */ 215 tridiagonal_solver(x, a, b, c, d, sim->tdma_c_prime, sim->tdma_d_prime, 216 sim->nz); 217 218 /* Store result in p_f_dot (rate) */ 219 for (i = 0; i < sim->nz; ++i) { 220 sim->p_f_dot[i] = (x[i] - sim->p_f_ghost[i + 1]) / sim->dt; 221 /* Store in impl array too for consistency/debug */ 222 sim->p_f_dot_impl[i] = sim->p_f_dot[i]; 223 } 224 225 add_darcy_iters(1); 226 solved = 1; 227 } 228 229 for (i = 0; i < sim->nz; ++i) 230 sim->p_f_next_ghost[i + 1] = 231 sim->p_f_dot[i] * sim->dt + sim->p_f_ghost[i + 1]; 232 233 set_fluid_bcs(sim->p_f_ghost, sim, p_f_top); 234 set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top); 235 #ifdef DEBUG 236 puts(".. p_f_dot_expl:"); 237 print_array(sim->p_f_dot_expl, sim->nz); 238 puts(".. p_f_dot_impl:"); 239 print_array(sim->p_f_dot_impl, sim->nz); 240 #endif 241 242 for (i = 0; i < sim->nz; ++i) 243 if (isnan(sim->p_f_dot_expl[i]) || isinf(sim->p_f_dot_expl[i])) 244 errx(1, "invalid: sim->p_f_dot_expl[%d] = %g (t = %g s)", i, 245 sim->p_f_dot_expl[i], sim->t); 246 247 for (i = 0; i < sim->nz; ++i) 248 if (isnan(sim->p_f_dot_impl[i]) || isinf(sim->p_f_dot_impl[i])) 249 errx(1, "invalid: sim->p_f_dot_impl[%d] = %g (t = %g s)", i, 250 sim->p_f_dot_impl[i], sim->t); 251 252 return solved - 1; 253 }