# Tutorial 2 of 3: Digging Deeper into OpenPNM¶

This tutorial will follow the same outline as the Tutorial 1 of 3: Getting Started with OpenPNM, but will dig a little bit deeper at each step to reveal the important features of OpenPNM that were glossed over previously.

Topics Covered in this Tutorial

- Tutorial 2 of 3: Digging Deeper into OpenPNM
- Building a Cubic Network
- Initialize and Build
*Multiple*Geometry Objects - Create a Phase Object and Assign Thermophysical Property Models
- Create Physics Objects for Each Geometry
- Pore-Scale Models: The Big Picture
- Determine Permeability Tensor by Changing Inlet and Outlet Boundary Conditions
- Using the Workspace Manager to Save (and reload) the simulation

**Learning Objectives**

- Explore different network topologies, and learn some handy topological query methods
- Create a
*heterogeneous*domain with different geometrical properties in different regions - Learn about data exchange between objects
- Utilize pore-scale models for calculating properties of all types
- Propagate changing geometrical and thermo-physical properties to all dependent properties
- Calculate the permeability tensor for the stratified media
- Use the Workspace Manager to save and load a simulation

## Building a Cubic Network¶

As usual, start by importing the OpenPNM and Scipy packages:

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

Let’s generate a cubic network again, but with a different connectivity:

```
>>> pn = OpenPNM.Network.Cubic(shape=[20, 20, 10], spacing=0.0001, connectivity=8)
```

- This
**Network**has pores distributed in a cubic lattice, but connected to diagonal neighbors due to the`connectivity`

being set to 8 (the default is 6 which is orthogonal neighbors). The various options are outlined in the**Cubic**class’s documentation which can be viewed with the Object Inspector in Spyder. - OpenPNM includes several other classes for generating networks including random topology based on Delaunay tessellations (
**Delaunay**). - It is also possible to import networks from external sources, such as networks extracted from tomographic images, or that networks generated by external code.

## Initialize and Build *Multiple* Geometry Objects¶

One of the main functionalities of OpenPNM is the ability to assign drastically different geometrical properties to different regions of the domain to create heterogeneous materials, such as layered structures. To demonstrate the motivation behind this feature, this tutorial will make a material that has different geometrical properties on the top and bottom surfaces compared to the internal pores. We need to create one **Geometry** object to manage the top and bottom pores, and a second to manage the remaining internal pores:

```
>>> Ps1 = pn.pores(['top', 'bottom'])
>>> Ts1 = pn.find_neighbor_throats(pores=Ps1, mode='union')
>>> geom1 = OpenPNM.Geometry.GenericGeometry(network=pn, pores=Ps1, throats=Ts1, name='surface')
>>> Ps2 = pn.pores(['top', 'bottom'], mode='not')
>>> Ts2 = pn.find_neighbor_throats(pores=Ps2, mode='intersection')
>>> geom2 = OpenPNM.Geometry.GenericGeometry(network=pn, pores=Ps2, throats=Ts2, name='core')
```

- The above statements result in two distinct
**Geometry**objects, each applying to different regions of the domain.`geom1`

applies to only the pores on the top and bottom surfaces (automatically labeled ‘top’ and ‘bottom’ during the network generation step), while`geom2`

applies to the pores ‘not’ on the top and bottom surfaces. - The assignment of throats is more complicated and illustrates the
`find_neighbor_throats`

method, which is one of the more useful topological query methods on the**Network**class. In both of these calls, all throats connected to the given set of pores (`Ps1`

or`Ps2`

) are found; however, the`mode`

argument alters which throats are returned. The terms`'union'`

and`'intersection'`

are used in the “set theory” sense, such that`'union'`

returns*all*throats connected to the pores in the supplied list, while`'intersection'`

returns the throats that are*only*connected to the supplied pores. More specifically, if pores 1 and 2 have throats [1, 2] and [2, 3] as neighbors, respectively, then the`'union'`

mode returns [1, 2, 3] and the`'intersection'`

mode returns [2]. A detailed description of this behavior is given in Representing Topology.

### Naming Objects¶

Each of the **Geometry** objects was assigned a `name`

during instantiation, and this is stored in the `name`

attribute:

```
>>> geom1.name # Inspect object's name
'surface'
>>> geom1.name = 'foobar' # Change object's name
>>> geom1.name # Ensure new name was set
'foobar'
>>> geom1.name = 'surface' # Replace original name
```

Naming objects in this way serves several purposes:

- It helps users keep track of which variable points to which object (i.e.
`geom1`

vs.`geom2`

). This is useful when interacting with the objects at the command line using`geom1.name`

, which will report`'surface'`

. - When any core object is instantiated, a
*label*is created in the**Network**based on the object’s name, indicating which pores and throats belong to which object. It this case, the pores assigned to`geom1`

can be quickly retrieved using`pn.pores('surface')`

or`pn.pores(geom1.name)`

. The use of*labels*is detailed in Data Storage. - Because the
*labels*are so integral to tracking which locations belong to which objects, all**Core**objects are automatically assigned a randomly generated name if none is specified during instantiation. - When an object is renamed, OpenPNM takes care of changing the names of the
*labels*throughout the simulation. Of course, no two objects can have the same name. In fact, an object cannot be given a name if it is already in use for another*label*.

### Assign Static Seed values to *Each* Geometry¶

In Tutorial 1 of 3: Getting Started with OpenPNM we only assigned ‘static’ values to the **Geometry** object, which we calculated explicitly. In this tutorial we will use the *pore-scale models* that are provided with OpenPNM. To get started, however, we’ll assign static random seed values between 0 and 1 to each pore on both **Geometry** objects, by assigning random numbers to each Geometry’s `'pore.seed'`

property:

```
>>> geom1['pore.seed'] = sp.rand(geom1.Np)
>>> geom2['pore.seed'] = sp.rand(geom2.Np)
```

- Each of the above lines produced an array of different length, corresponding to the number of pores assigned to each
**Geometry**object. This is accomplished by the calls to`geom1.Np`

and`geom2.Np`

, which return the number of pores on each object. - Every Core object in OpenPNM possesses the same set of methods for managing their data, such as counting the number of pore and throat values they represent; thus,
`pn.Np`

returns 1000 while`geom1.Np`

and`geom2.Np`

return 200 and 800 respectively.

### Accessing Data Distributed Between Geometries¶

The segmentation of the data between separate Geometry objects is essential to the management of pore-scale models, although it does create a complication: it’s not easy to obtain a single array containing *all* the values of a given property for the whole network. It is technically possible to piece this data together manually since we know the locations where each **Geometry** object applies, but this is tedious so OpenPNM provides a shortcut. First, let’s illustrate the manual approach using the `'pore.seed'`

values we have defined:

```
>>> seeds = sp.zeros_like(pn.Ps, dtype=float)
>>> seeds[pn.pores(geom1.name)] = geom1['pore.seed']
>>> seeds[pn.pores(geom2.name)] = geom2['pore.seed']
>>> assert sp.all(seeds > 0) # Ensure all zeros are overwritten
```

The following code illustrates the shortcut approach, which accomplishes the same result as above in a single line:

```
>>> seeds = pn['pore.seed']
```

- This shortcut works because the
`pn`

dictionary does not contain an array called`'pore.seed'`

, so all associated**Geometry**objects are then checked for the requested array(s). If it is found, then OpenPNM essentially performs the*interleaving*of the data as demonstrated by the manual approach and returns all the values together in a single full-size array. If it is not found, then a standard*KeyError*message is received. - This exchange of data between
**Network**and**Geometry**makes sense if you consider that**Network**objects act as a sort of master object relative**Geometry**objects.**Networks**apply to*all*pores and throats in the domain, while**Geometries**apply to subsets of the domain, so if the**Network**needs some values from all pores it has direct access.

### Add Pore Size Distribution Models to Each Geometry¶

Pore-scale models are mathematical functions that are applied to each pore (or throat) in the network to produce some local property value. Each of the modules in OpenPNM (Network, Geometry, Phase and Physics) have a “library” of pre-written models located under “models” (i.e. *Geometry.models*). Below this level, the models are further categorized according to what property they calculate, and there are typical 2-3 models for each. For instance, under `Geometry.models.pore_diameter`

you will see `random`

, `normal`

and `weibull`

among others.

Pore size distribution models are assigned to each Geometry object as follows:

```
>>> geom1.models.add(propname='pore.diameter',
... model=OpenPNM.Geometry.models.pore_diameter.normal,
... scale=0.00002, loc=0.000001,
... seeds='pore.seed')
>>> geom2.models.add(propname='pore.diameter',
... model=OpenPNM.Geometry.models.pore_diameter.weibull,
... shape=1.2, scale=0.00004, loc=0.000001,
... seeds='pore.seed')
```

Pore-scale models tend to be the most complex (i.e. confusing) aspects of OpenPNM, so it’s worth dwelling on the important points of the above two commands:

- Both
`geom1`

and`geom2`

have a`models`

attribute where the parameters specified in the`add`

command are stored for future use if/when needed. The`models`

attribute actually contains a**ModelsDict**object which is a customized dictionary for storing and managing this type of information. - The
`propname`

argument specifies which property the model calculates. This means that the numerical results of the model calculation will be saved in their respective**Geometry**objects as`geom1['pore.diameter']`

and`geom2['pore.diameter']`

. - Each model stores it’s result under the same
`propname`

but these values do not conflict since each**Geometry**object presides over a unique subset of pores and throats. - The
`model`

argument contains a*handle*to the desired function, which is extracted from the*models*library of the relevant*Module*(**Geometry**in this case). Each**Geometry**object has been assigned a different statistical model,*normal*and*weibull*. This ability to apply different models to different regions of the domain is reason multiple**Geometry**objects are permitted. The added complexity is well worth the added flexibility. - The remaining arguments are those required by the chosen
*model*. In the above cases, these are the parameters that define the statistical distribution. Note that the mean pore size for`geom1`

will be 20 um (set by`scale`

) while for`geom2`

it will be 50 um, thus creating the smaller surface pores as intended. The pore-scale models are well documented regarding what arguments are required and their meaning; as usual these can be viewed with Object Inspector in Spyder.

Now that we’ve added pore diameter models the each **Geometry** we can visualize the network in Paraview to confirm that distinctly different pore sizes on the surface regions:

### Add Additional Pore-Scale Models to Each Geometry¶

In addition to pore diameter, there are several other geometrical properties needed to perform a permeability simulation. Let’s start with throat diameter:

```
>>> geom1.models.add(propname='throat.diameter',
... model=OpenPNM.Geometry.models.throat_misc.neighbor,
... pore_prop='pore.diameter',
... mode='min')
>>> geom2.models.add(propname='throat.diameter',
... model=OpenPNM.Geometry.models.throat_misc.neighbor,
... pore_prop='pore.diameter',
... mode='min')
```

Instead of using statistical distribution functions, the above lines use the `neighbor`

model which determines each throat value based on the values found `'pore_prop'`

from it’s neighboring pores. In this case, each throat is assigned the minimum pore diameter of it’s two neighboring pores. Other options for `mode`

include `'max'`

and `'mean'`

.

We’ll also need throat length as well as the cross-sectional area of pores and throats, for calculating the hydraulic conductance model later.

```
>>> geom1.models.add(propname='throat.length',
... model=OpenPNM.Geometry.models.throat_length.straight)
>>> geom2.models.add(propname='throat.length',
... model=OpenPNM.Geometry.models.throat_length.straight)
>>> geom1.models.add(propname='throat.area',
... model=OpenPNM.Geometry.models.throat_area.cylinder)
>>> geom2.models.add(propname='throat.area',
... model=OpenPNM.Geometry.models.throat_area.cylinder)
>>> geom1.models.add(propname='pore.area',
... model=OpenPNM.Geometry.models.pore_area.spherical)
>>> geom2.models.add(propname='pore.area',
... model=OpenPNM.Geometry.models.pore_area.spherical)
```

## Create a Phase Object and Assign Thermophysical Property Models¶

For this tutorial, we will create a generic **Phase** object for water, then assign some pore-scale models for calculating their properties. Alternatively, we could use the prewritten **Water** class included in OpenPNM, which comes complete with the necessary pore-scale models, but this would defeat the purpose of the tutorial.

```
>>> water = OpenPNM.Phases.GenericPhase(network=pn)
>>> air = OpenPNM.Phases.GenericPhase(network=pn)
```

Note that all **Phase** objects are automatically assigned standard temperature and pressure conditions when created. This can be adjusted:

```
>>> water['pore.temperature'] = 353 # K
```

A variety of pore-scale models are available for calculating **Phase** properties, generally taken from correlations in the literature. An empirical correlation specifically for the viscosity of water is available:

```
>>> water.models.add(propname='pore.viscosity',
... model=OpenPNM.Phases.models.viscosity.water)
```

## Create Physics Objects for Each Geometry¶

**Physics** objects are where geometric information and thermophysical properties are combined to produce the pore and throat scale transport parameters. Thus we need to create one **Physics** object for *EACH* **Phase** and *EACH* **Geometry**:

```
>>> phys1 = OpenPNM.Physics.GenericPhysics(network=pn, phase=water,
... geometry=geom1)
>>> phys2 = OpenPNM.Physics.GenericPhysics(network=pn, phase=water,
... geometry=geom2)
```

Next add the Hagan-Poiseuille model to both:

```
>>> mod = OpenPNM.Physics.models.hydraulic_conductance.hagen_poiseuille
>>> phys1.models.add(propname='throat.hydraulic_conductance', model=mod)
>>> phys2.models.add(propname='throat.hydraulic_conductance', model=mod)
```

- The same function (
`mod`

) was passed as the`model`

argument to both**Physics**objects. This means that both objects will calculate the hydraulic conductance using the same function. A model*must*be assigned to both objects in order for the`'throat.hydraulic_conductance'`

property be defined everywhere in the domain since each**Physics**applies to a unique selection of pores and throats. - The “pore-scale model” mechanism was specifically designed to allow for users to easily create their own custom models. Creating custom models is outlined in advanced_usage.

### Accessing Data Distributed Between Multiple Physics Objects¶

Just as **Network** objects can retrieve data from separate **Geometries** as a single array with values in the correct locations, **Phase** objects can retrieve data from **Physics** objects as follows:

```
>>> g = water['throat.hydraulic_conductance']
```

- Each
**Physics**applies to the same subset for pores and throats as the**Geometries**so its values are distributed spatially, but each**Physics**is also associated with a single**Phase**object. Consequently, only a**Phase**object can to request all of the values within the domain pertaining to itself. - In other words, a
**Network**object cannot aggregate the**Physics**data because it doesn’t know which**Phase**is referred to. For instance, when asking for`'throat.hydraulic_conductance'`

it could refer to water or air conductivity, so it can only be requested by water or air.

## Pore-Scale Models: The Big Picture¶

Having created all the necessary objects with pore-scale models, it is now time to demonstrate why the OpenPNM pore-scale model approach is so powerful. First, let’s inspect the current value of hydraulic conductance in throat 1 on `phys1`

and `phys2`

:

```
>>> g1 = phys1['throat.hydraulic_conductance'] # Save this for later
>>> g2 = phys2['throat.hydraulic_conductance'] # Save this for later
```

Now, let’s alter the **Geometry** objects by assigning new random seeds, and adjust the temperature of `water`

.

```
>>> geom1['pore.seed'] = sp.rand(geom1.Np)
>>> geom2['pore.seed'] = sp.rand(geom2.Np)
>>> water['pore.temperature'] = 370 # K
```

So far we have not run the `regenerate`

command on any of these objects, which means that the above changes have not yet been applied to all the dependent properties. Let’s do this and examine what occurs at each step:

```
>>> geom1.models.regenerate()
>>> geom2.models.regenerate()
```

These two lines trigger the re-calculation of all the size related models on each **Geometry** object.

```
>>> water.models.regenerate()
```

This line causes the viscosity to be recalculated at the new temperature. Let’s confirm that the hydraulic conductance has NOT yet changed since we have not yet regenerated the **Physics** objects’ models:

```
>>> sp.all(phys1['throat.hydraulic_conductance'] == g1) # g1 was saved above
True
>>> sp.all(phys2['throat.hydraulic_conductance'] == g2) # g2 was saved above
True
```

Finally, if we regenerate `phys1`

and `phys2`

we can see that the hydraulic conductance will be updated to reflect the new sizes on the **Geometries** *and* the new temperature on the **Phase**:

```
>>> phys1.models.regenerate()
>>> phys2.models.regenerate()
>>> sp.all(phys1['throat.hydraulic_conductance'] != g1)
True
>>> sp.all(phys2['throat.hydraulic_conductance'] != g2)
True
```

## Determine Permeability Tensor by Changing Inlet and Outlet Boundary Conditions¶

The getting started tutorial already demonstrated the process of performing a basic permeability simulation. In this tutorial, we’ll perform the simulation in all three perpendicular dimensions to obtain the permeability tensor of our heterogeneous anisotropic material.

```
>>> alg = OpenPNM.Algorithms.StokesFlow(network=pn, phase=water)
```

Set boundary conditions for flow in the X-direction:

```
>>> alg.set_boundary_conditions(bctype='Dirichlet', bcvalue=202650,
... pores=pn.pores('right'))
>>> alg.set_boundary_conditions(bctype='Dirichlet', bcvalue=101325,
... pores=pn.pores('left'))
>>> alg.run()
```

The resulting pressure field can be seen using Paraview:

To determine the permeability coefficient we must find the flow rate through the network to use in Darcy’s law. The **StokesFlow** class (and all analogous transport algorithms) possess a `rate`

method that calculates the net transport through a given set of pores:

```
>>> Q = alg.rate(pores=pn.pores('left'))
```

To find K, we need to solve Darcy’s law: Q = KA/(mu*L)(P_in - P_out). This requires knowing the viscosity and macroscopic network dimensions:

```
>>> mu = sp.mean(water['pore.viscosity'])
```

The dimensions of the network can be determined manually from the `shape`

and `spacing`

specified during its generation:

```
>>> L = 20 * 0.0001
>>> A = 20 * 10 * (0.0001**2)
```

The pressure drop was specified as 1 atm when setting boundary conditions, so `Kxx`

can be found as:

```
>>> Kxx = Q * mu * L / (A * 101325)
```

We can either create 2 new **Algorithm** objects to perform the simulations in the other two directions, or reuse `alg`

by adjusting the boundary conditions and re-running it.

```
>>> alg.set_boundary_conditions(bctype='Dirichlet', bcvalue=202650,
... pores=pn.pores('front'),
... mode='overwrite')
>>> alg.set_boundary_conditions(bctype='Dirichlet', bcvalue=101325,
... pores=pn.pores('back'),
... mode='merge')
>>> alg.run()
```

The first call to `set_boundary_conditions`

used the `overwrite`

mode, which replaces all existing boundary conditions on the `alg`

object with the specified values. The second call uses the `merge`

mode which adds new boundary conditions to any already present, which is the default behavior.

A new value for the flow rate must be recalculated, but all other parameters are equal to the X-direction:

```
>>> Q = alg.rate(pores=pn.pores('back'))
>>> Kyy = Q * mu * L / (A * 101325)
```

The values of `Kxx`

and `Kyy`

should be nearly identical since both these two directions are parallel to the small surface pores. For the Z-direction:

```
>>> alg.set_boundary_conditions(bctype='Dirichlet', bcvalue=202650,
... pores=pn.pores('top'),
... mode='overwrite')
>>> alg.set_boundary_conditions(bctype='Dirichlet', bcvalue=101325,
... pores=pn.pores('bottom'))
>>> alg.run()
>>> Q = alg.rate(pores=pn.pores('bottom'))
>>> L = 10 * 0.0001
>>> A = 20 * 20 * (0.0001**2)
>>> Kzz = Q * mu * L / (A * 101325)
```

The permeability in the Z-direction is about half that in the other two directions due to the constrictions caused by the small surface pores.

## Using the Workspace Manager to Save (and reload) the simulation¶

OpenPNM includes a **Workspace** manager that provides the type of functionality found on the *menu-bar* of a typical GUI-based application Specifically, this enables *saving* and *loading* of all active networks, or individual objects.

To use these feature it is necessary to instantiate an instance:

```
mgr = OpenPNM.Base.Workspace()
mgr.save('filename.pnm')
```

Some of the more common functions of the **Workspace** are available via short-cuts under the main package, such that `op.save`

is equivalent to calling `mgr.save`

.

The **Workspace** also offers a few other functions, such as `purge_object``

which removes an object from the simulation, including all traces of its labels and references on other objects. It is also possible to ``clear``

the entire workspace, which is useful for *clearing the slate* when importing a new network. Of course, there is also a ``load``

function to load saved *pnm* files:

```
mgr.clear()
mgr.load('filename.pnm')
```