fluid.c (7638B)
1 #include "fluid.h" 2 #include "arrays.h" 3 #include "simulation.h" 4 #include <math.h> 5 #include <stdlib.h> 6 7 void hydrostatic_fluid_pressure_distribution(struct simulation *sim) { 8 int i; 9 10 for (i = 0; i < sim->nz; ++i) 11 sim->p_f_ghost[i + 1] = 12 sim->p_f_top + sim->phi[i] * sim->rho_f * sim->G * 13 (sim->origo_z + 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 (uses top 86 porosity phi[nz-1] at the ice-till interface) */ 87 set_bc_dirichlet(p_f_ghost, sim->nz, +1, 88 p_f_top - 89 sim->phi[sim->nz - 1] * sim->rho_f * sim->G * sim->dz); 90 p_f_ghost[sim->nz] = p_f_top; /* pin top physical node */ 91 set_bc_neumann(p_f_ghost, sim->nz, -1, sim->phi[0] * sim->rho_f * sim->G, 92 sim->dz); 93 } 94 95 int darcy_solver_1d(struct simulation *sim) { 96 int i, solved = 0; 97 double p_f_top; 98 double *a = sim->tdma_a; 99 double *b = sim->tdma_b; 100 double *c = sim->tdma_c; 101 double *d = sim->tdma_d; 102 double *x = sim->tdma_x; 103 double *k_n = sim->darcy_k_n; 104 double *phi_n = sim->darcy_phi_n; 105 106 for (i = 0; i < sim->nz; ++i) 107 sim->p_f_dot_impl[i] = 0.0; 108 109 if (isnan(sim->p_f_mod_pulse_time)) 110 p_f_top = sim->p_f_top + sine_wave(sim->t, sim->p_f_mod_ampl, 111 sim->p_f_mod_freq, sim->p_f_mod_phase); 112 else if (sim->p_f_mod_pulse_shape == 1) 113 p_f_top = 114 sim->p_f_top + square_pulse(sim->t, sim->p_f_mod_ampl, 115 sim->p_f_mod_freq, sim->p_f_mod_pulse_time); 116 else 117 p_f_top = sim->p_f_top + triangular_pulse(sim->t, sim->p_f_mod_ampl, 118 sim->p_f_mod_freq, 119 sim->p_f_mod_pulse_time); 120 121 /* set fluid BCs (1 of 2) */ 122 set_fluid_bcs(sim->p_f_ghost, sim, p_f_top); 123 set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top); 124 125 /* implicit solution with TDMA */ 126 { 127 #ifdef DEBUG 128 printf("\nIMPLICIT SOLVER IN %s\n", __func__); 129 #endif 130 /* Predictor step for nonlinear coefficients */ 131 if (sim->transient) 132 for (i = 0; i < sim->nz; ++i) { 133 phi_n[i] = sim->phi[i] + sim->dt * sim->phi_dot[i]; 134 k_n[i] = kozeny_carman(sim->d, phi_n[i]); 135 } 136 else 137 for (i = 0; i < sim->nz; ++i) { 138 phi_n[i] = sim->phi[i]; 139 k_n[i] = sim->k[i]; 140 } 141 142 /* Build Tridiagonal System */ 143 for (i = 0; i < sim->nz; ++i) { 144 if (sim->D > 0.0) { 145 /* Constant diffusivity mode */ 146 double coeff = sim->D * sim->dt / (sim->dz * sim->dz); 147 a[i] = -coeff; 148 b[i] = 1.0 + 2.0 * coeff; 149 c[i] = -coeff; 150 /* RHS is just p_old (plus potential source terms if applicable, 151 but explicit function for D > 0 uses only diffusion term) */ 152 d[i] = sim->p_f_ghost[i + 1]; 153 } else { 154 /* Coefficients calculation matches the discretized equation */ 155 double k_i = k_n[i]; 156 double k_zn, k_zp; 157 158 if (i == 0) 159 k_zn = k_i; 160 else 161 k_zn = k_n[i - 1]; 162 163 if (i == sim->nz - 1) 164 k_zp = k_i; 165 else 166 k_zp = k_n[i + 1]; 167 168 /* Harmonic means for conductivity at faces */ 169 double k_harm_p = 2.0 * k_zp * k_i / fmax(k_zp + k_i, 1e-30); 170 double k_harm_n = 2.0 * k_zn * k_i / fmax(k_zn + k_i, 1e-30); 171 172 /* Diffusion terms */ 173 double coupling = 174 1.0 / ((sim->alpha + sim->beta_f * phi_n[i]) * sim->mu_f); 175 double porosity_term = 176 -1.0 / ((sim->alpha + sim->beta_f * phi_n[i]) * (1.0 - phi_n[i])); 177 178 /* Matrix coefficients (LHS) */ 179 /* term: - dt * coupling * (k_n * (p_i - p_{i-1}) / dz^2) */ 180 double alpha_i = -sim->dt * coupling * k_harm_n / 181 (sim->dz * sim->dz); // coeff for p_{i-1} 182 double gamma_i = -sim->dt * coupling * k_harm_p / 183 (sim->dz * sim->dz); // coeff for p_{i+1} 184 double beta_i = 185 1.0 - 186 (alpha_i + gamma_i); // coeff for p_i (sum of abs(off-diags) + 1) 187 188 a[i] = alpha_i; 189 b[i] = beta_i; 190 c[i] = gamma_i; 191 192 /* RHS: p_old + source term */ 193 d[i] = 194 sim->p_f_ghost[i + 1] + sim->dt * porosity_term * sim->phi_dot[i]; 195 } 196 } 197 198 /* Apply Boundary Conditions to Linear System */ 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 at the physical top node (ice-till interface): 207 p[nz-1] = p_f_top(t). Identity row so the solved top-node pressure 208 equals the imposed value exactly. */ 209 a[sim->nz - 1] = 0.0; 210 b[sim->nz - 1] = 1.0; 211 c[sim->nz - 1] = 0.0; 212 d[sim->nz - 1] = p_f_top; 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_impl:"); 237 print_array(sim->p_f_dot_impl, sim->nz); 238 #endif 239 240 for (i = 0; i < sim->nz; ++i) 241 if (isnan(sim->p_f_dot_impl[i]) || isinf(sim->p_f_dot_impl[i])) { 242 fprintf(stderr, "invalid: sim->p_f_dot_impl[%d] = %g (t = %g s)\n", i, 243 sim->p_f_dot_impl[i], sim->t); 244 return 1; 245 } 246 247 return solved - 1; 248 }