cngf-pf

commit cc5fef7b04e30a4b6c1d41c4e3d47dddd2e45173
parent a91a9d97e214821989e66263c3de85574dc6871e
Author: Ian Madden <iamadden@stanford.edu>
Date:   Thu, 11 Feb 2021 09:41:27 -0800

First try at a more appropriate time step when coupling the fluid and transient solvers.

Signed-off-by: Anders Damsgaard <anders@adamsgaard.dk>

Diffstat:
M
cngf-pf.c
 | 
5
+++++
M
simulation.c
 | 
24
+++++++++++++++++++++++-
M
simulation.h
 | 
2
++

3 files changed, 30 insertions(+), 1 deletion(-)
diff --git a/cngf-pf.c b/cngf-pf.c
@@ -254,8 +254,13 @@ main(int argc, char *argv[])
 				free_arrays(&sim);
 				return 20;
 			}
+			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,6 +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,7 +853,6 @@ coupled_shear_solver(struct simulation *sim,
 				}
 			}
 		} while (sim->transient);
-
 		if (!isnan(sim->v_x_limit) || !isnan(sim->v_x_fix)) {
 			if (!isnan(sim->v_x_limit)) {
 				vel_res_norm = (sim->v_x_limit - sim->v_x[sim->nz - 1])
@@ -895,3 +895,25 @@ 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;
+}
diff --git a/simulation.h b/simulation.h
@@ -153,6 +153,8 @@ void print_output(struct simulation *sim, FILE *fp, const int normalize);
 int coupled_shear_solver(struct simulation *sim,
                          const int max_iter,
                          const double rel_tol);
+                         
+int set_coupled_fluid_transient_timestep(struct simulation *sim, const double safety);
 
 double find_flux(const struct simulation *sim);