#!/usr/bin/env python
# coding: utf-8

# <center>
# <h1>OpenPNM Jupyter Notebook Demo</h1>
# <center>by</center>
# <center>Andreas Putz</center>
# </center>

# 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
get_ipython().run_line_magic('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)


# # 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()


# 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)


# 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)


# ## 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


# 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())


# 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')


# 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)


# In[18]:


delaunay['geo'].plot_histograms()
delaunay['geo'].plot_fibre_slice(plane=[0.5,0,0])
delaunay['geo'].plot_porosity_profile();


# In[19]:


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

print("="*50,'\nGeometry Properties')
print(delaunay['geo'].props())


# ## Setup Phases and Physics

# In[20]:


# Define phases
delaunay['Hg']  = pnm.Phases.Mercury(network=delaunay['network'],name='Delaunay:Phase:Hg')
delaunay['Air'] = pnm.Phases.Air(network=delaunay['network'],name='Delaunay:Phase:Air')
delaunay['H2O'] = pnm.Phases.Water(network=delaunay['network'],name='Delaunay:Phase:H2O')


# In[21]:


Ps = delaunay['network'].pores()
Ts = delaunay['network'].throats()


# In[22]:


delaunay['phys_Hg'] = pnm.Physics.Standard(network=delaunay['network'], phase=delaunay['Hg'],
                                           pores=Ps, throats=Ts,
                                           dynamic_data=True,
                                           name='standard_Hg_physics')
delaunay['phys_Hg'].props()


# In[23]:


delaunay['phys_H2O'] = pnm.Physics.Standard(network=delaunay['network'], phase=delaunay['H2O'],
                                            pores=Ps, throats=Ts,
                                            dynamic_data=True,
                                            name='standard_H2O_physics')
delaunay['phys_H2O'].props()


# In[24]:


delaunay['phys_Air'] = pnm.Physics.Standard(network=delaunay['network'], phase=delaunay['Air'],
                                            pores=Ps, throats=Ts,
                                            dynamic_data=True,
                                            name='standard_Air_physics')
delaunay['phys_Air'].props()


# ## MIP

# In[25]:


delaunay['MIP'] = pnm.Algorithms.Drainage(network=delaunay['network'],name='Delaunay:Alg:MIP')
delaunay['MIP'].setup(invading_phase=delaunay['phys_Hg'], defending_phase=delaunay['phys_Air'])
delaunay['MIP'].set_inlets(pores=delaunay['network'].pores(['boundary']))


# In[26]:


delaunay['MIP'].run()
fig = delaunay['MIP'].plot_drainage_curve();


# In[27]:


df = pd.DataFrame.from_dict(delaunay['MIP'].get_drainage_data())
df


# In[28]:


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()


# In[29]:


delaunay['MIP'].return_results(Pc=delaunay['MIP'].get_drainage_data()['capillary_pressure'][-1])
pnm.export_data(network=delaunay['network'], filename='Voronoy_MIP')
df.to_csv('Voronoy_MIP.csv')


# ## Ordinary Percolation

# In[30]:


inlets = delaunay['network'].pores('bottom_boundary')
# Using every other pore in the bottom and boundary as an inlet prevents
# extremely small diffusivity and permeability values in the z direction
used_inlets = [inlets[x] for x in range(0, len(inlets), 2)]

delaunay['Alg:OP'] = pnm.Algorithms.OrdinaryPercolation(network=delaunay['network'],
                                                        invading_phase=delaunay['H2O'],
                                                        defending_phase=delaunay['Air'],
                                                        name='Delaunay:Alg:OP')
delaunay['Alg:OP'].set_inlets(pores=used_inlets)
delaunay['Alg:OP'].run(npts=100)


# In[31]:


delaunay['Alg:OP'].plot_drainage_curve();


# In[32]:


delaunay['Alg:OP'].plot_primary_drainage_curve();


# In[33]:


delaunay['Alg:OP'].return_results(sat=1.0)


# In[34]:


fig=plt.figure(figsize=(9,7))
pnm.Postprocessing.Plots.drainage_curves(delaunay['Alg:OP'], timing=None, fig=fig);


# In[35]:


pnm.export_data(network=delaunay['network'], filename='Voronoy_OP')