class InvasionPercolation(*args, **kwargs)[source]

Bases: openpnm.algorithms.GenericAlgorithm.GenericAlgorithm

A classic/basic invasion percolation algorithm optimized for speed.


network (OpenPNM Network object) – The Network upon which the invasion will occur.


This algorithm uses a binary heap to store all a list of all accessible throats, sorted according to entry pressure. This means that item [0] in the heap is the most easily invaded throat, so looking up which throat to invade next is computationally trivial. In order to keep the list sorted new throats to the list takes more time, however, the heap data structure is very efficient at this. Interested users can consult the wikipedia page on binary heaps for more information.


Start by importing the usual packages:

>>> import openpnm as op
>>> import scipy as sp
>>> import matplotlib.pyplot as plt

Create 2D cubic network for easier visualizaiton:

>>> S = np.array([100, 100, 1])
>>> pn =, spacing=0.0001, name='pn11')

Add a basic geometry:

>>> geom = op.geometry.StickAndBall(network=pn, pores=pn.Ps, throats=pn.Ts)

Create an invading phase, and attach the capillary pressure model:

>>> water = op.phases.Water(network=pn)
>>> water.add_model(propname='throat.entry_pressure',
...                 model=op.models.physics.capillary_pressure.washburn)

Initialize an invasion percolation object and define inlets:

>>> ip = op.algorithms.InvasionPercolation(network=pn)
>>> ip.setup(phase=water)
>>> ip.set_inlets(pores=0)

After running the algorithm the invading phase configuration at a given saturation can be obtained and assigned to the phase object:

>>> water.update(ip.results(Snwp=0.5))

Because it was a 2D network it’s easy to quickly visualize the invasion pattern as an image for verification:


Because the network is 2D and cubic, an image can be generated with color corresponding to a value. The following plots the entire invasion sequence, and the water configuraiton at Snwp = 0.5.

plt.subplot(1, 2, 1)

plt.imshow(np.reshape(ip['pore.invasion_sequence'], newshape=S[S > 1]))

plt.subplot(1, 2, 2)

plt.imshow(np.reshape(water['pore.occupancy'], newshape=S[S > 1]))


Apply trapping based on algorithm described by Y. Masson [1]. It is applied as a post-process and runs the percolation algorithm in reverse assessing the occupancy of pore neighbors. Consider the following scenario when running standard IP without trapping, three situations can happen after each invasion step:

  • The number of defending clusters stays the same and clusters can shrink

  • A cluster of size one is suppressed

  • A cluster is split into multiple clusters

In reverse the following opposite situations can happen:

  • The number of defending clusters stays the same and clusters can grow

  • A cluster of size one is created

  • Mutliple clusters merge into one cluster

With trapping the reversed rules are adjusted so that only clusters that do not connect to a sink can grow and merge. At the point that a neighbor connected to a sink is touched the trapped cluster stops growing as this is the point of trapping in forward invasion time.

Logger info displays the invasion sequence and pore index and a message with condition number based on the modified trapping rules and the assignment of the pore to a given cluster.

Initially all invaded pores are given cluster label -1 Outlets / Sinks are given -2 New clusters that grow into fully trapped clusters are either identified at the point of breakthrough or grow from nothing if the full invasion sequence is run, they are assigned numbers from 0 up.

Ref: [1] Masson, Y., 2016. A fast two-step algorithm for invasion percolation with trapping. Computers & Geosciences, 90, pp.41-48

  • outlets (list or array of pore indices for defending fluid to escape) –

  • through


  • Creates a throat array called ‘pore.clusters’ in the Algorithm

  • dictionary. Any positive number is a trapped cluster

  • Also creates 2 boolean arrays Np and Nt long called ‘<element>.trapped’


Get the percolation data as the invader volume or number fraction vs the capillary capillary pressure.


Plot the percolation curve as the invader volume or number fraction vs the capillary capillary pressure.


Returns the phase configuration at the specified non-wetting phase (invading phase) saturation.


Snwp (scalar, between 0 and 1) – The network saturation for which the phase configuration is desired.


  • Two dictionary containing arrays that describe the pore and throat

  • distribution at the given saturation. Specifically, these are

  • **’pore.occupancy’** (1 indicates the pores is invaded and 0)

  • otherwise.

  • **’throat.occupancy’** (Same as described above but for throats.)


Perform the algorithm


n_steps (int) – The number of throats to invaded during this step

set_inlets(pores=[], overwrite=False)[source]

pores (array_like) – The list of inlet pores from which the Phase can enter the Network

setup(phase, entry_pressure='', pore_volume='', throat_volume='')[source]

Set up the required parameters for the algorithm

  • phase (OpenPNM Phase object) – The phase to be injected into the Network. The Phase must have the capillary entry pressure values for the system.

  • entry_pressure (string) – The dictionary key to the capillary entry pressure. If none is supplied then the current value is retained. The default is ‘throat.capillary_pressure’.

  • pore_volume (string) – The dictionary key to the pore volume. If none is supplied then the current value is retained. The default is ‘pore.volume’.

  • throat_volume (string) – The dictionary key to the throat volume. If none is supplied then the current value is retained. The default is ‘throat.volume’.