1d_fd_simple_shear

commit 53e7a92719f78be98796bb7a7427803051da8776
parent 676dd6649e534885d2f7266b8aca2a6135a6e140
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Wed, 11 Dec 2019 14:38:28 +0100

Begin implementing solver

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

1 file changed, 41 insertions(+), 16 deletions(-)
diff --git a/max_depth_simple_shear.c b/max_depth_simple_shear.c
@@ -14,7 +14,7 @@
 #include "parameter_defaults.h"
 
 /* relative tolerance criteria for solver */
-#define RTOL 1e-5
+#define TOL 1e-3
 #define MAX_ITER 10000
 
 /* uncomment to print time spent per time step to stdout */
@@ -44,7 +44,8 @@ usage(void)
 	        " -q, --fluid-pressure-freq VAL   frequency of sinusoidal pressure variations [s^-1]\n"
 	        "                                 (default %g)\n"
 	        "\n"
-	        "The output value is the max. deformation depth from the top in meters.\n"
+	        "The first column of the output is the max. deformation depth from the top in meters.\n"
+	        "The second column is the skin depth in meters.\n"
 	        ,
 	        __func__, PROGNAME,
 	        sim.G,
@@ -72,7 +73,6 @@ version(void)
 	, PROGNAME, VERSION);
 }
 
-
 double
 skin_depth(const struct simulation* sim)
 {
@@ -80,6 +80,20 @@ skin_depth(const struct simulation* sim)
 			    (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_)
+{
+	return -(sim->rho_s - sim->rho_f)*sim->G*d_s / sim->p_f_mod_ampl
+		   - exp(z_/d_s);
+}
+
 int
 main(int argc, char* argv[])
 {
@@ -88,6 +102,7 @@ main(int argc, char* argv[])
 	const char* optstring;
 	unsigned long iter;
 	double new_phi, new_k;
+	double diff, d_s, depth;
 #ifdef BENCHMARK_PERFORMANCE
 	clock_t t_begin, t_end;
 	double t_elapsed;
@@ -207,30 +222,40 @@ main(int argc, char* argv[])
 	check_simulation_parameters(&sim);
 
 	iter = 0;
-	double diff = INFINITY;
-	while (diff > RTOL) {
+	depth = 0.0;
+	d_s = skin_depth(&sim);
 
 #ifdef BENCHMARK_PERFORMANCE
-		t_begin = clock();
+	t_begin = clock();
 #endif
 
-		diff = 0.0;
+	if (fabs(sim.p_f_mod_ampl) > 1e-6) {
+		do {
+			diff = analytical_solution_rhs(&sim, d_s, depth)
+				 - analytical_solution_lhs(d_s, depth);
 
+			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);
 
+			if (diff > 0.0)
+				depth += 1e-3*d_s/depth;
+			else
+				depth -= 1e-3*d_s/depth;
 
-#ifdef BENCHMARK_PERFORMANCE
-		t_end = clock();
-		t_elapsed = (double)(t_end - t_begin)/CLOCKS_PER_SEC;
-		printf("time spent per time step = %g s\n", t_elapsed);
-#endif
-
-		/* printf("%s\n", max_depth); */
-		iter++;
+			iter++;
+		} while (diff > TOL);
 	}
 
-	printf("skin depth: %g m\n", skin_depth(&sim));
+#ifdef BENCHMARK_PERFORMANCE
+	t_end = clock();
+	t_elapsed = (double)(t_end - t_begin)/CLOCKS_PER_SEC;
+	printf("time spent = %g s\n", t_elapsed);
+#endif
 
+	printf("\nresult:\n");
+	printf("%g\t%g\n", depth, d_s);
 
 	free_arrays(&sim);
 	return 0;