1d_fd_simple_shear_transient

commit 6c48f55a2b5bc7078851062650f5e4e31f90d471
parent 599c08baed18df47fd5a6fac8b587e5d51553519
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Mon,  1 Jul 2019 21:42:16 +0200

Fix Jacobi method in iterative solver

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

1 file changed, 7 insertions(+), 12 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -124,9 +124,7 @@ darcy_solver_1d(struct simulation* sim,
 	 *     epsilon = 0.0: explicit
 	 *     epsilon = 0.5: Crank-Nicolson
 	 *     epsilon = 1.0: implicit */
-	/* epsilon = 0.5; */
-	/* epsilon = 0.0; */
-	epsilon = 1.0;
+	epsilon = 0.5;
 
 	/* choose relaxation factor, parameter in ]0.0; 1.0]
 	 *     theta in ]0.0; 1.0]: underrelaxation
@@ -192,23 +190,19 @@ darcy_solver_1d(struct simulation* sim,
 				       i, dp_f_dt_expl[i], i, dp_f_dt_impl[i]);
 #endif
 
-				p_f_old = sim->p_f_ghost[idx1g(i)];
-
-				/* 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
+				p_f_ghost_new[idx1g(i)] = p_f_ghost_old[idx1g(i)]
 				                          + epsilon*dp_f_dt_impl[i]*sim->dt;
 
 				if (epsilon < 1.0)
 					p_f_ghost_new[idx1g(i)] += (1.0 - epsilon)
 					                           *dp_f_dt_expl[i]*sim->dt;
 
-				p_f_ghost_new[idx1g(i)] = p_f_old*(1.0 - theta)
+				p_f_ghost_new[idx1g(i)] = p_f_ghost_old[idx1g(i)]*(1.0 - theta)
 				                          + p_f_ghost_new[idx1g(i)]*theta;
 
-				r_norm[i] = (p_f_ghost_new[idx1g(i)] - p_f_old)
-				            /(p_f_old + 1e-16);
+				r_norm[i] = (p_f_ghost_new[idx1g(i)]
+				            - sim->p_f_ghost[idx1g(i)])
+				            /(sim->p_f_ghost[idx1g(i)] + 1e-16);
 			}
 
 			r_norm_max = max(r_norm, sim->nz);
@@ -252,6 +246,7 @@ darcy_solver_1d(struct simulation* sim,
 		print_array(sim->p_f_ghost, sim->nz+2);
 #endif
 	}
+	set_fluid_bcs(sim, p_f_top);
 
 	if (epsilon < 1.0)
 		free(dp_f_dt_expl);