Tutorial 1 of 3: Getting Started with OpenPNM

This tutorial is intended to show the basic outline of how OpenPNM works, and necessarily skips many of the more useful and powerful features of the package. So if you find yourself asking “why is this step so labor intensive” it’s probably because this tutorial deliberately simplifies some features to provide a more gentle introduction. The second and third tutorials of this User-Guide dive into the package more deeply, but those features are best appreciated once the basics are understood.

Learning Objectives

  1. Introduce the main OpenPNM objects and their roles
  2. Explore the way OpenPNM stores data, including network topology
  3. Learn some handy tools for working with objects
  4. Generate a standard cubic Network topology
  5. Calculate geometrical properties and assign them to a Geometry object
  6. Calculate thermophysical properties and assign to a Phase object
  7. Define pore-scale physics and assign transport parameters to a Physics object
  8. Run a permeability simulation using the pre-defined Algorithm
  9. Use the package to calculate the permeability coefficient of a porous media


Python and Numpy Tutorials

  • OpenPNM is written in Python. One of the best guides to learning Python is the excellent interactive online tutorial called Learn Python the Hard Way, which was originally a website and is now a book and video series. The Python Module of the Week website is another excellent resource with nice examples illustrating various features of Python standard library, and has also been published in a book.
  • For information on using Numpy, Scipy and generally doing scientific computing in Python checkout the Scipy lecture notes. The Scipy website also offers as solid introduction to using Numpy arrays.
  • The Stackoverflow website is an incredible resource for all computing related questions, including simple usage of Python, Scipy and Numpy functions.
  • For users more familiar with Matlab, there is a Matlab-Numpy cheat sheet that explains how to translate familiar Matlab commands to Numpy.

Overview of Data Storage in OpenPNM

Before creating an OpenPNM simulation it is necessary to give a quick description of how data is stored in OpenPNM; after all, a significant part of OpenPNM is dedicated to data storage and handling.

Python Dictionaries or dicts

OpenPNM employs 5 main objects which each store and manage a different type of information or data:

  1. Network: Manages topological data such as pore spatial locations and pore-to-pore connections
  2. Geometry: Manages geometrical properties such as pore diameter and throat length
  3. Phase: Manages thermophysical properties such as temperature and viscosity
  4. Physics: Manages pore-scale transport parameters such as hydraulic conductance
  5. Algorithm: Contains algorithms that use the data from other objects to perform simulations, such as diffusion or drainage

We will encounter each of these objects in action before the end of this tutorial.

Each of the above objects is a subclass of the Python dictionary or dict, which is a very general storage container that allows values to be accessed by a name using syntax like:

>>> foo = dict()  # Create an empty dict
>>> foo['bar'] = 1  # Store an integer under the key 'bar'
>>> foo['bar']  # Retrieve the integer stored in 'bar'

A detailed tutorial on dictionaries can be found here. The dict does not offer much functionality aside from basic storage of arbitrary objects, and it is meant to be extended. OpenPNM extends the dict to have functionality specifically suited for dealing with OpenPNM data. More information about the functionality of OpenPNM’s subclassed dicts can be found in the Overall Design.

Numpy Arrays of Pore and Throat Data

All data are stored in arrays which can accessed using standard array syntax. More details on the data storage scheme are given in Data Storage, but the following gives a quick overview:

  1. All pore and throat properties are stored in Numpy arrays. All data will be automatically converted to a Numpy array if necessary.
  2. The data for pore i (or throat i) can be found in element of i of an array. This means that pores and throat have indices which are implied by their position in arrays. When we speak of retrieving pore locations, it refers to the indices in the Numpy arrays.
  3. Arrays that store pore data are Np-long, while arrays that store throat data are Nt-long, where Np is the number of pores and Nt is the number of throats in the network.
  4. Each property is stored in it’s own array, meaning that ‘pore diameter’ and ‘throat volume’ are each stored in a separate array.
  5. Arrays can be any size in the other dimensions. For instance, triplets of pore coordinates (i.e. [x, y, z]) can be stored for each pore creating an Np-by-3 array.
  6. The storage of topological connections is also very nicely accomplished with this ‘list-based’ format, by creating an array ('throat.conns') that stores which pore indices are found on either end of a throat. This leads to an Nt-by-2 array. The implications and advantages of this storage scheme are discussed further in Representing Topology.

OpenPNM Objects: Combining dicts and Numpy Arrays

OpenPNM objects combine the above two levels of data storage, meaning they are dicts that are filled with Numpy arrays. OpenPNM enforces several rules to help maintain data consistency:

  1. When storing arrays in an OpenPNM object, their name (or dictionary key) must be prefixed with 'pore.' or 'throat.'.
  2. OpenPNM uses the prefix of the dictionary key to infer how long the array must be.
  3. The specific property that is stored in each array is indicated by the suffix such as 'pore.diameter' or 'throat.length'.
  4. Writing scalar values to OpenPNM objects automatically results in conversion to a full length array filled with the scalar value.
  5. Arrays containing Boolean data are treated as labels, which are explained later in this tutorial.

The following code snippets give examples of how all these pieces fit together using an Empty network as an example:

>>> import OpenPNM
>>> import scipy as sp
>>> net = OpenPNM.Network.Empty(Np=10, Nt=10)  # Instantiate an empty network object with 10 pores and 10 throats
>>> net['pore.foo'] = sp.ones([net.Np, ])  # Assign an Np-long array of ones
>>> net['pore.bar'] = range(0, net.Np)  # Assign an Np-long array of increasing ints
>>> type(net['pore.bar'])  # The Python range iterator is converted to a proper Numpy array
<class 'numpy.ndarray'>
>>> net['pore.foo'][4] = 44.0  # Overwrite values in the array
>>> net['pore.foo'][4]  # Retrieve values from the array
>>> net['pore.foo'][2:6]  # Extract a slice of the array
array([  1.,   1.,  44.,   1.])
>>> net['pore.foo'][[2, 4, 6]]  # Extract specific locations
array([  1.,  44.,   1.])
>>> net['throat.foo'] = 2  # Assign a scalar
>>> len(net['throat.foo'])  # The scalar values is converted to an Nt-long array
>>> net['throat.foo'][4]  # The scalar value was placed into all locations

Generate a Cubic Network

Now that we have seen the rough outline of how OpenPNM objects store data, we can begin building a simulation. Start by importing OpenPNM and the Scipy package:

>>> import OpenPNM
>>> import scipy as sp

Next, generate a Network by choosing the Cubic class, then create an instance with the desired parameters:

>>> pn = OpenPNM.Network.Cubic(shape=[4, 3, 1], spacing=0.0001)

The Network object stored in pn contains pores at the correct spatial positions and connections between the pores according the cubic topology.

  • The shape argument specifies the number of pores in the [X, Y, Z] directions of the cube. Networks in OpenPNM are always 3D dimensional, meaning that a 2D or “flat” network is still 1 layer of pores “thick” so [X, Y, Z] = [20, 10, 1], thus pn in this tutorial is 2D which is easier for visualization.
  • The spacing argument controls the center-to-center distance between pores and it can be a scalar or vector (i.e. [0.0001, 0.0002, 0.0003]).

The resulting network looks like:


This image was creating using Paraview, using the instructions given here: Example in the OpenPNM-Example collection

Inspecting Object Properties

OpenPNM objects have additional methods for querying their relevant properties, like the number of pores or throats, which properties have been defined, and so on:

>>> pn.num_pores()
>>> pn.Np  # Shortcut to get number of pores
>>> pn.num_throats()
>>> pn.Nt
>>> pn.props()
['pore.coords', 'pore.index', 'throat.conns']

More information about these various functions is given in Overall Design. It is also convenient to type print(pn) at the command line to view a nicely formatted table showing the current state of pn.

Accessing Pores and Throats

One simple but important feature of OpenPNM is the ability to label pores and throats. When a Cubic network is created, several labels are automatically created: the pores on each face are labeled ‘left’, ‘right’, etc. These labels can be used as follows:

>>> pn.pores('left')
array([0, 3, 6, 9])

The ability to retrieve pore indices is handy for querying pore properties, such as retrieving the pore coordinates of all pores on the ‘left’ face:

>>> pn['pore.coords'][pn.pores('left')]
array([[  5.00000000e-05,   5.00000000e-05,   5.00000000e-05],
       [  1.50000000e-04,   5.00000000e-05,   5.00000000e-05],
       [  2.50000000e-04,   5.00000000e-05,   5.00000000e-05],
       [  3.50000000e-04,   5.00000000e-05,   5.00000000e-05]])

A list of all labels currently assigned to the network can be obtained with:

>>> pn.labels()
['pore.all', 'pore.back', 'pore.bottom', 'pore.front', 'pore.internal', 'pore.left', 'pore.right', 'pore.top', 'throat.all']

The existing labels are also listed when an object is printed using print(pn). Detailed use of labels is given in Data Storage.

Create a Geometry Object and Assign Geometric Properties to Pores and Throats

The Network pn does not contain any information about pore and throat sizes at this point. The next step is to create a Geometry object to manage the geometrical properties.

>>> geom = OpenPNM.Geometry.GenericGeometry(network=pn, pores=pn.Ps, throats=pn.Ts)

This statement contains three arguments:

  • network tells the Geometry object which Network it is associated with. There can be multiple networks defined in a given session, so all objects must be associated with a single network.
  • pores and throats indicate the locations in the Network where this Geometry object will apply. In this tutorial geom applies to all pores and throats, but there are many cases where different regions of the network have different geometrical properties, so OpenPNM allows multiple Geometry objects to be created for managing the data in each region, but this is a subject for Tutorial 2 of 3: Digging Deeper into OpenPNM.

Add Pore and Throat Size Information

This freshly instantiated Geometry object (geom) contains no geometric properties as yet. For this tutorial we’ll use the direct assignment of manually calculated values.

We’ll start by assigning diameters to each pore from a random distribution, spanning 0 um to 100 um. The upper limit matches the spacing of the Network which was set to 0.0001 m (i.e. 100 um), so pore diameters exceeding 100 um might overlap with their neighbors. Using the Scipy rand function creates an array of random numbers between 0 and 0.0001 that is Np-long, meaning each pore is assigned a unique random number

>>> geom['pore.diameter'] = sp.rand(pn.Np)*0.0001  # Units of meters

We usually want the throat diameters to always be smaller than the two pores which it connects to maintain physical consistency. This requires understanding a little bit about how OpenPNM stores network topology. Consider the following:

>>> P12 = pn['throat.conns']  # An Nt x 2 list of pores on the end of each throat
>>> D12 = geom['pore.diameter'][P12]  # An Nt x 2 list of pore diameters
>>> Dt = sp.amin(D12, axis=1)  # An Nt x 1 list of the smaller pore from each pair
>>> geom['throat.diameter'] = Dt

Let’s dissect the above lines.

  • Firstly, P12 is a direct copy of the Network’s 'throat.conns' array, which contains the indices of the pore-pair connected by each throat.
  • Next, this Nt-by-2 array is used to index into the 'pore.diameter' array, resulting in another Nt-by-2 array containing the diameters of the pores on each end of a throat.
  • Finally, the Scipy function amin is used to find the minimum diameter of each pore-pair by specifying the axis argument as 1, and the resulting Nt-by-1 array is assigned to geom['throat.diameter'].
  • This trick of using 'throat.conns' to index into a pore property array is commonly used in OpenPNM and you should have a second look at the above code to understand it fully. Refer to Representing Topology for a full discussion.

We must still specify the remaining geometrical properties of the pores and throats. Since we’re creating a “Stick-and-Ball” geometry, the sizes are calculated from the geometrical equations for spheres and cylinders. For pore volumes, assume a sphere:

>>> Rp = geom['pore.diameter']/2
>>> geom['pore.volume'] = (4/3)*3.14159*(Rp)**3

The length of each throat is the center-to-center distance between pores, minus the radius of each of two neighboring pores.

>>> C2C = 0.0001  # The center-to-center distance between pores
>>> Rp12 = Rp[pn['throat.conns']]
>>> geom['throat.length'] = C2C - sp.sum(Rp12, axis=1)

The volume of each throat is found assuming a cylinder:

>>> Rt = geom['throat.diameter']/2
>>> Lt = geom['throat.length']
>>> geom['throat.volume'] = 3.14159*(Rt)**2*Lt

The basic geometrical properties of the network are now defined. The Geometry class possesses a method called plot_histograms that produces a plot of the most pertinent geometrical properties. The following figure doesn’t look very good since the network in this example has only 12 pores, but the utility of the plot for quick inspection is apparent.


Create a Phase Object

The simulation is now topologically and geometrically defined. It has pore coordinates, pore and throat sizes and so on. In order to perform any simulations it is necessary to define a Phase object to manage all the thermophysical properties of the fluids in the simulation:

>>> water = OpenPNM.Phases.GenericPhase(network=pn)
  • pn is passed as an argument because Phases must know to which Network they belong.
  • Note that pores and throats are NOT specified; this is because Phases are mobile and can exist anywhere or everywhere in the domain, so providing specific locations does not make sense. Algorithms for dynamically determining actual phase distributions are discussed later.

Add Thermophysical Properties

Now it is necessary to fill this Phase object with the desired thermophysical properties. OpenPNM includes a framework for calculating thermophysical properties from models and correlations, but this is covered in Tutorial 2 of 3: Digging Deeper into OpenPNM. For this tutorial, we’ll use the basic approach of simply assigning static values as follows:

>>> water['pore.temperature'] = 298.0
>>> water['pore.viscosity'] = 0.001
  • The above lines utilize the fact that OpenPNM converts scalars to full length arrays, essentially setting the temperature in each pore to 298.0 K.

Create a Physics Object

We are still not ready to perform any simulations. The last step is to define the desired pore-scale physics models, which dictate how the phase and geometrical properties interact to give the transport parameters. A classic example of this is the Hagen-Poiseuille equation for fluid flow through a throat to predict the flow rate as a function of the pressure drop. The flow rate is proportional to the geometrical size of the throat (radius and length) as well as properties of the fluid (viscosity) and thus combines geometrical and thermophysical properties:

>>> phys_water = OpenPNM.Physics.GenericPhysics(network=pn, phase=water, geometry=geom)
  • As with all objects, the Network must be specified
  • Physics objects combine information from a Phase (i.e. viscosity) and a Geometry (i.e. throat diameter), so each of these must be specified.
  • Physics objects do not require the specification of which pores and throats where they apply, since this information is implied by the geometry argument which was already assigned to specific locations.

Specify Desired Pore-Scale Transport Parameters

We need to calculate the numerical values representing our chosen pore-scale physics. To continue with the Hagen-Poiseuille example lets calculate the hydraulic conductance of each throat in the network. The throat radius and length are easily accessed as:

>>> R = geom['throat.diameter']/2
>>> L = geom['throat.length']

The viscosity of the Phases was only defined in the pores; however, the hydraulic conductance must be calculated for each throat. There are several options, but to keep this tutorial simple we’ll create a scalar value:

>>> mu_w = 0.001
>>> phys_water['throat.hydraulic_conductance'] = 3.14159*R**4/(8*mu_w*L)

Numpy arrays support vectorization, so since both L and R are arrays of Nt-length, their multiplication in this way results in another array that is also Nt-long.

Create an Algorithm Object for Performing a Permeability Simulation

Finally, it is now possible to run some useful simulations. The code below estimates the permeability through the network by applying a pressure gradient across and calculating the flux. This starts by creating a StokesFlow algorithm, which is pre-defined in OpenPNM:

>>> alg = OpenPNM.Algorithms.StokesFlow(network=pn, phase=water)
  • Like all the above objects, Algorithms must be assigned to a Network via the network argument.
  • This algorithm is also associated with a Phase object, in this case water, which dictates which pore-scale Physics properties to use (recall that phys_water was associated with water).

Next the boundary conditions are applied using the set_boundary_conditions method on the Algorithm object. Let’s apply a 1 atm pressure gradient between the left and right sides of the domain:

>>> BC1_pores = pn.pores('front')
>>> alg.set_boundary_conditions(bctype='Dirichlet', bcvalue=202650,
...                             pores=BC1_pores)
>>> BC2_pores = pn.pores('back')
>>> alg.set_boundary_conditions(bctype='Dirichlet', bcvalue=101325,
...                             pores=BC2_pores)

To actually run the algorithm use the run method:

>>> alg.run()

This builds the coefficient matrix from the existing values of hydraulic conductance, and inverts the matrix to solve for pressure in each pore, and stores the results within the Algorithm’s dictionary under 'pore.pressure'.

To determine the permeability coefficient, we must invoke Darcy’s law: Q = KA/uL(Pin - Pout). Everything in this equation is known except for the volumetric flow rate Q. The StokesFlow algorithm possesses a rate method that calculates the rate of a quantity leaving a specified set of pores:

>>> Q = alg.rate(pores=pn.pores('top'))
>>> A = 0.0001*3*1  # Cross-sectional area for flow
>>> L = 0.0001*4  # Length of flow path
>>> del_P = 101325  # Specified pressure gradient
>>> K = Q*mu_w*L/(A*del_P)

The StokesFlow class was developed with permeability simulations in mind, so a specific method is available for determining the permeability coefficient that essentially applies to recipe from above. This method could struggle with non-uniform geometries though, so use with caution:

>>> K = alg.calc_eff_permeability()

The results ('pore.pressure') are held within the alg object and must be explicitly returned to the air object by the user if they wish to use these values in a subsequent calculation. The point of this data containment is to prevent unintentional overwriting of data. Each algorithm has a method called return_results which places the pertinent values back onto the appropriate Phase object.

>>> alg.return_results()

Using Paraview for Visualization, the resulting pressure gradient across the network can be seen: