Created Jupyter notebook for new code

Rather than working on new code in Python file I put the code to work on in a Jupyter notebook, Long-Time Flow Development.ipynb

.gitignore

@@ -145,3 +145,4 @@ dmypy.json
 # SPE files
 data/production_data.csv
 data/well_data.csv
+*.jupyterlab-workspace

docs/Long-Time Flow Development.ipynb

@@ -0,0 +1,193 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "4e0c34ef-53d0-4d78-964e-8c55dd8c8074",
+   "metadata": {},
+   "source": [
+    "# Long-time flow development space"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c05601c2-9bb1-4bee-914f-44f0005b5422",
+   "metadata": {},
+   "source": [
+    "This is where we can work on writing the code to implement long-time flow.\n",
+    "\n",
+    "Han-Ting was right. It's simpler using $r=1/s$. The equation becomes\n",
+    "$${\\partial \\tilde m\\over \\partial t}={\\alpha(\\tilde m)\\over\\alpha(\\tilde m_i)} \\left [ s^4\\,\\left({{\\partial^2 \\tilde m}\\over{\\partial\\,s^2}} \\right)+s^3\\,\\left({{\\partial}\\tilde m\\over{\\partial\\,s\n",
+    " }}\\,\\right)\n",
+    " \\right ]$$\n",
+    "\n",
+    " The boundary conditions are \n",
+    " $$ \\tilde m(s=1,t)=\\tilde m_b.$$\n",
+    " Here $\\tilde m_b=\\tilde m(p_b)$ is the pseudopressure due to the bottom-hole pressure. In the analytical work what comes out formally is that one should set this to the long-time limit of the bottom-hole pressure, so we should do that here and just use a constant.\n",
+    " $$\\tilde m(s=0,t)=\\tilde m_i.$$\n",
+    " This is the initial pseudopressure, $\\tilde m_i=\\tilde m(p_i)$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "bf618988-3730-46b9-b4d5-5102946486ee",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "The autoreload extension is already loaded. To reload it, use:\n",
+      "  %reload_ext autoreload\n"
+     ]
+    }
+   ],
+   "source": [
+    "%load_ext autoreload\n",
+    "%autoreload 2\n",
+    "\n",
+    "import numpy as np\n",
+    "import scipy as sp\n",
+    "from scipy.interpolate import interp1d\n",
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from bluebonnet.flow import (\n",
+    "    IdealReservoir,\n",
+    "    FlowProperties,\n",
+    "    SinglePhaseReservoir,\n",
+    "    RelPermParams,\n",
+    "    relative_permeabilities,\n",
+    ")\n",
+    "from bluebonnet.fluids.fluid import Fluid, pseudopressure\n",
+    "from bluebonnet.plotting import (\n",
+    "    plot_pseudopressure,\n",
+    "    plot_recovery_factor,\n",
+    "    plot_recovery_rate,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "82e8c67d-e7c7-42bd-9904-8767792ea26e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class LongTimeSinglePhaseReservoir(SinglePhaseReservoir):\n",
+    "    def simulate(\n",
+    "        self, time: ndarray[float], pressure_fracface: ndarray[float] | None = None\n",
+    "    ):\n",
+    "        \"\"\"Calculate simulation pressure over time.\n",
+    "\n",
+    "        Args\n",
+    "        ----\n",
+    "        time : ndarray\n",
+    "            times to solve for pressure\n",
+    "        pressure_fracface : Iterable[float] | None, optional\n",
+    "            pressure at frac-face over time. Defaults to None, which is no change\n",
+    "\n",
+    "        Raises\n",
+    "        ------\n",
+    "        ValueError: wrong length changing pressure at frac-face\n",
+    "        \"\"\"\n",
+    "        self.time = time\n",
+    "        dx_squared = (1 / self.nx) ** 2\n",
+    "        pseudopressure = np.empty((len(time), self.nx))\n",
+    "        if pressure_fracface is None:\n",
+    "            pressure_fracface = np.full(len(time), self.pressure_fracface)\n",
+    "        else:\n",
+    "            if len(pressure_fracface) != len(time):\n",
+    "                raise ValueError(\n",
+    "                    \"Pressure time series does not match time variable:\"\n",
+    "                    f\" {len(pressure_fracface)} versus {len(time)}\"\n",
+    "                )\n",
+    "            self.pressure_fracface = pressure_fracface\n",
+    "        m_i = self.fluid.m_i\n",
+    "        m_f = self.fluid.m_scaled_func(pressure_fracface)\n",
+    "        pseudopressure_initial = np.full(self.nx, m_i)\n",
+    "        pseudopressure_initial[0] = m_f[0]\n",
+    "        pseudopressure[0, :] = pseudopressure_initial\n",
+    "\n",
+    "        for i in range(len(time) - 1):\n",
+    "            mesh_ratio = (time[i + 1] - time[i]) / dx_squared\n",
+    "            b = np.minimum(pseudopressure[i].copy(), m_i)\n",
+    "            # Enforce the boundary condition at x=0\n",
+    "            b[0] = m_f[i] + self.alpha_scaled(m_f[i]) * m_f[i] * mesh_ratio\n",
+    "            try:\n",
+    "                alpha_scaled = self.alpha_scaled(b)\n",
+    "            except ValueError:\n",
+    "                raise ValueError(\n",
+    "                    f\"scaling failed where m_initial={m_i}  m_fracface={m_f}\"\n",
+    "                )\n",
+    "            kt_h2 = mesh_ratio * alpha_scaled\n",
+    "            \"\"\"\n",
+    "            This is basically where the code has to be modified for long-time behavior.\n",
+    "            We need a\n",
+    "            \"\"\"\n",
+    "            a_matrix = _build_long_time_matrix(kt_h2)\n",
+    "            pseudopressure[i + 1], _ = sparse.linalg.bicgstab(a_matrix, b)\n",
+    "        self.pseudopressure = pseudopressure"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7bedfcae-fa29-49dd-bd8f-a78df998e04e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def _build_long_time_matrix(kt_h2: ndarray) -> sparse.spmatrix:\n",
+    "    r\"\"\"\n",
+    "    This has to be modified for the new equation.\n",
+    "\n",
+    "    Set up A matrix for timestepping.\n",
+    "\n",
+    "    Follows :math: `A x = b` -> :math: `x = A \\ b`\n",
+    "\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    kt_h2: ndarray\n",
+    "        diffusivity * dt / dx^2\n",
+    "\n",
+    "    Returns\n",
+    "    -------\n",
+    "    a_matrix: sp.sparse.matrix\n",
+    "        The A matrix\n",
+    "    \"\"\"\n",
+    "    diagonal_long = 1.0 + 2 * kt_h2\n",
+    "    # diagonal_long[0] = -1.0\n",
+    "    # diagonal_low = np.concatenate([[0], -kt_h2[2:-1], [-2 * kt_h2[-1]]])\n",
+    "    # diagonal_upper = np.concatenate([[0, -kt_h2[1]], -kt_h2[2:-1]])\n",
+    "    diagonal_long[-1] = 1.0 + kt_h2[-1]\n",
+    "    diagonal_low = -kt_h2[1:]\n",
+    "    diagonal_upper = -kt_h2[0:-1]\n",
+    "    a_matrix = sparse.diags(\n",
+    "        [diagonal_low, diagonal_long, diagonal_upper], [-1, 0, 1], format=\"csr\"\n",
+    "    )\n",
+    "    return a_matrix"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.18"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/flow.ipynb

@@ -378,9 +378,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "bluebonnet",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
-   "name": "bluebonnet"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -392,7 +392,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.13"
+   "version": "3.8.18"
   },
   "vscode": {
    "interpreter": {

src/bluebonnet/flow/reservoir.py

@@ -209,62 +209,6 @@ def simulate(
 
 
 @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
-
-
 @dataclass
 class TwoPhaseReservoir(SinglePhaseReservoir):
     """Oil-gas reservoir simulation.
@@ -400,32 +344,3 @@ def _build_matrix(kt_h2: ndarray) -> sparse.spmatrix:
     )
     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