cngf-pf

commit 5ea71c9d6514081e1f48b4f781e0856304a941f7
parent 1b7a4d4eb534ef7b24357bb256d4bb7348047ce4
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Tue,  8 Dec 2020 10:53:19 +0100

Merge branch 'master' of src.adamsgaard.dk:cngf-pf

Diffstat:
M
Makefile
 | 
14
+++++++-------
M
cngf-pf.1
 | 
6
+++++-
M
cngf-pf.c
 | 
3
+--
M
simulation.c
 | 
26
++++++++++----------------

4 files changed, 23 insertions(+), 26 deletions(-)
diff --git a/Makefile b/Makefile
@@ -39,25 +39,25 @@ OBJ = ${SRC:.c=.o}
 ${OBJ}: ${HDR}
 
 .o:
-	${CC} ${CNGFPFLDFLAGS} -o $@ $<
+	${CC} -o $@ $< ${CNGFPFLDFLAGS}
 
 .c.o:
 	${CC} ${CNGFPFCFLAGS} ${CNGFPFCPPFLAGS} -o $@ -c $<
 
 cngf-pf: cngf-pf.o arrays.o fluid.o simulation.o
-	${CC} ${CNGFPFLDFLAGS}\
+	${CC}\
 		cngf-pf.o arrays.o fluid.o simulation.o\
-		-o $@
+		-o $@ ${CNGFPFLDFLAGS}
 
 max_depth_simple_shear: max_depth_simple_shear.o arrays.o fluid.o simulation.o
-	${CC} ${CNGFPFLDFLAGS}\
+	${CC}\
 		max_depth_simple_shear.o arrays.o fluid.o simulation.o\
-		-o $@
+		-o $@ ${CNGFPFLDFLAGS}
 
 shear_flux: shear_flux.o
-	${CC} ${CNGFPFLDFLAGS}\
+	${CC}\
 		shear_flux.o\
-		-o $@
+		-o $@ ${CNGFPFLDFLAGS}
 
 dist:
 	rm -rf "${NAME}-${VERSION}"
diff --git a/cngf-pf.1 b/cngf-pf.1
@@ -250,12 +250,16 @@ granular solver error
 .It 20
 time step error
 .El
-.\" .Sh EXAMPLES
+.Sh EXAMPLES
+Plot critical-state shear velocity as a function of depth:
+.Pp
+.Dl $ cngf-pf | gnuplot -e 'set term pdf; p \&"-\&" u 2:1 w lp' > out.pdf
 .\" .Sh DIAGNOSTICS
 .\" For sections 1, 4, 6, 7, 8, and 9 printf/stderr messages only.
 .\" .Sh ERRORS
 .\" For sections 2, 3, 4, and 9 errno settings only.
 .Sh SEE ALSO
+.Xr gnuplot 1
 .Xr max_depth_simple_shear 1
 .Xr shear_flux 1
 .Sh AUTHORS
diff --git a/cngf-pf.c b/cngf-pf.c
@@ -10,7 +10,7 @@
 
 #include "arg.h"
 
-/* relative tolerance criteria for granular fluidity solver */
+/* relative tolerance criteria for the solvers */
 #define RTOL 1e-5
 #define MAX_ITER_1D_FD_SIMPLE_SHEAR 100000
 
@@ -76,7 +76,6 @@ main(int argc, char *argv[])
 #ifdef BENCHMARK_PERFORMANCE
 	clock_t t_begin, t_end;
 	double t_elapsed;
-
 #endif
 
 #ifdef __OpenBSD__
diff --git a/simulation.c b/simulation.c
@@ -7,20 +7,14 @@
 #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
+/* iteration limits for solvers */
+#define MAX_ITER_GRANULAR 100000
 #define MAX_ITER_DARCY 10000
-
-/* relative tolerance criteria when shear velocity is restricted */
-#define RTOL_STRESS 1e-3
 #define MAX_ITER_STRESS 20000
 
+/* tolerance criteria when in velocity driven or velocity limited mode */
+#define RTOL_VELOCITY 1e-3
+
 
 /* Simulation settings */
 void
@@ -457,7 +451,7 @@ compute_friction(struct simulation *sim)
 
 	if (sim->transient)
 		for (i = 0; i < sim->nz; ++i)
-			sim->mu[i] = sim->mu_c[i] - sim->tan_psi[i];
+			sim->mu[i] = sim->mu_c[i] + sim->tan_psi[i];
 	else
 		for (i = 0; i < sim->nz; ++i)
 			sim->mu[i] = sim->mu_c[i];
@@ -640,7 +634,7 @@ poisson_solver_1d_cell_update(int i,
 	g_out[i + 1] = 1.0 / (1.0 + coorp_term)
 	               * (coorp_term * g_local[i] + g_in[i + 2] / 2.0
 	                  + g_in[i] / 2.0);
-	r_norm[i] = residual(g_out[i + 1], g_in[i + 1]);
+	r_norm[i] = fabs(residual(g_out[i + 1], g_in[i + 1]));
 
 #ifdef DEBUG
 	printf("-- %d --------------\n", i);
@@ -805,7 +799,7 @@ coupled_shear_solver(struct simulation *sim,
 
 			/* step 5, Eq. 13 */
 			if (sim->fluid)
-				if (darcy_solver_1d(sim, MAX_ITER_DARCY, RTOL_DARCY))
+				if (darcy_solver_1d(sim, MAX_ITER_DARCY, rel_tol))
 					exit(11);
 
 			/* step 6 */
@@ -818,7 +812,7 @@ coupled_shear_solver(struct simulation *sim,
 			/* step 8, Eq. 11 */
 			if (implicit_1d_jacobian_poisson_solver(sim,
 			                                        MAX_ITER_GRANULAR,
-			                                        RTOL_GRANULAR))
+			                                        rel_tol))
 				exit(12);
 
 			/* step 9 */
@@ -880,7 +874,7 @@ coupled_shear_solver(struct simulation *sim,
 			return 10;
 		}
 	} while ((!isnan(sim->v_x_fix) || !isnan(sim->v_x_limit))
-	        && fabs(vel_res_norm) > RTOL_STRESS);
+	        && fabs(vel_res_norm) > RTOL_VELOCITY);
 
 	if (!isnan(sim->v_x_limit))
 		sim->mu_wall = mu_wall_orig;