This tutorial will follow the same outline as the getting_started, 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
- 2 Tutorial 2 of 3: Digging Deeper into OpenPNM
- 2.1 Building a Cubic Network
- 2.2 Initialize and Build Multiple Geometry Objects
- 2.3 Create a Phase Object and Assign Thermophysical Property Models
- 2.4 Create Physics Objects for Each Geometry
- 2.5 Pore-Scale Models: The Big Picture
- 2.6 Determine Permeability Tensor by Changing Inlet and Outlet Boundary Conditions
- 2.7 Using the Workspace Manager to Save (and reload) the simulation
- 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
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
connectivitybeing 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.
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.
geom1applies to only the pores on the top and bottom surfaces (automatically labeled ‘top’ and ‘bottom’ during the network generation step), while
geom2applies to the pores ‘not’ on the top and bottom surfaces.
- The assignment of throats is more complicated and illustrates the
find_neighbor_throatsmethod, 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 (
Ps2) are found; however, the
modeargument alters which throats are returned. The terms
'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 . A detailed description of this behavior is given in Representing Topology.
Each of the Geometry objects was assigned a
name during instantiation, and this is stored in the
>>> 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.
geom2). This is useful when interacting with the objects at the command line using
geom1.name, which will report
- 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
geom1can be quickly retrieved using
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.
In getting_started 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
>>> 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
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.Npreturns 1000 while
geom2.Npreturn 200 and 800 respectively.
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
pndictionary 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.
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
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:
modelsattribute where the parameters specified in the
addcommand are stored for future use if/when needed. The
modelsattribute actually contains a ModelsDict object which is a customized dictionary for storing and managing this type of information.
propnameargument 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
- Each model stores it’s result under the same
propnamebut these values do not conflict since each Geometry object presides over a unique subset of pores and throats.
modelargument 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
geom1will be 20 um (set by
scale) while for
geom2it 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:
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
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)
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)
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
modelargument 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.
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.
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
>>> 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
>>> 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.
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
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
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
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
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.
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
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: