OpenPNM Jupyter Notebook Demo

by
Andreas Putz
In [1]:
import IPython
from IPython.display import HTML, display

Notebook Setup

This section deals with the module imports for this notebook. This section needs to complete for the notebook to work correctly.

My anaconda setup:

  • Python 3.5 environment
  • Anaconda notebook extension + community notebook extensions
  • conda-forge channel activated

Additional installation steps:

pip install openpnm or OpenPNM
In [2]:
import IPython
from IPython.display import HTML, display
import datetime

from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

import scipy as sp
import pandas as pd

import pylab as plt
import matplotlib.mlab as mlab
%matplotlib notebook

import urllib.request
import os
import sys
import gzip
import zipfile
import shutil

try:
    import OpenPNM as pnm
except Exception as ex:
    print('OpenPNM is not installed. Please us pip or conda to install.')
    print('='*50)
    print(ex)
    
In [3]:
print('Notebook Executed:\t ' + str(datetime.datetime.now()))
print('='*80)
print('Python Version:')
print('-'*80)
print(sys.version)
print('-'*80)
print('OpenPNM Version:', pnm.__version__)
print('='*80)
Notebook Executed:	 2017-09-08 09:37:38.187090
================================================================================
Python Version:
--------------------------------------------------------------------------------
3.5.3 |Anaconda custom (64-bit)| (default, Mar  6 2017, 11:58:13) 
[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
--------------------------------------------------------------------------------
OpenPNM Version: 1.6.2
================================================================================

OpenPNM Basics - MIP on a simple cubic network

Create basic network

In [4]:
# Generate a cubic network and define a simple sitik an ball geometry
pnm.clear()
cubic = {}
cubic['network'] = pnm.Network.Cubic(shape=[10, 10, 10], spacing=0.0001, name= 'cubic:TOP')
cubic['geo']     = pnm.Geometry.Stick_and_Ball(network=cubic['network'], 
                                               pores=cubic['network'].Ps, throats=cubic['network'].Ts,
                                               name='cubic:GEO')
In [5]:
cubic['geo'].plot_histograms()
In [6]:
cubic['geo'].labels()
Out[6]:
['pore.all', 'throat.all']
In [7]:
fig = pnm.Network.tools.plot_coordinates(network=cubic['geo'],pores=cubic['geo'].pores('internal'))
fig = pnm.Network.tools.plot_coordinates(network=cubic['geo'],pores=cubic['geo'].pores('left'),c='r',fig=fig)
fig = pnm.Network.tools.plot_coordinates(network=cubic['geo'],pores=cubic['geo'].pores('right'),c='g',fig=fig)
2017-09-08 09:37:39,812 | WARNING  | root._parse_labels | 'pore.internal' not found
2017-09-08 09:37:40,177 | WARNING  | root._parse_labels | 'pore.left' not found
2017-09-08 09:37:40,179 | WARNING  | OpenPNM.Network.tools._scale_3d_axes | Axes is already scaled to previously plotted data
2017-09-08 09:37:40,184 | WARNING  | root._parse_labels | 'pore.right' not found
2017-09-08 09:37:40,186 | WARNING  | OpenPNM.Network.tools._scale_3d_axes | Axes is already scaled to previously plotted data
In [8]:
fig = pnm.Network.tools.plot_connections(network=cubic['network'])
fig = pnm.Network.tools.plot_coordinates(network=cubic['network'],pores=cubic['network'].pores('internal'), fig = fig)
fig = pnm.Network.tools.plot_coordinates(network=cubic['network'],pores=cubic['network'].pores('left'),c='r',fig=fig)
fig = pnm.Network.tools.plot_coordinates(network=cubic['network'],pores=cubic['network'].pores('right'),c='g',fig=fig)
2017-09-08 09:37:40,256 | WARNING  | OpenPNM.Network.tools._scale_3d_axes | Axes is already scaled to previously plotted data
2017-09-08 09:37:40,263 | WARNING  | OpenPNM.Network.tools._scale_3d_axes | Axes is already scaled to previously plotted data
2017-09-08 09:37:40,267 | WARNING  | OpenPNM.Network.tools._scale_3d_axes | Axes is already scaled to previously plotted data

Setup and run MIP simulation

In [9]:
# Define phases
cubic['Hg']   = pnm.Phases.Mercury(network=cubic['network'],name='cubic:phase:Hg')
cubic['Air']  = pnm.Phases.Air(network=cubic['network'],name='cubic:phase:Air')
cubic['phys'] = pnm.Physics.Standard(network=cubic['network'], 
                                     phase=cubic['Hg'], pores=cubic['network'].Ps, throats=cubic['network'].Ts,
                                     name='cubic:PHYS')
In [10]:
# Setup MIP simulation
cubic['MIP'] = pnm.Algorithms.Drainage(network=cubic['network'],name='cubic:ALG:MIP')
cubic['MIP'].setup(invading_phase=cubic['Hg'], defending_phase=cubic['Air'])
cubic['MIP'].set_inlets(pores=cubic['network'].pores(['left']))
cubic['MIP'].run()
fig = cubic['MIP'].plot_drainage_curve();
In [11]:
df = pd.DataFrame.from_dict(cubic['MIP'].get_drainage_data())
df
Out[11]:
capillary_pressure defending_phase_saturation invading_phase_saturation
0 41177.466844 1.000000 0.000000
1 43625.292886 0.994956 0.005044
2 46218.631821 0.992188 0.007812
3 48966.133774 0.967121 0.032879
4 51876.963084 0.949058 0.050942
5 54960.828872 0.650863 0.349137
6 58228.017422 0.496649 0.503351
7 61689.426497 0.356912 0.643088
8 65356.601684 0.272725 0.727275
9 69241.774908 0.185201 0.814799
10 73357.905228 0.127435 0.872565
11 77718.722067 0.085269 0.914731
12 82338.771001 0.050703 0.949297
13 87233.462282 0.028109 0.971891
14 92419.122233 0.016376 0.983624
15 97913.047709 0.010057 0.989943
16 103733.563792 0.003991 0.996009
17 109900.084910 0.002000 0.998000
18 116433.179596 0.001099 0.998901
19 123354.639098 0.000699 0.999301
20 130687.550059 0.000699 0.999301
21 138456.371526 0.000699 0.999301
22 146687.016533 0.000443 0.999557
23 155406.938534 0.000204 0.999796
24 164645.222974 0.000000 1.000000
In [12]:
plt.figure()
plt.plot(df.capillary_pressure,df.invading_phase_saturation,'ro-',label='Invading Phase')
plt.plot(df.capillary_pressure,df.defending_phase_saturation,'gx--',label='Defending Phase')
plt.legend()

plt.xlabel('Capillary Pressure $[Pa]$')
plt.ylabel('Saturation $[1]$')
plt.grid()

Export results for PostProcessing in Paraview

In [13]:
cubic['MIP'].return_results(Pc=cubic['MIP'].get_drainage_data()['capillary_pressure'][-1])
pnm.export_data(network=cubic['network'], filename='cubic_MIP')
df.to_csv('cubic_MIP.csv')

Delaunay Networks

In [14]:
pnm.clear()
delaunay = {}

Create Topology and Geometry

In [15]:
# Generate Geometry
delaunay['network'] = pnm.Network.Delaunay(num_pores=200, domain_size=[1.5e-4, 1.5e-4, 1.5e-4],name='Delaunay:TOP')

print("="*50,'\nOriginal Network Lables')
print(delaunay['network'].labels())

delaunay['network'].add_boundaries()

fig = pnm.Network.tools.plot_connections(network = delaunay['network']);
fig = pnm.Network.tools.plot_coordinates(network = delaunay['network'],fig=fig)
fig = pnm.Network.tools.plot_coordinates(network = delaunay['network'],fig=fig,
                                         pores=delaunay['network'].pores('bottom_boundary'),c='r')
fig = pnm.Network.tools.plot_coordinates(network = delaunay['network'],fig=fig,
                                         pores=delaunay['network'].pores('top_boundary'),c='g')

print("="*50,'\nNetwork lables after boundaries were added')
print(delaunay['network'].labels())
================================================== 
Original Network Lables
------------------------------------------------------------
1	: pore.all
2	: throat.all
------------------------------------------------------------
2017-09-08 09:37:43,781 | WARNING  | OpenPNM.Network.tools._scale_3d_axes | Axes is already scaled to previously plotted data
2017-09-08 09:37:43,787 | WARNING  | OpenPNM.Network.tools._scale_3d_axes | Axes is already scaled to previously plotted data
2017-09-08 09:37:43,791 | WARNING  | OpenPNM.Network.tools._scale_3d_axes | Axes is already scaled to previously plotted data
================================================== 
Network lables after boundaries were added
------------------------------------------------------------
1	: pore.all
2	: pore.back_boundary
3	: pore.bottom_boundary
4	: pore.boundary
5	: pore.front_boundary
6	: pore.left_boundary
7	: pore.right_boundary
8	: pore.top_boundary
9	: throat.all
10	: throat.boundary
------------------------------------------------------------
In [16]:
# Generate Geometry (Pore radii + fibre radii)
fibre_rad = 5e-6

Ps = delaunay['network'].pores()
Ts = delaunay['network'].find_neighbor_throats(pores=Ps, mode='intersection', flatten=True)

delaunay['geo'] = pnm.Geometry.Voronoi(network=delaunay['network'], pores=Ps, throats=Ts,
                                       fibre_rad=fibre_rad, voxel_vol=True,
                                       vox_len=1e-6, name='Delaunay:GEO')
2017-09-08 09:37:54,183 | WARNING  | Delaunay:TOP.check_network_health | 38 pores have no neighbors
2017-09-08 09:37:54,185 | WARNING  | Delaunay:TOP.check_network_health | Isolated clusters exist in the network
2017-09-08 09:37:54,195 | WARNING  | OpenPNM.Network.tools.trim | Isolated pores exist!  Run check_network_health to ID                         which pores to remove.
2017-09-08 09:37:54,199 | WARNING  | Delaunay:TOP.check_network_health | 38 pores have no neighbors
2017-09-08 09:37:54,202 | WARNING  | Delaunay:TOP.check_network_health | Isolated clusters exist in the network
In [17]:
#plot all pores in the Voronoi geometry
import OpenPNM.Utilities.vertexops as vo

fig = vo.plot_pore(delaunay['geo'], delaunay['geo'].pores())
#fig = pnm.Network.tools.plot_connections(network = vor,fig=fig)
fig = pnm.Network.tools.plot_coordinates(network = delaunay['network'],
                                         fig=fig,pores=delaunay['network'].pores('boundary'),
                                         c='g')
throats = delaunay['network'].find_neighbor_throats(pores=[0])
#plot all throats connected to the first pore in the network
fig = vo.plot_throat(delaunay['geo'], throats)