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.14/x64/lib/python3.10/site-packages/numpy/lib/_npyio_impl.py:455, in load(file, mmap_mode, allow_pickle, fix_imports, encoding, max_header_size)
    453     own_fid = False
    454 else:
--> 455     fid = stack.enter_context(open(os.fspath(file), "rb"))
    456     own_fid = True
    458 # 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.