foo = dict()  # Create an empty dict
foo['bar'] = 1  # Store an integer under the key 'bar'
print(foo['bar'])  # Retrieve the integer stored in 'bar'

"""
Explanation: Tutorial 1 of 3: Getting Started with OpenPNM

This tutorial is intended to show the basic outline of how OpenPNM works, and necessarily skips many of the more useful and powerful features of the package.  So if you find yourself asking "why is this step so labor intensive" it's probably because this tutorial deliberately simplifies some features to provide a more smooth introduction.  The second and third tutorials of this User-Guide dive into the package more deeply, but those features are best appreciated once the basics are understood.

Learning Objectives

Introduce the main OpenPNM objects and their roles
Explore the way OpenPNM stores data, including network topology
Learn some handy tools for working with objects
Generate a standard cubic Network topology
Calculate geometrical properties and assign them to a Geometry object
Calculate thermophysical properties and assign to a Phase object
Define pore-scale physics and assign transport parameters to a Physics object
Run a permeability simulation using the pre-defined Algorithm
Use the package to calculate the permeability coefficient of a porous media


Python and Numpy Tutorials
Before diving into OpenPNM it is probably a good idea to become familar with Python and Numpy. The following resources should be helpful.
* OpenPNM is written in Python.  One of the best guides to learning Python is the set of Tutorials available on the official Python website. The web is literally overrun with excellent Python tutorials owing to the popularity and importance of the language.  The official Python website also provides a long list of resources
* For information on using Numpy, Scipy and generally doing scientific computing in Python checkout the Scipy lecture notes.  The Scipy website also offers as solid introduction to using Numpy arrays.
* The Stackoverflow website is an incredible resource for all computing related questions, including simple usage of Python, Scipy and Numpy functions.
* For users more familiar with Matlab, there is a Matlab-Numpy cheat sheet that explains how to translate familiar Matlab commands to Numpy.

Overview of Data Storage in OpenPNM
Before creating an OpenPNM simulation it is necessary to give a quick description of how data is stored in OpenPNM; after all, a significant part of OpenPNM is dedicated to data storage and handling.
Python Dictionaries or dicts
OpenPNM employs 5 main objects which each store and manage a different type of information or data:

Network: Manages topological data such as pore spatial locations and pore-to-pore connections
Geometry: Manages geometrical properties such as pore diameter and throat length
Phase: Manages thermophysical properties such as temperature and viscosity
Physics: Manages pore-scale transport parameters such as hydraulic conductance
Algorithm: Contains algorithms that use the data from other objects to perform simulations, such as diffusion or drainage

We will encounter each of these objects in action before the end of this tutorial.
Each of the above objects is a subclass of the Python dictionary or dict, which is a very general storage container that allows values to be accessed by a name using syntax like:
End of explanation
"""


import openpnm as op
import scipy as sp
net = op.network.GenericNetwork(Np=10, Nt=10)  # Instantiate an empty network object with 10 pores and 10 throats
net['pore.foo'] = sp.ones([net.Np, ])  # Assign an Np-long array of ones
net['pore.bar'] = range(0, net.Np)  # Assign an Np-long array of increasing ints
print(type(net['pore.bar']))  # The Python range iterator is converted to a proper Numpy array

net['pore.foo'][4] = 44.0  # Overwrite values in the array
print(net['pore.foo'][4])  # Retrieve values from the array

print(net['pore.foo'][2:6])  # Extract a slice of the array

print(net['pore.foo'][[2, 4, 6]])  # Extract specific locations

net['throat.foo'] = 2  # Assign a scalar
print(len(net['throat.foo']))  # The scalar values is converted to an Nt-long array

print(net['throat.foo'][4])  # The scalar value was placed into all locations

"""
Explanation: A detailed tutorial on dictionaries can be found here &lt;http://learnpythonthehardway.org/book/ex39.html&gt;_.  The dict does not offer much functionality aside from basic storage of arbitrary objects, and it is meant to be extended.  OpenPNM extends the dict to have functionality specifically suited for dealing with OpenPNM data.  More information about the functionality of OpenPNM's subclassed dicts can be found in the :ref:overall_design.
Numpy Arrays of Pore and Throat Data
All data are stored in arrays which can accessed using standard array syntax.  More details on the data storage scheme are given in :ref:data_storage, but the following gives a quick overview:

All pore and throat properties are stored in Numpy arrays.  All data will be automatically converted to a Numpy array if necessary.
The data for pore i (or throat i) can be found in element of i of an array.  This means that pores and throat have indices which are implied by their position in arrays.  When we speak of retrieving pore locations, it refers to the indices in the Numpy arrays.
Each property is stored in it's own array, meaning that 'pore diameter' and 'throat volume' are each stored in a separate array.
Arrays that store pore data are Np-long, while arrays that store throat data are Nt-long, where Np is the number of pores and Nt is the number of throats in the network.
Arrays can be any size in the other dimensions.  For instance, triplets of pore coordinates (i.e. [x, y, z]) can be stored for each pore creating an Np-by-3 array.
The storage of topological connections is also very nicely accomplished with this 'list-based' format, by creating an array ('throat.conns') that stores which pore indices are found on either end of a throat.  This leads to an Nt-by-2 array.  

OpenPNM Objects: Combining dicts and Numpy Arrays
OpenPNM objects combine the above two levels of data storage, meaning they are dicts that are filled with Numpy arrays.  OpenPNM enforces several rules to help maintain data consistency:

When storing arrays in an OpenPNM object, their name (or dictionary key) must be prefixed with 'pore.' or 'throat.'.
OpenPNM uses the prefix of the dictionary key to infer how long the array must be.
The specific property that is stored in each array is indicated by the suffix such as 'pore.diameter' or 'throat.length'.
Writing scalar values to OpenPNM objects automatically results in conversion to a full length array filled with the scalar value.
Arrays containing Boolean data are treated as labels, which are explained later in this tutorial.

The following code snippets give examples of how all these pieces fit together using an Empty network as an example:
End of explanation
"""


import openpnm as op
import scipy as sp

"""
Explanation: Now that we have a better idea of how data is stored and retrieved from OpenPNM's objects, we can begin an actual simulation.
Generate a Cubic Network
Now that we have seen the rough outline of how OpenPNM objects store data, we can begin building a simulation.  Start by importing OpenPNM and the Scipy package:
End of explanation
"""


pn = op.network.Cubic(shape=[4, 3, 1], spacing=0.0001)

"""
Explanation: Next, generate a Network by choosing the Cubic class, then create an instance with the desired parameters:
End of explanation
"""


print('The total number of pores on the network is:', pn.num_pores())
print('A short-cut to the total number of pores is:', pn.Np)
print('The total number of throats on the network is:', pn.num_throats())
print('A short-cut to the total number of throats is:', pn.Nt)
print('A list of all calculated properties is availble with:\n', pn.props())

"""
Explanation: The Network object stored in pn contains pores at the correct spatial positions and connections between the pores according the cubic topology.

The shape argument specifies the number of pores in the [X, Y, Z] directions of the cube.  Networks in OpenPNM are always 3D dimensional, meaning that a 2D or "flat" network is still 1 layer of pores "thick" so [X, Y, Z] = [20, 10, 1], thus pn in this tutorial is 2D which is easier for visualization.
The spacing argument controls the center-to-center distance between pores and it can be a scalar or vector (i.e. [0.0001, 0.0002, 0.0003]).

The resulting network looks like:

This image was creating using Paraview, using the instructions given here: Example in the OpenPNM-Example collection
Inspecting Object Properties
OpenPNM objects have additional methods for querying their relevant properties, like the number of pores or throats, which properties have been defined, and so on:
End of explanation
"""


print(pn.pores('left'))

"""
Explanation: Accessing Pores and Throats
One simple but important feature of OpenPNM is the ability to label pores and throats.  When a Cubic network is created, several labels are automatically created: the pores on each face are labeled 'left', 'right', etc.  These labels can be used as follows:
End of explanation
"""


print(pn['pore.coords'][pn.pores('left')])

"""
Explanation: The ability to retrieve pore indices is handy for querying pore properties, such as retrieving the pore coordinates of all pores on the 'left' face:
End of explanation
"""


print(pn.labels())

"""
Explanation: A list of all labels currently assigned to the network can be obtained with:
End of explanation
"""


geom = op.geometry.GenericGeometry(network=pn, pores=pn.Ps, throats=pn.Ts)

"""
Explanation: Create a Geometry Object and Assign Geometric Properties to Pores and Throats
The Network pn does not contain any information about pore and throat sizes at this point.  The next step is to create a Geometry object to manage the geometrical properties.
End of explanation
"""


geom['pore.diameter'] = sp.rand(pn.Np)*0.0001  # Units of meters

"""
Explanation: This statement contains three arguments:

network tells the Geometry object which Network it is associated with.  There can be multiple networks defined in a given session, so all objects must be associated with a single network.
pores and throats indicate the locations in the Network where this Geometry object will apply.  In this  tutorial geom applies to all pores and throats, but there are many cases where different regions of the network have different geometrical properties, so OpenPNM allows multiple Geometry objects to be created for managing the data in each region, but this is a subject for :ref:intermediate_usage.

Add Pore and Throat Size Information
This freshly instantiated Geometry object (geom) contains no geometric properties as yet.  For this tutorial we'll use the direct assignment of manually calculated values.
We'll start by assigning diameters to each pore from a random distribution, spanning 0 um to 100 um.  The upper limit matches the spacing of the Network which was set to 0.0001 m (i.e. 100 um), so pore diameters exceeding 100 um might overlap with their neighbors.  Using the Scipy rand function creates an array of random numbers between 0 and 0.0001 that is Np-long, meaning each pore is assigned a unique random number
End of explanation
"""


P12 = pn['throat.conns']  # An Nt x 2 list of pores on the end of each throat
D12 = geom['pore.diameter'][P12]  # An Nt x 2 list of pore diameters
Dt = sp.amin(D12, axis=1)  # An Nt x 1 list of the smaller pore from each pair
geom['throat.diameter'] = Dt

"""
Explanation: We usually want the throat diameters to always be smaller than the two pores which it connects to maintain physical consistency. This requires understanding a little bit about how OpenPNM stores network topology.  Consider the following:
End of explanation
"""


Rp = geom['pore.diameter']/2
geom['pore.volume'] = (4/3)*3.14159*(Rp)**3

"""
Explanation: Let's dissect the above lines.

Firstly, P12 is a direct copy of the Network's 'throat.conns' array, which contains the indices of the pore-pair connected by each throat.
Next, this Nt-by-2 array is used to index into the 'pore.diameter' array, resulting in another Nt-by-2 array containing the diameters of the pores on each end of a throat.
Finally, the Scipy function amin is used to find the minimum diameter of each pore-pair by specifying the axis argument as 1, and the resulting Nt-by-1 array is assigned to geom['throat.diameter'].
This trick of using 'throat.conns' to index into a pore property array is commonly used in OpenPNM and you should have a second look at the above code to understand it fully.

We must still specify the remaining geometrical properties of the pores and throats. Since we're creating a "Stick-and-Ball" geometry, the sizes are calculated from the geometrical equations for spheres and cylinders.
For pore volumes, assume a sphere:
End of explanation
"""


C2C = 0.0001  # The center-to-center distance between pores
Rp12 = Rp[pn['throat.conns']]
geom['throat.length'] = C2C - sp.sum(Rp12, axis=1)

"""
Explanation: The length of each throat is the center-to-center distance between pores, minus the radius of each of two neighboring pores.
End of explanation
"""


Rt = geom['throat.diameter']/2
Lt = geom['throat.length']
geom['throat.volume'] = 3.14159*(Rt)**2*Lt

"""
Explanation: The volume of each throat is found assuming a cylinder:
End of explanation
"""


water = op.phases.GenericPhase(network=pn)

"""
Explanation: The basic geometrical properties of the network are now defined.  The Geometry class possesses a method called plot_histograms that produces a plot of the most pertinent geometrical properties.  The following figure doesn't look very good since the network in this example has only 12 pores, but the utility of the plot for quick inspection is apparent.

Create a Phase Object
The simulation is now topologically and geometrically defined.  It has pore coordinates, pore and throat sizes and so on.  In order to perform any simulations it is necessary to define a Phase object to manage all the thermophysical properties of the fluids in the simulation:
End of explanation
"""


water['pore.temperature'] = 298.0
water['pore.viscosity'] = 0.001

"""
Explanation: Some notes on this line:
* pn is passed as an argument because Phases must know to which Network they belong.
* Note that pores and throats are NOT specified; this is because Phases are mobile and can exist anywhere or everywhere in the domain, so providing specific locations does not make sense.  Algorithms for dynamically determining actual phase distributions are discussed later.
Add Thermophysical Properties
Now it is necessary to fill this Phase object with the desired thermophysical properties.  OpenPNM includes a framework for calculating thermophysical properties from models and correlations, but this is covered in :ref:intermediate_usage.  For this tutorial, we'll use the basic approach of simply assigning static values as follows:
End of explanation
"""


phys_water = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom)

"""
Explanation: The above lines utilize the fact that OpenPNM converts scalars to full length arrays, essentially setting the temperature in each pore to 298.0 K.

Create a Physics Object
We are still not ready to perform any simulations.  The last step is to define the desired pore-scale physics models, which dictate how the phase and geometrical properties interact to give the transport parameters.  A classic example of this is the Hagen-Poiseuille equation for fluid flow through a throat to predict the flow rate as a function of the pressure drop.  The flow rate is proportional to the geometrical size of the throat (radius and length) as well as properties of the fluid (viscosity) and thus combines geometrical and thermophysical properties:
End of explanation
"""


R = geom['throat.diameter']/2
L = geom['throat.length']

"""
Explanation: As with all objects, the Network must be specified
Physics objects combine information from a Phase (i.e. viscosity) and a Geometry (i.e. throat diameter), so each of these must be specified.
Physics objects do not require the specification of which pores and throats where they apply, since this information is implied by the geometry argument which was already assigned to specific locations.

Specify Desired Pore-Scale Transport Parameters
We need to calculate the numerical values representing our chosen pore-scale physics.  To continue with the Hagen-Poiseuille example lets calculate the hydraulic conductance of each throat in the network.  The throat radius and length are easily accessed as:
End of explanation
"""


mu_w = 0.001
phys_water['throat.hydraulic_conductance'] = 3.14159*R**4/(8*mu_w*L)

"""
Explanation: The viscosity of the Phases was only defined in the pores; however, the hydraulic conductance must be calculated for each throat.  There are several options, but to keep this tutorial simple we'll create a scalar value:
End of explanation
"""


alg = op.algorithms.StokesFlow(network=pn)
alg.setup(phase=water)

"""
Explanation: Numpy arrays support vectorization, so since both L and R are arrays of Nt-length, their multiplication in this way results in another array that is also Nt-long.
Create an Algorithm Object for Performing a Permeability Simulation
Finally, it is now possible to run some useful simulations.  The code below estimates the permeability through the network by applying a pressure gradient across and calculating the flux.  This starts by creating a StokesFlow algorithm, which is pre-defined in OpenPNM:
End of explanation
"""


BC1_pores = pn.pores('front')
alg.set_value_BC(values=202650, pores=BC1_pores)
BC2_pores = pn.pores('back') 
alg.set_value_BC(values=101325, pores=BC2_pores)

"""
Explanation: Like all the above objects, Algorithms must be assigned to a Network via the network argument.
This algorithm is also associated with a Phase object, in this case water, which dictates which pore-scale Physics properties to use (recall that phys_water was associated with water).  This can be passed as an argument to the instantiation or to the setup function.

Next the boundary conditions are applied using the set_boundary_conditions method on the Algorithm object.  Let's apply a 1 atm pressure gradient between the left and right sides of the domain:
End of explanation
"""


alg.run()

"""
Explanation: To actually run the algorithm use the run method:
End of explanation
"""


Q = alg.rate(pores=pn.pores('front'))
A = 0.0001*3*1  # Cross-sectional area for flow
L = 0.0001*4  # Length of flow path
del_P = 101325  # Specified pressure gradient
K = Q*mu_w*L/(A*del_P)
print(K)

"""
Explanation: This builds the coefficient matrix from the existing values of hydraulic conductance, and inverts the matrix to solve for pressure in each pore, and stores the results within the Algorithm's dictionary under 'pore.pressure'.
To determine the permeability coefficient, we must invoke Darcy's law: Q = KA/uL(Pin - Pout).  Everything in this equation is known except for the volumetric flow rate Q.  The StokesFlow algorithm possesses a rate method that calculates the rate of a quantity leaving a specified set of pores:
End of explanation
"""


K = alg.calc_effective_permeability(domain_area=A, domain_length=L)
print(K)

"""
Explanation: The StokesFlow class was developed with permeability simulations in mind, so a specific method is available for determining the permeability coefficient that essentially applies the recipe from above.  This method could struggle with non-uniform geometries though, so use with caution:
End of explanation
"""


water.update(alg.results())

"""
Explanation: The results ('pore.pressure') are held within the alg object and must be explicitly returned to the Phase object by the user if they wish to use these values in a subsequent calculation.  The point of this data containment is to prevent unintentional overwriting of data.  Each algorithm has a method called results which returns a dictionary of the pertinent simulation results, which can be added to the phase of interest using the update method.
End of explanation
"""


