cngf-pf

commit b47478c82b40a6ea1967da7c5eb9161e51688815
parent d9168f9f815c36795fe12bac4cb8b147afcc59e6
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Thu, 16 Apr 2020 15:36:26 +0200

Use linear operation and common function for darcy and coupled solver

Diffstat:
M
fluid.c
 | 
4
++--
M
simulation.c
 | 
25
+++++++++++++------------
M
simulation.h
 | 
15
++++++++-------

3 files changed, 23 insertions(+), 21 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -248,8 +248,8 @@ darcy_solver_1d(struct simulation *sim,
 				p_f_ghost_new[i+1] = p_f_ghost_old[i+1]*(1.0 - theta)
 				                     + p_f_ghost_new[i+1]*theta;
 
-				r_norm[i] = fabs((p_f_ghost_new[i+1] - sim->p_f_ghost[i+1])
-				            /(sim->p_f_ghost[i+1] + 1e-16));
+				r_norm[i] = residual_normalized(p_f_ghost_new[i+1],
+												sim->p_f_ghost[i+1]);
 			}
 
 			r_norm_max = max(r_norm, sim->nz);
diff --git a/simulation.c b/simulation.c
@@ -598,10 +598,10 @@ set_bc_dirichlet(double* g_ghost,
 	}
 }
 
-static double
+double
 residual_normalized(double new, double old)
 {
-	return pow(new - old, 2.0)/(pow(new, 2.0) + 1e-16);
+	return fabs((new - old)/(new + 1e-16));
 }
 
 static void
@@ -619,7 +619,9 @@ 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_normalized(g_out[i+1], g_in[i+1]);
+	/* TODO: try out residual_normalized(..) */
+	r_norm[i] = pow(g_out[i+1] - g_in[i+1], 2.0)
+	            /(pow(g_out[i+1], 2.0) + 1e-16);
 
 #ifdef DEBUG
 	printf("-- %d --------------\n", i);
@@ -739,13 +741,13 @@ coupled_shear_solver(struct simulation *sim,
 {
 	int i, coupled_iter, stress_iter;
 	double vel_res_norm, mu_wall_orig;
-	double *r_norm, r_norm_max, *v_x0;
+	double *r_norm, r_norm_max, *oldval;
 
 	vel_res_norm = NAN;
 	mu_wall_orig = sim->mu_wall;
 	if (sim->transient) {
 		r_norm = empty(sim->nz);
-		v_x0 = empty(sim->nz);
+		oldval = empty(sim->nz);
 	}
 
 	stress_iter = 0;
@@ -755,7 +757,7 @@ coupled_shear_solver(struct simulation *sim,
 		do {  /* coupled iterations */
 
 			if (sim->transient) {
-				copy_values(sim->v_x, v_x0, sim->nz);
+				copy_values(sim->gamma_dot_p, oldval, sim->nz);
 
 				/* step 1 */
 				compute_inertia_number(sim);			/* Eq. 1 */
@@ -811,15 +813,14 @@ coupled_shear_solver(struct simulation *sim,
 			/* stable velocity field == coupled solution converged */
 			if (sim->transient) {
 				for (i=0; i<sim->nz; ++i)
-					r_norm[i] = residual_normalized(sim->v_x[i], v_x0[i]);
+					r_norm[i] = residual_normalized(sim->gamma_dot_p[i], oldval[i]);
 				r_norm_max = max(r_norm, sim->nz);
-
-				if (r_norm_max <= rel_tol) {
+				if (r_norm_max <= rel_tol)
 					break;
-				}
+
 				if (coupled_iter++ >= max_iter) {
 					free(r_norm);
-					free(v_x0);
+					free(oldval);
 					fprintf(stderr, "coupled_shear_solver: ");
 					fprintf(stderr, "Transient solution did not converge "
 					                "after %d iterations\n", coupled_iter);
@@ -860,7 +861,7 @@ coupled_shear_solver(struct simulation *sim,
 
 	if (sim->transient) {
 		free(r_norm);
-		free(v_x0);
+		free(oldval);
 	}
 
 	return 0;
diff --git a/simulation.h b/simulation.h
@@ -131,6 +131,13 @@ void compute_porosity_change(struct simulation *sim);
 void compute_permeability(struct simulation *sim);
 void compute_tan_dilatancy_angle(struct simulation *sim);
 
+void compute_cooperativity_length(struct simulation *sim);
+void compute_local_fluidity(struct simulation *sim);
+void compute_shear_strain_rate_plastic(struct simulation *sim);
+void compute_shear_velocity(struct simulation *sim);
+void compute_effective_stress(struct simulation *sim);
+void compute_friction(struct simulation *sim);
+
 void set_bc_neumann(double* g_ghost,
                     const int nz,
                     const int boundary,
@@ -142,13 +149,7 @@ void set_bc_dirichlet(double* g_ghost,
                       const int boundary,
                       const double value);
 
-void compute_cooperativity_length(struct simulation *sim);
-void compute_local_fluidity(struct simulation *sim);
-void compute_shear_strain_rate_plastic(struct simulation *sim);
-void compute_shear_velocity(struct simulation *sim);
-void compute_effective_stress(struct simulation *sim);
-void compute_friction(struct simulation *sim);
-
+double residual_normalized(double new, double old);
 int implicit_1d_jacobian_poisson_solver(struct simulation *sim,
                                         const int max_iter,
                                         const double rel_tol);