cngf-pf

commit a06d8b2272c0dd6273b7f1abe88e6d4921db6a5a
parent 4f4e70d16f79d863d20480649887c9c1431cf783
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Mon,  6 Apr 2020 12:47:31 +0200

Add debug output, fix porosity computation, change parameter defaults

Diffstat:
M
1d_fd_simple_shear.c
 | 
15
+++++++++++++++
M
parameter_defaults.h
 | 
4
++--
M
simulation.c
 | 
2
+-

3 files changed, 18 insertions(+), 3 deletions(-)
diff --git a/1d_fd_simple_shear.c b/1d_fd_simple_shear.c
@@ -264,23 +264,38 @@ main(int argc, char* argv[])
 
 		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)
diff --git a/parameter_defaults.h b/parameter_defaults.h
@@ -42,8 +42,8 @@ init_sim()
 	sim.phi = initval(0.25, 1);
 	sim.d = 0.04;        /* Damsgaard et al 2013 */
 
-	sim.phi_min = 0.2;
-	sim.phi_max = 0.8;
+	sim.phi_min = 0.25;
+	sim.phi_max = 0.55;
 	sim.dilatancy_angle = 1.0;
 
 	/* Iverson et al 1997, 1998: Storglaciaren till */
diff --git a/simulation.c b/simulation.c
@@ -323,7 +323,7 @@ compute_porosity_change()
 {
 	int i;
 	for (i=0; i<sim.nz; ++i) {
-		sim.phi_dot[i] = sim.tan_psi[i]*sim.gamma_dot_p[i];
+		sim.phi_dot[i] = sim.tan_psi[i]*sim.gamma_dot_p[i]*sim.phi[i];
 		sim.phi[i] += sim.phi_dot[i]*sim.dt;
 	}
 }