1d_fd_simple_shear

commit 182fc9c9dfc2e04105703ae653fffbc9d3cf09f9
parent 2ac6bc191520c22444978580b39a69749d6bb804
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Thu, 16 Apr 2020 13:05:04 +0200

Move struct simulation out of global scope for helper programs

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

1 file changed, 28 insertions(+), 25 deletions(-)
diff --git a/max_depth_simple_shear.c b/max_depth_simple_shear.c
@@ -23,7 +23,6 @@
 #endif
 
 char *argv0;
-struct simulation sim;
 
 static void
 usage(void)
@@ -45,16 +44,16 @@ usage(void)
 	exit(1);
 }
 
-double
-skin_depth()
+static 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));
+	return sqrt(sim->k[0]/
+	       (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
-eff_normal_stress_gradient(double d_s, double z_)
+static double
+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 "
@@ -63,15 +62,17 @@ eff_normal_stress_gradient(double d_s, double z_)
 	}
 
 	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);
+	       + (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(double d_s,
-              double (*f)(double, double),
-              double x1,
-              double x2,
-              double tol)
+static double
+zbrent(struct simulation *sim,
+       double d_s,
+       double (*f)(const 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;
@@ -79,8 +80,8 @@ double zbrent(double d_s,
 	a = x1;
 	b = x2;
 	c = x2;
-	fa = (*f)(d_s, a);
-	fb = (*f)(d_s, b);
+	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,
@@ -141,7 +142,7 @@ double zbrent(double d_s,
 			b += d;
 		else
 			b += ((xm) >= 0.0 ? fabs(tol1) : -fabs(tol1));
-		fb = (*f)(d_s, b);
+		fb = (*f)(sim, d_s, b);
 	}
 	fprintf(stderr,
 			"error: %s: exceeded maximum number of iterations",
@@ -156,6 +157,7 @@ main(int argc, char* argv[])
 	int i;
 	double new_phi, new_k;
 	double d_s, depth, depth_limit1, depth_limit2;
+	struct simulation sim;
 #ifdef BENCHMARK_PERFORMANCE
 	clock_t t_begin, t_end;
 	double t_elapsed;
@@ -167,9 +169,8 @@ main(int argc, char* argv[])
 		exit(1);
 	}
 #endif
-	atexit(free_arrays);
 
-	init_sim();
+	init_sim(&sim);
 	new_phi = sim.phi[0];
 	new_k = sim.k[0];
 
@@ -219,7 +220,7 @@ main(int argc, char* argv[])
 		usage();
 
 	sim.nz = 2;
-	prepare_arrays();
+	prepare_arrays(&sim);
 
 	if (!isnan(new_phi))
 		for (i=0; i<sim.nz; ++i)
@@ -228,10 +229,10 @@ main(int argc, char* argv[])
 		for (i=0; i<sim.nz; ++i)
 			sim.k[i] = new_k;
 
-	check_simulation_parameters();
+	check_simulation_parameters(&sim);
 
 	depth = 0.0;
-	d_s = skin_depth();
+	d_s = skin_depth(&sim);
 
 #ifdef BENCHMARK_PERFORMANCE
 	t_begin = clock();
@@ -241,14 +242,15 @@ main(int argc, char* argv[])
 	 * 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(d_s, 0.0) < 0.0) {
+		eff_normal_stress_gradient(&sim, d_s, 0.0) < 0.0) {
 
 		depth_limit1 = 0.0;
 		depth_limit2 = 5.0*d_s;
 
-		depth = zbrent(d_s,
-                       (double (*)(double, double))
-		                   eff_normal_stress_gradient,
+		depth = zbrent(&sim,
+					   d_s,
+                       (double (*)(const struct simulation*, double, double))
+                       eff_normal_stress_gradient,
                        depth_limit1,
                        depth_limit2,
                        TOL);
@@ -262,5 +264,6 @@ main(int argc, char* argv[])
 
 	printf("%.17g\t%.17g\n", depth, d_s);
 
+	free_arrays(&sim);
 	return 0;
 }