cngf-pf

commit ad92a3ee9ff0823547bdc761073f8bb3370751b4
parent cc5fef7b04e30a4b6c1d41c4e3d47dddd2e45173
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Tue, 16 Feb 2021 15:30:10 +0100

style fixes for cc5fef7

Diffstat:
M
cngf-pf.c
 | 
4
+---
M
simulation.c
 | 
59
++++++++++++++++++++++++++++++++++++++---------------------

2 files changed, 39 insertions(+), 24 deletions(-)
diff --git a/cngf-pf.c b/cngf-pf.c
@@ -254,13 +254,11 @@ main(int argc, char *argv[])
 				free_arrays(&sim);
 				return 20;
 			}
-			if (sim.transient) {
+			if (sim.transient)
 			  set_coupled_fluid_transient_timestep(&sim, 0.5);
-			}
 		}
 	}
 	fprintf(stderr, "t_val = %g \n", sim.dt);
-	
 
 	if (sim.dt > sim.file_dt)
 		sim.dt = sim.file_dt;
diff --git a/simulation.c b/simulation.c
@@ -834,7 +834,7 @@ coupled_shear_solver(struct simulation *sim,
 				printf("sim->mu[%d] = %g\n", i, sim->mu[i]);
 			}
 #endif
-      
+
 			/* stable porosity change field == coupled solution converged */
 			if (sim->transient) {
 				for (i = 0; i < sim->nz; ++i)
@@ -852,6 +852,7 @@ coupled_shear_solver(struct simulation *sim,
 					return 1;
 				}
 			}
+
 		} while (sim->transient);
 		if (!isnan(sim->v_x_limit) || !isnan(sim->v_x_fix)) {
 			if (!isnan(sim->v_x_limit)) {
@@ -896,24 +897,40 @@ find_flux(const struct simulation *sim)
 	return flux;
 }
 
-int set_coupled_fluid_transient_timestep(struct simulation *sim, const double safety) {
-  // Compute Maximum Strain Rate Expected
-  double max_gamma_dot = 1/sim->d;
-  if (!isnan(sim->v_x_fix)) max_gamma_dot = sim->v_x_fix / sim->dz;
-  if (!isnan(sim->v_x_limit)) max_gamma_dot = sim->v_x_limit / sim->dz;
-  // Compute estimate for mu for cooperativity length
-  double mu = (sim->mu_wall / ((sim->sigma_n[sim->nz-1]- sim->p_f_mod_ampl)/ (sim->P_wall - sim->p_f_top))) + sim->dilatancy_constant * sim->phi[sim->nz-1];
-  // Compute estimate for cooperativity length 
-  double xi = cooperativity_length(sim->A,
-		                                  sim->d,
-		                                  mu,
-		                                  (sim->sigma_n[sim->nz-1]- sim->p_f_mod_ampl),
-		                                  sim->mu_s,
-		                                  sim->C);
-  // Compute Maximum Expected Inertia Number
-  double max_I = max_gamma_dot * sim->d / sqrt((sim->sigma_n[sim->nz-1]- sim->p_f_mod_ampl) / sim->rho_s);
-	double t_val = xi  *(sim->alpha + sim->phi[sim->nz-1]*sim->beta_f) * sim->mu_f
-                / (sim->phi[sim->nz-1] * sim->phi[sim->nz-1]* sim->phi[sim->nz-1]*sim->L_z* max_I);
-  if (sim->dt > safety * t_val) sim->dt = safety * t_val;
-  return 0;
+int
+set_coupled_fluid_transient_timestep(struct simulation *sim, const double safety)
+{
+	double max_gamma_dot, mu, xi, max_i, t_val;
+
+	/* max expected strain rate */
+	max_gamma_dot = 1.0/sim->d;
+	if (!isnan(sim->v_x_fix))
+		max_gamma_dot = sim->v_x_fix / sim->dz;
+	if (!isnan(sim->v_x_limit))
+		max_gamma_dot = sim->v_x_limit / sim->dz;
+
+	/* estimate for shear friction */
+	mu = (sim->mu_wall / ((sim->sigma_n[sim->nz-1] - sim->p_f_mod_ampl)
+	     / (sim->P_wall - sim->p_f_top)))
+	     + sim->dilatancy_constant * sim->phi[sim->nz-1];
+
+	/* estimate for cooperativity length  */
+	xi = cooperativity_length(sim->A,
+	                          sim->d,
+	                          mu,
+	                          (sim->sigma_n[sim->nz - 1] - sim->p_f_mod_ampl),
+	                          sim->mu_s,
+	                          sim->C);
+
+	/* max expected inertia number */
+	max_I = max_gamma_dot * sim->d
+	        / sqrt((sim->sigma_n[sim->nz - 1] - sim->p_f_mod_ampl) / sim->rho_s);
+	t_val = xi * (sim->alpha + sim->phi[sim->nz - 1]*sim->beta_f) * sim->mu_f
+	        / (sim->phi[sim->nz - 1] * sim->phi[sim->nz - 1]
+	           * sim->phi[sim->nz - 1] * sim->L_z * max_I);
+
+	if (sim->dt > safety * t_val)
+		sim->dt = safety * t_val;
+
+	return 0;
 }