1d_fd_simple_shear

commit 057284afb4fe90fd16e123b3415f282b009ee37e
parent 53e7a92719f78be98796bb7a7427803051da8776
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Wed, 11 Dec 2019 16:09:14 +0100

Further modify solver

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

1 file changed, 54 insertions(+), 12 deletions(-)
diff --git a/max_depth_simple_shear.c b/max_depth_simple_shear.c
@@ -13,9 +13,12 @@
 
 #include "parameter_defaults.h"
 
-/* relative tolerance criteria for solver */
-#define TOL 1e-3
-#define MAX_ITER 10000
+/* tolerance criteria in meter for solver */
+#define RTOL 1e-5
+/* #define MAX_ITER 10000 */
+#define MAX_ITER 1000
+
+#define VERBOSE
 
 /* uncomment to print time spent per time step to stdout */
 /*#define BENCHMARK_PERFORMANCE*/
@@ -90,6 +93,12 @@ analytical_solution_lhs(double d_s, double z_)
 double
 analytical_solution_rhs(const struct simulation* sim, double d_s, double z_)
 {
+	if (z_/d_s > 10.0) {
+		fprintf(stderr, "error: %s: unrealistic depth: %g m "
+				"relative to skin depth %g m\n", __func__, z_, d_s);
+		exit(1);
+	}
+
 	return -(sim->rho_s - sim->rho_f)*sim->G*d_s / sim->p_f_mod_ampl
 		   - exp(z_/d_s);
 }
@@ -229,23 +238,56 @@ 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
+			}
 
-			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);
+			res_norm = (depth - depth_prev)/depth;
+			printf("%g\t%g\n", diff, res_norm);
 
-			if (diff > 0.0)
-				depth += 1e-3*d_s/depth;
-			else
-				depth -= 1e-3*d_s/depth;
+#ifdef VERBOSE
+			/* printf("adjusting depth to %g\n", depth); */
+#endif
 
 			iter++;
-		} while (diff > TOL);
+			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);
+			}
+
+#ifdef VERBOSE
+			/* puts(""); */
+#endif
+		} while (fabs(diff) > RTOL);
 	}
 
 #ifdef BENCHMARK_PERFORMANCE