Tutorial 3 - Advanced

Tutorial 3 of 3: Advanced Topics and Usage

Learning Outcomes

• Use different methods to add boundary pores to a network

• Manipulate network topology by adding and removing pores and throats

• Explore the ModelsDict design, including copying models between objects, and changing model parameters

• Write a custom pore-scale model and a custom Phase

• Access and manipulate objects associated with the network

• Combine multiple algorithms to predict relative permeability

Build and Manipulate Network Topology

For the present tutorial, we’ll keep the topology simple to help keep the focus on other aspects of OpenPNM.

import warnings
import numpy as np
import scipy as sp
import openpnm as op
%config InlineBackend.figure_formats = ['svg']
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(10)
ws = op.Workspace()
ws.settings['loglevel'] = 40
np.set_printoptions(precision=4)
pn = op.network.Cubic(shape=[10, 10, 10], spacing=0.00006, name='net')

Adding Boundary Pores

When performing transport simulations it is often useful to have ‘boundary’ pores attached to the surface(s) of the network where boundary conditions can be applied. When using the Cubic class, two methods are available for doing this: add_boundaries, which is specific for the Cubic class, and add_boundary_pores, which is a generic method that can also be used on other network types and which is inherited from GenericNetwork. The first method automatically adds boundaries to ALL six faces of the network and offsets them from the network by 1/2 of the value provided as the network spacing. The second method provides total control over which boundary pores are created and where they are positioned, but requires the user to specify to which pores the boundary pores should be attached to. Let’s explore these two options:

pn.add_boundary_pores(labels=['top', 'bottom'])

Let’s quickly visualize this network with the added boundaries:

fig, ax = plt.subplots(figsize=(5, 5))
op.topotools.plot_connections(pn, c='r', ax=ax)
op.topotools.plot_coordinates(pn, c='b', ax=ax)

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f19762ec970>

Adding and Removing Pores and Throats

OpenPNM uses a list-based data storage scheme for all properties, including topological connections. One of the benefits of this approach is that adding and removing pores and throats from the network is essentially as simple as adding or removing rows from the data arrays. The one exception to this ‘simplicity’ is that the 'throat.conns' array must be treated carefully when trimming pores, so OpenPNM provides the extend and trim functions for adding and removing, respectively. To demonstrate, let’s reduce the coordination number of the network to create a more random structure:

Ts = np.random.rand(pn.Nt) < 0.1  # Create a mask with ~10% of throats labeled True
op.topotools.trim(network=pn, throats=Ts)  # Use mask to indicate which throats to trim

When the trim function is called, it automatically checks the health of the network afterwards, so logger messages might appear on the command line if problems were found such as isolated clusters of pores or pores with no throats. This health check is performed by calling the Network’s check_network_health method which returns a HealthDict containing the results of the checks:

a = pn.check_network_health()
print(a)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Key                                 Value
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
headless_throats                    []
looped_throats                      []
isolated_pores                      [1010, 1012, 1015, 1040, 1059, 1061, 1067, 1071, 1075, 1080, 1101, 1105, 1114, 1120, 1136, 1141, 1146, 1152, 1153, 1159, 1170, 1184]
disconnected_pores                  [1010, 1012, 1015, 1040, 1059, 1061, 1067, 1071, 1075, 1080, 1101, 1105, 1114, 1120, 1136, 1141, 1146, 1152, 1153, 1159, 1170, 1184]
duplicate_throats                   []
bidirectional_throats               []
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

The HealthDict contains several lists including things like duplicate throats and isolated pores, but also a suggestion of which pores to trim to return the network to a healthy state. Also, the HealthDict has a health attribute that is False is any checks fail.

op.topotools.trim(network=pn, pores=a['trim_pores'])

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Input In [6], in <cell line: 1>()
----> 1 op.topotools.trim(network=pn, pores=a['trim_pores'])

KeyError: 'trim_pores'

Let’s take another look at the network to see the trimmed pores and throats:

fig, ax = plt.subplots(figsize=(5, 5))
op.topotools.plot_connections(pn, c='r', ax=ax)
op.topotools.plot_coordinates(pn, c='b', ax=ax)

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f19715098b0>

Define Geometry Objects

The boundary pores we’ve added to the network should be treated a little bit differently. Specifically, they should have no volume or length (as they are not physically representative of real pores). To do this, we create two separate Geometry objects, one for internal pores and one for the boundaries:

Ps = pn.pores('*boundary', mode='not')
Ts = pn.throats('*boundary', mode='not')
geom = op.geometry.SpheresAndCylinders(network=pn, pores=Ps, throats=Ts, name='intern')
Ps = pn.pores('*boundary')
Ts = pn.throats('*boundary')
boun = op.geometry.Boundary(network=pn, pores=Ps, throats=Ts, name='boun')

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [8], in <cell line: 6>()
      4 Ps = pn.pores('*boundary')
      5 Ts = pn.throats('*boundary')
----> 6 boun = op.geometry.Boundary(network=pn, pores=Ps, throats=Ts, name='boun')

AttributeError: module 'openpnm.geometry' has no attribute 'Boundary'

The SpheresAndCylinders class is preloaded with the pore-scale models to calculate all the necessary size information (pore diameter, pore.volume, throat lengths, throat.diameter, etc). The Boundary class is speciall and is only used for the boundary pores. In this class, geometrical properties are set to small fixed values such that they don’t affect the simulation results.

Define Multiple Phase Objects

In order to simulate relative permeability of air through a partially water-filled network, we need to create each Phase object. OpenPNM includes pre-defined classes for each of these common fluids:

air = op.phases.Air(network=pn)
water = op.phases.Water(network=pn)
water['throat.contact_angle'] = 110
water['throat.surface_tension'] = 0.072

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [9], in <cell line: 1>()
----> 1 air = op.phases.Air(network=pn)
      2 water = op.phases.Water(network=pn)
      3 water['throat.contact_angle'] = 110

AttributeError: module 'openpnm' has no attribute 'phases'

Aside: Creating a Custom Phase Class

In many cases you will want to create your own fluid, such as an oil or brine, which may be commonly used in your research. OpenPNM cannot predict all the possible scenarios, but luckily it is easy to create a custom Phase class as follows:

from openpnm.phases import GenericPhase

class Oil(GenericPhase):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        self.add_model(propname='pore.viscosity',
                       model=op.models.misc.polynomial,
                       prop='pore.temperature',
                       a=[1.82082e-2, 6.51E-04, -3.48E-7, 1.11E-10])
        self['pore.molecular_weight'] = 116  # g/mol

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Input In [10], in <cell line: 1>()
----> 1 from openpnm.phases import GenericPhase
      3 class Oil(GenericPhase):
      4     def __init__(self, **kwargs):

ModuleNotFoundError: No module named 'openpnm.phases'

• Creating a Phase class basically involves placing a series of self.add_model commands within the __init__ section of the class definition. This means that when the class is instantiated, all the models are added to itself (i.e. self).

• **kwargs is a Python trick that captures all arguments in a dict called kwargs and passes them to another function that may need them. In this case they are passed to the __init__ method of Oil’s parent by the super function. Specifically, things like name and network are expected.

• The above code block also stores the molecular weight of the oil as a constant value

• Adding models and constant values in this way could just as easily be done in a run script, but the advantage of defining a class is that it can be saved in a file (i.e. ‘my_custom_phases’) and reused in any project.

oil = Oil(network=pn)
print(oil)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [11], in <cell line: 1>()
----> 1 oil = Oil(network=pn)
      2 print(oil)

NameError: name 'Oil' is not defined

Define Physics Objects for Each Geometry and Each Phase

In the tutorial #2 we created two Physics object, one for each of the two Geometry objects used to handle the stratified layers. In this tutorial, the internal pores and the boundary pores each have their own Geometry, but there are two Phases, which also each require a unique Physics:

phys_water_internal = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom)
phys_air_internal = op.physics.GenericPhysics(network=pn, phase=air, geometry=geom)
phys_water_boundary = op.physics.GenericPhysics(network=pn, phase=water, geometry=boun)
phys_air_boundary = op.physics.GenericPhysics(network=pn, phase=air, geometry=boun)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [12], in <cell line: 1>()
----> 1 phys_water_internal = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom)
      2 phys_air_internal = op.physics.GenericPhysics(network=pn, phase=air, geometry=geom)
      3 phys_water_boundary = op.physics.GenericPhysics(network=pn, phase=water, geometry=boun)

NameError: name 'water' is not defined

To reiterate, one Physics object is required for each Geometry AND each Phase, so the number can grow to become annoying very quickly Some useful tips for easing this situation are given below.

Create a Custom Pore-Scale Physics Model

Perhaps the most distinguishing feature between pore-network modeling papers is the pore-scale physics models employed. Accordingly, OpenPNM was designed to allow for easy customization in this regard, so that you can create your own models to augment or replace the ones included in the OpenPNM models libraries. For demonstration, let’s implement the capillary pressure model proposed by Mason and Morrow in 1994. They studied the entry pressure of non-wetting fluid into a throat formed by spheres, and found that the converging-diverging geometry increased the capillary pressure required to penetrate the throat. As a simple approximation they proposed \(P_c = -2 \sigma \cdot cos(2/3 \theta) / R_t\)

Pore-scale models are written as basic function definitions:

def mason_model(target, diameter='throat.diameter', theta='throat.contact_angle',
                sigma='throat.surface_tension', f=0.6667):
    proj = target.project
    network = proj.network
    phase = proj.find_phase(target)
    Dt = network[diameter]
    theta = phase[theta]
    sigma = phase[sigma]
    Pc = 4*sigma*np.cos(f*np.deg2rad(theta))/Dt
    return Pc[phase.throats(target.name)]

Let’s examine the components of above code:

• The function receives a target object as an argument. This indicates which object the results will be returned to.

• The f value is a scale factor that is applied to the contact angle. Mason and Morrow suggested a value of 2/3 as a decent fit to the data, but we’ll make this an adjustable parameter with 2/3 as the default.

• Note the pore.diameter is actually a Geometry property, but it is retrieved via the network using the data exchange rules outlined in the second tutorial.

• All of the calculations are done for every throat in the network, but this pore-scale model may be assigned to a target like a Physics object, that is a subset of the full domain. As such, the last line extracts values from the Pc array for the location of target and returns just the subset.

• The actual values of the contact angle, surface tension, and throat diameter are NOT sent in as numerical arrays, but rather as dictionary keys to the arrays. There is one very important reason for this: if arrays had been sent, then re-running the model would use the same arrays and hence not use any updated values. By having access to dictionary keys, the model actually looks up the current values in each of the arrays whenever it is run.

• It is good practice to include the dictionary keys as arguments, such as sigma = 'throat.contact_angle'. This way the user can control where the contact angle could be stored on the target object.

Copy Models Between Physics Objects

As mentioned above, the need to specify a separate Physics object for each Geometry and Phase can become tedious. It is possible to copy the pore-scale models assigned to one object onto another object. First, let’s assign the models we need to phys_water_internal:

mod = op.models.physics.hydraulic_conductance.hagen_poiseuille
phys_water_internal.add_model(propname='throat.hydraulic_conductance',
                              model=mod)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [14], in <cell line: 2>()
      1 mod = op.models.physics.hydraulic_conductance.hagen_poiseuille
----> 2 phys_water_internal.add_model(propname='throat.hydraulic_conductance',
      3                               model=mod)

NameError: name 'phys_water_internal' is not defined

phys_water_internal.add_model(propname='throat.entry_pressure',
                              model=mason_model)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [15], in <cell line: 1>()
----> 1 phys_water_internal.add_model(propname='throat.entry_pressure',
      2                               model=mason_model)

NameError: name 'phys_water_internal' is not defined

Now make a copy of the models on phys_water_internal and apply it all the other water Physics objects:

phys_water_boundary.models = phys_water_internal.models

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [16], in <cell line: 1>()
----> 1 phys_water_boundary.models = phys_water_internal.models

NameError: name 'phys_water_internal' is not defined

The only ‘gotcha’ with this approach is that each of the Physics objects must be regenerated in order to place numerical values for all the properties into the data arrays:

phys_water_boundary.regenerate_models()
phys_air_internal.regenerate_models()
phys_air_internal.regenerate_models()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [17], in <cell line: 1>()
----> 1 phys_water_boundary.regenerate_models()
      2 phys_air_internal.regenerate_models()
      3 phys_air_internal.regenerate_models()

NameError: name 'phys_water_boundary' is not defined

Adjust Pore-Scale Model Parameters

The pore-scale models are stored in a ModelsDict object that is itself stored under the models attribute of each object. This arrangement is somewhat convoluted, but it enables integrated storage of models on the object’s wo which they apply. The models on an object can be inspected with print(phys_water_internal), which shows a list of all the pore-scale properties that are computed by a model, and some information about the model’s regeneration mode.

Each model in the ModelsDict can be individually inspected by accessing it using the dictionary key corresponding to pore-property that it calculates, i.e. print(phys_water_internal)['throat.capillary_pressure']). This shows a list of all the parameters associated with that model. It is possible to edit these parameters directly:

phys_water_internal.models['throat.entry_pressure']['f'] = 0.75  # Change value
phys_water_internal.regenerate_models()  # Regenerate model with new 'f' value

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [18], in <cell line: 1>()
----> 1 phys_water_internal.models['throat.entry_pressure']['f'] = 0.75  # Change value
      2 phys_water_internal.regenerate_models()

NameError: name 'phys_water_internal' is not defined

More details about the ModelsDict and ModelWrapper classes can be found in :ref:models.

Perform Multiphase Transport Simulations

Use the Built-In Drainage Algorithm to Generate an Invading Phase Configuration

inv = op.algorithms.Porosimetry(network=pn, phase=water)
inv.set_inlets(pores=pn.pores(['top', 'bottom']))
inv.run()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [19], in <cell line: 1>()
----> 1 inv = op.algorithms.Porosimetry(network=pn, phase=water)
      2 inv.set_inlets(pores=pn.pores(['top', 'bottom']))
      3 inv.run()

AttributeError: module 'openpnm.algorithms' has no attribute 'Porosimetry'

• The inlet pores were set to both 'top' and 'bottom' using the pn.pores method. The algorithm applies to the entire network so the mapping of network pores to the algorithm pores is 1-to-1.

• The run method automatically generates a list of 25 capillary pressure points to test, but you can also specify more pores, or which specific points to tests. See the methods documentation for the details.

• Once the algorithm has been run, the resulting capillary pressure curve can be viewed with plot_drainage_curve. If you’d prefer a table of data for plotting in your software of choice you can use get_drainage_data which prints a table in the console.

Set Pores and Throats to Invaded

After running, the mip object possesses an array containing the pressure at which each pore and throat was invaded, stored as 'pore.inv_Pc' and 'throat.inv_Pc'. These arrays can be used to obtain a list of which pores and throats are invaded by water, using Boolean logic:

Pi = inv['pore.invasion_pressure'] < 5000
Ti = inv['throat.invasion_pressure'] < 5000

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [20], in <cell line: 1>()
----> 1 Pi = inv['pore.invasion_pressure'] < 5000
      2 Ti = inv['throat.invasion_pressure'] < 5000

NameError: name 'inv' is not defined

The resulting Boolean masks can be used to manually adjust the hydraulic conductivity of pores and throats based on their phase occupancy. The following lines set the water filled throats to near-zero conductivity for air flow:

Ts = phys_water_internal.map_throats(~Ti, origin=water)
phys_water_internal['throat.hydraulic_conductance'][Ts] = 1e-20

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [21], in <cell line: 1>()
----> 1 Ts = phys_water_internal.map_throats(~Ti, origin=water)
      2 phys_water_internal['throat.hydraulic_conductance'][Ts] = 1e-20

NameError: name 'phys_water_internal' is not defined

• The logic of these statements implicitly assumes that transport between two pores is only blocked if the throat is filled with the other phase, meaning that both pores could be filled and transport is still permitted. Another option would be to set the transport to near-zero if either or both of the pores are filled as well.

• The above approach can get complicated if there are several Geometry objects, and it is also a bit laborious. There is a pore-scale model for this under Physics.models.multiphase called conduit_conductance. The term conduit refers to the path between two pores that includes 1/2 of each pores plus the connecting throat.

Calculate Relative Permeability of Each Phase

We are now ready to calculate the relative permeability of the domain under partially flooded conditions. Instantiate an StokesFlow object:

water_flow = op.algorithms.StokesFlow(network=pn, phase=water)
water_flow.set_value_BC(pores=pn.pores('left'), values=200000)
water_flow.set_value_BC(pores=pn.pores('right'), values=100000)
water_flow.run()
Q_partial, = water_flow.rate(pores=pn.pores('right'))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [22], in <cell line: 1>()
----> 1 water_flow = op.algorithms.StokesFlow(network=pn, phase=water)
      2 water_flow.set_value_BC(pores=pn.pores('left'), values=200000)
      3 water_flow.set_value_BC(pores=pn.pores('right'), values=100000)

NameError: name 'water' is not defined

The relative permeability is the ratio of the water flow through the partially water saturated media versus through fully water saturated media; hence we need to find the absolute permeability of water. This can be accomplished by regenerating the phys_water_internal object, which will recalculate the 'throat.hydraulic_conductance' values and overwrite our manually entered near-zero values from the inv simulation using phys_water_internal.models.regenerate(). We can then re-use the water_flow algorithm:

phys_water_internal.regenerate_models()
water_flow.run()
Q_full, = water_flow.rate(pores=pn.pores('right'))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [23], in <cell line: 1>()
----> 1 phys_water_internal.regenerate_models()
      2 water_flow.run()
      3 Q_full, = water_flow.rate(pores=pn.pores('right'))

NameError: name 'phys_water_internal' is not defined

And finally, the relative permeability can be found from:

K_rel = Q_partial/Q_full
print(f"Relative permeability: {K_rel:.5f}")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [24], in <cell line: 1>()
----> 1 K_rel = Q_partial/Q_full
      2 print(f"Relative permeability: {K_rel:.5f}")

NameError: name 'Q_partial' is not defined

Notes: * The ratio of the flow rates gives the normalized relative permeability since all the domain size, viscosity and pressure differential terms cancel each other. * To generate a full relative permeability curve the above logic would be placed inside a for loop, with each loop increasing the pressure threshold used to obtain the list of invaded throats (Ti). * The saturation at each capillary pressure can be found be summing the pore and throat volume of all the invaded pores and throats using Vp = geom['pore.volume'][Pi] and Vt = geom['throat.volume'][Ti].