diff --git a/src/bluebonnet/flow/reservoir.py b/src/bluebonnet/flow/reservoir.py index 1f95dfa..9951fd6 100644 --- a/src/bluebonnet/flow/reservoir.py +++ b/src/bluebonnet/flow/reservoir.py @@ -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 @@ -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 @@ -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): @@ -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 @@ -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 @@ -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: @@ -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 @@ -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 @@ -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) @@ -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 @@ -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. @@ -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. @@ -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 @@ -319,7 +344,7 @@ def alpha_scaled( ---------- pseudopressure : ndarray scaled pseudopressure - saturation : ndarray + saturaion : ndarray record array with So, Sg, Sw records Returns @@ -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 @@ -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