drainage#
Simulate pressure-controlled injection of non-wetting fluid into an image using image-based sphere insertion [1].
import matplotlib.pyplot as plt
import numpy as np
import porespy as ps
ps.visualization.set_mpl_style()
im#
Can be a 2D or 3D image. To use the drainage function properly a capillary transform should be computed and provided to pc argument. If pc is not given, however, the function will still work as pc is computed assuming \(pc = 2/dt\).
im = ps.generators.blobs(shape=[200, 200], porosity=0.7, blobiness=1, seed=0)
drn = ps.simulations.drainage(im=im)
The function returns a Results object with a images containing the invasion sequence map , the saturation map, and capillary pressure map. These are shown below:
def plot_results(res):
fig, ax = plt.subplots(1, 3, figsize=[15, 4])
h = ax[0].imshow(res.im_seq / im, origin="lower")
cbar = fig.colorbar(h, aspect=10, format="{x:.0f}")
cbar.ax.set_title("Invasion Sequence", rotation=90, y=0, c="w")
ax[0].axis(False)
ax[0].set_title("(a) Sequence Map")
h = ax[1].imshow(res.im_snwp / im, origin="lower")
cbar = fig.colorbar(h, aspect=10, format="{x:.1f}")
cbar.ax.set_title("Non-Wetting Phase Saturation", rotation=90, y=0, c="w")
ax[1].axis(False)
ax[1].set_title("(b) Saturation Map")
h = ax[2].imshow(np.log10(res.im_pc) / im, origin="lower")
cbar = fig.colorbar(h, aspect=10, format="{x:.1f}")
cbar.ax.set_title("log(Capillary Pressure)", rotation=90, y=0, c="w")
ax[2].axis(False)
ax[2].set_title("(c) Capillary Pressure Map")
plot_results(drn)
inlets#
By default the drainage process begins from all faces/edges, but the inlets can be specified, for instance to inject from y=0:
inlets = ps.generators.faces(im.shape, inlet=1)
drn = ps.simulations.drainage(
im=im,
inlets=inlets,
)
plot_results(drn)
pc#
Supplying a capillary pressure transform as pc is how you tell drainage what the capillary entry pressure of each voxel is. This is calculated separately (e.g. using porespy.filters.capillary_transform) so that you have maximum control over how this is calculated, for instance to include gravity. In the following the effect of gravity can be noticed as the invasion preferentially occurs from the bottom left portion of the image (recall that the left edge is set as inlets)
pc = ps.filters.capillary_transform(
im=im,
sigma=0.01,
theta=180,
g=9.81,
rho_nwp=1000,
rho_wp=0,
voxel_size=1e-4,
)
drn = ps.simulations.drainage(
im=im,
pc=pc,
inlets=inlets,
)
plot_results(drn)
outlets#
Of outlets are specified then trapping of the wetting phase is included in the returned displacement maps.
outlets = ps.generators.faces(im.shape, outlet=1)
drn = ps.simulations.drainage(
im=im,
pc=pc,
inlets=inlets,
outlets=outlets,
)
plot_results(drn)
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
File ~/work/porespy/porespy/.venv/lib/python3.13/site-packages/numba/core/serialize.py:30, in _numba_unpickle(address, bytedata, hashed)
27 _unpickled_memo = {}
---> 30 def _numba_unpickle(address, bytedata, hashed):
31 """Used by `numba_unpickle` from _helperlib.c
32
33 Parameters
(...) 42 unpickled object
43 """
KeyboardInterrupt:
The above exception was the direct cause of the following exception:
SystemError Traceback (most recent call last)
Cell In[6], line 3
1 outlets = ps.generators.faces(im.shape, outlet=1)
----> 3 drn = ps.simulations.drainage(
4 im=im,
5 pc=pc,
6 inlets=inlets,
7 outlets=outlets,
8 )
10 plot_results(drn)
File ~/work/porespy/porespy/src/porespy/simulations/_drainage.py:621, in drainage(im, pc, dt, inlets, outlets, residual, steps, conn, min_size, smooth)
619 coords = np.where(edges) # Find (i, j, k) coordinates of edges
620 radii = dt[coords] # Extract sphere sizes to insert at each new location
--> 621 nwp_mask = _insert_disks_at_points_parallel(
622 im=nwp_mask,
623 coords=np.vstack(coords),
624 radii=radii.astype(int),
625 v=True,
626 smooth=smooth,
627 overwrite=False,
628 )
629 nwp_mask[seeds] = True # Fill in center in case spheres did not reach
630 # Connect residual to invasion front
SystemError: CPUDispatcher(<function _insert_disks_at_points_parallel at 0x7fc180ad7d80>) returned a result with an exception set
residual#
It is possible to specify the locations of any residual non-wetting phase which may be present prior to the fluid injection. These residual blobs impact the connecting of the invading fluid as well as the trapping of wetting phase:
residual = ps.filters.local_thickness(im) > 15
mask = ps.filters.trim_disconnected_voxels(
residual, inlets=ps.generators.borders(im.shape, mode="faces")
)
residual = residual * ~mask
fig, ax = plt.subplots(figsize=[4, 4])
ax.imshow(residual / im)
ax.axis(False)
ax.set_title("Locations of residual non-wetting phase")
drn = ps.simulations.drainage(
im=im,
pc=pc,
inlets=inlets,
residual=residual,
)
plot_results(drn)
min_size#
The min_size argument is passed to the find_trapped_regions function. It limits the minimum size of a cluster that can be considered trapped. For instance, in the above image there can be seen many small isolated trapped pixels next to the walls. These occur because of the digitized nature of the images, and are probably not realistic. Supplying min_size > 0 removes these pixels:
drn = ps.simulations.drainage(
im=im,
pc=pc,
inlets=inlets,
outlets=outlets,
min_size=5,
)
plot_results(drn)
Other Arguments#
dt#
Like many of the functions in PoreSpy, invasion optionally accepts a distance transform. This is computed automatically if dt is not provided, but providing one can save a small amount of time if the distance transform has already been computed (which is often the case).
return_sizes#
This flag indicates whether or not to compute the size map. Setting this to False can save some time if this result is not required.
steps#
Controls the pressure steps that are applied. If an integer is supplied than that number of pressure steps between the highest and lowest values in pc are used. If a list is supplied, then those values are used (in ascending order). If None is supplied, then every unique value in pc is used. This last option is necessary when comparing to the results of invasion for instance.