Invasion Percolation#

The next percolation algorithm to be demonstrated is known as Invasion Percolation. Instead of identifying connected clusters and invading them all in one go, as Ordinary Percolation does, this algorithm progresses one invasion step at a time. This is a more dynamic process and better simulates scenarios where instead of controlling the pressure at the network boundaries something else such as mass flow rate is controlled as the pressure is allowed to fluctuate up and down in order to meet the lowest available entry pressure for the growing cluster(s).

[1]:
import sys
import openpnm as op
import numpy as np
np.random.seed(10)
import matplotlib.pyplot as plt
import porespy as ps
from ipywidgets import interact, IntSlider
from openpnm.topotools import trim
%matplotlib inline
ws = op.Workspace()
ws.settings["loglevel"] = 50

In order to also showcase some other network generation options we first start with a small 2D network with SpheresAndCylinders geometry.

[2]:
spacing=2.5e-5
net = op.network.Cubic([20, 20, 1], spacing=spacing)
geo = op.geometry.SpheresAndCylinders(network=net, pores=net.Ps, throats=net.Ts)

We then trim all the surface pores to obtain disctint sets of boundary edge pores.

[3]:
net.labels()
net.num_throats('surface')
trim(network=net, throats=net.throats('surface'))
h = net.check_network_health()
trim(network=net, pores=h['trim_pores'])
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Input In [3], in <cell line: 5>()
      3 trim(network=net, throats=net.throats('surface'))
      4 h = net.check_network_health()
----> 5 trim(network=net, pores=h['trim_pores'])

KeyError: 'trim_pores'

Then we use a function from our porespy package to generate a tomography style image of the abstract network providing the number of pixels in each dimension.

[4]:
im = op.topotools.generate_voxel_image(net, max_dim=1000)
100%|██████████| 400/400 [00:00<00:00, 2458.12it/s]
100%|██████████| 684/684 [01:20<00:00,  8.47it/s]
[5]:
print(im.shape)
(1000, 1000, 34)

This creates a 3D image but we can crop it to get the central slice in 2D for visualization.

[6]:
fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(im[25:-25, 25:-25, 18].T)
[6]:
<matplotlib.image.AxesImage at 0x7ff39f22bdc0>
../../../_images/examples_simulations_percolation_B_invasion_percolation_11_1.png

Next the snow algorithm is used to do network extraction on the tomography style image. Of course if you have your own tomogrpahy image this can be used instead.

[7]:
crop = im[25:-25, 25:-25, :]
snow_out = ps.networks.snow2(crop > 0, voxel_size=4e-7)
[8]:
print(snow_out.regions.shape)
(956, 956, 23)

The SNOW algorithm provides a labelled region image containing the pore index. As zero is used for the background it is actually the pore index + 1 because python references arrays with first element as zero and we do not explicitly store the pore index.

[9]:
# NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots(figsize=(5, 5))
reg = snow_out.regions.astype(float) - 1
reg[reg == -1] = np.nan
region_slice = snow_out.regions[:, :, 18] - 1
mask = region_slice >= 0
ax.imshow(region_slice.T);
../../../_images/examples_simulations_percolation_B_invasion_percolation_16_0.png

Now our new network is extracted we can fill a network object with all the properties and begin simulation.

[10]:
wrk = op.Workspace()
wrk.clear()
[11]:
net, geo = op.io.PoreSpy.import_data(snow_out.network)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [11], in <cell line: 1>()
----> 1 net, geo = op.io.PoreSpy.import_data(snow_out.network)

AttributeError: module 'openpnm.io' has no attribute 'PoreSpy'

A helper function is defined for plotting a particular data set.

[12]:
def update_image(data):
    data = data.astype(float)
    out_im = np.ones(region_slice.shape, dtype=float)*-1
    out_im[mask] = data[region_slice[mask]]
    out_im[~mask] = np.nan
    return out_im
[13]:
# NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots(figsize=(5, 5))
out = update_image(net['pore.diameter'])
ax.imshow(out.T);
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [13], in <cell line: 3>()
      1 # NBVAL_IGNORE_OUTPUT
      2 fig, ax = plt.subplots(figsize=(5, 5))
----> 3 out = update_image(net['pore.diameter'])
      4 ax.imshow(out.T)

File ~/work/OpenPNM/OpenPNM/openpnm/network/_generic.py:138, in GenericNetwork.__getitem__(self, key)
    136     self._gen_ids()
    137     return self.get(f"{element}._id")
--> 138 vals = super().__getitem__(key)
    139 return vals

File ~/work/OpenPNM/OpenPNM/openpnm/core/_mixins.py:29, in ParamMixin.__getitem__(self, key)
     27         vals = self.network._params[key]
     28 else:
---> 29     vals = super().__getitem__(key)
     30 return vals

File ~/work/OpenPNM/OpenPNM/openpnm/core/_base.py:196, in Base.__getitem__(self, key)
    193 if key in self.keys():
    194     # Get values if present on self
    195     vals = super().__getitem__(key)
--> 196 elif key in self.keys(mode='all', deep=True):
    197     # Interleave values from geom if found there
    198     vals = self.interleave_data(key)
    199 elif any([k.startswith(key + '.') for k in self.keys()]):
    200     # Create a subdict of values present on self

File ~/work/OpenPNM/OpenPNM/openpnm/core/_base.py:450, in Base.keys(self, element, mode, deep)
    448             temp += item.keys(element=element, mode=mode, deep=False)
    449     if self._isa('network'):
--> 450         for item in self.project.geometries().values():
    451             temp += item.keys(element=element, mode=mode, deep=False)
    453 return temp

AttributeError: 'NoneType' object has no attribute 'geometries'
../../../_images/examples_simulations_percolation_B_invasion_percolation_22_1.png

Again, stadard physics is used to define the capillary entry pressures. And these are shown as a histogram for all the throats in the network.

[14]:
water = op.phases.Water(network=net)
phys = op.physics.Basic(network=net, geometry=geo, phase=water)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [14], in <cell line: 1>()
----> 1 water = op.phases.Water(network=net)
      2 phys = op.physics.Basic(network=net, geometry=geo, phase=water)

AttributeError: module 'openpnm' has no attribute 'phases'
[15]:
# NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots(figsize=[5, 5])
ax.hist(phys['throat.entry_pressure'], bins=10);
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [15], in <cell line: 3>()
      1 # NBVAL_IGNORE_OUTPUT
      2 fig, ax = plt.subplots(figsize=[5, 5])
----> 3 ax.hist(phys['throat.entry_pressure'], bins=10)

NameError: name 'phys' is not defined
../../../_images/examples_simulations_percolation_B_invasion_percolation_25_1.png

Next, the algorithm is defined and run with no arguments or outlets defined. This will proceed step by step assessing which pores are currently invaded (i.e. inlets first), which throats connect to an uninvaded pore and of these, which throat has the lowest capillary entry pressure for invasion. Invasion then proceeds along the path of least capillary resistance.

[16]:
# NBVAL_IGNORE_OUTPUT
alg_ip = op.algorithms.InvasionPercolation(network=net, phase=water)
alg_ip.set_inlets(pores=net.pores('xmin'))
alg_ip.run()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [16], in <cell line: 2>()
      1 # NBVAL_IGNORE_OUTPUT
----> 2 alg_ip = op.algorithms.InvasionPercolation(network=net, phase=water)
      3 alg_ip.set_inlets(pores=net.pores('xmin'))
      4 alg_ip.run()

NameError: name 'water' is not defined
[17]:
# NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots(figsize=(5, 5))
out = update_image(alg_ip['pore.invasion_sequence'])
plt.imshow(out.T);
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [17], in <cell line: 3>()
      1 # NBVAL_IGNORE_OUTPUT
      2 fig, ax = plt.subplots(figsize=(5, 5))
----> 3 out = update_image(alg_ip['pore.invasion_sequence'])
      4 plt.imshow(out.T)

NameError: name 'alg_ip' is not defined
../../../_images/examples_simulations_percolation_B_invasion_percolation_28_1.png
[18]:
def plot_invasion(seq):
    data = alg_ip['pore.invasion_sequence'] < seq
    fig, ax = plt.subplots(figsize=(5, 5))
    out = update_image(data)
    plt.imshow(out.T);

Using the slider below we can interactively plot the saturation at each invasion step (this works best using the left and right arrow keys).

[19]:
max_seq = alg_ip['pore.invasion_sequence'].max()
interact(plot_invasion, seq=IntSlider(min=0, max=max_seq, step=1, value=200));
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [19], in <cell line: 1>()
----> 1 max_seq = alg_ip['pore.invasion_sequence'].max()
      2 interact(plot_invasion, seq=IntSlider(min=0, max=max_seq, step=1, value=200))

NameError: name 'alg_ip' is not defined

As with Ordinary Percolation we can plot a drainage or intrusion curve but this time the capillary pressure is plotted from one step to the next as a continuous process with dynamic pressure boundary conditions and so is allowed to increase and decrease to meet the next lowest entry pressure for the invading cluster.

[20]:
fig, ax = plt.subplots(figsize=(5, 5))
alg_ip.plot_intrusion_curve(ax)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [20], in <cell line: 2>()
      1 fig, ax = plt.subplots(figsize=(5, 5))
----> 2 alg_ip.plot_intrusion_curve(ax)

NameError: name 'alg_ip' is not defined
../../../_images/examples_simulations_percolation_B_invasion_percolation_33_1.png

We can compare the results of the two algorithms and see that the pressure envelope, i.e. maximum pressure reached historically by the invasion process is the same as the ordinary percolation value.

[21]:
fig, ax = plt.subplots(figsize=(5, 5))
alg_op = op.algorithms.OrdinaryPercolation(network=net, phase=water)
alg_op.set_inlets(net.pores('xmin'))
alg_op.settings._update({'pore_volume': 'pore.volume',
                         'throat_volume': 'throat.volume'})
alg_op.run(points=1000)
alg_op.plot_intrusion_curve(ax)
alg_ip.plot_intrusion_curve(ax)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [21], in <cell line: 2>()
      1 fig, ax = plt.subplots(figsize=(5, 5))
----> 2 alg_op = op.algorithms.OrdinaryPercolation(network=net, phase=water)
      3 alg_op.set_inlets(net.pores('xmin'))
      4 alg_op.settings._update({'pore_volume': 'pore.volume',
      5                          'throat_volume': 'throat.volume'})

NameError: name 'water' is not defined
../../../_images/examples_simulations_percolation_B_invasion_percolation_35_1.png

An additional feature of the algorithm is the ability to identify where the defending phase becomes trapped. Whether this happens in reality in-fact relies on the connectivity of the defending phase and whether it can reside in the invaded pores as thin wetting films. If not then the defending phase is completely pushed out of a pore when invaded and it can become isolated and trapped when encircled by the invading phase. OpenPNM actually calculates this trapping as a post-process, employing some clever logic described by Masson 2016.

[22]:
alg_ip_t = op.algorithms.InvasionPercolation(network=net, phase=water)
alg_ip_t.set_inlets(pores=net.pores('xmin'))
alg_ip_t.run()
alg_ip_t.apply_trapping(outlets=net.pores(['boundary']))
fig, ax = plt.subplots(figsize=(5, 5))
out = update_image(alg_ip_t['pore.trapped'])
ax.imshow(out.T);
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [22], in <cell line: 1>()
----> 1 alg_ip_t = op.algorithms.InvasionPercolation(network=net, phase=water)
      2 alg_ip_t.set_inlets(pores=net.pores('xmin'))
      3 alg_ip_t.run()

NameError: name 'water' is not defined

Here a reasonable fraction of the pore space is not invaded due to trapping of the defending phase. Generally this fraction will be lower in truly 3D networks as there are more routes out of the network because pores have higher connectivity. Also, typically if a defending phase is considered to be wetting then film flow is assumed to allow residual defending phase to escape. However, we can show the differences on one plot with and without trapping below.

[23]:
fig, ax = plt.subplots(figsize=(5, 5))
alg_ip.plot_intrusion_curve(ax)
alg_ip_t.plot_intrusion_curve(ax)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [23], in <cell line: 2>()
      1 fig, ax = plt.subplots(figsize=(5, 5))
----> 2 alg_ip.plot_intrusion_curve(ax)
      3 alg_ip_t.plot_intrusion_curve(ax)

NameError: name 'alg_ip' is not defined
../../../_images/examples_simulations_percolation_B_invasion_percolation_39_1.png