<jupyter_start><jupyter_text>1D Reservoir Simulation using PyReSim<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
!git clone https://github.com/yohanesnuwara/pyresim
import sys
sys.path.append('/content/pyresim/functions')
from input_output import read_input, read_depth
from boundary import boundary_floweq1d
from gridding import source1d
from solver import lhs_coeffs1d_welltype, rhs_constant1d_welltype
from transmissibility import transmissibility1d, transmissibility1d_boundary
from wellblock import fraction_wellblock_geometric_factor
from wellsimulation import solution_well1d
from potential import potential_term1d
from synthetics import constant_depth1d<jupyter_output><empty_output><jupyter_text>Additional input for Slightly Compressible simulation<jupyter_code>from wellsimulation import well1d_report_slicomp<jupyter_output><empty_output><jupyter_text>Incompressible SimulationFor `incompressible` simulation, `CFLUID` (fluid compressibility) and `CPORE` (pore compressibility) are set as 0 microsip. `B` (FVF) is 1.   Non-elevated<jupyter_code>import warnings
warnings.filterwarnings("ignore")

"""""""""""
INPUT TXT DATA
"""""""""""

filepath = '/content/pyresim/input/intermediate/intermediate1d_input.txt'
depthpath = '/content/pyresim/input/basic/basic1d_depth.txt'

# call the read_input function to read input data
reservoir_input, wells, west, east = read_input(filepath)

"""""""""""
RESERVOIR DEPTH

* If you have depth data (reservoir is ELEVATED), read the data using 
  'read_depth' function
* If you don't have and you want to simulate ELEVATED reservoir, 
  create your own synthetic data
* If you don't have and you want to simulate NON-ELEVATED reservoir, 
  create constant depth using 'constant_depth1d'
"""""""""""

z = 3000 # reservoir is NON-ELEVATED, constant depth
z_array = constant_depth1d(z, reservoir_input['xi']) # use function

"""""""""""
BOUNDARY DEPTH

Manually input the boundary depth
* If the reservoir is ELEVATED, input the depth
* If the reservoir is NON-ELEVATED, boundary depth must equal with the 
  reservoir depth
"""""""""""

west_depth = east_depth = z # NON-ELEVATED. Boundary depth = reservoir depth

"""""""""""
INPUT PROCESSING
"""""""""""

# number of blocks in x and y
xi = reservoir_input['xi']

# depth input processing
depth = z_array[1] # reservoir depth
west['depth'] = west_depth # west boundary depth
east['depth'] = east_depth # east boundary depth

# define parameters, in this case blocks are homogeneous and same in size
dx = reservoir_input['dx'] 
dy = reservoir_input['dy']
dz = reservoir_input['dz']
kx = reservoir_input['kx']
B = reservoir_input['B']
mu = reservoir_input['mu']
rho = reservoir_input['rho']

# well information (dataframe)
well_name = wells['well_name']
well_rw = wells['well_rw'] / 12  # wellbore radius, inch to ft
well_skin = wells['well_skin']
well_loc = wells['well_loc']
well_condition = wells['well_rw']
well_value = wells['well_value']
well_config = wells['well_config']

"""""""""""
WELL INFORMATION PROCESSING
"""""""""""

# source block (production or injection well)
well_loc = np.array(well_loc).astype(int)
qsc = source1d(well_value, well_loc, xi)  # call function HERE

# Calculate wellblock geometric factor
kh, r_eq, Gw = [], [], []

for i in range(len(well_config)):
    kh_, r_eq_, Gw_ = fraction_wellblock_geometric_factor(dx, dy, kx, kx, well_skin[i],
                                                          well_rw[i], dz, well_config[i])
    kh.append(kh_);
    r_eq.append(r_eq_);
    Gw.append(Gw_)

# create dataframe from the created dictionary
well_df = pd.DataFrame.from_dict(wells)

well_mu = np.full(len(well_loc), mu)
well_B = np.full(len(well_loc), B)
well_dz = np.full(len(well_loc), dz)

well_df['well_loc'] = well_loc
well_df['well_dz'] = dz
well_df['well_kh'] = kh
well_df['well_mu'] = mu
well_df['well_B'] = B
well_df['well_Gw'] = Gw

"""""""""""
SIMULATION
"""""""""""

lhs_mat = np.array([[0.] * xi] * xi)
rhs_mat = np.array([[0.] * 1] * xi)

for i in range(xi):

    T_min, T_plus = transmissibility1d(dx, dy, dz, kx, mu, B)

    if i == 0:
        # left boundary

        ## flow rate from boundary
        qsc_b = boundary_floweq1d(west['type'], dx, dy, dz, kx, mu, B, value=west['value'], no_block=(i + 1))

        ## LHS coeffs
        T_plus = T_plus
        T_min = transmissibility1d_boundary(west['type'], dx, dy, dz, kx, mu, B)

        p_plus, p_min, p = lhs_coeffs1d_welltype('west', (i + 1), well_df, qsc, T_plus, T_min, solver='incompressible')

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min, T_plus, west['depth'], depth[i], depth[i + 1])

        ## RHS constants
        rhs = rhs_constant1d_welltype(west['type'], (i + 1), well_df, potential_term, qsc_b, T_pb=T_min,
                                    p_b=west['value'], solver='incompressible')

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i] = p
        lhs_mat[i, i + 1] = p_plus

        ## fill in RHS matrix (rhs_mat)   
        rhs_mat[i, 0] = rhs

    elif i == xi - 1:
        # right boundary

        ## flow rate from boundary
        qsc_b = boundary_floweq1d(east['type'], dx, dy, dz, kx, mu, B, value=east['value'], no_block=(i + 1))

        ## LHS coeffs
        T_plus = transmissibility1d_boundary(east['type'], dx, dy, dz, kx, mu, B)
        T_min = T_min
        p_plus, p_min, p = lhs_coeffs1d_welltype('east', (i + 1), well_df, qsc, T_plus, T_min, solver='incompressible')

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min, T_plus, depth[i - 1], depth[i], east['depth'])

        ## RHS constants        
        rhs = rhs_constant1d_welltype(east['type'], (i + 1), well_df, potential_term, qsc_b,
                                    T_pb=T_plus, p_b=east['value'], solver='incompressible')

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i - 1] = p_min
        lhs_mat[i, i] = p

        ## fill in RHS matrix (rhs_mat)    
        rhs_mat[i, 0] = rhs

    else:

        ## calculate potential term for the RHS constant

        T_plus, T_min = T_plus, T_min
        p_plus, p_min, p = lhs_coeffs1d_welltype(None, (i+1), well_df, qsc, T_plus, T_min, solver='incompressible') 

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min, T_plus, depth[i - 1], depth[i], depth[i + 1])

        ## RHS constants      
        rhs = -qsc[i] + potential_term

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i - 1] = p_min
        lhs_mat[i, i] = p
        lhs_mat[i, i + 1] = p_plus

        ## fill in RHS matrix (rhs_mat)    
        rhs_mat[i, 0] = rhs  

"""""""""""
PRESSURE SOLVER
"""""""""""

p_sol = np.linalg.solve(lhs_mat, rhs_mat)

"""""""""""
DISPLAY PRESSURE
"""""""""""

p_sol = p_sol.T.reshape(-1)
p_sol_ = np.expand_dims(p_sol, axis=0)  # or axis=1
plt.imshow(p_sol_)
plt.colorbar(orientation='horizontal')
plt.show() 
p_sol
"""""""""""
WELL SIMULATION
"""""""""""

well_report = solution_well1d(well_df, p_sol)
well_report = well_report.drop(['well_dz', 'well_kh', 'well_mu', 'well_B', 'well_Gw'], 1)
well_report<jupyter_output><empty_output><jupyter_text>Elevated<jupyter_code>import warnings
warnings.filterwarnings("ignore")

"""""""""""
INPUT TXT DATA
"""""""""""

filepath = '/content/pyresim/input/intermediate/intermediate1d_input.txt'
depthpath = '/content/pyresim/input/advanced/advanced1d_depth.txt'

# call the read_input function to read input data
reservoir_input, wells, west, east = read_input(filepath)

"""""""""""
RESERVOIR DEPTH

* If you have depth data (reservoir is ELEVATED), read the data using 
  'read_depth' function
* If you don't have and you want to simulate ELEVATED reservoir, 
  create your own synthetic data
* If you don't have and you want to simulate NON-ELEVATED reservoir, 
  create constant depth using 'constant_depth1d'
"""""""""""

z_array = read_depth(depthpath)

"""""""""""
BOUNDARY DEPTH

Manually input the boundary depth
* If the reservoir is ELEVATED, input the depth
* If the reservoir is NON-ELEVATED, boundary depth must equal with the 
  reservoir depth
"""""""""""

west_depth = 4000
east_depth = 2900

"""""""""""
INPUT PROCESSING
"""""""""""

# number of blocks in x and y
xi = reservoir_input['xi']

# depth input processing
depth = z_array[1] # reservoir depth
west['depth'] = west_depth # west boundary depth
east['depth'] = east_depth # east boundary depth

# define parameters, in this case blocks are homogeneous and same in size
dx = reservoir_input['dx'] 
dy = reservoir_input['dy']
dz = reservoir_input['dz']
kx = reservoir_input['kx']
B = reservoir_input['B']
mu = reservoir_input['mu']
rho = reservoir_input['rho']

# well information (dataframe)
well_name = wells['well_name']
well_rw = wells['well_rw'] / 12  # wellbore radius, inch to ft
well_skin = wells['well_skin']
well_loc = wells['well_loc']
well_condition = wells['well_rw']
well_value = wells['well_value']
well_config = wells['well_config']

"""""""""""
WELL INFORMATION PROCESSING
"""""""""""

# source block (production or injection well)
well_loc = np.array(well_loc).astype(int)
qsc = source1d(well_value, well_loc, xi)  # call function HERE

# Calculate wellblock geometric factor
kh, r_eq, Gw = [], [], []

for i in range(len(well_config)):
    kh_, r_eq_, Gw_ = fraction_wellblock_geometric_factor(dx, dy, kx, kx, well_skin[i],
                                                          well_rw[i], dz, well_config[i])
    kh.append(kh_);
    r_eq.append(r_eq_);
    Gw.append(Gw_)

# create dataframe from the created dictionary
well_df = pd.DataFrame.from_dict(wells)

well_mu = np.full(len(well_loc), mu)
well_B = np.full(len(well_loc), B)
well_dz = np.full(len(well_loc), dz)

well_df['well_loc'] = well_loc
well_df['well_dz'] = dz
well_df['well_kh'] = kh
well_df['well_mu'] = mu
well_df['well_B'] = B
well_df['well_Gw'] = Gw

"""""""""""
SIMULATION
"""""""""""

lhs_mat = np.array([[0.] * xi] * xi)
rhs_mat = np.array([[0.] * 1] * xi)

for i in range(xi):

    T_min, T_plus = transmissibility1d(dx, dy, dz, kx, mu, B)

    if i == 0:
        # left boundary

        ## flow rate from boundary
        qsc_b = boundary_floweq1d(west['type'], dx, dy, dz, kx, mu, B, value=west['value'], no_block=(i + 1))

        ## LHS coeffs
        T_plus = T_plus
        T_min = transmissibility1d_boundary(west['type'], dx, dy, dz, kx, mu, B)

        p_plus, p_min, p = lhs_coeffs1d_welltype('west', (i + 1), well_df, qsc, T_plus, T_min, solver='incompressible')

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min, T_plus, west['depth'], depth[i], depth[i + 1])

        ## RHS constants
        rhs = rhs_constant1d_welltype(west['type'], (i + 1), well_df, potential_term, qsc_b, T_pb=T_min,
                                    p_b=west['value'], solver='incompressible')

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i] = p
        lhs_mat[i, i + 1] = p_plus

        ## fill in RHS matrix (rhs_mat)   
        rhs_mat[i, 0] = rhs

    elif i == xi - 1:
        # right boundary

        ## flow rate from boundary
        qsc_b = boundary_floweq1d(east['type'], dx, dy, dz, kx, mu, B, value=east['value'], no_block=(i + 1))

        ## LHS coeffs
        T_plus = transmissibility1d_boundary(east['type'], dx, dy, dz, kx, mu, B)
        T_min = T_min
        p_plus, p_min, p = lhs_coeffs1d_welltype('east', (i + 1), well_df, qsc, T_plus, T_min, solver='incompressible')

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min, T_plus, depth[i - 1], depth[i], east['depth'])

        ## RHS constants        
        rhs = rhs_constant1d_welltype(east['type'], (i + 1), well_df, potential_term, qsc_b,
                                    T_pb=T_plus, p_b=east['value'], solver='incompressible')

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i - 1] = p_min
        lhs_mat[i, i] = p

        ## fill in RHS matrix (rhs_mat)    
        rhs_mat[i, 0] = rhs

    else:

        ## calculate potential term for the RHS constant

        T_plus, T_min = T_plus, T_min
        p_plus, p_min, p = lhs_coeffs1d_welltype(None, (i+1), well_df, qsc, T_plus, T_min, solver='incompressible') 

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min, T_plus, depth[i - 1], depth[i], depth[i + 1])

        ## RHS constants      
        rhs = -qsc[i] + potential_term

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i - 1] = p_min
        lhs_mat[i, i] = p
        lhs_mat[i, i + 1] = p_plus

        ## fill in RHS matrix (rhs_mat)    
        rhs_mat[i, 0] = rhs  

"""""""""""
PRESSURE SOLVER
"""""""""""

p_sol = np.linalg.solve(lhs_mat, rhs_mat)

"""""""""""
DISPLAY PRESSURE
"""""""""""

p_sol = p_sol.T.reshape(-1)
p_sol_ = np.expand_dims(p_sol, axis=0)  # or axis=1
plt.imshow(p_sol_)
plt.colorbar(orientation='horizontal')
plt.show() 
p_sol
"""""""""""
WELL SIMULATION
"""""""""""

well_report = solution_well1d(well_df, p_sol)
well_report = well_report.drop(['well_dz', 'well_kh', 'well_mu', 'well_B', 'well_Gw'], 1)
well_report<jupyter_output><empty_output><jupyter_text>Slightly Compressible SimulationFor `slicomp` simulation, `CFLUID` is increased to 10 microsip, `CPORE` increased to 1 microsip, and `B` increased to 1.5.  Non-elevated<jupyter_code>import warnings
warnings.filterwarnings("ignore")

"""""""""""
PUT THE INPUT TXT DATA
"""""""""""

filepath = '/content/pyresim/input/intermediate/intermediate1d_input.txt'
depthpath = '/content/pyresim/input/basic/basic1d_depth.txt'

# call the read_input function to read input data
reservoir_input, wells, west, east = read_input(filepath)

"""""""""""
SLIGHTLY INCOMPRESSIBLE INPUTS
"""""""""""

p_initial = np.full(reservoir_input['xi'], 4000) # initial pressure
timestep = 1 # time increment, day
schedule = 20 # how long will be simulated, days

"""""""""""
RESERVOIR DEPTH

* If you have depth data (reservoir is ELEVATED), read the data using 
  'read_depth' function
* If you don't have and you want to simulate ELEVATED reservoir, 
  create your own synthetic data
* If you don't have and you want to simulate NON-ELEVATED reservoir, 
  create constant depth using 'constant_depth1d'
"""""""""""

z = 3000 # reservoir is NON-ELEVATED, constant depth
z_array = constant_depth1d(z, reservoir_input['xi']) # use function

"""""""""""
BOUNDARY DEPTH

Manually input the boundary depth
* If the reservoir is ELEVATED, input the depth
* If the reservoir is NON-ELEVATED, boundary depth must equal with the 
  reservoir depth
"""""""""""

west_depth = east_depth = z # NON-ELEVATED. Boundary depth = reservoir depth

"""""""""""
INPUT PROCESSING
"""""""""""

# number of blocks in x and y
xi = reservoir_input['xi']

# depth input processing
depth = z_array[1] # reservoir depth
west['depth'] = west_depth # west boundary depth
east['depth'] = east_depth # east boundary depth

# define parameters, in this case blocks are homogeneous and same in size
dx = reservoir_input['dx'] 
dy = reservoir_input['dy']
dz = reservoir_input['dz']
kx = reservoir_input['kx']
B = reservoir_input['B']
mu = reservoir_input['mu']
rho = reservoir_input['rho']

# well information (dataframe)
well_name = wells['well_name']
well_rw = wells['well_rw'] / 12  # wellbore radius, inch to ft
well_skin = wells['well_skin']
well_loc = wells['well_loc']
well_condition = wells['well_rw']
well_value = wells['well_value']
well_config = wells['well_config']

" For Slightly Incompressible sim, there are information of CFLUID and CPORE "

# convert unit microsip to sip (* 10^-6)
reservoir_input['cfluid'] = reservoir_input['cfluid'] * 1E-6
reservoir_input['cpore'] = reservoir_input['cpore'] * 1E-6

"""""""""""
WELL INFORMATION PROCESSING
"""""""""""

# source block (production or injection well)
well_loc = np.array(well_loc).astype(int)
qsc = source1d(well_value, well_loc, xi)  # call function HERE

# Calculate wellblock geometric factor
kh, r_eq, Gw = [], [], []

for i in range(len(well_config)):
    kh_, r_eq_, Gw_ = fraction_wellblock_geometric_factor(dx, dy, kx, kx, well_skin[i],
                                                          well_rw[i], dz, well_config[i])
    kh.append(kh_);
    r_eq.append(r_eq_);
    Gw.append(Gw_)

# create dataframe from the created dictionary
well_df = pd.DataFrame.from_dict(wells)

well_mu = np.full(len(well_loc), mu)
well_B = np.full(len(well_loc), B)
well_dz = np.full(len(well_loc), dz)

well_df['well_loc'] = well_loc
well_df['well_dz'] = dz
well_df['well_kh'] = kh
well_df['well_mu'] = mu
well_df['well_B'] = B
well_df['well_Gw'] = Gw

"""""""""""
SIMULATION
"""""""""""

" Produce LHS matrix. In Slightly Compressible simulation, LHS can be assumed CONSTANT "

lhs_mat = np.array([[0.] * xi] * xi)

# set up empty array for boundary flow rate, and two transmissibilities 
# for later computation of RHS and pressure solve at each timestep

qsc_b_array = []; T_min_array = []; T_plus_array = []

for i in range(xi):

    T_min, T_plus = transmissibility1d(dx, dy, dz, kx, mu, B)

    if i == 0:
        # left boundary

        ## flow rate from boundary
        qsc_b = boundary_floweq1d(west['type'], dx, dy, dz, kx, mu, B, value=west['value'], no_block=(i + 1))

        ## LHS coeffs
        T_plus = T_plus
        T_min = transmissibility1d_boundary(west['type'], dx, dy, dz, kx, mu, B)

        p_plus, p_min, p = lhs_coeffs1d_welltype('west', (i + 1), well_df, qsc, T_plus, T_min,
                                                 solver='slicomp', reservoir_input=reservoir_input, 
                                                 timestep=timestep)
        qsc_b_array.append(qsc_b)
        T_min_array.append(T_min)
        T_plus_array.append(T_plus)        

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i] = p
        lhs_mat[i, i + 1] = p_plus

    elif i == xi - 1:
        # right boundary

        ## flow rate from boundary
        qsc_b = boundary_floweq1d(east['type'], dx, dy, dz, kx, mu, B, value=east['value'], no_block=(i + 1))

        ## LHS coeffs
        T_plus = transmissibility1d_boundary(east['type'], dx, dy, dz, kx, mu, B)
        T_min = T_min

        p_plus, p_min, p = lhs_coeffs1d_welltype('east', (i + 1), well_df, qsc, T_plus, T_min,
                                                 solver='slicomp', reservoir_input=reservoir_input, 
                                                 timestep=timestep)  

        qsc_b_array.append(qsc_b)
        T_min_array.append(T_min)
        T_plus_array.append(T_plus)  

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i - 1] = p_min
        lhs_mat[i, i] = p


    else:
        # interior block
        T_plus, T_min = T_plus, T_min
        p_plus, p_min, p = lhs_coeffs1d_welltype(None, (i+1), well_df, qsc, T_plus, T_min,
                                                 solver='slicomp', reservoir_input=reservoir_input, 
                                                 timestep=timestep)

        qsc_b_array.append(qsc_b)
        T_min_array.append(T_min)
        T_plus_array.append(T_plus)                                                           

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i - 1] = p_min
        lhs_mat[i, i] = p
        lhs_mat[i, i + 1] = p_plus

" Timestep evolution of computing RHS and solving the pressure "

# initiate solution pressure with the initial pressure array (p_initial)
p_sol = p_initial

p_sol_record = []
for k in range(schedule):

  rhs_mat = np.array([[0.] * 1] * xi)  
  for i in range(xi):
    # the pressure in grid block. FIRST timestep, it equals 'p_initial'
    # the NEXT timesteps, it equals the solved and updated pressure   
    p_sol = p_sol

    if i == 0:
        # left boundary

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min_array[i], T_plus_array[i], west['depth'], depth[i], depth[i + 1])

        ## RHS constants
        rhs = rhs_constant1d_welltype(west['type'], (i + 1), well_df, potential_term, qsc_b_array[i], T_pb=T_min_array[i],
                                    p_b=west['value'], solver='slicomp', p_initial=p_sol[i],
                                    reservoir_input=reservoir_input, timestep=timestep)

        ## fill in RHS matrix (rhs_mat)   
        rhs_mat[i, 0] = rhs

    elif i == xi - 1:
        # right boundary                                               

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min_array[i], T_plus_array[i], depth[i - 1], depth[i], east['depth'])

        ## RHS constants       
        rhs = rhs_constant1d_welltype(east['type'], (i + 1), well_df, potential_term, qsc_b_array[i],
                                    T_pb=T_plus_array[i], p_b=east['value'], solver='slicomp', p_initial=p_sol[i],
                                    reservoir_input=reservoir_input, timestep=timestep) 

        ## fill in RHS matrix (rhs_mat)    
        rhs_mat[i, 0] = rhs

    else:   
        # interior blocks

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min_array[i], T_plus_array[i], depth[i - 1], depth[i], depth[i + 1])

        # RHS constants
        # add term 
        Vb = reservoir_input['dx'] * reservoir_input['dy'] * reservoir_input['dz']
        ct = reservoir_input['cpore'] + reservoir_input['cfluid']
        rhs_term = (Vb * reservoir_input['poro'] * ct) / (5.614583 * reservoir_input['B'] * timestep) 

        rhs = (-qsc[i] + potential_term) - (rhs_term * p_sol[i])

        ## fill in RHS matrix (rhs_mat)    
        rhs_mat[i, 0] = rhs 

  """""""""""
  PRESSURE SOLVER
  """""""""""
  p_sol = np.linalg.solve(lhs_mat, rhs_mat)  
  p_sol_record.append(p_sol) 

"""""""""""
OUTPUT
"""""""""""

p_sol_ = []

for i in range(len(p_sol_record)):
  p_sol = p_sol_record[i].T.reshape(-1)
  p_sol = np.expand_dims(p_sol, axis=0)
  p_sol_.append(p_sol)  

# merge the initial pressure to the solved pressure afterwards
p_initial_ = p_initial.reshape((1,1,xi))
p_sol_ = np.concatenate((p_initial_, p_sol_), axis=0)
"""""""""""
INTERACTIVE ANIMATION
"""""""""""

from __future__ import print_function
from ipywidgets import interact, interactive, fixed, interact_manual, ToggleButtons
import ipywidgets as widgets

min, max = np.round(np.amin(p_sol_)), np.round(np.amax(p_sol_))

@interact

def display_pressure(day=(0, schedule)):
  fig, ax = plt.subplots()  
  im = ax.imshow(p_sol_[day], vmin=min, vmax=max)
  ax.set_title('Pressure at day {}'.format(day))
  cax = fig.add_axes([0.1, 0.2, 0.8, 0.05])
  fig.colorbar(im, cax=cax, orientation='horizontal')

  # # add values on the grid blocks
  # for (j,i),label in np.ndenumerate(np.round(p_sol_[day], 3)):
  #     ax.text(i,j,label,ha='center',va='center')  

  plt.show() 
"""""""""""
WELL SIMULATION
"""""""""""

well_report = well1d_report_slicomp('B', wells, well_df, p_sol_record, p_initial, schedule)

plt.plot(well_report['time (day)'], well_report['grid block pressure (psi)'], label='Grid Block Pressure (psi)')
plt.plot(well_report['time (day)'], well_report['fbhp_sol (psi)'], label='Flowing Borehole Pressure (psi)')
plt.xlabel('Days'); plt.ylabel('Pressure (psi)')
plt.xlim(0, schedule)
plt.legend()
plt.show()

well_report<jupyter_output><empty_output><jupyter_text>Elevated<jupyter_code>import warnings
warnings.filterwarnings("ignore")

"""""""""""
PUT THE INPUT TXT DATA
"""""""""""

filepath = '/content/pyresim/input/intermediate/intermediate1d_input.txt'
depthpath = '/content/pyresim/input/advanced/advanced1d_depth.txt'

# call the read_input function to read input data
reservoir_input, wells, west, east = read_input(filepath)

"""""""""""
SLIGHTLY INCOMPRESSIBLE INPUTS
"""""""""""

p_initial = np.full(reservoir_input['xi'], 4000) # initial pressure
timestep = 1 # time increment, day
schedule = 20 # how long will be simulated, days

"""""""""""
RESERVOIR DEPTH

* If you have depth data (reservoir is ELEVATED), read the data using 
  'read_depth' function
* If you don't have and you want to simulate ELEVATED reservoir, 
  create your own synthetic data
* If you don't have and you want to simulate NON-ELEVATED reservoir, 
  create constant depth using 'constant_depth1d'
"""""""""""

z_array = read_depth(depthpath)

"""""""""""
BOUNDARY DEPTH

Manually input the boundary depth
* If the reservoir is ELEVATED, input the depth
* If the reservoir is NON-ELEVATED, boundary depth must equal with the 
  reservoir depth
"""""""""""

west_depth = 4000
east_depth = 2900

"""""""""""
INPUT PROCESSING
"""""""""""

# number of blocks in x and y
xi = reservoir_input['xi']

# depth input processing
depth = z_array[1] # reservoir depth
west['depth'] = west_depth # west boundary depth
east['depth'] = east_depth # east boundary depth

# define parameters, in this case blocks are homogeneous and same in size
dx = reservoir_input['dx'] 
dy = reservoir_input['dy']
dz = reservoir_input['dz']
kx = reservoir_input['kx']
B = reservoir_input['B']
mu = reservoir_input['mu']
rho = reservoir_input['rho']

# well information (dataframe)
well_name = wells['well_name']
well_rw = wells['well_rw'] / 12  # wellbore radius, inch to ft
well_skin = wells['well_skin']
well_loc = wells['well_loc']
well_condition = wells['well_rw']
well_value = wells['well_value']
well_config = wells['well_config']

" For Slightly Incompressible sim, there are information of CFLUID and CPORE "

# convert unit microsip to sip (* 10^-6)
reservoir_input['cfluid'] = reservoir_input['cfluid'] * 1E-6
reservoir_input['cpore'] = reservoir_input['cpore'] * 1E-6

"""""""""""
WELL INFORMATION PROCESSING
"""""""""""

# source block (production or injection well)
well_loc = np.array(well_loc).astype(int)
qsc = source1d(well_value, well_loc, xi)  # call function HERE

# Calculate wellblock geometric factor
kh, r_eq, Gw = [], [], []

for i in range(len(well_config)):
    kh_, r_eq_, Gw_ = fraction_wellblock_geometric_factor(dx, dy, kx, kx, well_skin[i],
                                                          well_rw[i], dz, well_config[i])
    kh.append(kh_);
    r_eq.append(r_eq_);
    Gw.append(Gw_)

# create dataframe from the created dictionary
well_df = pd.DataFrame.from_dict(wells)

well_mu = np.full(len(well_loc), mu)
well_B = np.full(len(well_loc), B)
well_dz = np.full(len(well_loc), dz)

well_df['well_loc'] = well_loc
well_df['well_dz'] = dz
well_df['well_kh'] = kh
well_df['well_mu'] = mu
well_df['well_B'] = B
well_df['well_Gw'] = Gw

"""""""""""
SIMULATION
"""""""""""

" Produce LHS matrix. In Slightly Compressible simulation, LHS can be assumed CONSTANT "

lhs_mat = np.array([[0.] * xi] * xi)

# set up empty array for boundary flow rate, and two transmissibilities 
# for later computation of RHS and pressure solve at each timestep

qsc_b_array = []; T_min_array = []; T_plus_array = []

for i in range(xi):

    T_min, T_plus = transmissibility1d(dx, dy, dz, kx, mu, B)

    if i == 0:
        # left boundary

        ## flow rate from boundary
        qsc_b = boundary_floweq1d(west['type'], dx, dy, dz, kx, mu, B, value=west['value'], no_block=(i + 1))

        ## LHS coeffs
        T_plus = T_plus
        T_min = transmissibility1d_boundary(west['type'], dx, dy, dz, kx, mu, B)

        p_plus, p_min, p = lhs_coeffs1d_welltype('west', (i + 1), well_df, qsc, T_plus, T_min,
                                                 solver='slicomp', reservoir_input=reservoir_input, 
                                                 timestep=timestep)
        qsc_b_array.append(qsc_b)
        T_min_array.append(T_min)
        T_plus_array.append(T_plus)        

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i] = p
        lhs_mat[i, i + 1] = p_plus

    elif i == xi - 1:
        # right boundary

        ## flow rate from boundary
        qsc_b = boundary_floweq1d(east['type'], dx, dy, dz, kx, mu, B, value=east['value'], no_block=(i + 1))

        ## LHS coeffs
        T_plus = transmissibility1d_boundary(east['type'], dx, dy, dz, kx, mu, B)
        T_min = T_min

        p_plus, p_min, p = lhs_coeffs1d_welltype('east', (i + 1), well_df, qsc, T_plus, T_min,
                                                 solver='slicomp', reservoir_input=reservoir_input, 
                                                 timestep=timestep)  

        qsc_b_array.append(qsc_b)
        T_min_array.append(T_min)
        T_plus_array.append(T_plus)  

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i - 1] = p_min
        lhs_mat[i, i] = p


    else:
        # interior block
        T_plus, T_min = T_plus, T_min
        p_plus, p_min, p = lhs_coeffs1d_welltype(None, (i+1), well_df, qsc, T_plus, T_min,
                                                 solver='slicomp', reservoir_input=reservoir_input, 
                                                 timestep=timestep)

        qsc_b_array.append(qsc_b)
        T_min_array.append(T_min)
        T_plus_array.append(T_plus)                                                           

        ## fill in LHS matrix (lhs_mat)
        lhs_mat[i, i - 1] = p_min
        lhs_mat[i, i] = p
        lhs_mat[i, i + 1] = p_plus

" Timestep evolution of computing RHS and solving the pressure "

# initiate solution pressure with the initial pressure array (p_initial)
p_sol = p_initial

p_sol_record = []
for k in range(schedule):

  rhs_mat = np.array([[0.] * 1] * xi)  
  for i in range(xi):
    # the pressure in grid block. FIRST timestep, it equals 'p_initial'
    # the NEXT timesteps, it equals the solved and updated pressure   
    p_sol = p_sol

    if i == 0:
        # left boundary

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min_array[i], T_plus_array[i], west['depth'], depth[i], depth[i + 1])

        ## RHS constants
        rhs = rhs_constant1d_welltype(west['type'], (i + 1), well_df, potential_term, qsc_b_array[i], T_pb=T_min_array[i],
                                    p_b=west['value'], solver='slicomp', p_initial=p_sol[i],
                                    reservoir_input=reservoir_input, timestep=timestep)

        ## fill in RHS matrix (rhs_mat)   
        rhs_mat[i, 0] = rhs

    elif i == xi - 1:
        # right boundary                                               

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min_array[i], T_plus_array[i], depth[i - 1], depth[i], east['depth'])

        ## RHS constants       
        rhs = rhs_constant1d_welltype(east['type'], (i + 1), well_df, potential_term, qsc_b_array[i],
                                    T_pb=T_plus_array[i], p_b=east['value'], solver='slicomp', p_initial=p_sol[i],
                                    reservoir_input=reservoir_input, timestep=timestep) 

        ## fill in RHS matrix (rhs_mat)    
        rhs_mat[i, 0] = rhs

    else:   
        # interior blocks

        ## calculate potential term for the RHS constant
        potential_term = potential_term1d(rho, T_min_array[i], T_plus_array[i], depth[i - 1], depth[i], depth[i + 1])

        # RHS constants
        # add term 
        Vb = reservoir_input['dx'] * reservoir_input['dy'] * reservoir_input['dz']
        ct = reservoir_input['cpore'] + reservoir_input['cfluid']
        rhs_term = (Vb * reservoir_input['poro'] * ct) / (5.614583 * reservoir_input['B'] * timestep) 

        rhs = (-qsc[i] + potential_term) - (rhs_term * p_sol[i])

        ## fill in RHS matrix (rhs_mat)    
        rhs_mat[i, 0] = rhs 

  """""""""""
  PRESSURE SOLVER
  """""""""""
  p_sol = np.linalg.solve(lhs_mat, rhs_mat)  
  p_sol_record.append(p_sol) 

"""""""""""
OUTPUT
"""""""""""

p_sol_ = []

for i in range(len(p_sol_record)):
  p_sol = p_sol_record[i].T.reshape(-1)
  p_sol = np.expand_dims(p_sol, axis=0)
  p_sol_.append(p_sol)  

# merge the initial pressure to the solved pressure afterwards
p_initial_ = p_initial.reshape((1,1,xi))
p_sol_ = np.concatenate((p_initial_, p_sol_), axis=0)
"""""""""""
INTERACTIVE ANIMATION
"""""""""""

from __future__ import print_function
from ipywidgets import interact, interactive, fixed, interact_manual, ToggleButtons
import ipywidgets as widgets

min, max = np.round(np.amin(p_sol_)), np.round(np.amax(p_sol_))

@interact

def display_pressure(day=(0, schedule)):
  fig, ax = plt.subplots()  
  im = ax.imshow(p_sol_[day], vmin=min, vmax=max)
  ax.set_title('Pressure at day {}'.format(day))
  cax = fig.add_axes([0.1, 0.2, 0.8, 0.05])
  fig.colorbar(im, cax=cax, orientation='horizontal')

  # # add values on the grid blocks
  # for (j,i),label in np.ndenumerate(np.round(p_sol_[day], 3)):
  #     ax.text(i,j,label,ha='center',va='center')  

  plt.show() 
"""""""""""
WELL SIMULATION
"""""""""""

well_report = well1d_report_slicomp('B', wells, well_df, p_sol_record, p_initial, schedule)

plt.plot(well_report['time (day)'], well_report['grid block pressure (psi)'], label='Grid Block Pressure (psi)')
plt.plot(well_report['time (day)'], well_report['fbhp_sol (psi)'], label='Flowing Borehole Pressure (psi)')
plt.xlabel('Days'); plt.ylabel('Pressure (psi)')
plt.xlim(0, schedule)
plt.legend()
plt.show()

well_report<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1391908483813973, 'fraction_numerical': 0.009223489952984158, 'mean_word_length': 3.239969527679025, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 72, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '153943', 'license': 'Apache-2.0', 'n_tokens_mistral': 12291, 'n_tokens_neox': 11707, 'n_words': 2904, 'path': 'simulators/PyReSim1d.ipynb', 'repo_name': 'yohanesnuwara/2016-ml-contest', 'size': '33393'}