Working with Extracted Networks#

Extracting a pore network using PoreSpy and loading into OpenPNM#

[1]:
import warnings
import scipy as sp
import numpy as np
import porespy as ps
import openpnm as op
import matplotlib.pyplot as plt
ws = op.Workspace()
ws.settings["loglevel"] = 40
warnings.filterwarnings('ignore')
%matplotlib inline
np.random.seed(10)
[2]:
# NBVAL_IGNORE_OUTPUT
im = ps.generators.overlapping_spheres(shape=[200, 200, 200], r=10, porosity=0.5, maxiter=0)
plt.imshow(im[:, :, 50]);
../../../_images/examples_reference_uncategorized_working_with_extracted_networks_2_0.png

Let’s check out the porosity of the generated image!

[3]:
eps = ps.metrics.porosity(im)
print(f"Porosity: {eps*100:.1f}%")
Porosity: 62.0%

Let’s visualize the image using porespy’s 3D visualizer: (this might take several seconds)

[4]:
# NBVAL_IGNORE_OUTPUT
im_3d = ps.visualization.show_3D(im)
plt.imshow(im_3d, cmap=plt.cm.magma);
../../../_images/examples_reference_uncategorized_working_with_extracted_networks_6_0.png
[5]:
# NBVAL_IGNORE_OUTPUT
snow = ps.networks.snow2(im, boundary_width=[[0, 3], 0, 0])

OpenPNM has an IO class specifically for importing the output from PoreSpy. The import_data method can either accept a handle to a dictionary (as output from the snow algorithm above), or it can accept a filename to a saved dctionary (saved using Python’s pickle library). All IO methods in OpenPNM return a project which is a list, in this case containing a network and a geometry object.

[6]:
# NBVAL_IGNORE_OUTPUT
proj = op.io.PoreSpy.import_data(snow.network)
print(proj)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [6], in <cell line: 2>()
      1 # NBVAL_IGNORE_OUTPUT
----> 2 proj = op.io.PoreSpy.import_data(snow.network)
      3 print(proj)

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

We can unpack the network and geometry objects from the project using the indices in the list as follows:

[7]:
# NBVAL_IGNORE_OUTPUT
net = proj[0]
geo = proj[1]
print(net)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [7], in <cell line: 2>()
      1 # NBVAL_IGNORE_OUTPUT
----> 2 net = proj[0]
      3 geo = proj[1]
      4 print(net)

NameError: name 'proj' is not defined

It is important to note that the net object only has topological information and labels. The geo object was created by the openpnm.io.PoreSpy import class to extract all geometric information from the supplied snow dict and put in on a geometry object. We can print geo to confirm:

[8]:
# NBVAL_IGNORE_OUTPUT
print(geo)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [8], in <cell line: 2>()
      1 # NBVAL_IGNORE_OUTPUT
----> 2 print(geo)

NameError: name 'geo' is not defined

Now let’s plot things to see what we have:

[9]:
# NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots(1, 1, figsize=[8, 8])
ax = op.topotools.plot_connections(network=net, alpha=0.8, color='grey', ax=ax)
ax = op.topotools.plot_coordinates(network=net, ax=ax, color='b', markersize=50)
ax = op.topotools.plot_coordinates(network=net, pores=net.pores('xmax'), ax=ax, color='r', markersize=50)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [9], in <cell line: 3>()
      1 # NBVAL_IGNORE_OUTPUT
      2 fig, ax = plt.subplots(1, 1, figsize=[8, 8])
----> 3 ax = op.topotools.plot_connections(network=net, alpha=0.8, color='grey', ax=ax)
      4 ax = op.topotools.plot_coordinates(network=net, ax=ax, color='b', markersize=50)
      5 ax = op.topotools.plot_coordinates(network=net, pores=net.pores('xmax'), ax=ax, color='r', markersize=50)

NameError: name 'net' is not defined
../../../_images/examples_reference_uncategorized_working_with_extracted_networks_15_1.png

This looks pretty good, but it only has boundary pores on the right face, indicated by the red dots. When we ran the snow algorithm we specifically told it to only put boundary pores the "right". We could have added them to all faces during the extraction, but for the sake of demonstration we can add them after the fact, although the result is slightly different, as you’ll see.

We’ll use the find_surface_pores function in the topotools module. This function applies a Delaunay tessellation between the interior pores and some fictitious “marker” nodes. Only pores that are on the surface will be connected to these marker nodes. To get the best result from the find_surface_pores function, it’s a good idea to supply your own markers, so let’s make a 2D plane of points, positioned outside the left face of the domain:

[10]:
m = np.meshgrid(range(50, 195, 10), range(50, 195, 10))
m = np.vstack([-10*np.ones_like(m[0].flatten()), m[0].flatten(), m[1].flatten()]).T

Now we pass these points in as markers to the find_surface_pores function:

[11]:
op.topotools.find_surface_pores(network=net, markers=m, label='left')
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [11], in <cell line: 1>()
----> 1 op.topotools.find_surface_pores(network=net, markers=m, label='left')

NameError: name 'net' is not defined

Lastly we want to “clone” these pores and translate them to domain edge:

[12]:
op.topotools.clone_pores(network=net, pores=net.pores('left'), labels='left_boundary')
net['pore.coords'][net.pores('left_boundary')] *= [0, 1, 1]
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [12], in <cell line: 1>()
----> 1 op.topotools.clone_pores(network=net, pores=net.pores('left'), labels='left_boundary')
      2 net['pore.coords'][net.pores('left_boundary')] *= [0, 1, 1]

NameError: name 'net' is not defined

Now let’s inspect the result using the quick plotting tools in the topotools module. First we’ll add a new label called 'right_boundary' to match the 'left_boundary' we added above, then we’ll plot the throats that connect to ther 'right_boundary' or 'left_boundary':

[13]:
# NBVAL_IGNORE_OUTPUT
Ps = net.pores('xmax')
net.set_label('right_boundary', pores=Ps)
Ts = net.find_neighbor_throats(pores=net.pores('right_boundary'), mode='or')
net.set_label('right_boundary', throats=Ts)

fig, ax = plt.subplots(1, 1, figsize=[5, 5])
ax = op.topotools.plot_coordinates(network=net, color='g', alpha=0.2, ax=ax)
ax = op.topotools.plot_connections(network=net, throats=net.throats('right_boundary'), color='r', ax=ax)
ax = op.topotools.plot_connections(network=net, throats=net.throats('left_boundary'), color='b', ax=ax)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [13], in <cell line: 2>()
      1 # NBVAL_IGNORE_OUTPUT
----> 2 Ps = net.pores('xmax')
      3 net.set_label('right_boundary', pores=Ps)
      4 Ts = net.find_neighbor_throats(pores=net.pores('right_boundary'), mode='or')

NameError: name 'net' is not defined

This result shows that the boundary pores added during the snow extraction (red) are randomly oriented while those added by the clone_pores_method are aligned with their internal counter-parts. It also seems like longer connections are made with the clone_pores_method which may be the result of the Delanauy tessellation identifying pores that are too deep into the domain.