cngf-pf

commit 4f4e70d16f79d863d20480649887c9c1431cf783
parent 8b7ade10afeebea5e616ab0a0f3efffec2116279
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Mon,  6 Apr 2020 12:20:17 +0200

Fix calculation of transient porosity, add porosity forcing to fluid equation

Diffstat:
M
fluid.c
 | 
7
++++++-
M
simulation.c
 | 
6
++++--

2 files changed, 10 insertions(+), 3 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -102,6 +102,7 @@ darcy_pressure_change_1d(const int i,
                          const int nz,
                          const double* p_f_ghost_in,
                          const double* phi,
+                         const double* phi_dot,
                          const double* k,
                          const double dz,
                          const double beta_f,
@@ -132,7 +133,9 @@ darcy_pressure_change_1d(const int i,
 	        i, phi[i], i, div_k_grad_p);
 #endif
 
-	return 1.0/(beta_f*phi[i]*mu_f)*div_k_grad_p; 
+	/* TODO: add advective term */
+	return 1.0/(beta_f*phi[i]*mu_f)*div_k_grad_p
+		   - 1.0/(beta_f*(1.0 - phi[i]))*phi_dot[i];
 }
 
 int
@@ -192,6 +195,7 @@ darcy_solver_1d(const int max_iter,
 			                                           sim.nz,
 			                                           sim.p_f_ghost,
 			                                           sim.phi,
+			                                           sim.phi_dot,
 			                                           sim.k,
 			                                           sim.dz,
 			                                           sim.beta_f,
@@ -220,6 +224,7 @@ darcy_solver_1d(const int max_iter,
 				                                           sim.nz,
 				                                           sim.p_f_ghost,
 				                                           sim.phi,
+				                                           sim.phi_dot,
 				                                           sim.k,
 				                                           sim.dz,
 				                                           sim.beta_f,
diff --git a/simulation.c b/simulation.c
@@ -322,8 +322,10 @@ void
 compute_porosity_change()
 {
 	int i;
-	for (i=0; i<sim.nz; ++i)
-		sim.phi[i] += sim.tan_psi[i]*sim.gamma_dot_p[i]*sim.phi[i]*sim.dt;
+	for (i=0; i<sim.nz; ++i) {
+		sim.phi_dot[i] = sim.tan_psi[i]*sim.gamma_dot_p[i];
+		sim.phi[i] += sim.phi_dot[i]*sim.dt;
+	}
 }
 
 void