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
import os
from pathlib import Path
Load a small image of Berea sandstone:
path = Path(os.getcwd(),
'../../tests/fixtures/berea_100_to_300.npz')
data = np.load(path.resolve())
im = data['im']
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[2], line 3
1 path = Path(os.getcwd(),
2 '../../tests/fixtures/berea_100_to_300.npz')
----> 3 data = np.load(path.resolve())
4 im = data['im']
File /opt/hostedtoolcache/Python/3.10.13/x64/lib/python3.10/site-packages/numpy/lib/npyio.py:427, in load(file, mmap_mode, allow_pickle, fix_imports, encoding, max_header_size)
425 own_fid = False
426 else:
--> 427 fid = stack.enter_context(open(os_fspath(file), "rb"))
428 own_fid = True
430 # Code to distinguish from NumPy binary files and pickles.
FileNotFoundError: [Errno 2] No such file or directory: '/home/runner/work/OpenPNM/OpenPNM/docs/tests/fixtures/berea_100_to_300.npz'
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',
parallelization=None,
)
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 Wed Nov 22 17:21:54 2023
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
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 0x183be2404f0>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# Properties Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
2 throat.conns 2510 / 2510
3 pore.coords 1420 / 1420
4 pore.region_label 1420 / 1420
5 pore.phase 1420 / 1420
6 throat.phases 2510 / 2510
7 pore.region_volume 1420 / 1420
8 pore.equivalent_diameter 1420 / 1420
9 pore.local_peak 1420 / 1420
10 pore.global_peak 1420 / 1420
11 pore.geometric_centroid 1420 / 1420
12 throat.global_peak 2510 / 2510
13 pore.inscribed_diameter 1420 / 1420
14 pore.extended_diameter 1420 / 1420
15 throat.inscribed_diameter 2510 / 2510
16 throat.total_length 2510 / 2510
17 throat.direct_length 2510 / 2510
18 throat.perimeter 2510 / 2510
19 pore.volume 1420 / 1420
20 pore.surface_area 1420 / 1420
21 throat.cross_sectional_area 2510 / 2510
22 throat.equivalent_diameter 2510 / 2510
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# Labels Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
2 pore.all 1420
3 throat.all 2510
4 pore.boundary 180
5 pore.xmin 103
6 pore.xmax 77
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
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, 783, 1023, 1115, 1158, 1179]
disconnected_pores [227, 246, 247, 319, 362, 456, 574, 783, 1023, 1115, 1158, 1179, 1181, 1189, 1204, 1345, 1346, 1347]
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()
[17:21:55] 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.09)
The Formation factor of the extracted network is 25.098142684548733
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.05)
Permeability coefficient is 1.3181940680110902 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.