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()

[17:49:25] 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

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

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

out = ps.filters.ibip(im=im, inlets=bd, maxiter=15000)
inv_seq, inv_size = out.inv_sequence, out.inv_sizes

/tmp/ipykernel_9267/3299454374.py:1: DeprecationWarning: Call to deprecated function (or staticmethod) ibip. (The ibip function will be moved to the ``simulations`` module in a future version)
  out = ps.filters.ibip(im=im, inlets=bd, maxiter=15000)

Do some post-processing

Convert invasion sequence values to saturation, and identify trapped regions:

inv_satn = ps.filters.seq_to_satn(seq=inv_seq)
inv_seq_trapping = ps.filters.find_trapped_regions(seq=inv_seq, bins=None,
                                                   return_mask=False)
inv_satn_trapping = ps.filters.seq_to_satn(seq=inv_seq_trapping)

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.

sizes = np.arange(int(dt.max())+1, 0, -1)
mio = ps.filters.porosimetry(im=im, inlets=bd, sizes=sizes, mode='mio')
mio_seq = ps.filters.size_to_seq(mio)
mio_seq[im*(mio_seq == 0)] = -1  # Adjust to set uninvaded to -1
mio_satn = ps.filters.seq_to_satn(mio_seq)
mio_seq_trapping = ps.filters.find_trapped_regions(seq=mio_seq, bins=None,
                                                   return_mask=False)

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 fields for both IBIP and MIO.

Saturation field 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')
ax[1].plot(satns, err, 'o-');

Comparing Capillary Pressure Curves

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

d = ps.metrics.pc_curve_from_ibip(im=im, seq=inv_seq, sizes=inv_size, voxel_size=1e-5);
e = ps.metrics.pc_curve(im=im, sizes=mio, voxel_size=1e-5);

/tmp/ipykernel_9267/1817964060.py:1: DeprecationWarning: Call to deprecated function (or staticmethod) pc_curve_from_ibip. (This function is deprecated, use pc_curve instead)
  d = ps.metrics.pc_curve_from_ibip(im=im, seq=inv_seq, sizes=inv_size, voxel_size=1e-5);

fig, ax = plt.subplots()
ax.plot(np.log10(d.pc), d.snwp, 'g-', linewidth=2, label='IBIP')
ax.step(np.array(np.log10(e.pc)), e.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')

The above capillary pressure curves are quite interesting. The MIO curve in red is the envelop of peak pressures obtained by the IPIP curve 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.

Compare results with trapping

One of the key benefits of the IBIP algorithm is that the fluid invades incrememtally instead of in steps. This means the appearance of trapped defending phase is immediately noted. By contrast, the MIO approach computes large changes in saturation, so any trapping that may occur between these large jumps is not resolved. Therefore, only the IBIP method can correctly identify trapped defending phase. This is illustrated below:

cmap = copy(plt.cm.turbo)
cmap.set_over(color='white')
cmap.set_under(color='grey')
fig, ax = plt.subplots(1, 3, figsize=[8, 4])
ax[0].imshow((inv_seq_trapping > 0)/im, vmin=1e-3, vmax=2,
             cmap=cmap, origin='lower', interpolation='none')
ax[0].axis('off')
ax[1].imshow((mio_seq_trapping > 0)/im, vmin=1e-3, vmax=1,
             cmap=cmap, origin='lower', interpolation='none')
ax[1].axis('off')
temp = (mio_seq_trapping > 0)*2.0 - (inv_seq_trapping > 0)*1.0
ax[2].imshow(temp/im, vmin=1e-3, vmax=2.0, 
             cmap=cmap, origin='lower', interpolation='none')
ax[2].axis('off');

The red color in the right image indicates regions that were invaded by the MIO algorithm, but not by the IBIP algorithm. The MIO will always underestimate trapping.

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)