Uploading all the required files

Class_problems/Additional1_petrophysics_differentials/differentials/assignment2_finitediffapprox.py

@@ -0,0 +1,195 @@
+"""
+reservoir simulation assignment 2
+Pressure values and derivative approximation
+Author: Mohammad Afzal Shadab
+Email: [email protected]
+Date modified: 9/22/2020
+"""
+
+#Importing required libraries
+from IPython import get_ipython
+get_ipython().magic('reset -sf')#for clearing everything
+get_ipython().run_line_magic('matplotlib', 'qt') #for plotting in separate window
+# Import curve fitting package from scipy
+from scipy.optimize import curve_fit
+
+import numpy as np #import numpy library
+import matplotlib.pyplot as plt #library for plotting
+plt.rcParams.update({'font.size': 22})
+
+# Function to calculate the power-law with constants a and b
+def power_law(x, a, b):
+ return a*np.power(x, b)
+
+#Making grid class
+class grid:
+ def __init__(self):
+ self.xmin = []
+ self.xmax = []
+ self.Nx = []
+ self.N = []
+
+#Input parameters
+pi = 1000 # initial reservior pressure [psi] 
+qw = 100 # well flow rate [ft^3/day] 
+k = 50 # permeability [Darcy]
+h = 20 # thickness of the reservior [feet]
+alpha = 10.0 # diffusivity contant [ft^2/day] 
+rw = 0.25 # well radius [feet] 
+mu = 1.0 # fluid viscosity [centipoise] 
+conv_fac = 1.0/(6.33E-3) #conversion factor for mD per Darcy to ft^2/psi-day
+tf = 1 # final time [day]
+dt_a = 0.01 # time step for approximation part a [day]
+dr_b = 0.001 # radial step for approximation part b [feet]
+dr_c = 0.001 # radial step for approximation part c [feet]
+
+#Calculate the constant beta
+beta = qw * mu / (4 * np.pi* k * h) *conv_fac 
+
+#Pressure as anonymous function
+P = lambda r,t: pi - beta*(np.log(r**2.0/(alpha*t)) + 0.57721)
+
+#Evaluating the analytical values
+dPbydt_analy = lambda t: beta / t
+dPbydr_analy = lambda r: -2.0 * beta/r
+d2Pbydr2_analy = lambda r: 2.0 * beta/r**2.0
+
+#PART B =========================================
+#subpart (a)
+dPbydt_frwd = (P(rw,tf+dt_a)-P(rw,tf))/dt_a #forward difference
+dPbydt_bkwrd = (P(rw,tf)-P(rw,tf-dt_a))/dt_a #backward difference
+dPbydt_cntr = (P(rw,tf+dt_a)-P(rw,tf-dt_a))/(2*dt_a) #centered difference
+
+#subpart (b)
+dPbydr_frwd = (P(rw+dr_b,tf)-P(rw,tf))/dr_b #forward difference
+dPbydr_bkwrd = (P(rw,tf)-P(rw-dr_b,tf))/dr_b #backward difference
+dPbydr_cntr = (P(rw+dr_b,tf)-P(rw-dr_b,tf))/(2*dr_b) #centered difference
+
+#subpart (c)
+d2Pbydr2_cntr = ( P(rw+dr_c,tf) -2.0 * P(rw,tf) + P(rw-dr_c,tf) ) / dr_c**2 #centered difference
+
+
+#PART C =============================================
+dr = np.logspace(-3,-1,100)
+dt = np.logspace(-4,-2,100)
+
+#subpart (a)
+VdPbydt_frwd = (P(rw,tf+dt)-P(rw,tf))/dt #forward difference V stands for vectorized
+VdPbydt_bkwrd = (P(rw,tf)-P(rw,tf-dt))/dt #backward difference
+VdPbydt_cntr = (P(rw,tf+dt)-P(rw,tf-dt))/(2*dt) #centered difference
+#Absolute errors
+EdPbydt_frwd = np.abs(VdPbydt_frwd - dPbydt_analy(tf)) #errors in forward difference E stands for vectorized errors
+EdPbydt_bkwrd = np.abs(VdPbydt_bkwrd - dPbydt_analy(tf)) #errors in backward difference E stands for vectorized errors
+EdPbydt_cntr = np.abs(VdPbydt_cntr - dPbydt_analy(tf)) #errors in centered difference E stands for vectorized errors
+#Power law fit to find the order of accuracy
+ZdPbydt_frwd, cov = curve_fit(f=power_law, xdata=dt, ydata=EdPbydt_frwd, bounds=(-np.inf, np.inf))
+ZdPbydt_bkwrd, cov = curve_fit(f=power_law, xdata=dt, ydata=EdPbydt_bkwrd, bounds=(-np.inf, np.inf))
+ZdPbydt_cntr, cov = curve_fit(f=power_law, xdata=dt, ydata=EdPbydt_cntr, bounds=(-np.inf, np.inf))
+
+#subpart (b)
+VdPbydr_frwd = (P(rw+dr,tf)-P(rw,tf))/dr #forward difference
+VdPbydr_bkwrd = (P(rw,tf)-P(rw-dr,tf))/dr #backward difference
+VdPbydr_cntr = (P(rw+dr,tf)-P(rw-dr,tf))/(2*dr) #centered difference
+#Absolute errors
+EdPbydr_frwd = np.abs(VdPbydr_frwd - dPbydr_analy(rw)) #errors in forward difference E stands for vectorized errors
+EdPbydr_bkwrd = np.abs(VdPbydr_bkwrd - dPbydr_analy(rw)) #errors in backward difference E stands for vectorized errors
+EdPbydr_cntr = np.abs(VdPbydr_cntr - dPbydr_analy(rw)) #errors in centered difference E stands for vectorized errors
+
+#Power law fit to find the order of accuracy
+ZdPbydr_frwd, cov = curve_fit(f=power_law, xdata=dr, ydata=EdPbydr_frwd, bounds=(-np.inf, np.inf))
+ZdPbydr_bkwrd, cov = curve_fit(f=power_law, xdata=dr, ydata=EdPbydr_bkwrd, bounds=(-np.inf, np.inf))
+ZdPbydr_cntr, cov = curve_fit(f=power_law, xdata=dr, ydata=EdPbydr_cntr, bounds=(-np.inf, np.inf))
+
+#subpart (c)
+Vd2Pbydr2_cntr = (P(rw+dr,tf) -2.0 * P(rw,tf) + P(rw-dr,tf) ) / dr**2 #centered difference
+#Absolute errors
+Ed2Pbydr2_cntr = np.abs(Vd2Pbydr2_cntr - d2Pbydr2_analy(rw)) #errors in centered difference E stands for vectorized errors
+#Polyfit to find the order of accuracy
+Zd2Pbydr2_cntr, cov = curve_fit(f=power_law, xdata=dr, ydata=Ed2Pbydr2_cntr, bounds=(-np.inf, np.inf))
+
+#PLOTTING
+#Higher order derivative approximations ===========================================================
+#(a) Plotting the forward, backward and center difference for dP/dt
+fig = plt.figure(figsize=(15,7.5) , dpi=100)
+plot = plt.plot(dt,VdPbydt_frwd,'ro',label=r'$Forward$')
+plot = plt.plot(dt,VdPbydt_bkwrd,'bd',label=r'$Backward$')
+plot = plt.plot(dt,VdPbydt_cntr,'gs',label=r'$Centred$')
+plt.hlines(dPbydt_analy(tf),dt[1],dt[99],label=r'$Analytical$')
+manager = plt.get_current_fig_manager()
+manager.window.showMaximized()
+plt.ylabel(r'$dP/dt$ [psi/day]')
+plt.xlabel(r'$\Delta t$ [day]')
+legend = plt.legend(loc='best', shadow=False, fontsize='x-large')
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.savefig(f'dpbydt.png',bbox_inches='tight', dpi = 600)
+
+#(b) Plotting the forward, backward and center difference for dP/dr
+fig = plt.figure(figsize=(15,7.5) , dpi=100)
+plot = plt.plot(dr,VdPbydr_frwd,'ro',label=r'$Forward$')
+plot = plt.plot(dr,VdPbydr_bkwrd,'bd',label=r'$Backward$')
+plot = plt.plot(dr,VdPbydr_cntr,'gs',label=r'$Centred$')
+plt.hlines(dPbydr_analy(rw),dr[1],dr[99],label=r'$Analytical$')
+manager = plt.get_current_fig_manager()
+manager.window.showMaximized()
+plt.ylabel(r'$dP/dr$ [psi/feet]')
+plt.xlabel(r'$\Delta r$ [feet]')
+legend = plt.legend(loc='best', shadow=False, fontsize='x-large')
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.savefig(f'dpbydr.png',bbox_inches='tight', dpi = 600)
+
+#(b) Plotting the center difference for d^2P/dr^2
+fig = plt.figure(figsize=(15,7.5) , dpi=100)
+plot = plt.plot(dr,Vd2Pbydr2_cntr,'ro',label=r'$Centred$')
+plt.hlines(d2Pbydr2_analy(rw),dr[1],dr[99],label=r'$Analytical$')
+manager = plt.get_current_fig_manager()
+manager.window.showMaximized()
+plt.ylabel(r'$d^2P/dr^2$ [psi/feet^2]')
+plt.xlabel(r'$\Delta r$ [feet]')
+legend = plt.legend(loc='best', shadow=False, fontsize='x-large')
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.savefig(f'd2pbydr2.png',bbox_inches='tight', dpi = 600)
+
+#ABSOLUTE ERRORS===========================================================
+#(a) Plotting the absolute errors in forward, backward and center difference for dP/dt
+ax = plt.figure(figsize=(15,7.5) , dpi=100)
+plot = plt.loglog(dt,EdPbydt_frwd,'ro',label=r'$Forward$')
+plot = plt.loglog(dt,power_law(dt, *ZdPbydt_frwd),'r--',label=r'$fit: O( \Delta t^{%5.3f})$' %ZdPbydt_frwd[1] )
+plot = plt.loglog(dt,EdPbydt_bkwrd,'bd',label=r'$Backward$')
+plot = plt.loglog(dt,power_law(dt, *ZdPbydt_bkwrd),'b--',label=r'$fit: O(\Delta t^{%5.3f})$' %ZdPbydt_bkwrd[1] )
+plot = plt.loglog(dt,EdPbydt_cntr,'gs',label=r'$Centred$')
+plot = plt.loglog(dt,power_law(dt, *ZdPbydt_cntr),'g--',label=r'$fit: O(\Delta t^{%5.3f})$' %ZdPbydt_cntr[1] )
+manager = plt.get_current_fig_manager()
+manager.window.showMaximized()
+plt.ylabel(r'$|dP/dt-dP/dt_{analytical}|$ [psi/day]')
+plt.xlabel(r'$\Delta t$ [day]')
+legend = plt.legend(loc='best', shadow=False, fontsize='x-large')
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.savefig(f'AbsErr_dpbydt.png',bbox_inches='tight', dpi = 600)
+
+#(b) Plotting the absolute errors in forward, backward and center difference for dP/dr
+fig = plt.figure(figsize=(15,7.5) , dpi=100)
+plot = plt.loglog(dr,EdPbydr_frwd,'ro',label=r'$Forward$')
+plot = plt.loglog(dr,power_law(dr, *ZdPbydr_frwd),'r--',label=r'$fit: O( \Delta r^{%5.3f})$' %ZdPbydr_frwd[1] )
+plot = plt.loglog(dr,EdPbydr_bkwrd,'bd',label=r'$Backward$')
+plot = plt.loglog(dr,power_law(dr, *ZdPbydr_bkwrd),'b--',label=r'$fit: O( \Delta r^{%5.3f})$' %ZdPbydr_bkwrd[1] )
+plot = plt.loglog(dr,EdPbydr_cntr,'gs',label=r'$Centred$')
+plot = plt.loglog(dr,power_law(dr, *ZdPbydr_cntr),'g--',label=r'$fit: O(\Delta r^{%5.3f})$' %ZdPbydr_cntr[1] )
+manager = plt.get_current_fig_manager()
+manager.window.showMaximized()
+plt.xlabel(r'$\Delta r$ [feet]')
+plt.ylabel(r'$|dP/dr-dP/dr_{analytical}|$ [psi/feet]')
+legend = plt.legend(loc='best', shadow=False, fontsize='x-large')
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.savefig(f'AbsErr_dpbydr.png',bbox_inches='tight', dpi = 600)
+
+#(c) Plotting the absolute errors in center difference for d^2P/dr^2
+fig = plt.figure(figsize=(15,7.5) , dpi=100)
+plot = plt.loglog(dr,Ed2Pbydr2_cntr,'gs',label=r'$Centred$')
+plot = plt.loglog(dr,power_law(dr, *Zd2Pbydr2_cntr),'g--',label=r'$fit: O(\Delta r^{%5.3f})$' %Zd2Pbydr2_cntr[1] )
+manager = plt.get_current_fig_manager()
+manager.window.showMaximized()
+plt.ylabel(r'$|d^2P/dr^2-d^2P/dr^2_{analytical}|$ [psi/feet^2]')
+plt.xlabel(r'$\Delta r$ [feet]')
+legend = plt.legend(loc='best', shadow=False, fontsize='x-large')
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.savefig(f'AbsErr_d2pbydr2.png',bbox_inches='tight', dpi = 600)

Class_problems/Problem_0_1D_reservoir/Analytical python.txt

@@ -0,0 +1,19 @@
+def analytical(time,blocks):
+ tol = 0.0001
+ x = np.arange(blocks)*dx + dx/2
+ P1 = 500-((D-Length*np.sin(dip*np.pi/180))*(rho/144))- 14.7; #based on boundary 500 psi at L P1 at pooint n check
+ P0 = 0
+ alpha = 6.33E-3*k/(phi*u*ct)
+ n = 0
+ SUM = 0
+ for i in range(10000000):
+ new = (-1)**(n+1)*np.cos((n+1/2)*np.pi*x/Length)*np.exp(-(2*n+1)**2/Length**2*alpha*np.pi**2*time/4)/(2*n+1)
+ SUM = 4/np.pi*(P1-P0)*new + SUM
+ error=max(abs(new));
+ if error > tol:
+ n = n +1
+ else:
+ break
+
+ P = SUM + p + P1
+ return P;

Class_problems/Problem_0_1D_reservoir/Analytical_with_gravity.m

@@ -0,0 +1,36 @@
+function [x,P] = Analytical_with_gravity (t,N,theta)
+%input t:time (days),Scalar
+%input N: number of grid blocks
+%input theta: reservoir dip angle above the horizon in radians (positive..
+%..for dipping up, negative for dipping down)
+%output x: reservoir x-location (ft) 1*N vector 
+%output P: pressure (psi) 1*N vector 
+tol=eps;
+rho_w=62.4;
+L=4000; %Reservoir Length (ft)
+za=2309; %Depth of left boundry of the reservoir (ft)
+grid_vec=((1:N));
+dx=L/N;
+x=[(dx/2)*(2*grid_vec-1)];
+phi = 0.25;
+k = 15;
+mu = 1;
+cf = 10^-5;
+dpg1=((za-L*sin(theta))*(rho_w/144))+14.6959502543;
+P1 = 500-dpg1 ;
+P0=0;
+alpha = 6.33*10^-3*k/(phi*mu*cf);
+n = 0;
+SUM = 0;
+for i=1:10000000
+ new = (-1)^(n+1)*cos((n+1/2)*pi*x/L)*exp(-(2*n+1)^2/L^2*alpha*pi^2*t/4)/(2*n+1);
+ SUM = 4/pi*(P1-P0).*new + SUM;
+ error=max(abs(new));
+ if error > tol
+ n = n+1;
+ else
+ break
+ end
+end
+dpg=((za-x*sin(theta))*(rho_w/144))+14.6959502543;
+P= SUM +P1+dpg;

Class_problems/Problem_0_1D_reservoir/analytical_python_new.py

@@ -0,0 +1,31 @@
+import numpy as np
+
+def analytical(time,blocks):
+ tol = 0.00000001
+ rho = 62.4
+ Length = 10000
+ za = 1000/(rho/144.0)
+ dip = 0.0*np.pi/6.0
+ k = 50.0
+ phi = 0.2
+ ct = 1E-6
+ mu = 1.0
+ dx = Length / blocks
+ x = np.arange(blocks)*dx + dx/2
+ P1 = 2000-((za-Length*np.sin(dip))*(rho/144)) - 14.6959502543; #based on boundary 500 psi at L P1 at pooint n check
+ P0 = 0
+ alpha = 6.33E-3*k/(phi*mu*ct)
+ n = 0
+ SUM = 0
+ for i in range(10000000):
+ new = (-1)**(n+1)*np.cos((n+1/2)*np.pi*x/Length)*np.exp(-(2*n+1)**2/Length**2*alpha*np.pi**2*time/4)/(2*n+1)
+ SUM = 4/np.pi*(P1-P0)*new + SUM
+ error=max(abs(new));
+ if error > tol:
+ n = n + 1
+ else:
+ break
+ dpg=((za-x*np.sin(dip))*(rho/144.0)) + 14.6959502543
+ P = SUM + dpg + P1
+ P = np.flip(P)
+ return P, x

Class_problems/Problem_0_1D_reservoir/input_file.py

@@ -0,0 +1,64 @@
+"""
+reservoir simulation assignment 2
+1D reservoir simulation Q3.4: Input files
+Author: Mohammad Afzal Shadab
+Email: [email protected]
+Date modified: 9/24/2020
+"""
+import numpy as np
+
+class fluid:
+ def __init__(self):
+ self.mu = []
+class numerical:
+ def __init__(self):
+ self.Bw = []
+class reservoir:
+ def __init__(self):
+ self.dt = [] 
+class grid:
+ def __init__(self):
+ self.xmin = []
+class BC:
+ def __init__(self):
+ self.xmin = []
+class IC:
+ def __init__(self):
+ self.xmin = []
+
+#fluid, reservoir and simulation parameters 
+def inputfile(fluid,reservoir,numerical,BC,IC):
+ fluid.mu = 1.0 #fluid viscosity [centipoise] 
+ fluid.Bw = 1.0 #formation volume factor of water [rb/stb] 
+ fluid.ct = 1E-6 #total compressibility [1/psi] 
+ fluid.rho = 62.4 #density of the fluid [lbm/ft^3]
+
+ numerical.dt = 1.0 #time step (days)
+ numerical.tfinal = 200 #final time [days]
+ numerical.N = 40000 #number of grid blocks
+ numerical.theta = 0 #type of method theta = 1-Explicit, 0-Implicit, 0.5-CN
+
+ reservoir.L = 10000.0 #length of the reservoir [feet] 
+ reservoir.A = 2E5 #cross sectinal area of the reservoir [feet^2] 
+ reservoir.phi= 0.2*np.ones((numerical.N, 1)) #porosity of the reservior vector
+ reservoir.k = 50.0*np.ones((numerical.N, 1)) #fluid permeability vector [mDarcy] 
+ reservoir.Dref=1000/(fluid.rho/144.0) #reference depth [feet]
+ reservoir.alpha=0.0*np.pi/6.0#dip angle [in radians]
+
+ numerical.dx = (reservoir.L/numerical.N)*np.ones((numerical.N, 1)) #block thickness vector
+ numerical.xc = np.empty((numerical.N, 1))
+ numerical.xc[0,0] = 0.5 * numerical.dx[0,0]
+
+ #position of the block centres
+ for i in range(0,numerical.N):
+ numerical.xc[i,0] = numerical.xc[i-1,0] + 0.5*(numerical.dx[i-1,0] + numerical.dx[i,0])
+
+ #depth vector 
+ numerical.D = reservoir.Dref-numerical.xc*np.sin(reservoir.alpha)
+
+ BC.type = [['Dirichlet'],['Neumann']] #type of BC as first cell, last cell
+ BC.value = [[2000],[0]] #value of the boundary condition: psi or ft^3/day
+ BC.depth = [[reservoir.Dref],[reservoir.Dref-reservoir.L*np.sin(reservoir.alpha)]]
+
+ #IC.P = 1000.0*np.ones((numerical.N, 1)) #initial condition: pressure in the domain
+ IC.P = 1000.0 + (fluid.rho/144.0)*(numerical.D - reservoir.Dref)