Import an Extracted Network and Predict Transport Properties#
This example illustrates the process of both extracting a pore network from an image (that has already been binarized), then opening this image with OpenPNM to perform some simulations. The results of the simulations are compared to the known transport properties of the image and good agreement is obtained.
import openpnm as op
import porespy as ps
import numpy as np
from pathlib import Path
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 2
1 import openpnm as op
----> 2 import porespy as ps
3 import numpy as np
4 from pathlib import Path
ModuleNotFoundError: No module named 'porespy'
Load a small image of Berea sandstone:
pkg_root = Path(op.__file__).resolve().parents[2]
path = pkg_root / "tests/fixtures/berea_100_to_300.npz"
data = np.load(path.resolve())
im = data['im']
Note meta data for this image, which we’ll use to compare the network predictions too later:
data = {
'shape': {
'x': im.shape[0],
'y': im.shape[1],
'z': im.shape[2],
},
'resolution': 5.345e-6,
'porosity': 19.6,
'permeability': {
'Kx': 1360,
'Ky': 1304,
'Kz': 1193,
'Kave': 1286,
},
'formation factor': {
'Fx': 23.12,
'Fy': 23.99,
'Fz': 25.22,
'Fave': 24.08,
},
}
Perform extraction using snow2 from PoreSpy.
snow = ps.networks.snow2(
phases=im,
voxel_size=data['resolution'],
boundary_width=[3, 0, 0],
accuracy='standard',
parallel_kw=None,
)
/Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/filters/_snows.py:446: FutureWarning: `cube` is deprecated since version 0.25 and will be removed in version 0.27. Use `skimage.morphology.footprint_rectangle` instead.
peaks_dil = spim.binary_dilation(input=peaks_i, structure=cube(3))
/Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/filters/_snows.py:594: FutureWarning: `cube` is deprecated since version 0.25 and will be removed in version 0.27. Use `skimage.morphology.footprint_rectangle` instead.
labels, N = spim.label(peaks > 0, structure=cube(3))
/Users/amin/Code/OpenPNM/.venv/lib/python3.13/site-packages/porespy/tools/_morphology.py:174: FutureWarning: `cube` is deprecated since version 0.25 and will be removed in version 0.27. Use `skimage.morphology.footprint_rectangle` instead.
strel = cube(w)
For more details on the use of the snow2 function see this example. The snow2 function returns a Results object with results attaches as attributes. We can see these attributes if we print it:
print(snow)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Results of snow2 generated at Fri Dec 5 00:33:54 2025
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
network Dictionary with 26 items
regions Array of size (206, 200, 200)
phases Array of size (206, 200, 200)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
We can open the network in OpenPNM easily as follows:
pn = op.io.network_from_porespy(snow.network)
It’s usually helpful to inspect the network by printing it:
print(pn)
══════════════════════════════════════════════════════════════════════════════
net : <openpnm.network.Network at 0x11af6f6b0>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# Properties Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
2 throat.conns 2638 / 2638
3 pore.coords 1422 / 1422
4 pore.region_label 1422 / 1422
5 pore.phase 1422 / 1422
6 throat.phases 2638 / 2638
7 pore.region_volume 1422 / 1422
8 pore.equivalent_diameter 1422 / 1422
9 pore.local_peak 1422 / 1422
10 pore.global_peak 1422 / 1422
11 pore.geometric_centroid 1422 / 1422
12 throat.global_peak 2638 / 2638
13 pore.inscribed_diameter 1422 / 1422
14 pore.extended_diameter 1422 / 1422
15 throat.inscribed_diameter 2638 / 2638
16 throat.total_length 2638 / 2638
17 throat.direct_length 2638 / 2638
18 throat.perimeter 2638 / 2638
19 pore.volume 1422 / 1422
20 pore.surface_area 1422 / 1422
21 throat.cross_sectional_area 2638 / 2638
22 throat.equivalent_diameter 2638 / 2638
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# Labels Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
2 pore.all 1422
3 throat.all 2638
4 pore.boundary 182
5 pore.xmin 104
6 pore.xmax 78
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
When we update the SNOW algorithm to snow2 we removed many of the opinionated decision that the original version made. For instance, snow2 does not pick which diameter should be the definitive one, so this decision needs to be made by the user once they can the network into OpenPNM. This is illustrated below:
pn['pore.diameter'] = pn['pore.equivalent_diameter']
pn['throat.diameter'] = pn['throat.inscribed_diameter']
pn['throat.spacing'] = pn['throat.total_length']
The user also needs to decide which ‘shape’ to assume for the pores and throats, which impacts how the transport conductance values are computed. Here we use the pyramids_and_cuboid models:
pn.add_model(propname='throat.hydraulic_size_factors',
model=op.models.geometry.hydraulic_size_factors.pyramids_and_cuboids)
pn.add_model(propname='throat.diffusive_size_factors',
model=op.models.geometry.diffusive_size_factors.pyramids_and_cuboids)
pn.regenerate_models()
More information about the “size factors” can be found in this example.
One issue that can happen in extracted networks is that some pores can be isolated from the remainder of the network. The problem with ‘disconnected networks’ is discussed in detail here. The end result is that we need to first check the network health, then deal with any problems:
h = op.utils.check_network_health(pn)
print(h)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Key Value
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
headless_throats []
looped_throats []
isolated_pores [456, 574, 1115, 1158, 1179]
disconnected_pores [456, 574, 1115, 1158, 1179, 1204, 1346, 1347, 1348]
duplicate_throats []
bidirectional_throats []
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
The above output shows there are both ‘single’ isolated pores as well as ‘groups’ of pores that are disconnected. These all need to be “trimmed”. We can use the list provided by h['disconnected_pores'] directly in the trim function:
op.topotools.trim(network=pn, pores=h['disconnected_pores'])
h = op.utils.check_network_health(pn)
print(h)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Key Value
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
headless_throats []
looped_throats []
isolated_pores []
disconnected_pores []
duplicate_throats []
bidirectional_throats []
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Now we are ready to perform some simulations, so let’s create a phase object to compute the thermophysical properties and transport conductances:
gas = op.phase.Phase(network=pn)
gas['pore.diffusivity'] = 1.0
gas['pore.viscosity'] = 1.0
gas.add_model_collection(op.models.collections.physics.basic)
gas.regenerate_models()
[00:33:54] WARNING throat.entry_pressure was not run since the following property is missing: _models.py:480 'throat.surface_tension'
Finally we can do a Fickian diffusion simulation to find the formation factor:
fd = op.algorithms.FickianDiffusion(network=pn, phase=gas)
fd.set_value_BC(pores=pn.pores('xmin'), values=1.0)
fd.set_value_BC(pores=pn.pores('xmax'), values=0.0)
fd.run()
dC = 1.0
L = (data['shape']['x'] + 6)*data['resolution']
A = data['shape']['y']*data['shape']['z']*data['resolution']**2
Deff = fd.rate(pores=pn.pores('xmin'))*(L/A)/dC
F = 1/Deff
print(f"The Formation factor of the extracted network is {F[0]}")
print(f"The compares to a value of {data['formation factor']['Fx']} from DNS")
np.testing.assert_allclose(F, data['formation factor']['Fx'], rtol=0.1)
The Formation factor of the extracted network is 22.33394375429123
The compares to a value of 23.12 from DNS
And a Stokes flow simulation to find Permeability coefficient:
sf = op.algorithms.StokesFlow(network=pn, phase=gas)
sf.set_value_BC(pores=pn.pores('xmin'), values=1.0)
sf.set_value_BC(pores=pn.pores('xmax'), values=0.0)
sf.run()
dP = 1.0
L = (data['shape']['x'] + 6)*data['resolution']
A = data['shape']['y']*data['shape']['z']*data['resolution']**2
K = sf.rate(pores=pn.pores('xmin'))*(L/A)/dP*1e12
print(f'Permeability coefficient is {K[0]} Darcy')
print(f"The compares to a value of {data['permeability']['Kx']/1000} from DNS")
np.testing.assert_allclose(K, data['permeability']['Kx']/1000, rtol=0.1)
Permeability coefficient is 1.4732149604234128 Darcy
The compares to a value of 1.36 from DNS
Both of the above simulations agree quite well with the known values for this sample. This is not always the case because network extraction is not always perfect. One problem that can occur is that the pore sizes are much too small due to anisotropic materials. In other cases the pores are too large and overlap each other too much. Basically, the user really need to double, then triple check that their extraction is ‘sane’. It is almost mandatory to compare the extraction to some known values as we have done above. It’s also a good idea visualize the network in Paraview, as explained here, which can reveal problems.