cngf-pf

continuum model for granular flows with pore-pressure dynamics
git clone git://src.adamsgaard.dk/1d_fd_simple_shear
Log | Files | Refs | README | LICENSE Back to index

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 }