1d_fd_simple_shear_transient

commit 599c08baed18df47fd5a6fac8b587e5d51553519
parent 802c24eb1aed07014ed6c8f0ba825298b24a8325
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Mon,  1 Jul 2019 21:12:56 +0200

Rename p_f_ghost_out to p_f_ghost_new, and add p_f_ghost_old with old vals saved

Diffstat:
M
fluid.c
 | 
33
++++++++++++++++++++-------------

1 file changed, 20 insertions(+), 13 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -117,7 +117,8 @@ darcy_solver_1d(struct simulation* sim,
 {
 	int i, iter, solved;
 	double epsilon, theta, p_f, p_f_top, r_norm_max;
-	double *dp_f_dt_expl, *dp_f_dt_impl, *p_f_ghost_out, *r_norm;
+	double *dp_f_dt_expl;
+	double *p_f_ghost_old, *dp_f_dt_impl, *p_f_ghost_new, *r_norm;
 
 	/* choose integration method, parameter in [0.0; 1.0]
 	 *     epsilon = 0.0: explicit
@@ -162,8 +163,10 @@ darcy_solver_1d(struct simulation* sim,
 	if (epsilon > 0.0) {
 		/* compute implicit solution with Jacobian iterations */
 		dp_f_dt_impl = zeros(sim->nz);
-		p_f_ghost_out = zeros(sim->nz+2);
+		p_f_ghost_new = zeros(sim->nz+2);
 		r_norm = zeros(sim->nz);
+		p_f_ghost_old = empty(sim->nz+2);
+		copy_values(sim->p_f_ghost, p_f_ghost_old, sim->nz+2);
 
 		for (iter=0; iter<max_iter; ++iter) {
 
@@ -189,29 +192,32 @@ darcy_solver_1d(struct simulation* sim,
 				       i, dp_f_dt_expl[i], i, dp_f_dt_impl[i]);
 #endif
 
-				p_f = sim->p_f_ghost[idx1g(i)];
+				p_f_old = sim->p_f_ghost[idx1g(i)];
 
-				p_f_ghost_out[idx1g(i)] = p_f
+				/* I keep adding dp/dt with every iteration
+				 * p_f_old is not the value from the previous time step,
+				 * but the value from the previous iteration. */
+				p_f_ghost_new[idx1g(i)] = p_f_old
 				                          + epsilon*dp_f_dt_impl[i]*sim->dt;
 
 				if (epsilon < 1.0)
-					p_f_ghost_out[idx1g(i)] += (1.0 - epsilon)
+					p_f_ghost_new[idx1g(i)] += (1.0 - epsilon)
 					                           *dp_f_dt_expl[i]*sim->dt;
 
-				/* apply relaxation */
-				p_f_ghost_out[idx1g(i)] = p_f*(1.0 - theta)
-				                          + p_f_ghost_out[idx1g(i)]*theta;
+				p_f_ghost_new[idx1g(i)] = p_f_old*(1.0 - theta)
+				                          + p_f_ghost_new[idx1g(i)]*theta;
 
-				r_norm[i] = (p_f_ghost_out[idx1g(i)] - p_f)/(p_f + 1e-16);
+				r_norm[i] = (p_f_ghost_new[idx1g(i)] - p_f_old)
+				            /(p_f_old + 1e-16);
 			}
 
 			r_norm_max = max(r_norm, sim->nz);
 #ifdef DEBUG
-			puts(".. p_f_ghost_out:");
-			print_array(p_f_ghost_out, sim->nz+2);
+			puts(".. p_f_ghost_new:");
+			print_array(p_f_ghost_new, sim->nz+2);
 #endif
 
-			copy_values(p_f_ghost_out, sim->p_f_ghost, sim->nz+2);
+			copy_values(p_f_ghost_new, sim->p_f_ghost, sim->nz+2);
 #ifdef DEBUG
 			puts(".. p_f_ghost after update:");
 			print_array(sim->p_f_ghost, sim->nz+2);
@@ -226,8 +232,9 @@ darcy_solver_1d(struct simulation* sim,
 			}
 		}
 		free(dp_f_dt_impl);
-		free(p_f_ghost_out);
+		free(p_f_ghost_new);
 		free(r_norm);
+		free(p_f_ghost_old);
 		if (!solved) {
 			fprintf(stderr, "darcy_solver_1d: ");
 			fprintf(stderr, "Solution did not converge after %d iterations\n",