cngf-pf

commit f74ed98b4c7ea4bde059ced6673690a20f85e016
parent 699306346b7a372e537e8b3edd7598c7ec419411
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Mon,  6 Apr 2020 13:50:25 +0200

Move main loop into simulation.c

Diffstat:
M
1d_fd_simple_shear.c
 | 
102
++++++-------------------------------------------------------------------------
M
fluid.c
 | 
2
+-
M
simulation.c
 | 
109
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
M
simulation.h
 | 
2
++

4 files changed, 119 insertions(+), 96 deletions(-)
diff --git a/1d_fd_simple_shear.c b/1d_fd_simple_shear.c
@@ -12,16 +12,8 @@
 #include "parameter_defaults.h"
 
 /* relative tolerance criteria for granular fluidity solver */
-#define RTOL_GRANULAR 1e-5
-#define MAX_ITER_GRANULAR 10000
-
-/* relative tolerance criteria for fluid-pressure solver */
-#define RTOL_DARCY 1e-5
-#define MAX_ITER_DARCY 10000
-
-/* relative tolerance criteria when shear velocity is restricted */
-#define RTOL_STRESS 1e-3
-#define MAX_ITER_STRESS 20000
+#define RTOL 1e-5
+#define MAX_ITER_1D_FD_SIMPLE_SHEAR 10000
 
 /* uncomment to print time spent per time step to stdout */
 /*#define BENCHMARK_PERFORMANCE*/
@@ -77,8 +69,8 @@ int
 main(int argc, char* argv[])
 {
 	int i, normalize;
-	unsigned long iter, stressiter;
-	double new_phi, new_k, filetimeclock, res_norm, mu_wall_orig;
+	unsigned long iter;
+	double new_phi, new_k, filetimeclock, mu_wall_orig;
 #ifdef BENCHMARK_PERFORMANCE
 	clock_t t_begin, t_end;
 	double t_elapsed;
@@ -87,7 +79,7 @@ main(int argc, char* argv[])
 #ifdef __OpenBSD__
 	if (pledge("stdio wpath cpath", NULL) == -1) {
 		fprintf(stderr, "error: pledge failed");
-		exit(1);
+		exit(2);
 	}
 #endif
 
@@ -255,94 +247,14 @@ main(int argc, char* argv[])
 	filetimeclock = 0.0;
 	iter = 0;
 	mu_wall_orig = sim.mu_wall;
-	res_norm = NAN;
 	while (sim.t <= sim.t_end) {
 
 #ifdef BENCHMARK_PERFORMANCE
 		t_begin = clock();
 #endif
 
-		stressiter = 0;
-		do {
-			printf("\niter %lu, sim.t = %g s\n", iter, sim.t);
-
-			if (sim.transient) {
-				/* step 1 */
-				compute_inertia_number();			/* Eq. 1 */
-				for (i=0; i<sim.nz; i++)
-					printf("sim.I[%d] = %g\n", i, sim.I[i]);
-				/* step 2 */
-				compute_critical_state_porosity();	/* Eq. 2 */
-				for (i=0; i<sim.nz; i++)
-					printf("sim.phi_c[%d] = %g\n", i, sim.phi_c[i]);
-				/* step 3 */
-				compute_tan_dilatancy_angle();		/* Eq. 5 */
-				for (i=0; i<sim.nz; i++)
-					printf("sim.tan_psi[%d] = %g\n", i, sim.tan_psi[i]);
-			}
-			compute_critical_state_friction();	/* Eq. 7 */
-			for (i=0; i<sim.nz; i++)
-				printf("sim.mu_c[%d] = %g\n", i, sim.mu_c[i]);
-
-			/* step 4 */
-			if (sim.transient) {
-				compute_porosity_change();			/* Eq. 3 */
-				for (i=0; i<sim.nz; i++)
-					printf("sim.phi_dot[%d] = %g\n", i, sim.phi_dot[i]);
-				compute_permeability();				/* Eq. 6 */
-				for (i=0; i<sim.nz; i++)
-					printf("sim.k[%d] = %g\n", i, sim.k[i]);
-			}
-			compute_friction();						/* Eq. 4 */
-			for (i=0; i<sim.nz; i++)
-				printf("sim.mu[%d] = %g\n", i, sim.mu[i]);
-
-			/* step 5, Eq. 13 */
-			if (sim.fluid)
-				if (darcy_solver_1d(MAX_ITER_DARCY, RTOL_DARCY))
-					exit(1);
-
-			/* step 6 */
-			compute_effective_stress();				/* Eq. 9 */
-
-			/* step 7 */
-			compute_local_fluidity();				/* Eq. 10 */
-			compute_cooperativity_length();			/* Eq. 12 */
-
-			/* step 8, Eq. 11 */
-			if (implicit_1d_jacobian_poisson_solver(MAX_ITER_GRANULAR,
-			                                        RTOL_GRANULAR))
-				exit(1);
-
-			/* step 9 */
-			compute_shear_strain_rate_plastic();	/* Eq. 8 */
-			compute_inertia_number();
-			compute_shear_velocity();
-
-			if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
-				if (!isnan(sim.v_x_limit)) {
-					res_norm = (sim.v_x_limit - sim.v_x[sim.nz-1])
-					           /(sim.v_x[sim.nz-1] + 1e-12);
-					if (res_norm > 0.0)
-						res_norm = 0.0;
-				} else {
-					res_norm = (sim.v_x_fix - sim.v_x[sim.nz-1])
-					           /(sim.v_x[sim.nz-1] + 1e-12);
-				}
-				sim.mu_wall *= 1.0 + (res_norm*1e-2);
-			}
-
-			if (++stressiter > MAX_ITER_STRESS) {
-				fprintf(stderr, "error: stress solution did not converge:\n");
-				fprintf(stderr, "v_x=%g, v_x_fix=%g, v_x_limit=%g, "
-				        "res_norm=%g, mu_wall=%g\n",
-				        sim.v_x[sim.nz-1], sim.v_x_fix, sim.v_x_limit,
-				        res_norm, sim.mu_wall);
-				return 10;
-			}
-
-		} while ((!isnan(sim.v_x_fix) || !isnan(sim.v_x_limit))
-		         && fabs(res_norm) > RTOL_STRESS);
+		if (coupled_shear_solver(MAX_ITER_1D_FD_SIMPLE_SHEAR, RTOL))
+			exit(10);
 
 #ifdef BENCHMARK_PERFORMANCE
 		t_end = clock();
diff --git a/fluid.c b/fluid.c
@@ -32,7 +32,7 @@ set_largest_fluid_timestep(const double safety)
 			fprintf(stderr, "error: cell spacing negative (%g) in cell %d\n",
 					dx[i], i);
 			free(dx);
-			exit(10);
+			exit(1);
 		}
 		if (dx[i] < dx_min) dx_min = dx[i];
 	}
diff --git a/simulation.c b/simulation.c
@@ -3,9 +3,22 @@
 #include <math.h>
 #include "arrays.h"
 #include "simulation.h"
+#include "fluid.h"
 
 /* #define SHOW_PARAMETERS */
 
+/* relative tolerance criteria for granular fluidity solver */
+#define RTOL_GRANULAR 1e-5
+#define MAX_ITER_GRANULAR 10000
+
+/* relative tolerance criteria for fluid-pressure solver */
+#define RTOL_DARCY 1e-5
+#define MAX_ITER_DARCY 10000
+
+/* relative tolerance criteria when shear velocity is restricted */
+#define RTOL_STRESS 1e-3
+#define MAX_ITER_STRESS 20000
+
 void
 prepare_arrays()
 {
@@ -621,3 +634,99 @@ print_output(FILE* fp, const int norm)
 
 	free(v_x_out);
 }
+
+int
+coupled_shear_solver(const int max_iter, 
+                     const double rel_tol)
+{
+	int coupled_iter, stress_iter;
+
+	double res_norm;
+	res_norm = NAN;
+
+	stress_iter = 0;
+	do {  /* stress iterations */
+
+		coupled_iter = 0;
+		do {  /* coupled iterations */
+
+			if (sim.transient) {
+				/* step 1 */
+				compute_inertia_number();			/* Eq. 1 */
+				/* step 2 */
+				compute_critical_state_porosity();	/* Eq. 2 */
+				/* step 3 */
+				compute_tan_dilatancy_angle();		/* Eq. 5 */
+			}
+			compute_critical_state_friction();		/* Eq. 7 */
+
+			/* step 4 */
+			if (sim.transient) {
+				compute_porosity_change();			/* Eq. 3 */
+				compute_permeability();				/* Eq. 6 */
+			}
+			compute_friction();						/* Eq. 4 */
+
+			/* step 5, Eq. 13 */
+			if (sim.fluid)
+				if (darcy_solver_1d(MAX_ITER_DARCY, RTOL_DARCY))
+					exit(11);
+
+			/* step 6 */
+			compute_effective_stress();				/* Eq. 9 */
+
+			/* step 7 */
+			compute_local_fluidity();				/* Eq. 10 */
+			compute_cooperativity_length();			/* Eq. 12 */
+
+			/* step 8, Eq. 11 */
+			if (implicit_1d_jacobian_poisson_solver(MAX_ITER_GRANULAR,
+													RTOL_GRANULAR))
+				exit(12);
+
+			/* step 9 */
+			compute_shear_strain_rate_plastic();	/* Eq. 8 */
+			compute_inertia_number();
+			compute_shear_velocity();
+
+#ifdef DEBUG
+			for (i=0; i<sim.nz) {
+				printf("\nsim.t = %g s\n", sim.t);
+				printf("sim.I[%d] = %g\n", i, sim.I[i]);
+				printf("sim.phi_c[%d] = %g\n", i, sim.phi_c[i]);
+				printf("sim.tan_psi[%d] = %g\n", i, sim.tan_psi[i]);
+				printf("sim.mu_c[%d] = %g\n", i, sim.mu_c[i]);
+				printf("sim.phi_dot[%d] = %g\n", i, sim.phi_dot[i]);
+				printf("sim.k[%d] = %g\n", i, sim.k[i]);
+				printf("sim.mu[%d] = %g\n", i, sim.mu[i]);
+			}
+#endif
+
+			if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
+				if (!isnan(sim.v_x_limit)) {
+					res_norm = (sim.v_x_limit - sim.v_x[sim.nz-1])
+							   /(sim.v_x[sim.nz-1] + 1e-12);
+					if (res_norm > 0.0)
+						res_norm = 0.0;
+				} else {
+					res_norm = (sim.v_x_fix - sim.v_x[sim.nz-1])
+							   /(sim.v_x[sim.nz-1] + 1e-12);
+				}
+				sim.mu_wall *= 1.0 + (res_norm*1e-2);
+			}
+
+			if (++stress_iter > MAX_ITER_STRESS) {
+				fprintf(stderr, "error: stress solution did not converge:\n");
+				fprintf(stderr, "v_x=%g, v_x_fix=%g, v_x_limit=%g, "
+						"res_norm=%g, mu_wall=%g\n",
+						sim.v_x[sim.nz-1], sim.v_x_fix, sim.v_x_limit,
+						res_norm, sim.mu_wall);
+				return 10;
+			}
+		} while (0);
+
+	} while ((!isnan(sim.v_x_fix) || !isnan(sim.v_x_limit))
+			 && fabs(res_norm) > RTOL_STRESS);
+	return 0;
+}
+
diff --git a/simulation.h b/simulation.h
@@ -153,4 +153,6 @@ int implicit_1d_jacobian_poisson_solver(const int max_iter,
 void write_output_file(const int normalize);
 void print_output(FILE* fp, const int normalize);
 
+int coupled_shear_solver(const int max_iter, const double rel_tol);
+
 #endif