Image-based invasion percoalation

This notebook illustrates the use of the image-based invasion percolation algorithm decribed here.

Import necessary packages and functions

import porespy as ps
import numpy as np
import matplotlib.pyplot as plt
from edt import edt
import imageio
from copy import copy
ps.visualization.set_mpl_style()

[01:07:55] ERROR    PARDISO solver not installed, run `pip install pypardiso`. Otherwise,          _workspace.py:56
                    simulations will be slow. Apple M chips not supported.                                         

Generate test image

im = ps.generators.blobs(shape=[500, 500], porosity=0.6, blobiness=2, seed=5)

Also generate an image defining the inlets

bd = np.zeros_like(im, dtype=bool)
bd[:, 0] = 1
bd *= im
im = ps.filters.trim_disconnected_blobs(im=im, inlets=bd)
dt = edt(im)

Apply IBIP to image

ibip = ps.simulations.qbip(im=im, inlets=bd, maxiter=25000, return_sizes=True)
inv_satn = ps.filters.seq_to_satn(seq=ibip.im_seq)

Exiting after 10355 steps

Run standard MIO for comparison:

Compute the MIO result, convert to sequation and saturation using the functions in porespy.tools, and also find trapped regions.

pc = ps.filters.capillary_transform(im=im)
outlets = ps.generators.faces(shape=im.shape, outlet=0)
sizes = np.unique(pc[pc != 0])
mio = ps.simulations.drainage(
    im=im, 
    inlets=bd, 
    bins=sizes, 
    return_sizes=True,
)
mio_seq = ps.filters.satn_to_seq(mio.im_satn)
mio_seq[im*(mio_seq == 0)] = -1  # Adjust to set uninvaded to -1
mio_satn = ps.filters.seq_to_satn(mio_seq)

Comparing the invasion configurations between the MIO and IBIP algorithms below suggests that the methods produce similar results when trapping is turned off. In fact, it can be shown that the MIO result is a subset of the IBIP result.

Compare MIO and IBIP at Equal Saturations

The comparison between the two methods can only be done when trapping is neglected. The difference in the trapped defending phase from each method is discussed further below.

The images below show the saturation maps for both IBIP and MIO.

Saturation map means that each voxel in the image contains the value of the saturation at the point is was invaded. In other words, the operation satn <= s will yeild and image containing the fluid configuration at saturation equal to s.

fig, ax = plt.subplots(1, 2, figsize=[8, 4])
ax[0].imshow(inv_satn/im, origin='lower', interpolation='none')
ax[1].imshow(mio_satn/im, origin='lower', interpolation='none');

The saturation fields can be compared directly by seeing which voxels are filled at each saturation in both methods. We must first find all the values of saturation that were obtained in the MIO result, which consists of discrete jumps (i.e. from s=0.33 to s=0.49) due to the pressure based nature of the simulations. Once we have a list of all unique saturation values, we can threshold the saturation fields from both simulations at these values and compare the configurations.

inv_satn_t = np.around(inv_satn, decimals=4)
mio_satn_t = np.around(mio_satn, decimals=4)
satns = np.unique(mio_satn_t)[1:]
err = []
diff = np.zeros_like(im, dtype=float)
for s in satns:
    ip_mask = (inv_satn_t <= s) * (inv_satn_t > 0)
    mio_mask = (mio_satn_t <= s) * (mio_satn_t > 0)
    diff[(mio_mask == 1)*(ip_mask == 0)*(im == 1)] = 1
    diff[(mio_mask == 0)*(ip_mask == 1)*(im == 1)] = -1
    err.append((mio_mask != ip_mask).sum())

In the figure below we plot the differences between the two results at equal saturation. On the left is the fluid configurations of each method overlaid on top of each other. Any voxels that are not green indicate locations where the two results disagree. There are very few of them and almost impossible to see (e.g. near [x, y]=[250, 50]), so the plot on the right shows the count of mis-matching voxels as a function of saturation. The total count never exceeds 70, which is a very small fraction of a 500 x 500 image. The differences are due to the numerical precision of the saturation values.

fig, ax = plt.subplots(1, 2, figsize=[8, 4])
ax[0].imshow(diff/im, origin='lower', vmin=-1, vmax=1)
ax[1].plot(satns, err, 'o-')
ax[1].set_xlabel('Saturation')
ax[1].set_ylabel('Number of Pixels');

Comparing Capillary Pressure Curves

Converting the invasion configurations to capillary pressure curves is possible using functions provided in the metrics module.

ibip.im_pc = ps.filters.size_to_pc(
    im=im, 
    size=ibip.im_size,
    sigma=1, 
    theta=180,
    voxel_size=1.0,
)
pc_curve_ibip = ps.metrics.pc_map_to_pc_curve(
    pc=ibip.im_pc,
    im=im,
    seq=ibip.im_seq,
)

fig, ax = plt.subplots()
ax.plot(np.log10(pc_curve_ibip.pc), pc_curve_ibip.snwp, 'g-', linewidth=2, label='IBIP')
ax.step(np.array(np.log10(pc_curve_mio.pc)), pc_curve_mio.snwp, 'r--', where='post', markersize=20, linewidth=3, alpha=0.6, label='MIO')
plt.xlabel('log(Capillary Pressure [Pa])')
plt.ylabel('Non-wetting Phase Saturation')
plt.legend()
ax.xaxis.grid(True, which='both')

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[10], line 3
      1 fig, ax = plt.subplots()
      2 ax.plot(np.log10(pc_curve_ibip.pc), pc_curve_ibip.snwp, 'g-', linewidth=2, label='IBIP')
----> 3 ax.step(np.array(np.log10(pc_curve_mio.pc)), pc_curve_mio.snwp, 'r--', where='post', markersize=20, linewidth=3, alpha=0.6, label='MIO')
      4 plt.xlabel('log(Capillary Pressure [Pa])')
      5 plt.ylabel('Non-wetting Phase Saturation')

NameError: name 'pc_curve_mio' is not defined

The above capillary pressure curves are quite interesting. The MIO curve in red is the envelop of peak pressures obtained by IBIP (in green).

The red line corresponds to the normal ‘porosimetry’ experiment where increasingly high pressures are applied (x-axis) and larger fractions of the void space are filled (y-axis). Each verticle jump in the red curve corresponds to a new region of the image being penetrated because the breakthrough pressure of a constriction was exceed. All the void space beyond this constriction is then filled at once.

The green line is the result of IBIP and corresponds to volume controlled injection, where fluid is injected at some constant rate. Every data point has a slightly higher saturation than the previous so there are no plateaus of constant saturation. However, the capillary pressure (x-axis) fluctuates up and down seemlingly randomly. This is because the advancing menisci squeeze through small constrictions (higher pressure), then enter into a larger region that was sheilded by the constriction to they expand (lower pressure). This process repeats itself through the entire simulation.

Make a movie of the invasion

The image below is a sample of the animated output:

Obtaining a movie such as the one above can be accomplished with something like the following, which creates an interactive widget in a jupyter notebook.

from IPython.display import HTML
%matplotlib inline
plt.rcParams["animation.html"] = "jshtml"
plt.rcParams["animation.embed_limit"] = 100.0
plt.rcParams["figure.dpi"] = 150
mov = ps.visualization.satn_to_movie(im=im, satn=inv_satn)