cngf-pf

commit df8b703a4833a956612cd7e1710087e7c2d844ba
parent 1d27bea446f7b34d61d3126d701dc9367a1b0ea6
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Mon,  1 Jul 2019 15:19:49 +0200

Allow any gradient with von Neumann BCs, fix diurnal example

Diffstat:
M
fluid.c
 | 
10
+++++++---
M
simulation.c
 | 
12
++++++++----
M
simulation.h
 | 
9
+++++----

3 files changed, 20 insertions(+), 11 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -65,7 +65,11 @@ set_fluid_bcs(struct simulation* sim, const double p_f_top)
 {
 	set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top);
 	sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; /* Include top node in BC */
-	set_bc_neumann(sim->p_f_ghost, sim->nz, -1);
+	set_bc_neumann(sim->p_f_ghost,
+	               sim->nz,
+	               -1,
+	               sim->phi[0]*sim->rho_f*sim->G,
+	               sim->z[1] - sim->z[0]);
 }
 
 static double
@@ -121,7 +125,7 @@ darcy_solver_1d(struct simulation* sim,
 	 *     epsilon = 1.0: implicit */
 	/* epsilon = 0.5; */
 	/* epsilon = 0.0; */
-	epsilon = 1.0;
+	epsilon = 0.0;
 
 	/* choose relaxation factor, parameter in ]0.0; 1.0]
 	 *     theta in ]0.0; 1.0]: underrelaxation
@@ -156,7 +160,7 @@ darcy_solver_1d(struct simulation* sim,
 	}
 
 	if (epsilon > 0.0) {
-		/* compute implicit solution to pressure change */
+		/* compute implicit solution with Jacobian iterations */
 		dp_f_dt_impl = zeros(sim->nz);
 		p_f_ghost_out = zeros(sim->nz+2);
 		r_norm = zeros(sim->nz);
diff --git a/simulation.c b/simulation.c
@@ -307,12 +307,16 @@ compute_local_fluidity(struct simulation* sim)
 }
 
 void
-set_bc_neumann(double* g_ghost, const int nz, const int boundary)
+set_bc_neumann(double* g_ghost,
+               const int nz,
+               const int boundary,
+               const double df,
+               const double dx)
 {
 	if (boundary == -1)
-		g_ghost[0] = g_ghost[1];
+		g_ghost[0] = g_ghost[1] + df*dx;
 	else if (boundary == +1)
-		g_ghost[nz+1] = g_ghost[nz];
+		g_ghost[nz+1] = g_ghost[nz] - df*dx;
 	else {
 		fprintf(stderr, "set_bc_neumann: Unknown boundary %d\n", boundary);
 		exit(1);
@@ -393,7 +397,7 @@ implicit_1d_jacobian_poisson_solver(struct simulation* sim,
 		set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0);
 
 		/* Neumann BCs resemble free surfaces */
-		/* set_bc_neumann(sim->g_ghost, sim->nz, +1); */
+		/* set_bc_neumann(sim->g_ghost, sim->nz, +1, 0.0, 0.0); */
 
 		for (i=0; i<sim->nz; ++i)
 			poisson_solver_1d_cell_update(i,
diff --git a/simulation.h b/simulation.h
@@ -100,10 +100,11 @@ void check_simulation_parameters(const struct simulation* sim);
 
 void lithostatic_pressure_distribution(struct simulation* sim);
 
-void set_bc_neumann(
-        double* g_ghost,
-        const int nz,
-        const int boundary);
+void set_bc_neumann(double* g_ghost,
+                    const int nz,
+                    const int boundary,
+                    const double df,
+                    const double dx);
 
 void set_bc_dirichlet(
         double* g_ghost,