1d_fd_simple_shear

commit 5e65e6215183d5f05212487071b69408d72bbe28
parent 77f3027af56d193b84b7e97d29a82fa9e685756c
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Thu, 16 Apr 2020 15:04:00 +0200

Check norm from velocity field, use common function for residual

Diffstat:
M
simulation.c
 | 
73
+++++++++++++++++++++++++++++++++++++++++++++++--------------------------

1 file changed, 47 insertions(+), 26 deletions(-)
diff --git a/simulation.c b/simulation.c
@@ -598,6 +598,12 @@ set_bc_dirichlet(double* g_ghost,
 	}
 }
 
+static double
+residual_normalized(double new, double old)
+{
+	return pow(new - old, 2.0)/(pow(new, 2.0) + 1e-16);
+}
+
 static void
 poisson_solver_1d_cell_update(int i,
                               const double* g_in,
@@ -613,8 +619,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] = pow(g_out[i+1] - g_in[i+1], 2.0)
-	            /(pow(g_out[i+1], 2.0) + 1e-16);
+	r_norm[i] = residual_normalized(g_out[i+1], g_in[i+1]);
 
 #ifdef DEBUG
 	printf("-- %d --------------\n", i);
@@ -732,13 +737,16 @@ coupled_shear_solver(struct simulation *sim,
                      const int max_iter, 
                      const double rel_tol)
 {
-	int coupled_iter, stress_iter;
+	int i, coupled_iter, stress_iter;
 	double vel_res_norm, mu_wall_orig;
-	double *r_norm;
+	double *r_norm, r_norm_max, *v_x0;
 
 	vel_res_norm = NAN;
-	r_norm = empty(sim->nz);
 	mu_wall_orig = sim->mu_wall;
+	if (sim->transient) {
+		r_norm = empty(sim->nz);
+		v_x0 = empty(sim->nz);
+	}
 
 	stress_iter = 0;
 	do {  /* stress iterations */
@@ -747,6 +755,8 @@ coupled_shear_solver(struct simulation *sim,
 		do {  /* coupled iterations */
 
 			if (sim->transient) {
+				copy_values(sim->v_x, v_x0, sim->nz);
+
 				/* step 1 */
 				compute_inertia_number(sim);			/* Eq. 1 */
 				/* step 2 */
@@ -798,37 +808,48 @@ coupled_shear_solver(struct simulation *sim,
 				printf("sim->mu[%d] = %g\n", i, sim->mu[i]);
 			}
 #endif
-
-			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])
-					               /(sim->v_x[sim->nz-1] + 1e-12);
-					if (vel_res_norm > 0.0)
-						vel_res_norm = 0.0;
-				} else {
-					vel_res_norm = (sim->v_x_fix - sim->v_x[sim->nz-1])
-					               /(sim->v_x[sim->nz-1] + 1e-12);
-				}
-				sim->mu_wall *= 1.0 + (vel_res_norm*1e-2);
+			/* 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_max = max(r_norm, sim->nz);
 			}
 
-			if (++stress_iter > MAX_ITER_STRESS) {
-				fprintf(stderr, "error: stress solution did not converge:\n");
-				fprintf(stderr, "v_x=%g, v_x_fix=%g, v_x_limit=%g, "
-				        "vel_res_norm=%g, mu_wall=%g\n",
-				        sim->v_x[sim->nz-1], sim->v_x_fix, sim->v_x_limit,
-				        vel_res_norm, sim->mu_wall);
-				return 10;
-			}
+
 		} while (0);
 
+		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])
+							   /(sim->v_x[sim->nz-1] + 1e-12);
+				if (vel_res_norm > 0.0)
+					vel_res_norm = 0.0;
+			} else {
+				vel_res_norm = (sim->v_x_fix - sim->v_x[sim->nz-1])
+							   /(sim->v_x[sim->nz-1] + 1e-12);
+			}
+			sim->mu_wall *= 1.0 + (vel_res_norm*1e-2);
+		}
+
+		if (++stress_iter > MAX_ITER_STRESS) {
+			fprintf(stderr, "error: stress solution did not converge:\n");
+			fprintf(stderr, "v_x=%g, v_x_fix=%g, v_x_limit=%g, "
+					"vel_res_norm=%g, mu_wall=%g\n",
+					sim->v_x[sim->nz-1], sim->v_x_fix, sim->v_x_limit,
+					vel_res_norm, sim->mu_wall);
+			return 10;
+		}
+
 	} while ((!isnan(sim->v_x_fix) || !isnan(sim->v_x_limit))
 			 && fabs(vel_res_norm) > RTOL_STRESS);
 
 	if (!isnan(sim->v_x_limit))
 		sim->mu_wall = mu_wall_orig;
 
-	free(r_norm);
+	if (sim->transient) {
+		free(r_norm);
+		free(v_x0);
+	}
 
 	return 0;
 }