single_phase_utils.py

from __future__ import division

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.sparse.linalg import gmres

"""
Tita Ristanto - 2018 

This file is an implementation of a discretized single-phase homogeneous isothermal 2D reservoir 
equation: T*p^(n+1)=D(p^(n+1)-p^(n))+Q, where T is a transmissibility matrix, D is an accumulation
matrix, and Q is a source/sink matrix (which represents well(s)). p^(n+1) is a vector of unknown 
(in this case, pressure) at time level n and p^(n) is the same variable at current time level n.
This case assumes no capillary pressure, isotropic permeability, and constant viscosity. 
Well treatment is not included. This code treats the unknown(pressure) implicitly and transmissibilities 
explicitly.

Generalized Minimal Residual iteration (GMRES) is used to solve the linear matrix system, 
which is an iterative method for the numerical solution of a nonsymmetric system of 
linear equations. This function solves for p at time level n+1. This new pressure vector 
updates formation volume factor B_o in every grid, which also updates transmissibility matrix, 
so then we can solve for pressure at the next time level. This process continues until 
it reaches the end of simulation time.

The main simulation loop in this file calls the `run_simulation()` function for
simulating the reservoir dynamics. The `main()` function runs three different example cases: 
- 1 producer
- 3 producers
- 1 producer and 1 injector

Variable of interest demonstrated here includes block pressure (in the middle of the reservoir)
and spatial pressure.

"""


class prop_rock(object):
    """
    This is a class that captures rock physical properties, including permeability, porosity, and
    compressibility.
    """

    def __init__(self, kx=0, ky=0, por=0, cr=0):
        self.kx = kx
        self.ky = ky
        self.por = por
        self.cr = cr


class prop_fluid(object):
    """
    This class captures fluid properties. Basic properties are assumed constant.
    Formation volume factor (b) is a function of pressure.
    """

    def __init__(self, c_o=0, mu_o=0, rho_o=0):
        self.c_o = c_o
        self.mu_o = mu_o
        self.rho_o = rho_o

    def calc_b(self, p):
        return 1 / (1 + self.c_o * (p - 14.7))


class prop_grid(object):
    """This describes grid dimension and numbers."""

    def __init__(self, Nx=0, Ny=0, Nz=0):
        self.Nx = Nx
        self.Ny = Ny
        self.Nz = Nz

    def grid_dimension_x(self, Lx):
        return self.Nx / Lx

    def grid_dimension_y(self, Ly):
        return self.Ny / Ly

    def grid_dimension_z(self, Lz):
        return self.Nz / Lz


class prop_res(object):
    """A class that captures reservoir dimension and initial pressure."""

    def __init__(self, Lx=0, Ly=0, Lz=0, p_init=0):
        self.Lx = Lx
        self.Ly = Ly
        self.Lz = Lz
        self.p_init = p_init


class prop_well(object):
    """Describes well location a flow rate. Also provides conversion from
    cartesian i,j coordinate to grid number"""

    def __init__(self, loc=0, q=0):
        self.loc = loc
        self.q = q

    def index_to_grid(self, Nx):
        return self.loc[1] * Nx + self.loc[0]


class prop_time(object):
    """Describes time-step (assumed constant) and time interval"""

    def __init__(self, tstep=0, timeint=0):
        self.tstep = tstep
        self.timeint = timeint


def load_data(filename):
    """Loads ECLIPSE simulation block pressure data as a comparison"""
    df = pd.read_csv(filename)

    t = df.loc[:, ['TIME']]  # Time in simulation: DAY
    p = df.loc[:, ['BPR:(18,18,1)']]
    print('Data loaded!')
    return t, p


def calc_transmissibility_x(k_x, mu, B_o, props, i, j):
    """Calculates transmissibility in x-direction"""
    dx = props['res'].Lx / props['grid'].Nx
    dy = props['res'].Ly / props['grid'].Ny
    dz = props['res'].Lz / props['grid'].Nz

    k_x = (k_x[j, i] + k_x[j, i + 1]) / 2
    mu = (mu[j, i] + mu[j, i + 1]) / 2
    B_o = (B_o[j, i] + B_o[j, i + 1]) / 2
    return k_x * dy * dz / mu / B_o / dx


def calc_transmissibility_y(k_y, mu, B_o, props, i, j):
    """Calculates transmissibility in y-direction"""
    dx = props['res'].Lx / props['grid'].Nx
    dy = props['res'].Ly / props['grid'].Ny
    dz = props['res'].Lz / props['grid'].Nz

    k_y = (k_y[j, i] + k_y[j + 1, i]) / 2
    mu = (mu[j, i] + mu[j + 1, i]) / 2
    B_o = (B_o[j, i] + B_o[j + 1, i]) / 2
    return k_y * dx * dz / mu / B_o / dy


def ij_to_grid(i, j, Nx):
    """Converts i,j coordinate to grid number"""
    return (i) + Nx * j


def construct_T(mat, params, props):
    """
    Given various rock and fluid properties and grid geometry, this function constructs
    transmissibility matrix T. Matrix T is a tri-diagonal matrix.
    """
    k_x = params['k_x']
    k_y = params['k_y']
    B_o = params['B_o']
    mu = params['mu']

    m = mat.shape[0]
    n = mat.shape[1]
    A = np.zeros((m * n, m * n))
    for j in range(m):
        for i in range(n):
            # 2 neighbors in x direction
            if i < n - 1:
                A[mat[j, i] - 1, mat[j, i + 1] - 1] = calc_transmissibility_x(k_x, mu, B_o, props, i, j)
                A[mat[j, i + 1] - 1, mat[j, i] - 1] = A[mat[j, i] - 1, mat[j, i + 1] - 1]
            # 2 neighbors in y direction
            if j < m - 1:
                A[mat[j, i] - 1, mat[j + 1, i] - 1] = calc_transmissibility_y(k_y, mu, B_o, props, i, j)
                A[mat[j + 1, i] - 1, mat[j, i] - 1] = A[mat[j, i] - 1, mat[j + 1, i] - 1]

    for k in range(A.shape[0]):
        p = np.sum(A[k, :]) * -1
        A[k, k] = p
    return A / 887.5


def run_simulation(props):
    """
    This function controls the simulation loop. Basically at every time-step, it
    solves the linear system, captures variable of interest, and update parameters.

    :param props: a dictionary containing rock, fluid, and well properties.
    :return: p_well_block: time-series block pressure
             p_grids: 2D pressure matrix representing spatial distribution of pressure
    """
    rock = props['rock']
    fluid = props['fluid']
    grid = props['grid']
    res = props['res']
    wells = props['well']
    sim_time = props['time']

    # Distribute properties (and their values) to the grid
    k_x = np.full((grid.Ny, grid.Nx), rock.kx)
    k_y = np.full((grid.Ny, grid.Nx), rock.ky)
    B_o = np.full((grid.Ny, grid.Nx), fluid.calc_b(res.p_init))
    mu = np.full((grid.Ny, grid.Nx), fluid.mu_o)
    p_grids = np.full((grid.Ny * grid.Nx, 1), res.p_init)
    params = {'k_x': k_x, 'k_y': k_y, 'B_o': B_o, 'mu': mu}

    # Construct transmissibility matrix T
    mat = np.reshape(np.arange(1, grid.Ny * grid.Nx + 1), (grid.Ny, grid.Nx))
    T = construct_T(mat, params, props)

    # Create matrix A = transmissibility matrix - accumulation matrix
    dx = res.Lx / grid.Nx
    dy = res.Lx / grid.Nx
    dz = res.Lx / grid.Nx
    V = dx * dy * dz
    accumulation = V * rock.por * fluid.c_o / 5.615 / (sim_time.tstep)
    A = T - np.eye(T.shape[0]) * accumulation

    # Assign well flow rate to Q matrix
    Q = np.zeros((T.shape[0], 1))
    for well in wells:
        Q[well.index_to_grid(grid.Nx)] = -well.q

    # Calculate right hand side of the equation
    p_n = np.full((grid.Ny * grid.Nx, 1), -accumulation * res.p_init)
    b = p_n - Q

    # Initiate variable of interest: pressure in block (18,18)
    p_well_block = []

    # Time-loop
    for t in sim_time.timeint:
        print('evaluating t = %1.1f (days)' % t)
        p_well_block.append(p_grids[wells[0].index_to_grid(grid.Nx)])

        # Calculate pressure at time level n+1
        p_grids = (gmres(A, b))[0]
        p_grids = np.reshape(p_grids, (len(p_grids), 1))

        # Update B, b, and transmissibility matrix
        for i in range(grid.Nx):
            for j in range(grid.Ny):
                B_o[i, j] = fluid.calc_b(p_grids[ij_to_grid(i, j, grid.Nx)])
        params['B_o'] = B_o
        A = construct_T(mat, params, props)
        A = A - np.eye(A.shape[0]) * accumulation

        b = -accumulation * p_grids - Q
    return p_well_block, p_grids


def plot_pressure(t, p_pred, label, color):
    """Plots pressure v time"""
    plt.plot(t, p_pred, color=color, markeredgecolor=color, label=label)
    plt.xlabel("Time (days)")
    plt.ylabel("Block Pressure Cell (18,18)", fontsize=9)
    plt.legend(loc="best", prop=dict(size=8))
    plt.xlim(0, max(t))
    plt.ylim(0, max(p_pred))
    plt.grid(True)
    plt.draw()


def spatial_map(p_2D, title):
    """Plots variable of interest on a 2D spatial map"""
    plt.matshow(p_2D)
    plt.colorbar()
    plt.xlabel('grid in x-direction')
    plt.ylabel('grid in y-direction')
    plt.title(title)
    plt.draw()


def derivatives(t, p_pred, p_act, title):
    """Computes the derivatives of pressure"""
    t_act = np.linspace(0, 400, len(p_act))
    dp_act = abs(p_act[0] - p_act)
    dp_pred = abs(p_pred[0] - p_pred)

    # Calculate Derivatives
    p_der_act = np.zeros(len(p_act) - 2)
    p_der_pred = np.zeros(len(p_pred) - 2)
    for i in range(1, len(p_act) - 1):
        p_der_act[i - 1] = t_act[i] / (t_act[i + 1] - t_act[i - 1]) * abs(dp_act[i + 1] - dp_act[i - 1])
    for i in range(1, len(p_pred) - 1):
        p_der_pred[i - 1] = t[i] / (t[i + 1] - t[i - 1]) * abs(dp_pred[i + 1] - dp_pred[i - 1])
    plot_derivatives(t_act, t, dp_act, dp_pred, p_der_act, p_der_pred, title)


def plot_derivatives(t_act, t, dp_act, dp_pred, p_der_act, p_der_pred, title):
    """Computes the derivatives of pressure"""
    plt.figure()
    plt.loglog(t_act, dp_act, 'k-', linewidth=3, label='Actual dP')
    plt.loglog(t, dp_pred, 'ro', label='Predicted dP')
    plt.loglog(t_act[1:-1], p_der_act, 'k*', linewidth=3, label='Actual Derivative')
    plt.loglog(t[1:-1], p_der_pred, 'gx', label='Predicted Derivative')
    plt.title(title)
    plt.xlabel('dt (hours)')
    plt.ylabel('dP & Derivative')
    plt.legend(loc="best", prop=dict(size=12))
    plt.grid()
    plt.draw()