import openpnm as op
import numpy as np
pn = op.network.Cubic([3, 3])
pn.add_model_collection(op.models.collections.geometry.spheres_and_cylinders)
pn.regenerate_models()Benefits of the RC Format
At the beginning of this section a few of the ‘high level’ benefits of the RC format were mentioned, specifically that it is human readable, self-explanatory, and compatible with common file formats. To end this section we will explore the numerical and computational benefits of using the “RC” format.
Vectorized Calculations are Easy
Having data in columns makes it very easy to perform “vectorized” calculations. For example, we can find the volume of a throat as follows:
Vt = np.pi/4*(pn['throat.diameter'])**2 * pn['throat.length']
pn['throat.diameter'] = VtThe 'throat.diameter' and 'throat.length' arrays are both the same length, and this calculation produces a 'throat.volume' array that is the same length again. In some programming languages, matrix multiplication is assumed by default, but Numpy arrays in Python assume “elementwise” calcuations. We can also see that the scalar values are applied to each element in calculation.
Vectorized Calculations are Fast
Vectorized calculations are fast since Numpy is a highly optimized package. Multiplying two very large Numpy arrays together takes just milliseconds. There are some packages that attempt to speed this up by parallelizing or using the GPU, but these are not really useful for pore-network simulations since the performance bottleneck, by a long way, is solving the systems of transport equations.
Fancy Indexing and Boolean Masking
Fancy indexing and masking avoid the need for almost all ‘for loops’. For instance, we can find all pores with a size smaller than some threshold, and set the volume of those pores to 0 using:
Ps = pn['pore.diameter'] < 0.5
pn['pore.volume'][Ps] = 0.0Or we can use so-called “fancy indexing” which works like a mask but using numerical array indices:
pn['pore.diameter'][[0, 5, 3, 4]] = 0.0Fancy indexing is particularly useful for determining throat properties based on the values in its 2 neighboring pores, but this will be covered elsewhere.
Data Storage is Easy
Data storage is simple since it is trivial to read and write columnar data to various file formats. The only caveat is that ‘multidimensional’ data like 'pore.coords' needs to be divided up into separate columns if writing to a CSV or other ‘Excel-type’ file. There are plenty of other formats that don’t require this though, like zarr and hdf5.
Adding Pores and Throats is Simple
Adding pores/throats to the network is as simple as adding rows to the end of the table. There is no benefit to adding rows to the middle of the table since this only changes the implied index number, which is immaterial. The only challenge when adding new elements (i.e. pore or throats) is that all the arrays must be increased in length, so we need to know what values to put into the new locations. For boolean arrays it is safe to assume False values, but for numerical arrays it is more challenging. The common options like 0 or -1 may have a physical meaning, so the only safe approach is to put nan in all locations, but this requires converting any int arrays to float.
Removing Pores and Throats is Simple
Removing pores/throats is almost as simple as removing rows from the table, but there is one major exception: when deleting pores, the values in the 'throat.conns' array need to be updated since they are essentially hard-coded to point to specific rows. For instance, if pore 0 is deleted, then pore 1 becomes pore 0, pore 2 becomes for 1, and so forth. So any entries in 'throat.conns' must be re-index accordingly. Despite this complication, removing elements is very straight-forward.
Avoiding for-loops by using conns to Scan Throats
Python is often accused of being slow, which it is if you try to do a for-loop on a list. However, there are many ways to avoid this. By combining vectorization with fancy indexing it is possible to “loop” through all throats in the network. For instance to set the throat diameter to be equal to the minimum of it’s two neighboring pores we can do:
conns = pn['throat.conns']
neighboring_pore_diameters = pn['pore.diameter'][conns]
neighboring_pore_diametersarray([[0. , 0.60829207],
[0.60829207, 0.28534076],
[0. , 0. ],
[0. , 0. ],
[0.34605188, 0.46065186],
[0.46065186, 0.35152799],
[0. , 0. ],
[0.60829207, 0. ],
[0.28534076, 0. ],
[0. , 0.34605188],
[0. , 0.46065186],
[0. , 0.35152799]])throat_diameter = np.amin(neighboring_pore_diameters, axis=1)
throat_diameter
pn['throat.diameter'] = throat_diameter*0.5 # Use 1/2 the found diameterAlthough it is (initially) more intuitive to scan through the network on a ‘pore-by-pore’ basis, it is actually much faster and more concise to scan through each throat as shown above. This small change in workflow results in large speed-ups so it’s worth it.
Avoiding for-loops Using “unbuffered” Numpy Operations
Although scanning ‘throat-by-throat’ works far more often than you might expect, sometimes it is not possible to accomplish the desired result by scanning the throats. In these cases, numpy has a pre-built but little known solution called ‘unbuffered’ operations.
Note to self: Using the incidence matrix with unbuffered operations is super easy! np.add.at(values, im.row, data[im.col])
First let’s illustrate the general idea by attempting to subtract 1 from a given location in an array 3 times.
Using the standard vectorization approach would look like this:
a = np.array([1, 2, 3, 4, 5])
b = np.array([1, 1, 1])
ind = [2, 2, 2]
a[ind] = a[ind] - b
print(a)[1 2 2 4 5]We can see that 1 was only subtracted from a once even though our index suggested that we wanted this to happen three times. This occurred because numpy buffers the operation for computational reasons, so only the last operation has an effect.
The correct way to do this is as follows:
a = np.array([1, 2, 3, 4, 5])
b = np.array([1, 1, 1])
ind = [2, 2, 2]
np.subtract.at(a, ind, b)
print(a)[1 2 0 4 5]Now we can see that the value of 1 has been subtracted the desired number of times.
The practical example of this is the determination of pore surface area. We start by assuming that all pores have the surface area of a sphere, then subtract the cross-sectional area of each throat attached to that pore. Because pore have many throats, we want to subtract several values from the pore surface area, hence need to use ‘unbuffered’ operations.
Here is an example of using the wrong approach with standard vectorization:
pore_surface_area = 4*np.pi*(pn['pore.diameter']/2)**2
throat_cross_section = np.pi*(pn['throat.diameter'])**2
conns = pn['throat.conns']
print("Original surface area:", pore_surface_area)
pore_surface_area[conns] = (pore_surface_area[conns].T - throat_cross_section).T
print("Incorrectly adjusted surface area:", pore_surface_area)Original surface area: [0. 1.16244974 0.25578643 0. 0. 0.
0.3762117 0.66664638 0.38821265]
Incorrectly adjusted surface area: [0. 1.16244974 0.19183982 0. 0. 0.
0.3762117 0.66664638 0.38821265]pore_surface_area = 4*np.pi*(pn['pore.diameter']/2)**2
throat_cross_section = np.pi*(pn['throat.diameter'])**2
conns = pn['throat.conns']
print("Original surface area:", pore_surface_area)
np.subtract.at(pore_surface_area, conns[:, 0], throat_cross_section)
np.subtract.at(pore_surface_area, conns[:, 1], throat_cross_section)
print("Correctly adjusted surface area:", pore_surface_area)Original surface area: [0. 1.16244974 0.25578643 0. 0. 0.
0.3762117 0.66664638 0.38821265]
Correctly adjusted surface area: [0. 1.09850313 0.19183982 0. 0. 0.
0.28215877 0.47554029 0.29115949]We can see that in every case the correctly computed pore surface area is smaller than the incorrect one since the incorrectly calculated value only had 1 throat cross-section subtracted from it, while the correct approach resulted in the cross-sectional area of multiple throats being subtracted as expected.