cngf-pf

commit 13337da5aff9712f2b6e986de9100d93ec5372f7
parent 4b66b1f82965a92a91d4cc3a8a93a14da04c576e
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Thu, 12 Dec 2019 12:42:57 +0100

Find roots to equation system using Brent's method

Diffstat:
M
max_depth_simple_shear.c
 | 
189
++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------

1 file changed, 115 insertions(+), 74 deletions(-)
diff --git a/max_depth_simple_shear.c b/max_depth_simple_shear.c
@@ -3,25 +3,25 @@
 #include <stdlib.h>
 #include <math.h>
 #include <getopt.h>
-#include <string.h>
 #include <time.h>
 
 #include "simulation.h"
-#include "fluid.h"
+/* #include "fluid.h" */
 
 #define PROGNAME "max_depth_simple_shear"
 
 #include "parameter_defaults.h"
 
-/* tolerance criteria in meter for solver */
-#define RTOL 1e-5
-/* #define MAX_ITER 10000 */
-#define MAX_ITER 1000
+#define EPS 1e-12
+
+/* depth accuracy criteria in meter for solver */
+#define TOL 1e-10
+#define MAX_ITER 100
 
 #define VERBOSE
 
 /* uncomment to print time spent per time step to stdout */
-/*#define BENCHMARK_PERFORMANCE*/
+/* #define BENCHMARK_PERFORMANCE */
 
 static void
 usage(void)
@@ -80,18 +80,12 @@ double
 skin_depth(const struct simulation* sim)
 {
 	return sqrt(sim->k[0]/
-			    (sim->phi[0]*sim->mu_f*sim->beta_f*M_PI*sim->p_f_mod_freq));
+	       (sim->phi[0]*sim->mu_f*sim->beta_f*M_PI*sim->p_f_mod_freq));
 }
 
 /* using alternate form: sin(x) + cos(x) = sqrt(2)*sin(x + pi/4) */
 double
-analytical_solution_lhs(double d_s, double z_)
-{
-	return sqrt(2.0) * sin((3.0*M_PI/2.0 - z_/d_s) + M_PI/4.0);
-}
-
-double
-analytical_solution_rhs(const struct simulation* sim, double d_s, double z_)
+eff_normal_stress_gradient(const struct simulation* sim, double d_s, double z_)
 {
 	if (z_/d_s > 10.0) {
 		fprintf(stderr, "error: %s: unrealistic depth: %g m "
@@ -99,19 +93,102 @@ analytical_solution_rhs(const struct simulation* sim, double d_s, double z_)
 		exit(1);
 	}
 
-	return -(sim->rho_s - sim->rho_f)*sim->G*d_s / sim->p_f_mod_ampl
-		   - exp(z_/d_s);
+	return sqrt(2.0)*sin((3.0*M_PI/2.0 - z_/d_s) + M_PI/4.0)
+	       + (sim->rho_s - sim->rho_f)*sim->G*d_s/sim->p_f_mod_ampl*exp(z_/d_s);
+}
+
+/* use Brent's method for finding roots (Press et al., 2007) */
+double zbrent(struct simulation* sim,
+              double d_s,
+              double (*f)(struct simulation* sim, double, double),
+              double x1,
+              double x2,
+              double tol)
+
+{
+	int iter;
+	double a, b, c, d, e, min1, min2, fa, fb, fc, p, q, r, s, tol1, xm;
+
+	a = x1; b=x2; c=x2;
+	fa = (*f)(sim, d_s, a);
+	fb = (*f)(sim, d_s, b);
+
+	if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
+		fprintf(stderr,
+				"error: %s: no root of stress gradient in range %g,%g m\n",
+				__func__, x1, x2);
+		exit(1);
+	}
+	fc = fb;
+	for (iter=0; iter<MAX_ITER; iter++) {
+		if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
+			c = a;
+			fc = fa;
+			d = b - a;
+			e = d;
+		}
+		if (fabs(fc) < fabs(fb)) {
+			a = b;
+			b = c;
+			c = a;
+			fa = fb;
+			fb = fc;
+			fc = fa;
+		}
+		tol1 = 2.0*EPS*fabs(b) + 0.5*tol;
+		xm = 0.5*(c - b);
+		if (fabs(xm) <= tol1 || fb == 0.0)
+			return b;
+		if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
+			s = fb/fa;
+			if (a == c) {
+				p = 2.0*xm*s;
+				q = 1.0 - s;
+			} else {
+				q = fa/fc;
+				r = fb/fc;
+				p = s*(2.0*xm*q*(q - r) - (b - a)*(r - 1.0));
+				q = (q - 1.0)*(r - 1.0)*(s - 1.0);
+			}
+			if (p > 0.0)
+				q = -q;
+			p = fabs(p);
+			min1 = 3.0*xm*q - fabs(tol1*q);
+			min2 = fabs(e*q);
+			if (2.0*p < (min1 < min2 ? min1 : min2)) {
+				e = d;
+				d = p/q;
+			} else {
+				d = xm;
+				e = d;
+			}
+		} else {
+			d = xm;
+			e = d;
+		}
+		a = b;
+		fa = fb;
+		if (fabs(d) > tol1)
+			b += d;
+		else
+			b += ((xm) >= 0.0 ? fabs(tol1) : -fabs(tol1));
+		fb = (*f)(sim, d_s, b);
+	}
+	fprintf(stderr,
+			"error: %s: exceeded maximum number of iterations",
+			__func__);
+	exit(10);
+	return NAN;
 }
 
 int
 main(int argc, char* argv[])
 {
-	int norm, opt, i;
+	int opt, i;
 	struct simulation sim;
 	const char* optstring;
-	unsigned long iter;
 	double new_phi, new_k;
-	double diff, d_s, depth;
+	double d_s, depth, depth_limit1, depth_limit2;
 #ifdef BENCHMARK_PERFORMANCE
 	clock_t t_begin, t_end;
 	double t_elapsed;
@@ -127,8 +204,6 @@ main(int argc, char* argv[])
 	/* load with default values */
 	sim = init_sim();
 
-	norm = 0;
-
 	optstring = "hvNn:G:P:m:s:l:V:A:b:f:C:Fp:d:r:o:L:c:i:R:k:O:a:q:H:T:S:t:T:D:I:";
 	const struct option longopts[] = {
 		{"help",                      no_argument,       NULL, 'h'},
@@ -221,16 +296,15 @@ main(int argc, char* argv[])
 		for (i=0; i<sim.nz; ++i)
 			sim.k[i] = new_k;
 
-	lithostatic_pressure_distribution(&sim);
+	/*lithostatic_pressure_distribution(&sim);*/
 
-	if (sim.fluid) {
+	/*if (sim.fluid) {
 		hydrostatic_fluid_pressure_distribution(&sim);
 		set_largest_fluid_timestep(&sim, 0.5);
-	}
+	}*/
 
 	check_simulation_parameters(&sim);
 
-	iter = 0;
 	depth = 0.0;
 	d_s = skin_depth(&sim);
 
@@ -238,56 +312,24 @@ main(int argc, char* argv[])
 	t_begin = clock();
 #endif
 
-	double depth_prev;
-	double res_norm;
-
-	if (fabs(sim.p_f_mod_ampl) > 1e-6) {
-		do {
-			depth_prev = depth;
-			diff = analytical_solution_rhs(&sim, d_s, depth)
-				 - analytical_solution_lhs(d_s, depth);
-			/* printf("%g\n", diff); */
-#ifdef VERBOSE
-			/* printf("rhs = %g\n", analytical_solution_rhs(&sim, d_s, depth)); */
-			/* printf("lhs = %g\n", analytical_solution_lhs(d_s, depth)); */
-			/* printf("%ld\t%g\t%g\n", iter, depth, diff); */
-#endif
-			depth -= diff*d_s;
-
-			/*if (diff > 0.0) {
-				depth += 1e-5*d_s;
-			} else {
-				depth -= 1e-5*d_s;
-			}*/
-
-			if (depth < 0.0) {
-				depth = fabs(depth);
-#ifdef VERBOSE
-				printf("bouncing back from zero depth\n");
-#endif
-			}
+	/* deep deformatin does not occur with approximately zero amplitude in 
+	 * water pressure forcing, or if the stress gradient is positive at
+	 * zero depth */
+	if (fabs(sim.p_f_mod_ampl) > 1e-6 && 
+		eff_normal_stress_gradient(&sim, d_s, 0.0) < 0.0) {
 
-			res_norm = fabs((depth - depth_prev)/(depth_prev + 1e-16));
-			printf("%g\t%g\n", diff, res_norm);
+		depth_limit1 = 0.0;
+		depth_limit2 = 5.0*d_s;
 
-#ifdef VERBOSE
-			/* printf("adjusting depth to %g\n", depth); */
-#endif
-
-			iter++;
-			if (iter > MAX_ITER) {
-				fprintf(stderr, "error: %s: solution did not converge\n"
-						"    iteration %ld\n"
-						"    diff = %g\n"
-						"    depth = %g\n",
-						__func__, iter, diff, depth);
-				exit(1);
-			}
+		depth = zbrent(&sim,
+                       d_s,
+                       (double (*)(struct simulation*, double, double))
+		                   eff_normal_stress_gradient,
+                       depth_limit1,
+                       depth_limit2,
+                       TOL);
 
-#ifdef VERBOSE
-			/* puts(""); */
-#endif
-		} while (fabs(diff) > RTOL);
+		/* printf("%g\t%g\n", diff, res_norm); */
 	}
 
 #ifdef BENCHMARK_PERFORMANCE
@@ -296,7 +338,6 @@ main(int argc, char* argv[])
 	printf("time spent = %g s\n", t_elapsed);
 #endif
 
-	printf("\nresult:\n");
 	printf("%g\t%g\n", depth, d_s);
 
 	free_arrays(&sim);