1d_fd_simple_shear

commit 3c032ad517b519a8b500dc877e3a6ce2067b770b
parent 84ff27ebe390512189687c4fc28adddef22a3d0a
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Wed, 11 Dec 2019 12:46:29 +0100

Add input size checking to array functions

Diffstat:
M
Makefile
 | 
2
+-
M
arrays.c
 | 
28
+++++++++++++++++++++++++++-
M
max_depth_simple_shear.c
 | 
43
++++++++++++++++++++++++++-----------------

3 files changed, 54 insertions(+), 19 deletions(-)
diff --git a/Makefile b/Makefile
@@ -1,4 +1,4 @@
-CFLAGS = -std=c99 -pedantic -Wall -O3
+CFLAGS = -std=c99 -pedantic -Wall -O2 -g
 LDFLAGS = -lm
 HDR = arrays.h fluid.h parameter_defaults.h simulation.h
 BIN = 1d_fd_simple_shear max_depth_simple_shear
diff --git a/arrays.c b/arrays.c
@@ -4,6 +4,16 @@
 
 #define DEG2RAD(x) (x*M_PI/180.0)
 
+void
+check_magnitude(const char* func_name, int limit, int value)
+{
+	if (value < limit) {
+		fprintf(stderr, "error: %s: input size %d is less than %d\n",
+						func_name, value, limit);
+		exit(1);
+	}
+}
+
 /* Translate a i,j,k index in grid with dimensions nx, ny, nz into a linear
  * index */
 unsigned int
@@ -58,6 +68,7 @@ linspace(const double lower, const double upper, const int n)
 	double *x;
 	double dx;
 
+	check_magnitude(__func__, 1, n);
 	x = malloc(n*sizeof(double));
 	dx = (upper - lower)/(double)(n-1);
 	for (i=0; i<n; ++i)
@@ -72,6 +83,7 @@ spacing(const double* x, const int n)
 	int i;
 	double *dx;
 
+	check_magnitude(__func__, 2, n);
 	dx = malloc((n-1)*sizeof(double));
 	for (i=0; i<n-1; ++i)
 		dx[i] = x[i+1] - x[i];
@@ -85,6 +97,7 @@ zeros(const int n)
 	int i;
 	double *x;
 
+	check_magnitude(__func__, 1, n);
 	x = malloc(n*sizeof(double));
 	for (i=0; i<n; ++i)
 		x[i] = 0.0;
@@ -98,6 +111,7 @@ ones(const int n)
 	int i;
 	double *x;
 
+	check_magnitude(__func__, 1, n);
 	x = malloc(n*sizeof(double));
 	for (i=0; i<n; ++i)
 		x[i] = 1.0;
@@ -111,6 +125,7 @@ initval(const double value, const int n)
 	int i;
 	double *x;
 
+	check_magnitude(__func__, 1, n);
 	x = malloc(n*sizeof(double));
 	for (i=0; i<n; ++i)
 		x[i] = value;
@@ -121,6 +136,7 @@ initval(const double value, const int n)
 double*
 empty(const int n)
 {
+	check_magnitude(__func__, 1, n);
 	return malloc(n*sizeof(double));
 }
 
@@ -131,6 +147,7 @@ max(const double* a, const int n)
 	int i;
 	double maxval;
 
+	check_magnitude(__func__, 1, n);
 	maxval = -INFINITY;
 	for (i=0; i<n; ++i)
 		if (a[i] > maxval)
@@ -145,6 +162,7 @@ min(const double* a, const int n)
 	int i;
 	double minval;
 
+	check_magnitude(__func__, 1, n);
 	minval = +INFINITY;
 	for (i=0; i<n; ++i)
 		if (a[i] < minval)
@@ -156,6 +174,7 @@ void
 print_array(const double* a, const int n)
 {
 	int i;
+	check_magnitude(__func__, 1, n);
 	for (i=0; i<n; ++i)
 		printf("%.17g\n", a[i]);
 }
@@ -164,6 +183,7 @@ void
 print_arrays(const double* a, const double* b, const int n)
 {
 	int i;
+	check_magnitude(__func__, 1, n);
 	for (i=0; i<n; ++i)
 		printf("%.17g\t%.17g\n", a[i], b[i]);
 }
@@ -173,7 +193,7 @@ print_arrays_2nd_normalized(const double* a, const double* b, const int n)
 {
 	int i;
 	double max_b;
-
+	check_magnitude(__func__, 1, n);
 	max_b = max(b, n);
 	for (i=0; i<n; ++i)
 		printf("%.17g\t%.17g\n", a[i], b[i]/max_b);
@@ -186,6 +206,7 @@ print_three_arrays(const double* a,
                    const int n)
 {
 	int i;
+	check_magnitude(__func__, 1, n);
 	for (i=0; i<n; ++i)
 		printf("%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]);
 }
@@ -194,6 +215,7 @@ void
 fprint_arrays(FILE* fp, const double* a, const double* b, const int n)
 {
 	int i;
+	check_magnitude(__func__, 1, n);
 	for (i=0; i<n; ++i)
 		fprintf(fp, "%.17g\t%.17g\n", a[i], b[i]);
 }
@@ -206,6 +228,7 @@ fprint_three_arrays(FILE* fp,
                     const int n)
 {
 	int i;
+	check_magnitude(__func__, 1, n);
 	for (i=0; i<n; ++i)
 		fprintf(fp, "%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]);
 }
@@ -214,6 +237,7 @@ void
 copy_values(const double* in, double* out, const int n)
 {
 	int i;
+	check_magnitude(__func__, 1, n);
 	for (i=0; i<n; ++i)
 		out[i] = in[i];
 }
@@ -222,6 +246,7 @@ double*
 copy(const double* in, const int n)
 {
 	double *out;
+	check_magnitude(__func__, 1, n);
 	out = empty(n);
 	copy_values(in, out, n);
 	return out;
@@ -234,6 +259,7 @@ normalize(const double* in, const int n)
 	double max_val;
 	double *out;
 
+	check_magnitude(__func__, 1, n);
 	out = malloc(n*sizeof(double));
 	copy_values(in, out, n);
 	max_val = max(out, n);
diff --git a/max_depth_simple_shear.c b/max_depth_simple_shear.c
@@ -27,12 +27,13 @@ usage(void)
 	printf("%s: %s [OPTIONS]\n"
 	        "outputs the maximum deformation depth in a shear system with sinusoidal\n"
 	        "fluid-pressure perturbations from the top, assuming infinite layer thickness.\n"
-			"Assumes infinite layer thickness.\n
+			"Assumes infinite layer thickness.\n"
 	        "\nOptional arguments:\n"
 	        " -v, --version                   show version information\n"
 	        " -h, --help                      show this message\n"
 	        " -G, --gravity VAL               gravity magnitude [m/s^2] (default %g)\n"
 	        " -P, --normal-stress VAL         normal stress on top [Pa] (default %g)\n"
+	        " -p, --porosity VAL              porosity fraction [-] (default %g)\n"
 	        " -r, --density VAL               grain material density [kg/m^3] (default %g)\n"
 	        " -c, --fluid-compressibility VAL fluid compressibility [Pa^-1] (default %g)\n"
 	        " -i, --fluid-viscosity VAL       fluid viscosity [Pa*s] (default %g)\n"
@@ -48,6 +49,7 @@ usage(void)
 	        __func__, PROGNAME,
 	        sim.G,
 	        sim.P_wall,
+			sim.phi[0],
 	        sim.rho_s,
 	        sim.beta_f,
 	        sim.mu_f,
@@ -70,6 +72,14 @@ version(void)
 	, PROGNAME, VERSION);
 }
 
+
+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));
+}
+
 int
 main(int argc, char* argv[])
 {
@@ -77,6 +87,7 @@ main(int argc, char* argv[])
 	struct simulation sim;
 	const char* optstring;
 	unsigned long iter;
+	double new_phi, new_k;
 #ifdef BENCHMARK_PERFORMANCE
 	clock_t t_begin, t_end;
 	double t_elapsed;
@@ -100,6 +111,7 @@ main(int argc, char* argv[])
 		{"version",                   no_argument,       NULL, 'v'},
 		{"gravity",                   required_argument, NULL, 'G'},
 		{"normal-stress",             required_argument, NULL, 'P'},
+		{"porosity",                  required_argument, NULL, 'p'},
 		{"density",                   required_argument, NULL, 'r'},
 		{"fluid-compressiblity",      required_argument, NULL, 'c'},
 		{"fluid-viscosity",           required_argument, NULL, 'i'},
@@ -135,6 +147,9 @@ main(int argc, char* argv[])
 			case 'P':
 				sim.P_wall = atof(optarg);
 				break;
+			case 'p':
+				new_phi = atof(optarg);
+				break;
 			case 'r':
 				sim.rho_s = atof(optarg);
 				break;
@@ -167,10 +182,10 @@ main(int argc, char* argv[])
 		}
 	}
 
-	if (optind>=argc) {
+	/*if (optind < argc) {
 		fprintf(stderr, "error: unknown argument specified\n");
 		return 1;
-	}
+	}*/
 
 	prepare_arrays(&sim);
 
@@ -190,12 +205,16 @@ main(int argc, char* argv[])
 
 	check_simulation_parameters(&sim);
 
-	while (sim.t <= sim.t_end) {
+	iter = 0;
+	double diff = INFINITY;
+	while (diff > RTOL) {
 
 #ifdef BENCHMARK_PERFORMANCE
 		t_begin = clock();
 #endif
 
+		diff = 0.0;
+
 
 
 
@@ -205,22 +224,12 @@ main(int argc, char* argv[])
 		printf("time spent per time step = %g s\n", t_elapsed);
 #endif
 
-		if (!isnan(sim.v_x_limit))
-			sim.mu_wall = mu_wall_orig;
-		sim.t += sim.dt;
-		filetimeclock += sim.dt;
+		/* printf("%s\n", max_depth); */
 		iter++;
-
-		if ((filetimeclock - sim.dt >= sim.file_dt || iter == 1) &&
-		    strcmp(sim.name, DEFAULT_SIMULATION_NAME) != 0) {
-			write_output_file(&sim, norm);
-			filetimeclock = 0.0;
-			if (iter == 1)
-				sim.t = 0.0;
-		}
 	}
 
-	printf("%s\n", max_depth);
+	printf("skin depth: %g m", skin_depth(&sim));
+
 
 	free_arrays(&sim);
 	return 0;