WIP: Add porosity

examples/Bourdet_1983_1/Bourdet_example_1.ini

@@ -61,6 +61,7 @@ rp_k	= 1;
 rp_C	= 1;
 rp_pi	= 1;
 rp_re   = 0;
+rp_phi  = 0;
 
 
 ###########################################################################
@@ -75,6 +76,7 @@ jac_k	=  0;
 jac_C	=  0;
 jac_pi	=  0;
 jac_re	=  0;
+jac_phi =  2;
 
 
 ###########################################################################

examples/Bourdet_1983_1/Bourdet_example_1b.ini

@@ -60,6 +60,7 @@ rp_k	= 1;
 rp_C	= 1;
 rp_pi	= 1;
 rp_re   = 0;
+rp_phi  = 0;
 
 
 ###########################################################################
@@ -74,6 +75,7 @@ jac_k	=  0;
 jac_C	=  0;
 jac_pi	=  0;
 jac_re	=  2;
+jac_phi = 2;
 
 
 ###########################################################################

examples/Bourdet_1983_2/Bourdet_example_2.ini

@@ -94,6 +94,7 @@ rp_pi	= 1;
 rp_re	= 0;
 rp_L	= 0;
 rp_xf	= 0;
+rp_phi  = 0;
 
 ###############################################################################
 #  Derivatives:
@@ -110,6 +111,7 @@ jac_pi	=  0;
 jac_re	=  2;
 jac_L	=  2;
 jac_xf	=  2;
+jac_phi =  2;
 
 
 [Smoothing parameters]

examples/Bourdet_1983b/Bourdet_1983b.ini

@@ -91,6 +91,7 @@ rp_C	  = 1;
 rp_pi	  = 1;
 rp_omega  = 1;
 rp_lambda = 1;
+rp_phi    = 0;
 
 
 ###############################################################################
@@ -107,6 +108,7 @@ jac_C	   =  2;
 jac_pi	   =  2;
 jac_omega  =  2;
 jac_lambda =  2;
+jac_phi    =  2;
 
 
 [Smoothing parameters]

examples/Bourdet_1983b/Bourdet_1983b_a.ini

@@ -88,6 +88,7 @@ rp_k	  = 1;
 rp_C	  = 1;
 rp_pi	  = 1;
 rp_L	  = 1;
+rp_phi    = 0;
 
 
 ###############################################################################
@@ -103,6 +104,7 @@ jac_k	   =  0;
 jac_C	   =  0;
 jac_pi	   =  0;
 jac_L	   =  0;
+jac_phi    =  2;
 
 
 [Smoothing parameters]

examples/Bourdet_1983b/Bourdet_1983b_b.ini

@@ -85,6 +85,7 @@ rp_C	  = 1;
 rp_pi	  = 0;
 rp_omega  = 1;
 rp_lambda = 1;
+rp_phi    = 0;
 
 
 ###############################################################################
@@ -101,6 +102,7 @@ jac_C	   =  2;
 jac_pi	   =  2;
 jac_omega  =  2;
 jac_lambda =  2;
+jac_phi    =  2;
 
 
 [Smoothing parameters]

examples/Bourdet_1984/Bourdet_1984.ini

@@ -83,6 +83,7 @@ rp_C	  = 1;
 rp_pi	  = 1;
 rp_omega  = 1;
 rp_lambda = 1;
+rp_phi    = 0;
 
 
 ###############################################################################
@@ -99,6 +100,7 @@ jac_C	   =  0;
 jac_pi	   =  0;
 jac_omega  =  2;
 jac_lambda =  2;
+jac_phi    =  2;
 
 
 [Smoothing parameters]

examples/DaPrat_1984_BU/DaPrat_1984_BU.ini

@@ -83,6 +83,7 @@ rp_C	  = 1;
 rp_pi	  = 1;
 rp_omega  = 1;
 rp_lambda = 1;
+rp_phi    = 0;
 
 
 ###############################################################################
@@ -99,6 +100,7 @@ jac_C	   =  2;
 jac_pi	   =  2;
 jac_omega  =  2;
 jac_lambda =  2;
+jac_phi    =  2;
 
 
 [Smoothing parameters]

examples/DaPrat_1984_DD/DaPrat_1984_DD.ini

@@ -83,6 +83,7 @@ rp_C	  = 1;
 rp_pi	  = 0;
 rp_omega  = 1;
 rp_lambda = 1;
+rp_phi    = 0;
 
 
 ###############################################################################
@@ -99,6 +100,7 @@ jac_C	   =  2;
 jac_pi	   =  2;
 jac_omega  =  2;
 jac_lambda =  2;
+jac_phi    =  2;
 
 
 [Smoothing parameters]

examples/Gringarten_1979/Gringarten_example_1.ini

@@ -60,6 +60,7 @@ rp_k	= 1;
 rp_C	= 1;
 rp_pi	= 1;
 rp_re	= 0;
+rp_phi  = 0;
 
 
 ###########################################################################
@@ -74,6 +75,7 @@ jac_k	=  0;
 jac_C	=  0;
 jac_pi	=  0;
 jac_re	=  0;
+jac_phi  = 2;
 
 ###########################################################################
 #  Model:

examples/Lee_2003_2_12/Lee_example_2_12.ini

@@ -7,7 +7,7 @@
 #     ____________________________________________________________________
 #
 #	phi = 0,039 		p_i = 4412 psia
-#	�   = 0,8 cp 		c_t = 17�10^-6 psi^-1
+#	�   = 0,8 cp 		c_t = 17�10^-6 psi^-1
 #	r_w = 0,198 ft 		B 	= 1,136
 #	q   = 250 STB/d 	h 	= 69 ft
 #     ____________________________________________________________________
@@ -95,6 +95,7 @@ rp_k	= 1;
 rp_C	= 1;
 rp_pi	= 1;
 rp_re   = 1;
+rp_phi  = 0;
 
 
 ###########################################################################
@@ -109,6 +110,8 @@ jac_k	=  0;
 jac_C	=  0;
 jac_pi	=  0;
 jac_re	=  0;
+jac_phi  = 2;
+
 
 
 [Smoothing parameters]

examples/Lee_2003_2_12/Lee_example_2_12b.ini

@@ -7,7 +7,7 @@
 #     ____________________________________________________________________
 #
 #	phi = 0,039 		p_i = 4412 psia
-#	�   = 0,8 cp 		c_t = 17�10^-6 psi^-1
+#	�   = 0,8 cp 		c_t = 17�10^-6 psi^-1
 #	r_w = 0,198 ft 		B 	= 1,136
 #	q   = 250 STB/d 	h 	= 69 ft
 #     ____________________________________________________________________
@@ -100,6 +100,7 @@ rp_w1x  = 1;
 rp_w2x  = 1;
 rp_w1y  = 1;
 rp_w2y  = 1;
+rp_phi  = 0;
 
 
 ###########################################################################
@@ -117,6 +118,7 @@ jac_w1x	=  2;
 jac_w2x	=  2;
 jac_w1y	=  2;
 jac_w2y	=  2;
+jac_phi =  2;
 
 
 [Smoothing parameters]

examples/Lee_2003_2_12/Lee_example_2_12c.ini

@@ -7,7 +7,7 @@
 #     ____________________________________________________________________
 #
 #	phi = 0,039 		p_i = 4412 psia
-#	�   = 0,8 cp 		c_t = 17�10^-6 psi^-1
+#	�   = 0,8 cp 		c_t = 17�10^-6 psi^-1
 #	r_w = 0,198 ft 		B 	= 1,136
 #	q   = 250 STB/d 	h 	= 69 ft
 #     ____________________________________________________________________
@@ -100,6 +100,7 @@ rp_w1x  = 1;
 rp_w2x  = 1;
 rp_w1y  = 1;
 rp_w2y  = 1;
+rp_phi  = 0;
 
 
 ###########################################################################
@@ -117,7 +118,7 @@ jac_w1x	=  2;
 jac_w2x	=  2;
 jac_w1y	=  2;
 jac_w2y	=  2;
-
+jac_phi =  2;
 
 [Smoothing parameters]
 Lder	= 0.1;

examples/Samad_1994_7/Samad_1994_7.ini

@@ -81,6 +81,7 @@ rp_k	  = 1;
 rp_C	  = 1;
 rp_pi	  = 0;
 rp_xf     = 1;
+rp_phi    = 0;
 
 
 ###############################################################################
@@ -96,6 +97,7 @@ jac_k	   =  2;
 jac_C	   =  2;
 jac_pi	   =  2;
 jac_xf     =  2;
+jac_phi    =  2;
 
 dx_xf      =  1.0e-3;
 

examples/Samad_1994_9/Samad_1994_9.ini

@@ -83,6 +83,7 @@ rp_C	  = 0;
 rp_pi	  = 0;
 rp_xf     = 1;
 rp_fc     = 1;
+rp_phi    = 0;
 
 
 ###############################################################################
@@ -99,6 +100,7 @@ jac_C	   =  2;
 jac_pi	   =  2;
 jac_xf     =  2;
 jac_fc     =  2;
+jac_phi    =  2;
 
 
 [Smoothing parameters]

src/captas.c

@@ -1,6 +1,6 @@
 /*
 
- Computer Aided Pressure Transient Analysis (and Simulation?) � CAPTA(S?)
+ Computer Aided Pressure Transient Analysis (and Simulation?) � CAPTA(S?)
      Copyright (C) 2013  Carlos E. Pico
 
     This program is free software: you can redistribute it and/or modify
@@ -534,7 +534,9 @@ void init_parameters(modelparameters *par)
 	par->parnames[PENETRATION_RATIO]       = "b";
 	par->parnames[MIDPOINT_ELEVATION]      = "zw",
 	par->parnames[EFFECTIVE_HEAT_CAPACITY] = "cpt",
-	par->parnames[INITIAL_TEMPERATURE]     = "Ti";
+	par->parnames[INITIAL_TEMPERATURE]     = "Ti",
+	par->parnames[COEFFICIENT_JOULE_THOMSON]  = "ejt",
+	par->parnames[POROSITY]                = "phi";
 
 	/******************** pointers to delta_pwf functions ********************/
 	for(i = 0; i <  NPMODELS; i++){
@@ -677,13 +679,11 @@ void read_inifile(modelparameters *par)
 	par->units    = iniparser_getint(ini, "Units system:units",        OILFIELD);
 	par->testtype = iniparser_getint(ini, "Test description:testtype", DRAWDOWN);
 
-	par->phi     = iniparser_getdouble(ini, "Test description:phi",    0.0);
 	par->B       = iniparser_getdouble(ini, "Test description:B",      0.0);
 	par->mu      = iniparser_getdouble(ini, "Test description:mu",     0.0);
 	par->h       = iniparser_getdouble(ini, "Test description:h",      0.0);
 	par->rw      = iniparser_getdouble(ini, "Test description:rw",     0.0);
 	par->ct      = iniparser_getdouble(ini, "Test description:ct",     0.0);
-	par->ejt     = iniparser_getdouble(ini, "Test description:ejt",    0.0);
 	par->rhosc   = iniparser_getdouble(ini, "Test description:rhosc",  0.0);
 	par->cp      = iniparser_getdouble(ini, "Test description:cp",     0.0);
 

src/models/dTwf.c

@@ -16,7 +16,7 @@
 * Transient temperature behavior and analysis of single-phase liquid-water
 * geothermal reservoirs during drawdown and buildup tests: Part I.
 * Theory, new analytical and approximate solutions.
-* Journal of Petroleum Science and Engineering, 146, 637�656.
+* Journal of Petroleum Science and Engineering, 146, 637�656.
 * https://doi.org/10.1016/j.petrol.2016.08.003
 */
 
@@ -28,21 +28,23 @@
 double dTwf(const modelparameters *p, double t)
 {
     double a, b, cpR, f, d;
-    double k, S, cpt;
+    double k, S, cpt, ejt;
 
     gsl_set_error_handler_off();
 
     k   = p->rpval[PERMEABILITY];
     S   = p->rpval[SKIN_FACTOR];
     cpt = p->rpval[EFFECTIVE_HEAT_CAPACITY];
+    ejt = p->rpval[COEFFICIENT_JOULE_THOMSON];
+    phi = p->rpval[POROSITY];
 
     a   = (p->C2 * p->qB * p->mu) / (k * p->h);
-    b   = (p->rw * p->rw * p->phi * p->ct * p->mu ) / (4.0 * k * p->C1);
+    b   = (p->rw * p->rw * phi * p->ct * p->mu ) / (4.0 * k * p->C1);
     cpR = (p->rhosc / p->B) * p->cp / cpt;
-    f   = p->phi * cpR * (p->ejt + 1.0 / (p->rhosc * p->cp / p->B)) - p->ejt;
-    d   = (p->phi * cpR * p->ct)  * a * 0.5;
+    f   = phi * cpR * (ejt + 1.0 / (p->rhosc * p->cp / p->B)) - ejt;
+    d   = (phi * cpR * p->ct)  * a * 0.5;
 
-    return -0.5*a*( -p->ejt*(gsl_sf_expint_E1(b/t) + 2.0*S) - f*gsl_sf_expint_E1(b/t + d) );
+    return -0.5*a*( -ejt*(gsl_sf_expint_E1(b/t) + 2.0*S) - f*gsl_sf_expint_E1(b/t + d) );
 
 }
 /*****************************************************************************/

src/models/dpwf.c

@@ -38,11 +38,12 @@ double dpwfbar(const void *parameters, double u)
 	k   = p->rpval[PERMEABILITY];
 	S   = p->rpval[SKIN_FACTOR];
 	C   = p->rpval[WELLBORE_STORAGE];
+	phi = p->rpval[POROSITY];
 
 	rws = p->rw*exp(-S);
 	a   = (p->qB * p->mu * p->C2)/ (p->h * k);
-	b   = (p->phi * p->mu * p->ct * rws * rws) / (k * p->C1);
-	CD  = (C * p->C3) / (p->phi * p->h  * p->ct * rws * rws);
+	b   = (phi * p->mu * p->ct * rws * rws) / (k * p->C1);
+	CD  = (C * p->C3) / (phi * p->h  * p->ct * rws * rws);
 	uD  = (u * b);
 
 	uDsqrt = sqrt(uD);
@@ -70,12 +71,13 @@ double ddpwf_dkbar(const void *parameters, double u)
 	k   = p->rpval[PERMEABILITY];
 	S   = p->rpval[SKIN_FACTOR];
 	C   = p->rpval[WELLBORE_STORAGE];
+	phi = p->rpval[POROSITY];
 
 	rws = p->rw*exp(-S);
 	a   = (p->qB * p->mu * p->C2) / (p->h * k);
-	y   = (p->phi * p->mu * p->ct * rws * rws)/(p->C1);
+	y   = (phi * p->mu * p->ct * rws * rws)/(p->C1);
 	b   = y / k;
-	CD  = (C * p->C3) / (p->phi * p->h  * p->ct * rws * rws);
+	CD  = (C * p->C3) / (phi * p->h  * p->ct * rws * rws);
 	uD  = (u * b);
 
 	uDsqrt = sqrt(uD);
@@ -104,12 +106,13 @@ double ddpwf_dSbar(const void *parameters, double u)
 	k   = p->rpval[PERMEABILITY];
 	S   = p->rpval[SKIN_FACTOR];
 	C   = p->rpval[WELLBORE_STORAGE];
+	phi = p->rpval[POROSITY];
 
 	rws = p->rw*exp(-S);
 	a   = (p->qB * p->mu * p->C2) / (p->h * k);
-	y   = (p->phi * p->mu * p->ct * p->rw * p->rw) / (k * p->C1);
-	b   = (p->phi * p->mu * p->ct *   rws *   rws) / (k * p->C1);
-	z   = (C * p->C3) / (p->phi * p->h  * p->ct * p->rw * p->rw);
+	y   = (phi * p->mu * p->ct * p->rw * p->rw) / (k * p->C1);
+	b   = (phi * p->mu * p->ct *   rws *   rws) / (k * p->C1);
+	z   = (C * p->C3) / (phi * p->h  * p->ct * p->rw * p->rw);
 	uD  = (u * b);
 
 	uDsqrt = sqrt(uD);