Skip to content

Commit

Permalink
Starting to modify code.
Browse files Browse the repository at this point in the history 
Haven't done anything yet but highlight where the work needs to be done.
  • Loading branch information
@chaosmarder
chaosmarder committed Oct 28, 2023
1 parent 1059713 commit b70656c
Showing 1 changed file with 140 additions and 93 deletions.
233 changes: 140 additions & 93 deletions src/bluebonnet/flow/reservoir.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@

from __future__ import annotations

from collections.abc import Callable
from dataclasses import dataclass
from typing import Callable

import numpy as np
import numpy.typing as npt
Expand All @@ -16,16 +16,12 @@

from bluebonnet.flow.flowproperties import FlowProperties

_ATOL = 1e-12


@dataclass
class IdealReservoir:
"""
Class for building scaling solutions of production from hydrofractured wells.
This maintains simulation parameters and fluid properties.
Parameters
----------
nx : int
Expand All @@ -36,18 +32,19 @@ class IdealReservoir:
reservoir pressure before production (psi)
fluid : FlowProperties
reservoir fluid PVT/flow properties
Methods
----------
simulate : calculate pressure over time
recovery_factor : calculate recovery factor over time
"""

nx: int
"""number of spatial nodes"""
pressure_fracface: float | npt.NDArray
"""drawdown pressure at :math:`x=0` (psi)"""
pressure_initial: float
"""reservoir pressure before production (psi)"""
fluid: FlowProperties | None = None
"""reservoir fluid PVT/flow properties"""
fluid: FlowProperties

def __post_init__(self):
def __post_init___(self):
"""Last initialization steps."""

def simulate(self, time: ndarray):
Expand All @@ -70,10 +67,10 @@ def simulate(self, time: ndarray):
alpha_scaled = self.alpha_scaled(b)
kt_h2 = mesh_ratio * alpha_scaled
a_matrix = _build_matrix(kt_h2)
pseudopressure[i + 1], _ = sparse.linalg.bicgstab(a_matrix, b, atol=_ATOL)
pseudopressure[i + 1], _ = sparse.linalg.bicgstab(a_matrix, b)
self.pseudopressure = pseudopressure

def recovery_factor(self, time: ndarray | None = None, density=False) -> ndarray:
def recovery_factor(self, time: ndarray | None = None) -> ndarray:
"""Calculate recovery factor over time.
If time has is not specified, requires that `simulate` has been run
Expand All @@ -91,23 +88,14 @@ def recovery_factor(self, time: ndarray | None = None, density=False) -> ndarray
if time is None:
try:
time = self.time
except AttributeError as e:
msg = "Need to run simulate before calculating recovery factor"
raise RuntimeError(msg) from e
except AttributeError:
raise RuntimeError(
"Need to run simulate before calculating recovery factor",
)
h_inv = self.nx - 1.0
if density:
pseudopressure_to_mass = interpolate.interp1d(
self.fluid.pvt_props["m-scaled"],
self.fluid.pvt_props["density"],
fill_value="extrapolate",
)
mass = pseudopressure_to_mass(self.pseudopressure)
mass_over_time = np.sum(mass, 1)
cumulative = 1.0 - mass_over_time / mass_over_time[0]
else:
pp = self.pseudopressure[:, :3]
rate = (-pp[:, 2] + 4 * pp[:, 1] - 3 * pp[:, 0]) * h_inv * 0.5 # dp_dx
cumulative = integrate.cumulative_trapezoid(rate, self.time, initial=0)
pp = self.pseudopressure[:, :3]
dp_dx = (-pp[:, 2] + 4 * pp[:, 1] - 3 * pp[:, 0]) * h_inv * 0.5
cumulative = integrate.cumulative_trapezoid(dp_dx, self.time, initial=0)
self.recovery = cumulative * self.fvf_scale()
return self.recovery

Expand All @@ -122,9 +110,10 @@ def recovery_factor_interpolator(self) -> Callable:
"""
try:
time = self.time
except AttributeError as e:
msg = "Need to run simulate"
raise RuntimeError(msg) from e
except AttributeError:
raise RuntimeError(
"Need to run simulate",
)
try:
recovery = self.recovery
except AttributeError:
Expand All @@ -142,9 +131,8 @@ def alpha_scaled(self, pseudopressure: ndarray) -> ndarray:
def fvf_scale(self) -> float:
"""Scaling for formation volume factor.
Returns
-------
FVF : float
Returns:
float: FVF
"""
return 1 - self.pressure_fracface / self.pressure_initial

Expand All @@ -160,9 +148,8 @@ def alpha_scaled(self, pseudopressure: ndarray) -> ndarray:
def fvf_scale(self) -> float:
"""Scaling for formation volume factor.
Returns
-------
FVF : float
Returns:
float: FVF
"""
return 1

Expand All @@ -189,11 +176,10 @@ def simulate(
pressure_fracface = np.full(len(time), self.pressure_fracface)
else:
if len(pressure_fracface) != len(time):
msg = (
raise ValueError(
"Pressure time series does not match time variable:"
f" {len(pressure_fracface)} versus {len(time)}"
)
raise ValueError(msg)
self.pressure_fracface = pressure_fracface
m_i = self.fluid.m_i
m_f = self.fluid.m_scaled_func(pressure_fracface)
Expand All @@ -208,12 +194,74 @@ def simulate(
b[0] = m_f[i] + self.alpha_scaled(m_f[i]) * m_f[i] * mesh_ratio
try:
alpha_scaled = self.alpha_scaled(b)
except ValueError as e:
msg = f"scaling failed where m_initial={m_i} m_fracface={m_f}"
raise ValueError(msg) from e
except ValueError:
raise ValueError(
f"scaling failed where m_initial={m_i} m_fracface={m_f}"
)
kt_h2 = mesh_ratio * alpha_scaled
"""
This is basically where the code has to be modified for long-time behavior.
We need a
"""
a_matrix = _build_matrix(kt_h2)
pseudopressure[i + 1], _ = sparse.linalg.bicgstab(a_matrix, b, atol=_ATOL)
pseudopressure[i + 1], _ = sparse.linalg.bicgstab(a_matrix, b)
self.pseudopressure = pseudopressure


@dataclass
class LongTimeSinglePhaseReservoir(SinglePhaseReservoir):
def simulate(
self, time: ndarray[float], pressure_fracface: ndarray[float] | None = None
):
"""Calculate simulation pressure over time.
Args
----
time : ndarray
times to solve for pressure
pressure_fracface : Iterable[float] | None, optional
pressure at frac-face over time. Defaults to None, which is no change
Raises
------
ValueError: wrong length changing pressure at frac-face
"""
self.time = time
dx_squared = (1 / self.nx) ** 2
pseudopressure = np.empty((len(time), self.nx))
if pressure_fracface is None:
pressure_fracface = np.full(len(time), self.pressure_fracface)
else:
if len(pressure_fracface) != len(time):
raise ValueError(
"Pressure time series does not match time variable:"
f" {len(pressure_fracface)} versus {len(time)}"
)
self.pressure_fracface = pressure_fracface
m_i = self.fluid.m_i
m_f = self.fluid.m_scaled_func(pressure_fracface)
pseudopressure_initial = np.full(self.nx, m_i)
pseudopressure_initial[0] = m_f[0]
pseudopressure[0, :] = pseudopressure_initial

for i in range(len(time) - 1):
mesh_ratio = (time[i + 1] - time[i]) / dx_squared
b = np.minimum(pseudopressure[i].copy(), m_i)
# Enforce the boundary condition at x=0
b[0] = m_f[i] + self.alpha_scaled(m_f[i]) * m_f[i] * mesh_ratio
try:
alpha_scaled = self.alpha_scaled(b)
except ValueError:
raise ValueError(
f"scaling failed where m_initial={m_i} m_fracface={m_f}"
)
kt_h2 = mesh_ratio * alpha_scaled
"""
This is basically where the code has to be modified for long-time behavior.
We need a
"""
a_matrix = _build_long_time_matrix(kt_h2)
pseudopressure[i + 1], _ = sparse.linalg.bicgstab(a_matrix, b)
self.pseudopressure = pseudopressure


Expand All @@ -224,13 +272,13 @@ class TwoPhaseReservoir(SinglePhaseReservoir):
References
----------
Ruiz Maraggi, L.M., Lake, L.W. and Walsh, M.P., 2020. "A Two-Phase Non-Linear One-
Dimensional Flow Model for Reserves Estimation in Tight Oil and Gas
Condensate Reservoirs Using Scaling Principles." In SPE Latin American and
Caribbean Petroleum Engineering Conference. OnePetro.
https://doi.org/10.2118/199032-MS
Dimensional Flow Model for Reserves Estimation in Tight Oil and Gas
Condensate Reservoirs Using Scaling Principles." In SPE Latin American and
Caribbean Petroleum Engineering Conference. OnePetro.
https://doi.org/10.2118/199032-MS
"""

Sw_init: float = 0.0
Sw_init: float

def simulate(self, time: ndarray):
"""Calculate simulation pressure over time.
Expand All @@ -245,32 +293,11 @@ def simulate(self, time: ndarray):

@dataclass
class MultiPhaseReservoir(SinglePhaseReservoir):
"""Reservoir with three phases: oil, gas, and water all flowing.
"""Reservoir with three phases: oil, gas, and water all flowing."""

Parameters
----------
nx : int
number of spatial nodes
pressure_fracface : float | NDArray
drawdown pressure at x=0 (psi)
pressure_initial : float
reservoir pressure before production (psi)
fluid : FlowProperties
reservoir fluid PVT/flow properties
So_init : float
oil saturation at the beginning
Sw_init : float
water saturation
Sg_init : float
gas saturation
"""

So_init: float = 1.0
"""Initial oil saturation."""
Sw_init: float = 0.0
"""Initial water saturation."""
Sg_init: float = 0.0
"""Initial gas saturation."""
So_init: float
Sg_init: float
Sw_init: float

def simulate(self, time: ndarray):
"""Calculate simulation pressure over time.
Expand Down Expand Up @@ -300,9 +327,7 @@ def simulate(self, time: ndarray):
alpha_scaled = self.alpha_scaled(b, sat)
kt_h2 = mesh_ratio * alpha_scaled
a_matrix = _build_matrix(kt_h2)
pseudopressure[i + 1], info = sparse.linalg.bicgstab(
a_matrix, b, atol=_ATOL
)
pseudopressure[i + 1], info = sparse.linalg.bicgstab(a_matrix, b)
saturations[i + 1] = self._step_saturation(sat, b, pseudopressure[i + 1])
self.time = time
self.pseudopressure = pseudopressure
Expand All @@ -319,7 +344,7 @@ def alpha_scaled(
----------
pseudopressure : ndarray
scaled pseudopressure
saturation : ndarray
saturaion : ndarray
record array with So, Sg, Sw records
Returns
Expand All @@ -332,25 +357,17 @@ def alpha_scaled(
return alpha(pseudopressure, s["So"], s["Sg"], s["Sw"]) / alpha(1)

def _step_saturation(
self,
saturation: ndarray, # noqa: ARG002
ppressure_old: ndarray, # noqa: ARG002
ppressure_new: ndarray, # noqa: ARG002
self, saturation: ndarray, ppressure_old: ndarray, ppressure_new: ndarray
) -> ndarray:
"""Calculate new saturation.
Args
----
saturation : ndarray
old saturation
ppressure_old : ndarray
old pressure
ppressure_new : ndarray
new pressure
Args:
saturation (ndarray): old saturation
ppressure_old (ndarray): old pressure
ppressure_new (ndarray): new pressure
Returns
-------
new saturation: ndarray
Returns:
ndarray: new saturation
"""
return NotImplementedError

Expand Down Expand Up @@ -382,3 +399,33 @@ def _build_matrix(kt_h2: ndarray) -> sparse.spmatrix:
[diagonal_low, diagonal_long, diagonal_upper], [-1, 0, 1], format="csr"
)
return a_matrix

def _build_long_time_matrix(kt_h2: ndarray) -> sparse.spmatrix:
r"""
This has to be modified for the new equation.
Set up A matrix for timestepping.
Follows :math: `A x = b` -> :math: `x = A \ b`
Parameters
----------
kt_h2: ndarray
diffusivity * dt / dx^2
Returns
-------
a_matrix: sp.sparse.matrix
The A matrix
"""
diagonal_long = 1.0 + 2 * kt_h2
# diagonal_long[0] = -1.0
# diagonal_low = np.concatenate([[0], -kt_h2[2:-1], [-2 * kt_h2[-1]]])
# diagonal_upper = np.concatenate([[0, -kt_h2[1]], -kt_h2[2:-1]])
diagonal_long[-1] = 1.0 + kt_h2[-1]
diagonal_low = -kt_h2[1:]
diagonal_upper = -kt_h2[0:-1]
a_matrix = sparse.diags(
[diagonal_low, diagonal_long, diagonal_upper], [-1, 0, 1], format="csr"
)
return a_matrix

0 comments on commit b70656c

Please sign in to comment.