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

max_depth_simple_shear.c (5777B)


      1 #include <unistd.h>
      2 #include <stdio.h>
      3 #include <stdlib.h>
      4 #include <math.h>
      5 #include <getopt.h>
      6 #include <time.h>
      7 
      8 #include "simulation.h"
      9 #include "arg.h"
     10 
     11 #define EPS 1e-12
     12 
     13 /* depth accuracy criteria in meter for solver */
     14 #define TOL 1e-10
     15 #define MAX_ITER 100
     16 
     17 /* uncomment to print time spent per time step to stdout */
     18 /* #define BENCHMARK_PERFORMANCE */
     19 
     20 #ifndef M_PI
     21 #define M_PI 3.1415926535897932384626433832795028841971693993751058209749445923078164062
     22 #endif
     23 
     24 char *argv0;
     25 
     26 static void
     27 usage(void)
     28 {
     29 	fprintf(stderr, "usage: %s "
     30 		"[-a fluid-pressure-ampl] "
     31 		"[-C fluid-compressibility] "
     32 		"[-D fluid-diffusivity] "
     33 		"[-g gravity-accel] "
     34 		"[-h] "
     35 		"[-i fluid-viscosity] "
     36 		"[-k fluid-permeability] "
     37 		"[-O fluid-pressure-top] "
     38 		"[-P grain-compressibility] "
     39 		"[-p grain-porosity] "
     40 		"[-q fluid-pressure-freq] "
     41 		"[-R fluid-density] "
     42 		"[-r grain-density] "
     43 		"[-v] "
     44 		"\n", argv0);
     45 	exit(1);
     46 }
     47 
     48 static double
     49 skin_depth(const struct simulation *sim)
     50 {
     51 	if (sim->D > 0.0)
     52 		return sqrt(sim->D / (M_PI * sim->p_f_mod_freq));
     53 	else
     54 		return sqrt(sim->k[0] / (sim->mu_f * (sim->alpha + sim->phi[0] * sim->beta_f) * M_PI * sim->p_f_mod_freq));
     55 }
     56 
     57 /* using alternate form: sin(x) + cos(x) = sqrt(2)*sin(x + pi/4) */
     58 static double
     59 eff_normal_stress_gradient(struct simulation * sim, double d_s, double z_)
     60 {
     61 	if (z_ / d_s > 10.0) {
     62 		fprintf(stderr, "error: %s: unrealistic depth: %g m "
     63 		        "relative to skin depth %g m\n", __func__, z_, d_s);
     64 		free_arrays(sim);
     65 		exit(10);
     66 	}
     67 	return sqrt(2.0) * sin((3.0 * M_PI / 2.0 - z_ / d_s) + M_PI / 4.0)
     68 		+ (sim->rho_s - sim->rho_f) * sim->G * d_s 
     69 		/ sim->p_f_mod_ampl * exp(z_ / d_s);
     70 }
     71 
     72 /* use Brent's method for finding roots (Press et al., 2007) */
     73 static double
     74 zbrent(struct simulation *sim,
     75        double d_s,
     76        double (*f) (struct simulation *sim, double, double),
     77        double x1,
     78        double x2,
     79        double tol)
     80 {
     81 	int iter;
     82 	double a, b, c, d, e, min1, min2, fa, fb, fc, p, q, r, s, tol1, xm;
     83 
     84 	a = x1;
     85 	b = x2;
     86 	c = x2;
     87 	fa = (*f) (sim, d_s, a);
     88 	fb = (*f) (sim, d_s, b);
     89 
     90 	if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
     91 		fprintf(stderr,
     92 			"error: %s: no root in range %g,%g m\n",
     93 			__func__, x1, x2);
     94 		free_arrays(sim);
     95 		exit(11);
     96 	}
     97 	fc = fb;
     98 	for (iter = 0; iter < MAX_ITER; iter++) {
     99 		if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
    100 			c = a;
    101 			fc = fa;
    102 			d = b - a;
    103 			e = d;
    104 		}
    105 		if (fabs(fc) < fabs(fb)) {
    106 			a = b;
    107 			b = c;
    108 			c = a;
    109 			fa = fb;
    110 			fb = fc;
    111 			fc = fa;
    112 		}
    113 		tol1 = 2.0 * EPS * fabs(b) + 0.5 * tol;
    114 		xm = 0.5 * (c - b);
    115 		if (fabs(xm) <= tol1 || fb == 0.0)
    116 			return b;
    117 		if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
    118 			s = fb / fa;
    119 			if (a == c) {
    120 				p = 2.0 * xm * s;
    121 				q = 1.0 - s;
    122 			} else {
    123 				q = fa / fc;
    124 				r = fb / fc;
    125 				p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0));
    126 				q = (q - 1.0) * (r - 1.0) * (s - 1.0);
    127 			}
    128 			if (p > 0.0)
    129 				q = -q;
    130 			p = fabs(p);
    131 			min1 = 3.0 * xm * q - fabs(tol1 * q);
    132 			min2 = fabs(e * q);
    133 			if (2.0 * p < (min1 < min2 ? min1 : min2)) {
    134 				e = d;
    135 				d = p / q;
    136 			} else {
    137 				d = xm;
    138 				e = d;
    139 			}
    140 		} else {
    141 			d = xm;
    142 			e = d;
    143 		}
    144 		a = b;
    145 		fa = fb;
    146 		if (fabs(d) > tol1)
    147 			b += d;
    148 		else
    149 			b += ((xm) >= 0.0 ? fabs(tol1) : -fabs(tol1));
    150 		fb = (*f) (sim, d_s, b);
    151 	}
    152 	fprintf(stderr,
    153 		"error: %s: exceeded maximum number of iterations",
    154 		__func__);
    155 	free_arrays(sim);
    156 	exit(12);
    157 	return NAN;
    158 }
    159 
    160 int
    161 main(int argc, char *argv[])
    162 {
    163 	int i;
    164 	double new_phi, new_k;
    165 	double d_s, depth, depth_limit1, depth_limit2;
    166 	struct simulation sim;
    167 
    168 #ifdef BENCHMARK_PERFORMANCE
    169 	clock_t t_begin, t_end;
    170 	double t_elapsed;
    171 
    172 #endif
    173 
    174 #ifdef __OpenBSD__
    175 	if (pledge("stdio", NULL) == -1) {
    176 		fprintf(stderr, "error: pledge failed");
    177 		exit(2);
    178 	}
    179 #endif
    180 
    181 	init_sim(&sim);
    182 	new_phi = sim.phi[0];
    183 	new_k = sim.k[0];
    184 
    185 	ARGBEGIN {
    186 	case 'a':
    187 		sim.p_f_mod_ampl = atof(EARGF(usage()));
    188 		break;
    189 	case 'C':
    190 		sim.beta_f = atof(EARGF(usage()));
    191 		break;
    192 	case 'D':
    193 		sim.D = atof(EARGF(usage()));
    194 		break;
    195 	case 'g':
    196 		sim.G = atof(EARGF(usage()));
    197 		break;
    198 	case 'h':
    199 		usage();
    200 		break;
    201 	case 'i':
    202 		sim.mu_f = atof(EARGF(usage()));
    203 		break;
    204 	case 'k':
    205 		new_k = atof(EARGF(usage()));
    206 		break;
    207 	case 'O':
    208 		sim.p_f_top = atof(EARGF(usage()));
    209 		break;
    210 	case 'P':
    211 		sim.alpha = atof(EARGF(usage()));
    212 		break;
    213 	case 'p':
    214 		new_phi = atof(EARGF(usage()));
    215 		break;
    216 	case 'q':
    217 		sim.p_f_mod_freq = atof(EARGF(usage()));
    218 		break;
    219 	case 'R':
    220 		sim.rho_f = atof(EARGF(usage()));
    221 		break;
    222 	case 'r':
    223 		sim.rho_s = atof(EARGF(usage()));
    224 		break;
    225 	case 'v':
    226 		printf("%s-" VERSION "\n", argv0);
    227 		return 0;
    228 		break;
    229 	default:
    230 		usage();
    231 	} ARGEND;
    232 
    233 	if (argc > 0)
    234 		usage();
    235 
    236 	sim.nz = 2;
    237 	prepare_arrays(&sim);
    238 
    239 	if (!isnan(new_phi))
    240 		for (i = 0; i < sim.nz; ++i)
    241 			sim.phi[i] = new_phi;
    242 	if (!isnan(new_k))
    243 		for (i = 0; i < sim.nz; ++i)
    244 			sim.k[i] = new_k;
    245 
    246 	check_simulation_parameters(&sim);
    247 
    248 	depth = 0.0;
    249 	d_s = skin_depth(&sim);
    250 
    251 #ifdef BENCHMARK_PERFORMANCE
    252 	t_begin = clock();
    253 #endif
    254 
    255 	/* deep deformatin does not occur with approximately zero
    256 	 * amplitude in water pressure forcing, or if the stress gradient
    257 	 * is positive at zero depth */
    258 	if (fabs(sim.p_f_mod_ampl) > 1e-6 &&
    259 	    eff_normal_stress_gradient(&sim, d_s, 0.0) < 0.0) {
    260 
    261 		depth_limit1 = 0.0;
    262 		depth_limit2 = 5.0 * d_s;
    263 
    264 		depth = zbrent(&sim,
    265 		               d_s,
    266 		               (double (*) (struct simulation *, double, double))
    267 		               eff_normal_stress_gradient,
    268 		               depth_limit1,
    269 		               depth_limit2,
    270 		               TOL);
    271 	}
    272 #ifdef BENCHMARK_PERFORMANCE
    273 	t_end = clock();
    274 	t_elapsed = (double) (t_end - t_begin) / CLOCKS_PER_SEC;
    275 	printf("time spent = %g s\n", t_elapsed);
    276 #endif
    277 
    278 	printf("%.17g\t%.17g\n", depth, d_s);
    279 
    280 	free_arrays(&sim);
    281 
    282 	return 0;
    283 }