examples/notebooks/networks/extraction/doddington_ICL_benchmark.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:37.100393Z",
     "iopub.status.busy": "2021-06-24T11:32:37.098780Z",
     "iopub.status.idle": "2021-06-24T11:32:37.808991Z",
     "shell.execute_reply": "2021-06-24T11:32:37.807549Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pickle as pk\n",
    "import openpnm as op\n",
    "%config InlineBackend.figure_formats = ['svg']\n",
    "import matplotlib.pyplot as plt\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "This example is from our series of permeability benchmarks, mostly for making sure that `OpenPNM` can accurately predict the permeability of known rock samples. In this example, we try to calculate the permeability of an extracted pore network model form a Bentheimer sandstone rock sample. The network has already been extracted using `PoreSpy`, and so here we only load the output generated from `PoreSpy` into `OpenPNM`.\n",
    "\n",
    "To load extracted networks from `PoreSpy` into `OpenPNM`, it's highly recommended to use `OpenPNM`'s exclusive IO class for `PoreSpy`. Here's how it works: (basically the IO class returns a project object)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:37.815817Z",
     "iopub.status.busy": "2021-06-24T11:32:37.814387Z",
     "iopub.status.idle": "2021-06-24T11:32:37.914524Z",
     "shell.execute_reply": "2021-06-24T11:32:37.915692Z"
    }
   },
   "outputs": [],
   "source": [
    "fname = '../../_fixtures/Extracted-Networks/Doddington'\n",
    "prj = op.io.PoreSpy.load(fname)\n",
    "net = prj.network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to assigning a network, the `PoreSpy` IO class automatically assigns a geometry using the loaded data as well. To fetch the geometry, we can do that by asking the project object, but we need to know the geometry name. If nothing is provided by the user, `OpenPNM` usually automatically assign names to objects in a sequential manner. For instance, the first assigned network or geometry is named \"net_01\", \"geo_01\", respectively. But let's say you didn't know about this. Let's ispect the project object and see what geometries are attached to it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:37.923603Z",
     "iopub.status.busy": "2021-06-24T11:32:37.922144Z",
     "iopub.status.idle": "2021-06-24T11:32:37.929229Z",
     "shell.execute_reply": "2021-06-24T11:32:37.930370Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'geo_01': <openpnm.geometry.Imported object at 0x7f982d3fe6d0>}\n"
     ]
    }
   ],
   "source": [
    
    "print(prj.geometries())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! Now, we know that there's only a singly geometry attached to the project object, which can be accessed via the key \"geo_01\". Let's fetch the geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:37.938513Z",
     "iopub.status.busy": "2021-06-24T11:32:37.937660Z",
     "iopub.status.idle": "2021-06-24T11:32:37.942466Z",
     "shell.execute_reply": "2021-06-24T11:32:37.941435Z"
    }
   },
   "outputs": [],
   "source": [
    "geo = prj.geometries()[\"geo_01\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's take a look at the network using `ParaView`:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://user-images.githubusercontent.com/14086031/79818906-dc9e2800-8356-11ea-9ab1-96e7055e15b5.png\" style=\"width: 65%\" align=\"left\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-processing\n",
    "Now, let's run `StokesFlow` algorithm on this network to calculate its permeability. Usually, in the process of network extraction, some isolated pores could exist in the final network. Having isolated pores in your network rightfully confuses `OpenPNM` as gives rise to `NaNs` in your solution. We need to get rid of those pores. To do that, first we call `check_network_health` on the network, which returns back a dictionary containing useful information regarding network health. For now, we're only interested in the the \"trim_pores\" key in that dictionary, which contains a list of pores that need to be trimmed off of the network in order for the pores to be fully connected. Here's how we do this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:37.951381Z",
     "iopub.status.busy": "2021-06-24T11:32:37.950112Z",
     "iopub.status.idle": "2021-06-24T11:32:37.980362Z",
     "shell.execute_reply": "2021-06-24T11:32:37.981442Z"
    }
   },
   "outputs": [],
   "source": [
    "# Look for disconnected pores and trim them\n",
    "net_health = net.check_network_health()\n",
    "op.topotools.trim(network=net, pores=net_health[\"trim_pores\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Permeability calculation\n",
    "Now, let's create our `StokesFlow` algorithm object, impose a constant pressure drop across the y-axis, and calculate the y-axis permeability based on the resulting flow rate.\n",
    "\n",
    "Note that first things first! For `OpenPNM` to be able to run any algorithm, we need to define a `phase` object and a `physics` object. More precisely, you need as many `physics` objects as your `geometry` objects. For instance, if your network consists of  3 geometries, you need to define 3 separate `physics` objects. Here, since there's only one `geometry` object, we only need to create 1 `physics` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:37.991757Z",
     "iopub.status.busy": "2021-06-24T11:32:37.990308Z",
     "iopub.status.idle": "2021-06-24T11:32:38.454607Z",
     "shell.execute_reply": "2021-06-24T11:32:38.455534Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Permeability (Doddington sandstone): 2.6072 Darcy\n"
     ]
    }
   ],
   "source": [
    "# Add phase and physics\n",
    "air = op.phases.Air(network=net)\n",
    "phys_air = op.physics.Standard(network=net, phase=air, geometry=geo)\n",
    "\n",
    "# Perform StokesFlow algorithm to calculate permeability\n",
    "perm = op.algorithms.StokesFlow(network=net)\n",
    "perm.setup(phase=air)\n",
    "perm.set_value_BC(pores=net.pores('left'), values=0)\n",
    "perm.set_value_BC(pores=net.pores('right'), values=101325)\n",
    "perm.run()\n",
    "air.update(perm.results())\n",
    "\n",
    "# Original image shape (tomography data) needed to calculate cross-sectional area and length\n",
    "shape = np.array([1000, 1000, 1000])\n",
    "voxel_size = 2.6969e-6\n",
    "area = shape[[0, 2]].prod() * voxel_size**2\n",
    "length = shape[1] * voxel_size\n",
    "K = perm.calc_effective_permeability(domain_area=area, domain_length=length)\n",
    "# 1e12 is the conversion factor to convert from SI to darcy units\n",
    "print(f\"Permeability (Doddington sandstone): {K[0]*1e12:.4f} Darcy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Effective diffusivity calculation\n",
    "Now, let's create our `FickianDiffusion` algorithm object, impose a constant concentration drop across the y-axis, and calculate the y-axis effective diffusivity based on the resulting mass flow rate.\n",
    "\n",
    "Since we've already defined `phase` and `physics` in the \"Permeability calculation\" section, we can proceed with creating the `FickianDiffusion` algorithm part. Also, we've already calculated the cross-sectional area and length, so no need to do that part again.\n",
    "\n",
    "**Important note**: Since in the previous step we solved Stokes algorithm, which returned back a pressure distribution, our phase object, i.e. `air`, contains pressure values. Since we eventually want to normalize the calculated effective diffusivity by that of our working fluid, i.e. air, we need to reset the pressure profile back to 1 atm. This is because the diffusivity model that is used for `air` is the Fuller diffusivity equation which is pressure-dependent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:38.460434Z",
     "iopub.status.busy": "2021-06-24T11:32:38.459074Z",
     "iopub.status.idle": "2021-06-24T11:32:38.600173Z",
     "shell.execute_reply": "2021-06-24T11:32:38.601148Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Effective diffusivity (Doddington sandstone): 1.1566e-06 m^2/s\n"
     ]
    }
   ],
   "source": [
    "# First things first! Let's reset the pressure profile back to 1 atm\n",
    "air[\"pore.pressure\"] = 101325\n",
    "# Now we need to regenearte the diffusivity model to reflect the pressure reset\n",
    "air.regenerate_models(propnames=\"pore.diffusivity\") \n",
    "\n",
    "# Perform FickianDiffusion algorithm to calculate effective diffusivity\n",
    "fd = op.algorithms.FickianDiffusion(network=net)\n",
    "fd.setup(phase=air)\n",
    "fd.set_value_BC(pores=net.pores('left'), values=0)\n",
    "fd.set_value_BC(pores=net.pores('right'), values=1)\n",
    "fd.run()\n",
    "air.update(fd.results())\n",
    "\n",
    "D_eff = fd.calc_effective_diffusivity(domain_area=area, domain_length=length)\n",
    "print(f\"Effective diffusivity (Doddington sandstone): {D_eff[0]:.4e} m^2/s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also normalize the effective diffusivity with respect to that of our working fluid, which is air for our case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:38.605591Z",
     "iopub.status.busy": "2021-06-24T11:32:38.604422Z",
     "iopub.status.idle": "2021-06-24T11:32:38.608306Z",
     "shell.execute_reply": "2021-06-24T11:32:38.609378Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative diffusivity (Doddington sandstone): 0.056\n"
     ]
    }
   ],
   "source": [
    "D_air = air[\"pore.diffusivity\"].mean()\n",
    "print(f\"Relative diffusivity (Doddington sandstone): {D_eff[0]/D_air:.3f}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}