Repository: PMEAL/OpenPNM
Pull request #368 (2015-05-19 13:46:23)

[WIP] 85 percent test coverage
The objective of this PR is to reach 85% test coverage before the release of v1.2. 

Please branch off of this branch and add tests as you see fit, then open a PR against this branch and we'll merge it. Once `85-percent-test-coverage` covered with 85% or greater, we'll merge into `develop`. 

Keep in mind if you are unsure what the test coverage is for each file, run `bin/test cov` and it will show the percentage break down. If you still need more information, please see [coveralls](https://coveralls.io/r/PMEAL/OpenPNM) for a line by line analysis of the source code.

/cc @PMEAL/openpnm-core 


Patch:
diff --git a/OpenPNM/Algorithms/__InvasionPercolationForImbibition__.py b/OpenPNM/Algorithms/__InvasionPercolationForImbibition__.py
index f0decc990..7aaaf222d 100644
--- a/OpenPNM/Algorithms/__InvasionPercolationForImbibition__.py
+++ b/OpenPNM/Algorithms/__InvasionPercolationForImbibition__.py
@@ -103,19 +103,6 @@ class InvasionPercolationForImbibition(InvasionPercolationTimed):
             IP_inv_seq      : 0 for uninvaded, simulation step for invaded
             IP_inv_time     : 0 for uninvaded, simulation time for invaded
 
-        Examples
-        --------
-        >>> import OpenPNM
-        >>> pn = OpenPNM.Network.TestNet()
-        >>> geo = OpenPNM.Geometry.TestGeometry(network=pn,pores=pn.pores(),throats=pn.throats())
-        >>> phase1 = OpenPNM.Phases.TestPhase(network=pn)
-        >>> phase1['pore.contact_angle'] = 20
-        >>> phase2 = OpenPNM.Phases.TestPhase(network=pn)
-        >>> phys1 = OpenPNM.Physics.TestPhysics(network=pn, phase=phase1,pores=pn.pores(),throats=pn.throats())
-        >>> phys2 = OpenPNM.Physics.TestPhysics(network=pn, phase=phase2,pores=pn.pores(),throats=pn.throats())
-        >>> from OpenPNM.Algorithms.__InvasionPercolationForImbibition__ import InvasionPercolationForImbibition
-        >>> IP = InvasionPercolationForImbibition(network=pn)
-
         Suggested Improvements ::
 
             a) Allow updating of cluster flow-rates (this will require a delta-t calculation at each step, instead of a total t calculation).
diff --git a/OpenPNM/Algorithms/__InvasionPercolationTimed__.py b/OpenPNM/Algorithms/__InvasionPercolationTimed__.py
index b8c062b49..41574a33e 100644
--- a/OpenPNM/Algorithms/__InvasionPercolationTimed__.py
+++ b/OpenPNM/Algorithms/__InvasionPercolationTimed__.py
@@ -25,23 +25,6 @@ class InvasionPercolationTimed(GenericAlgorithm):
     name : string
         The name this algorithm will go by
 
-    Examples
-    --------
-    >>> import OpenPNM
-    >>> pn = OpenPNM.Network.TestNet()
-    >>> geo = OpenPNM.Geometry.TestGeometry(network=pn,pores=pn.pores(),throats=pn.throats())
-    >>> phase1 = OpenPNM.Phases.TestPhase(network=pn)
-    >>> phase2 = OpenPNM.Phases.TestPhase(network=pn)
-    >>> phys1 = OpenPNM.Physics.TestPhysics(network=pn, phase=phase1,pores=pn.pores(),throats=pn.throats())
-    >>> phys2 = OpenPNM.Physics.TestPhysics(network=pn, phase=phase2,pores=pn.pores(),throats=pn.throats())
-    >>> IP = OpenPNM.Algorithms.InvasionPercolationTimed(network=pn)
-    >>> IP.run(invading_phase=phase1, defending_phase=phase2, inlets=pn.pores('top'), outlets=pn.pores('bottom'),report=0)
-         IP algorithm at 0 % completion at 0.0 seconds
-         IP algorithm at 100% completion at  0.0  seconds
-    >>> IP.return_results()
-    >>> max(phase1['pore.IP_inv_seq']) #unless something changed with our test objects, this should print "60"
-    60
-
     Suggested Improvements ::
 
         a) Allow updating of cluster flow-rates (this will require a delta-t calculation at each step, instead of a total t calculation).
diff --git a/OpenPNM/Base/__Core__.py b/OpenPNM/Base/__Core__.py
index b7eac0a6a..bdc369251 100644
--- a/OpenPNM/Base/__Core__.py
+++ b/OpenPNM/Base/__Core__.py
@@ -70,7 +70,7 @@ class Core(dict):
         # Enforce correct dict naming
         element = key.split('.')[0]
         if (element != 'pore') and (element != 'throat'):
-            print('Array name \''+key+'\' does not begin with \'pore\' or \'throat\'')
+            logger.error('Array name \''+key+'\' does not begin with \'pore\' or \'throat\'')
             return
         # Convert value to an ndarray
         value = sp.array(value,ndmin=1)
@@ -86,7 +86,7 @@ class Core(dict):
                     super(Core, self).__setitem__(key,value)
                 else:
                     logger.warning(key+' is already defined.')
-                    pass
+                    return
             else:
                 logger.debug(key+' is being defined.')
                 super(Core, self).__setitem__(key,value)
@@ -447,7 +447,7 @@ class Core(dict):
                     temp[i] = list(labels[arr[i,:]])
                 return temp
             else:
-                print('unrecognized mode')
+                logger.error('unrecognized mode:'+mode)
 
     def labels(self,element='',pores=[],throats=[],mode='union'):
         r"""
@@ -457,7 +457,12 @@ class Core(dict):
         ----------
         pores (or throats) : array_like
             The pores (or throats) whose labels are sought.  If left empty a
-            dictionary containing all pore and throat labels is returned.
+            list containing all pore and throat labels is returned.
+
+        element : string
+            Controls whether pore or throat labels are returned.  If empty then
+            both are returned.
+
         mode : string, optional
             Controls how the query should be performed
 
diff --git a/OpenPNM/Base/__ModelsDict__.py b/OpenPNM/Base/__ModelsDict__.py
index a8ca8bf62..d574cecea 100644
--- a/OpenPNM/Base/__ModelsDict__.py
+++ b/OpenPNM/Base/__ModelsDict__.py
@@ -4,11 +4,11 @@ ModelsDict:  Abstract Class for Containing Models
 ###############################################################################
 """
 import inspect
-import scipy as sp
 from collections import OrderedDict
 from OpenPNM.Base import logging, Controller
 logger = logging.getLogger()
 
+
 class ModelWrapper(dict):
     r"""
     Accepts a model from the OpenPNM model library, as well as all required
@@ -48,8 +48,8 @@ class ModelWrapper(dict):
                 if item not in defs.keys():
                     defs[item] = '---'
                 lines.append('{0:<20s} {1} / ({2})'.format(item,
-                                                          self[item],
-                                                          defs[item]))
+                                                           self[item],
+                                                           defs[item]))
         lines.append(horizontal_rule)
         return '\n'.join(lines)
 
@@ -91,6 +91,7 @@ class ModelWrapper(dict):
                             'method.')
         return master[0]
 
+
 class GenericModel(ModelWrapper):
     r"""
     This class was deprecated, and replaced by ModelWrapper.  Unfortunately,
@@ -100,6 +101,7 @@ class GenericModel(ModelWrapper):
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
 
+
 class ModelsDict(OrderedDict):
     r"""
     This custom dictionary stores the models that are associated with each
@@ -350,8 +352,7 @@ class ModelsDict(OrderedDict):
             raise Exception('More than one master found! This model dictionary '
                             'has been associated with multiple objects. To use the '
                             'same dictionary multiple times use the copy method.')
-        return master[0]
+        elif len(master) == 0:
+            raise Exception('ModelsDict has no master.')
 
-if __name__ == '__main__':
-    pn = OpenPNM.Network.TestNet()
-    a = ModelsDict()
+        return master[0]
diff --git a/OpenPNM/Geometry/models/pore_area.py b/OpenPNM/Geometry/models/pore_area.py
index 387942a2b..731ba908f 100644
--- a/OpenPNM/Geometry/models/pore_area.py
+++ b/OpenPNM/Geometry/models/pore_area.py
@@ -9,7 +9,19 @@ import scipy as _sp
 
 def spherical(geometry, pore_diameter='pore.diameter', **kwargs):
     r"""
-    Calculate cross-sectional area for a spherical pore
+    Calculate cross-sectional area assuming the pore body is a sphere
+
+    Parameters
+    ----------
+    geometry : OpenPNM Geometry Object
+        The Geometry object which this model is associated with. This controls
+        the length of the calculated array, and also provides access to other
+        necessary geometric properties.
+
+    pore_diameter : string
+        The dictionary key of the array on the Geometry object containing the
+        pore diameter values necessary to find the area.
+
     """
     diams = geometry[pore_diameter]
     value = _sp.constants.pi/4*(diams)**2
@@ -18,7 +30,20 @@ def spherical(geometry, pore_diameter='pore.diameter', **kwargs):
 
 def cubic(geometry, pore_diameter='pore.diameter', **kwargs):
     r"""
-    Calculate cross-sectional pore area for a cubic pore
+    Calculate cross-sectional area assuming the pore body is a cube
+
+    Parameters
+    ----------
+    geometry : OpenPNM Geometry Object
+        The Geometry object which this model is associated with. This controls
+        the length of the calculated array, and also provides access to other
+        necessary geometric properties.
+
+    pore_diameter : string
+        The dictionary key of the array on the Geometry object containing the
+        pore diameter values necessary to find the area.
+
+
     """
     diams = geometry[pore_diameter]
     value = diams**2
diff --git a/OpenPNM/Geometry/models/pore_centroid.py b/OpenPNM/Geometry/models/pore_centroid.py
index 6c9045670..0fab83a14 100644
--- a/OpenPNM/Geometry/models/pore_centroid.py
+++ b/OpenPNM/Geometry/models/pore_centroid.py
@@ -1,7 +1,7 @@
 r"""
 ===============================================================================
-pore_centroid -- N.B Neither method is accurate but is quicker than that implemented
-                 in pore.volume
+pore_centroid -- N.B Neither method is accurate but is quicker than that
+                 implemented in pore.volume
 ===============================================================================
 
 """
diff --git a/OpenPNM/Geometry/models/pore_diameter.py b/OpenPNM/Geometry/models/pore_diameter.py
index cd8e5c4eb..de8a77619 100644
--- a/OpenPNM/Geometry/models/pore_diameter.py
+++ b/OpenPNM/Geometry/models/pore_diameter.py
@@ -13,7 +13,30 @@ from scipy.spatial import ConvexHull
 def sphere(geometry, psd_name, psd_shape, psd_loc, psd_scale,
            pore_seed='pore.seed', psd_offset=0, **kwargs):
     r"""
-    Calculate pore diameter from seed values for a spherical pore body
+    Calculate pore diameter from given seed values.
+
+    Parameters
+    ----------
+    geometry : OpenPNM Geometry Object
+        The Geometry object which this model is associated with. This controls
+        the length of the calculated array, and also provides access to other
+        necessary geometric properties.
+
+    psd_name : string
+        The name of the statistical distribution to use. This model uses the
+        Scipy.stats module, so any of the distributions available there are
+        suitable options.
+
+    psd_shape, loc and scale : scalars
+        The parameters to send to the selected statistics model.  Most of the
+        Scipy.stats models accept the same keyword arguments.  Note that the
+        psd_ prefix is added by OpenPNM to indicate 'pore size distribution'.
+
+    psd_offset : scalar
+        Controls the minimum value in the pore size distribution by shifting
+        the entire set of values by the given offset.  This is useful for
+        avoiding pore sizes too close to zero.
+
     """
     import scipy.stats as spst
     prob_fn = getattr(spst, psd_name)
@@ -22,9 +45,51 @@ def sphere(geometry, psd_name, psd_shape, psd_loc, psd_scale,
     return value
 
 
+def equivalent_sphere(geometry, pore_volume='pore.volume', **kwargs):
+    r"""
+    Calculate pore diameter as the diameter of a sphere with an equivalent
+    volume.
+
+    Parameters
+    ----------
+    geometry : OpenPNM Geometry Object
+        The Geometry object which this model is associated with. This controls
+        the length of the calculated array, and also provides access to other
+        necessary geometric properties.
+
+    pore_volume : string
+        The dictionary key containing the pore volume values
+    """
+    from scipy.special import cbrt
+    pore_vols = geometry[pore_volume]
+    value = cbrt(6*pore_vols/_sp.pi)
+    return value
+
+
+def equivalent_cube(geometry, pore_volume='pore.volume', **kwargs):
+    r"""
+    Calculate pore diameter as the width of a cube with an equivalent volume.
+
+    Parameters
+    ----------
+    geometry : OpenPNM Geometry Object
+        The Geometry object which this model is associated with. This controls
+        the length of the calculated array, and also provides access to other
+        necessary geometric properties.
+
+    pore_volume : string
+        The dictionary key containing the pore volume values
+    """
+    from scipy.special import cbrt
+    pore_vols = geometry[pore_volume]
+    value = cbrt(pore_vols)
+    return value
+
+
 def voronoi(geometry, pore_volume='pore.volume', **kwargs):
     r"""
-    Calculate pore diameter from equivalent sphere - volumes must be calculated first
+    Calculate pore diameter from equivalent sphere - volumes must be calculated
+    first.
     """
     from scipy.special import cbrt
     pore_vols = geometry[pore_volume]
diff --git a/OpenPNM/Geometry/models/pore_misc.py b/OpenPNM/Geometry/models/pore_misc.py
index 57de99c36..df22a28d6 100644
--- a/OpenPNM/Geometry/models/pore_misc.py
+++ b/OpenPNM/Geometry/models/pore_misc.py
@@ -23,7 +23,7 @@ def random(geometry, seed=None, num_range=[0, 1], **kwargs):
     """
     range_size = num_range[1]-num_range[0]
     range_min = num_range[0]
-    _sp.random.seed(seed)
+    _sp.random.seed(seed=seed)
     value = _sp.random.rand(geometry.num_pores(),)
     value = value*range_size + range_min
     return value
diff --git a/OpenPNM/Geometry/models/pore_seed.py b/OpenPNM/Geometry/models/pore_seed.py
index 63c80142e..d4c1aa246 100644
--- a/OpenPNM/Geometry/models/pore_seed.py
+++ b/OpenPNM/Geometry/models/pore_seed.py
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 r"""
 ===============================================================================
-pore_seeds -- Methods for generating fields of values for use as seeds in
+pore_seed -- Methods for generating fields of values for use as seeds in
 statistical pore size distributions
 ===============================================================================
 
@@ -9,75 +9,29 @@ statistical pore size distributions
 import scipy as _sp
 
 
-def perlin_noise(geometry, freq=1, octaves=4, mode='classic', **kwargs):
+def random(geometry, seed=None, num_range=[0, 1], **kwargs):
     r"""
-    Generate pore seed values using the Perlin noise algorithm.  This approach
-    imparts some spatial clumpiness to the pore seeds.
+    Assign random number to pores, for use in statistical distributions that
+    return pore size
 
     Parameters
     ----------
-    freq, octaves : int
-        Parameters that control the qualities of the noise.  Lower frequency
-        results in more smaller clumps.  Higher octaves gives a more textured
-        noise.
-    mode : {'classic','simplex'}
-        The algorithm to use when generating the noise.
+    seed : int
+        The starting seed value to send to Scipy's random number generator.
+        The default is None, which means different distribution is returned
+        each time the model is run.
 
-        * 'classic' : (default) The original algorithm developed by Perlin
-        * 'simplex' : A newer algorithm that is supposedly faster and results
-        in a more natural texture.
-
-    Returns
-    -------
-    A list of pore seeds values between 0:1 that are spatially correlated (i.e.
-    similar values are clumped together)
-
-    Notes
-    -----
-    - This method uses image analysis type tools, so only works on Cubic
-    networks
-    - This method requires the 'noise' module is installed
-
-    Examples
-    --------
-    >>> import OpenPNM
-    >>> pn = OpenPNM.Network.Cubic(shape=[50, 50, 50])
-    >>> geom = OpenPNM.Geometry.GenericGeometry(network=pn,pores=pn.Ps,throats=pn.Ts)
-    >>> geom.add_model(propname='pore.seed',
-    ...                model=OpenPNM.Geometry.models.pore_seed.perlin_noise)
-    >>> im = pn.asarray(geom['pore.seed'])
-
-    Visualizing the end result can be done with:
-
-    ``matplotlib.pyplot.imshow(im[:, 25, :],interpolation='none')``
+    num_range : list
+        A two element list indicating the low and high end of the returned
+        numbers.
 
     """
-    from noise import pnoise3, snoise3
-    import scipy.stats as spst
-
-    net = geometry._net
-    if mode == 'classic':
-        model = pnoise3
-    elif mode == 'simplex':
-        model = snoise3
-    freq = freq * octaves
-    # The following will only work on Cubic networks
-    x = net._shape[0]
-    y = net._shape[1]
-    z = net._shape[2]
-    temp = _sp.ndarray((x, y, z))
-    for k in range(z):
-        for j in range(y):
-            for i in range(x):
-                temp[i, j, k] = model(i / freq, j / freq, k / freq, octaves) + 0.5
-    # Assuming that the noise is normally distributed, find seeds of that dist
-    temp = _sp.reshape(temp, (temp.size,))
-    x_mean = _sp.mean(temp)
-    x_sigma = _sp.sqrt(1/(temp.size-1)*_sp.sum((temp - x_mean)**2))
-    fn1 = spst.norm(loc=x_mean, scale=x_sigma)
-    values = fn1.cdf(temp)
-    values = values[geometry.map_pores(target=net, pores=geometry.Ps)]
-    return values.flatten()
+    range_size = num_range[1]-num_range[0]
+    range_min = num_range[0]
+    _sp.random.seed(seed=seed)
+    value = _sp.random.rand(geometry.num_pores(),)
+    value = value*range_size + range_min
+    return value
 
 
 def distance_from_inclusion(geometry, p, **kwargs):
@@ -115,7 +69,9 @@ def distance_from_inclusion(geometry, p, **kwargs):
 
     Visualizing the end result can be done with:
 
-    ``matplotlib.pyplot.imshow(im[:,25,:],interpolation='none')``
+    .. code-block:: python
+
+        matplotlib.pyplot.imshow(im[:,25,:],interpolation='none')
 
     """
     import scipy.ndimage as _spim
@@ -191,7 +147,9 @@ def spatially_correlated(geometry, network, weights=None, strel=None, **kwargs):
 
     Visualizing the end result can be done with:
 
-    ``matplotlib.pyplot.imshow(im[:, 25, :],interpolation='none')``
+    .. code-block:: python
+
+        matplotlib.pyplot.imshow(im[:, 25, :],interpolation='none')
 
     """
     import scipy.ndimage as spim
diff --git a/OpenPNM/Geometry/models/pore_surface_area.py b/OpenPNM/Geometry/models/pore_surface_area.py
index 50310d420..99e19144f 100644
--- a/OpenPNM/Geometry/models/pore_surface_area.py
+++ b/OpenPNM/Geometry/models/pore_surface_area.py
@@ -1,18 +1,47 @@
 r"""
 ===============================================================================
-Submodule -- throat_surface_area
+Submodule -- pore_surface_area
 ===============================================================================
 
 """
 import scipy as _sp
 
 
-def sphere(geometry, pore_diameter='pore.diameter', **kwargs):
+def sphere(geometry, network,
+           pore_diameter='pore.diameter',
+           throat_area='throat.area',
+           **kwargs):
     r"""
-    Calculate internal surface area for a spherical pore
+    Calculates internal surface area of pore bodies assuming they are spherical.
+    then subtracts the area of the neighboring throats in a crude way, by
+    simply considering the throat cross-sectional area, thus not accounting
+    for the actual curvature of the intersection.
+
+    Parameters
+    ----------
+    geometry : OpenPNM Geometry Object
+        The Geometry object which this model is associated with. This controls
+        the length of the calculated array, and also provides access to other
+        necessary geometric properties.
+
+    network : OpenPNM Network Object
+        The Network object associated with the Geometry.  This is needed to
+        provide some topological information such as throat connections, and
+        neighboring pores.
+
+    pore_diameter : string
+        The dictionary key to the pore diameter array.
+
+    throat_area : string
+        The dictioanry key to the throat area array.  Throat areas are needed
+        since their insection with the pore are removed from the computation.
+
     """
     R = geometry[pore_diameter]/2
-    value = 4*_sp.constants.pi*R**2
+    Asurf = 4*_sp.constants.pi*R**2
+    Tn = network.find_neighbor_throats(pores=geometry.Ps, flatten=False)
+    Tsurf = _sp.array([_sp.sum(network[throat_area][Ts]) for Ts in Tn])
+    value = Asurf - Tsurf
     return value
 
 
diff --git a/OpenPNM/Geometry/models/pore_volume.py b/OpenPNM/Geometry/models/pore_volume.py
index 474c05a6b..ae7843f8e 100644
--- a/OpenPNM/Geometry/models/pore_volume.py
+++ b/OpenPNM/Geometry/models/pore_volume.py
@@ -10,6 +10,46 @@ from scipy.spatial import Delaunay
 import OpenPNM.Utilities.misc as misc
 
 
+def sphere(geometry, pore_diameter='pore.diameter', **kwargs):
+    r"""
+    Calculate pore volume from diameter assuming a spherical pore body
+
+    Parameters
+    ----------
+    geometry : OpenPNM Geometry Object
+        The Geometry object which this model is associated with. This controls
+        the length of the calculated array, and also provides access to other
+        necessary geometric properties.
+
+    pore_diameter : string
+        The dictionary key to the pore diameter values
+
+    """
+    diams = geometry[pore_diameter]
+    value = _sp.pi/6*diams**3
+    return value
+
+
+def cube(geometry, pore_diameter='pore.diameter', **kwargs):
+    r"""
+    Calculate pore volume from diameter assuming a cubic pore body
+
+    Parameters
+    ----------
+    geometry : OpenPNM Geometry Object
+        The Geometry object which this model is associated with. This controls
+        the length of the calculated array, and also provides access to other
+        necessary geometric properties.
+
+    pore_diameter : string
+        The dictionary key to the pore diameter values
+
+    """
+    diams = geometry[pore_diameter]
+    value = diams**3
+    return value
+
+
 def _get_hull_volume(points):
     r"""
     Calculate the volume of a set of points by dividing the bounding surface
@@ -84,24 +124,6 @@ def _get_hull_volume(points):
     return hull_volume, hull_COM
 
 
-def sphere(geometry, pore_diameter='pore.diameter', **kwargs):
-    r"""
-    Calculate pore volume from diameter for a spherical pore body
-    """
-    diams = geometry[pore_diameter]
-    value = _sp.pi/6*diams**3
-    return value
-
-
-def cube(geometry, pore_diameter='pore.diameter', **kwargs):
-    r"""
-    Calculate pore volume from diameter for a cubic pore body
-    """
-    diams = geometry[pore_diameter]
-    value = diams**3
-    return value
-
-
 def voronoi(network, geometry, **kwargs):
     r"""
     Calculate volume from the convex hull of the offset vertices making the throats
diff --git a/OpenPNM/Geometry/models/throat_area.py b/OpenPNM/Geometry/models/throat_area.py
index 401f731b6..21deb26d6 100644
--- a/OpenPNM/Geometry/models/throat_area.py
+++ b/OpenPNM/Geometry/models/throat_area.py
@@ -11,7 +11,18 @@ import OpenPNM.Utilities.vertexops as vo
 
 def cylinder(geometry, throat_diameter='throat.diameter', **kwargs):
     r"""
-    Calculate throat area for a cylindrical throat
+    Calculate throat cross-sectional area for a cylindrical throat
+
+    Parameters
+    ----------
+    geometry : OpenPNM Geometry Object
+        The Geometry object which this model is associated with. This controls
+        the length of the calculated array, and also provides access to other
+        necessary geometric properties.
+
+    throat_diameter : string
+        Dictionary key to the throat diameter values
+
     """
     diams = geometry[throat_diameter]
     value = _sp.constants.pi/4*(diams)**2
@@ -20,7 +31,18 @@ def cylinder(geometry, throat_diameter='throat.diameter', **kwargs):
 
 def cuboid(geometry, throat_diameter='throat.diameter', **kwargs):
     r"""
-    Calculate throat area for a cuboid throat
+    Calculate throat cross-sectional area for a cuboid throat
+
+    Parameters
+    ----------
+    geometry : OpenPNM Geometry Object
+        The Geometry object which this model is associated with. This controls
+        the length of the calculated array, and also provides access to other
+        necessary geometric properties.
+
+    throat_diameter : string
+        Dictionary key to the throat diameter values
+
     """
     diams = geometry[throat_diameter]
     value = (diams)**2
diff --git a/OpenPNM/Geometry/models/throat_diameter.py b/OpenPNM/Geometry/models/throat_diameter.py
index b159b10ee..92a9f8a8b 100644
--- a/OpenPNM/Geometry/models/throat_diameter.py
+++ b/OpenPNM/Geometry/models/throat_diameter.py
@@ -9,7 +9,6 @@ import numpy as np
 
 import OpenPNM.Utilities.transformations as tr
 from scipy.spatial import ConvexHull
-import time
 
 
 def cylinder(geometry, tsd_name, tsd_shape, tsd_loc, tsd_scale,
diff --git a/OpenPNM/Geometry/models/throat_misc.py b/OpenPNM/Geometry/models/throat_misc.py
index c890a7953..a632ddb24 100644
--- a/OpenPNM/Geometry/models/throat_misc.py
+++ b/OpenPNM/Geometry/models/throat_misc.py
@@ -14,7 +14,7 @@ def random(geometry, seed=None, num_range=[0, 1], **kwargs):
     """
     range_size = num_range[1] - num_range[0]
     range_min = num_range[0]
-    _sp.random.seed(seed)
+    _sp.random.seed(seed=seed)
     value = _sp.random.rand(geometry.num_throats(),)
     value = value*range_size + range_min
     return value
diff --git a/OpenPNM/Geometry/models/throat_seed.py b/OpenPNM/Geometry/models/throat_seed.py
index 0dda5733c..098842ea0 100644
--- a/OpenPNM/Geometry/models/throat_seed.py
+++ b/OpenPNM/Geometry/models/throat_seed.py
@@ -7,5 +7,54 @@ Submodule -- throat_seeds
 import scipy as _sp
 
 
-def na():
-    pass
+def random(geometry, seed=None, num_range=[0, 1], **kwargs):
+    r"""
+    Assign random number to throats, for use in statistical distributions that
+    return pore size
+
+    Parameters
+    ----------
+    seed : int
+        The starting seed value to send to Scipy's random number generator.
+        The default is None, which means different distribution is returned
+        each time the model is run.
+
+    num_range : list
+        A two element list indicating the low and high end of the returned
+        numbers.
+    """
+    range_size = num_range[1] - num_range[0]
+    range_min = num_range[0]
+    _sp.random.seed(seed)
+    value = _sp.random.rand(geometry.num_throats(),)
+    value = value*range_size + range_min
+    return value
+
+
+def neighbor(geometry, network, pore_prop='pore.seed', mode='min', **kwargs):
+    r"""
+    Adopt a value based on the values in the neighboring pores
+
+    Parameters
+    ----------
+    mode : string
+        Indicates how to select the values from the neighboring pores.  The
+        options are:
+
+        - min : (Default) Uses the minimum of the value found in the neighbors
+        - max : Uses the maximum of the values found in the neighbors
+        - mean : Uses an average of the neighbor values
+
+    pore_prop : string
+        The dictionary key containing the pore property to be used.
+    """
+    throats = network.throats(geometry.name)
+    P12 = network.find_connected_pores(throats)
+    pvalues = network[pore_prop][P12]
+    if mode == 'min':
+        value = _sp.amin(pvalues, axis=1)
+    if mode == 'max':
+        value = _sp.amax(pvalues, axis=1)
+    if mode == 'mean':
+        value = _sp.mean(pvalues, axis=1)
+    return value
diff --git a/OpenPNM/Geometry/models/throat_surface_area.py b/OpenPNM/Geometry/models/throat_surface_area.py
index 493329ec4..c9da35eff 100644
--- a/OpenPNM/Geometry/models/throat_surface_area.py
+++ b/OpenPNM/Geometry/models/throat_surface_area.py
@@ -33,8 +33,8 @@ def extrusion(geometry, throat_perimeter='throat.perimeter',
               throat_length='throat.length', **kwargs):
     r"""
     Calculate surface area from perimeter and length -
-    perimeter calculated when throat area is calculated so must be run in correct
-    order
+    perimeter calculated when throat area is calculated so must be run in
+    correct order
     """
     P = geometry[throat_perimeter]
     L = geometry[throat_length]
diff --git a/OpenPNM/Geometry/models/throat_vector.py b/OpenPNM/Geometry/models/throat_vector.py
index 8a80ef112..11020a943 100644
--- a/OpenPNM/Geometry/models/throat_vector.py
+++ b/OpenPNM/Geometry/models/throat_vector.py
@@ -16,7 +16,7 @@ def pore_to_pore(geometry, network, **kwargs):
     There is an important impicit assumption here: the positive direction is
     taken as the direction from the pore with the lower index to the higher.
     This corresponds to the pores in the 1st and 2nd columns of the
-    'conns' array as stored on the network.
+    'throat.conns' array as stored on the etwork.
     """
     throats = network.throats(geometry.name)
     pores = network.find_connected_pores(throats, flatten=False)
diff --git a/OpenPNM/Geometry/models/throat_volume.py b/OpenPNM/Geometry/models/throat_volume.py
index 18559c19e..1af1f9cf2 100644
--- a/OpenPNM/Geometry/models/throat_volume.py
+++ b/OpenPNM/Geometry/models/throat_volume.py
@@ -11,7 +11,8 @@ def cylinder(geometry, throat_length='throat.length',
              throat_diameter='throat.diameter', **kwargs):
     r"""
     Calculate throat diameter from seeds for a cylindrical throat
-    - note: this will need to account for volume taken up by spherical pore bodies
+    - note: this will need to account for volume taken up by spherical pore
+    bodies
     """
     leng = geometry[throat_length]
     diam = geometry[throat_diameter]
diff --git a/OpenPNM/Network/__Cubic__.py b/OpenPNM/Network/__Cubic__.py
index 8fb9e984b..37bb0248a 100644
--- a/OpenPNM/Network/__Cubic__.py
+++ b/OpenPNM/Network/__Cubic__.py
@@ -127,7 +127,7 @@ class Cubic(GenericNetwork):
         elif connectivity == 8:
             joints = corner_joints
         elif connectivity == 12:
-            joints = corner_joints
+            joints = edge_joints
         elif connectivity == 14:
             joints = face_joints + corner_joints
         elif connectivity == 18:
diff --git a/OpenPNM/Network/__GenericNetwork__.py b/OpenPNM/Network/__GenericNetwork__.py
index 7fb9bc710..ba2d1b9c2 100644
--- a/OpenPNM/Network/__GenericNetwork__.py
+++ b/OpenPNM/Network/__GenericNetwork__.py
@@ -258,7 +258,7 @@ class GenericNetwork(Core):
         if Ts.dtype == bool:
             Ts = self.toindices(Ts)
         if sp.size(Ts) == 0:
-            return sp.ndarray([0,2], dtype=int)
+            return sp.ndarray([0, 2], dtype=int)
         Ps = self['throat.conns'][Ts]
         if flatten:
             Ps = sp.unique(sp.hstack(Ps))
@@ -486,11 +486,11 @@ class GenericNetwork(Core):
         >>> pn.num_neighbors(pores=[0, 2], flatten=True)
         6
         """
-        pores = sp.array(pores,ndmin=1)
+        pores = sp.array(pores, ndmin=1)
         if pores.dtype == bool:
             pores = self.toindices(pores)
         if sp.size(pores) == 0:
-            return sp.array([],ndmin=1,dtype=int)
+            return sp.array([], ndmin=1, dtype=int)
 
         # Count number of neighbors
         if flatten:
@@ -584,7 +584,9 @@ class GenericNetwork(Core):
             temp = sp.array(conns[:, 0]*conns[:, 1], dtype=bool)
         else:
             raise Exception('Mask received was neither Nt nor Np long')
-        temp = self.create_adjacency_matrix(data=temp, sprsfmt='csr', dropzeros=True)
+        temp = self.create_adjacency_matrix(data=temp,
+                                            sprsfmt='csr',
+                                            dropzeros=True)
         clusters = sprs.csgraph.connected_components(csgraph=temp,
                                                      directed=False)[1]
         return clusters
@@ -618,40 +620,40 @@ class GenericNetwork(Core):
         Examples
         --------
         >>> import OpenPNM
-        >>> import scipy as sp
-        >>> sp.random.seed = 0  # Set seed value for random number generator
-        >>> pn = OpenPNM.Network.Cubic(shape=[25,25,1])
+        >>> pn = OpenPNM.Network.Cubic(shape=[25, 25, 1])
         >>> geom = OpenPNM.Geometry.GenericGeometry(network=pn,
-                                                    pores=pn.Ps,
-                                                    throats=pn.Ts)
+        ...                                         pores=pn.Ps,
+        ...                                         throats=pn.Ts)
         >>> geom['pore.seed'] = sp.rand(pn.Np)
         >>> geom['throat.seed'] = sp.rand(pn.Nt)
 
         Bond percolation is achieved by sending a list of invaded throats:
 
         >>> (p_bond,t_bond) = pn.find_clusters2(mask=geom['throat.seed'] < 0.3,
-                                                t_labels=True)
+        ...                                     t_labels=True)
 
         Site percolation is achieved by sending a list of invaded pores:
 
         >>> (p_site,t_site) = pn.find_clusters2(mask=geom['pore.seed'] < 0.3,
-                                      t_labels=True)
+        ...                                     t_labels=True)
 
-        To visualize the invasion pattern, use matplotlib's matshow method:
+        To visualize the invasion pattern, use matplotlib's matshow method
+        along with the Cubic Network's asarray method which converts list based
+        data to square arrays:
 
         .. code-block:: python
 
             import matplotlib.pyplot as plt
-            im_bond = pn.asarray(p_bond)[:,:,0]
-            im_site = pn.asarray(p_site)[:,:,0]
-            plt.subplot(1,2,1)
-            plt.imshow(im_site,interpolation='none')
-            plt.subplot(1,2,2)
-            plt.imshow(im_bond,interpolation='none')
+            im_bond = pn.asarray(p_bond)[:, :, 0]
+            im_site = pn.asarray(p_site)[:, :, 0]
+            plt.subplot(1, 2, 1)
+            plt.imshow(im_site, interpolation='none')
+            plt.subplot(1, 2, 2)
+            plt.imshow(im_bond, interpolation='none')
 
         """
         # Parse the input arguments
-        mask = sp.array(mask,ndmin=1)
+        mask = sp.array(mask, ndmin=1)
         if mask.dtype != bool:
             raise Exception('Mask must be a boolean array of Np or Nt length')
 
@@ -668,7 +670,7 @@ class GenericNetwork(Core):
         else:
             return p_clusters
 
-    def _site_percolation(self,pmask):
+    def _site_percolation(self, pmask):
         r"""
         """
         # Find throats that produce site percolation
@@ -691,14 +693,14 @@ class GenericNetwork(Core):
         p_clusters = (clusters + 1)*(pmask) - 1
         # Label invaded throats with their neighboring pore's label
         t_clusters = clusters[self['throat.conns']]
-        ind = (t_clusters[:,0] == t_clusters[:,1])
-        t_clusters = t_clusters[:,0]
+        ind = (t_clusters[:, 0] == t_clusters[:, 1])
+        t_clusters = t_clusters[:, 0]
         # Label non-invaded throats with -1
         t_clusters[~ind] = -1
 
         return (p_clusters, t_clusters)
 
-    def _bond_percolation(self,tmask):
+    def _bond_percolation(self, tmask):
         r"""
         """
         # Perform the clustering using scipy.csgraph
@@ -710,11 +712,11 @@ class GenericNetwork(Core):
 
         # Convert clusters to a more usable output:
         # Find pores attached to each invaded throats
-        Ps = self.find_connected_pores(throats=tmask,flatten=True)
+        Ps = self.find_connected_pores(throats=tmask, flatten=True)
         # Adjust cluster numbers such that non-invaded pores are labelled -1
         p_clusters = (clusters + 1)*(self.tomask(pores=Ps).astype(int)) - 1
         # Label invaded throats with their neighboring pore's label
-        t_clusters = clusters[self['throat.conns']][:,0]
+        t_clusters = clusters[self['throat.conns']][:, 0]
         # Label non-invaded throats with -1
         t_clusters[~tmask] = -1
 
@@ -777,7 +779,10 @@ class GenericNetwork(Core):
     stitch.__doc__ = topo.stitch.__doc__
 
     def connect_pores(self, pores1, pores2, labels=[]):
-        topo.connect_pores(network=self, pores1=pores1, pores2=pores2, labels=labels)
+        topo.connect_pores(network=self,
+                           pores1=pores1,
+                           pores2=pores2,
+                           labels=labels)
     connect_pores.__doc__ = topo.connect_pores.__doc__
 
     def check_network_health(self):
@@ -838,21 +843,26 @@ class GenericNetwork(Core):
         v = sp.array(self['throat.all'], dtype=int)
         Np = self.num_pores()
         adjmat = sprs.coo_matrix((v, (i, j)), [Np, Np])
-        temp = adjmat.tocsr()  # Convert to CSR to combine duplicates
-        temp = adjmat.tocoo()  # And back to COO
-        mergedTs = sp.where(temp.data > 1)
-        Ps12 = sp.vstack((temp.row[mergedTs], temp.col[mergedTs])).T
-        dupTs = []
-        for i in range(0, sp.shape(Ps12)[0]):
-            dupTs.append(self.find_connecting_throat(Ps12[i, 0], Ps12[i, 1]).tolist)
-        health['duplicate_throats'] = dupTs
+        temp = adjmat.tolil()  # Convert to lil to combine duplicates
+        # Compile lists of which specfic throats are duplicates
+        # Be VERY careful here, as throats are not in order
+        mergeTs = []
+        for i in range(0, self.Np):
+            if sp.any(sp.array(temp.data[i]) > 1):
+                ind = sp.where(sp.array(temp.data[i]) > 1)[0]
+                P = sp.array(temp.rows[i])[ind]
+                Ts = self.find_connecting_throat(P1=i, P2=P)[0]
+                mergeTs.append(Ts)
+        health['duplicate_throats'] = mergeTs
 
         # Check for bidirectional throats
         num_full = adjmat.sum()
         temp = sprs.triu(adjmat, k=1)
         num_upper = temp.sum()
         if num_full > num_upper:
-            health['bidirectional_throats'] = str(num_full-num_upper) + ' detected!'
+            biTs = sp.where(self['throat.conns'][:, 0] >
+                            self['throat.conns'][:, 1])[0]
+            health['bidirectional_throats'] = biTs.tolist()
 
         return health
 
@@ -954,9 +964,9 @@ class GenericNetwork(Core):
             f1 = self['pore.coords'][face_1]
             f2 = self['pore.coords'][face_2]
             distavg = [0, 0, 0]
-            distavg[0] = sp.absolute(sp.average(f1[:, 0]) - sp.average(f2[:, 0]))
-            distavg[1] = sp.absolute(sp.average(f1[:, 1]) - sp.average(f2[:, 1]))
-            distavg[2] = sp.absolute(sp.average(f1[:, 2]) - sp.average(f2[:, 2]))
+            distavg[0] = sp.absolute(sp.average(f1[:, 0])-sp.average(f2[:, 0]))
+            distavg[1] = sp.absolute(sp.average(f1[:, 1])-sp.average(f2[:, 1]))
+            distavg[2] = sp.absolute(sp.average(f1[:, 2])-sp.average(f2[:, 2]))
             L = max(distavg)
         return L
 
diff --git a/OpenPNM/Phases/models/density.py b/OpenPNM/Phases/models/density.py
index ebd9ad6f8..0c6baf6b6 100644
--- a/OpenPNM/Phases/models/density.py
+++ b/OpenPNM/Phases/models/density.py
@@ -7,55 +7,83 @@ Submodule -- density
 import scipy as _sp
 
 
-def standard(phase, **kwargs):
+def standard(phase,
+             pore_MW='pore.molecular_weight',
+             pore_molar_density='pore.molar_density',
+             **kwargs):
     r"""
     Calculates the mass density from the molecular weight and molar density
+
+    Parameters
+    ----------
+    pore_MW : string
+        The dictionary key containing the molecular weight values
+
+    pore_molar_density : string
+        The dictionary key containing the molar density values
+
     """
-    MW = phase['pore.molecular_weight']
-    rho = phase['pore.molar_density']
+    MW = phase[pore_MW]
+    rho = phase[pore_molar_density]
     value = rho*MW
     return value
 
 
-def ideal_gas(phase, **kwargs):
+def ideal_gas(phase,
+              pore_pressure='pore.pressure',
+              pore_temperature='pore.temperature',
+              pore_MW='pore.molecular_weight',
+              **kwargs):
     r"""
     Uses ideal gas law to calculate the mass density of an ideal gas
 
     Parameters
     ----------
-    P, T, MW: float, array_like
-        P pressure of the gas in [Pa]
-        T temperature of the gas in [K]
-        MW molecular weight of the gas in [kg/kmole]
+    pore_pressure : string
+        The dictionary key containing the pressure values (Pa)
 
-    Returns
-    -------
-    rho, the density in [mol/m3]
+    pore_temperature : string
+        The dictionary key containing the temperature values (K)
+
+    pore_molecular_weight : string
+        The dictionary key containing the molecular weight values (kg/mol)
 
     """
 
-    P = phase['pore.pressure']
-    T = phase['pore.temperature']
-    MW = phase['pore.molecular_weight']
+    P = phase[pore_pressure]
+    T = phase[pore_temperature]
+    MW = phase[pore_MW]
     R = 8.31447
     value = P/(R*T)*MW
     return value
 
 
-def water(phase, **kwargs):
+def water(phase,
+          pore_temperature='pore.temperature',
+          pore_salinity='pore.salinity',
+          **kwargs):
     r"""
     Calculates density of pure water or seawater at atmospheric pressure
     using Eq. (8) given by Sharqawy et. al [1]_. Values at temperature higher
-    than the normal boiling temperature are calculated at the saturation pressure.
+    than the normal boiling temperature are calculated at the saturation
+    pressure.
 
     Parameters
     ----------
-    T, S: strings
-        Property names where phase temperature and salinity are located.
+    phase : OpenPNM Phase Object
+
+    pore_temperature : string
+        The dictionary key containing the temperature values.  Temperature must
+        be in Kelvin for this emperical equation to work
+
+    pore_salinity : string
+        The dictionary key containing the salinity values.  Salinity must be
+        expressed in g of salt per kg of solution (ppt).
+
 
     Returns
     -------
-    rho_sw, the density of water/seawater in [kg/m3]
+    The density of water/seawater in [kg/m3]
 
     Notes
     -----
@@ -66,8 +94,8 @@ def water(phase, **kwargs):
 
     References
     ----------
-    [1] Sharqawy M. H., Lienhard J. H., and Zubair, S. M., Desalination and Water
-        Treatment, 2010.
+    [1] Sharqawy M. H., Lienhard J. H., and Zubair, S. M., Desalination and
+        Water Treatment, 2010.
 
     """
     T = phase['pore.temperature']
diff --git a/OpenPNM/Phases/models/diffusivity.py b/OpenPNM/Phases/models/diffusivity.py
index fc1105803..2bfa10b11 100644
--- a/OpenPNM/Phases/models/diffusivity.py
+++ b/OpenPNM/Phases/models/diffusivity.py
@@ -4,10 +4,13 @@ Submodule -- diffusivity
 ===============================================================================
 
 """
-import scipy as sp
+import scipy as _sp
 
 
-def fuller(phase, MA, MB, vA, vB, **kwargs):
+def fuller(phase, MA, MB, vA, vB,
+           pore_temperature='pore.temperature',
+           pore_pressure='pore.pressure',
+           **kwargs):
     r"""
     Uses Fuller model to estimate diffusion coefficient for gases from first
     principles at conditions of interest
@@ -16,16 +19,25 @@ def fuller(phase, MA, MB, vA, vB, **kwargs):
     ----------
     MA : float, array_like
         Molecular weight of component A [kg/mol]
+
     MB : float, array_like
         Molecular weight of component B [kg/mol]
+
     vA:  float, array_like
         Sum of atomic diffusion volumes for component A
+
     vB:  float, array_like
         Sum of atomic diffusion volumes for component B
+
+    pore_pressure : string
+        The dictionary key containing the pressure values in Pascals (Pa)
+
+    pore_temperature : string
+        The dictionary key containing the temperature values in Kelvin (K)
     """
 
-    T = phase['pore.temperature']
-    P = phase['pore.pressure']
+    T = phase[pore_temperature]
+    P = phase[pore_pressure]
     MAB = 2*(1.0/MA+1.0/MB)**(-1)
     MAB = MAB*1e3
     P = P*1e-5
@@ -33,10 +45,13 @@ def fuller(phase, MA, MB, vA, vB, **kwargs):
     return value
 
 
-def fuller_scaling(phase, DABo, To, Po, **kwargs):
+def fuller_scaling(phase, DABo, To, Po,
+                   pore_temperature='pore.temperature',
+                   pore_pressure='pore.pressure',
+                   **kwargs):
     r"""
-    Uses Fuller model to adjust a diffusion coefficient for gases from reference
-    conditions to conditions of interest
+    Uses Fuller model to adjust a diffusion coefficient for gases from
+    reference conditions to conditions of interest
 
     Parameters
     ----------
@@ -47,42 +62,60 @@ def fuller_scaling(phase, DABo, To, Po, **kwargs):
 
     Po, To : float, array_like
         Pressure & temperature at reference conditions, respectively
+
+    pore_pressure : string
+        The dictionary key containing the pressure values in Pascals (Pa)
+
+    pore_temperature : string
+        The dictionary key containing the temperature values in Kelvin (K)
     """
-    Ti = phase['pore.temperature']
-    Pi = phase['pore.pressure']
+    Ti = phase[pore_temperature]
+    Pi = phase[pore_pressure]
     value = DABo*(Ti/To)**1.75*(Po/Pi)
     return value
 
 
-def tyn_calus(phase, VA, VB, sigma_A, sigma_B, **kwargs):
+def tyn_calus(phase, VA, VB, sigma_A, sigma_B,
+              pore_temperature='pore.temperature',
+              pore_viscosity='pore.viscosity',
+              **kwargs):
     r"""
     Uses Tyn_Calus model to estimate diffusion coefficient in a dilute liquid
     solution of A in B from first principles at conditions of interest
 
     Parameters
     ----------
-    T :  float, array_like
-        Temperature of interest (K)
-    viscosity :  float, array_like
-        Viscosity of solvent (Pa.s)
     VA : float, array_like
         Molar volume of component A at boiling temperature (m3/mol)
+
     VB : float, array_like
         Molar volume of component B at boiling temperature (m3/mol)
+
     sigmaA:  float, array_like
         Surface tension of component A at boiling temperature (N/m)
+
     sigmaB:  float, array_like
         Surface tension of component B at boiling temperature (N/m)
 
+    pore_pressure : string
+        The dictionary key containing the pressure values in Pascals (Pa)
+
+    pore_temperature : string
+        The dictionary key containing the temperature values in Kelvin (K)
+
     """
     T = phase['pore.temperature']
     mu = phase['pore.viscosity']
-    value = \
-        8.93e-8*(VB*1e6)**0.267/(VA*1e6)**0.433*T*(sigma_B/sigma_A)**0.15/(mu*1e3)
+    A = 8.93e-8*(VB*1e6)**0.267/(VA*1e6)**0.433*T
+    B = (sigma_B/sigma_A)**0.15/(mu*1e3)
+    value = A*B
     return value
 
 
-def tyn_calus_Scaling(phase, DABo, To, mu_o, **kwargs):
+def tyn_calus_scaling(phase, DABo, To, mu_o,
+                      pore_temperature='pore.temperature',
+                      pore_viscosity='pore.viscosity',
+                      **kwargs):
     r"""
     Uses Tyn_Calus model to adjust a diffusion coeffciient for liquids from
     reference conditions to conditions of interest
@@ -96,8 +129,14 @@ def tyn_calus_Scaling(phase, DABo, To, mu_o, **kwargs):
 
     mu_o, To : float, array_like
         Viscosity & temperature at reference conditions, respectively
+
+    pore_pressure : string
+        The dictionary key containing the pressure values in Pascals (Pa)
+
+    pore_temperature : string
+        The dictionary key containing the temperature values in Kelvin (K)
     """
-    Ti = phase['pore.temperature']
-    mu_i = phase['pore.viscosity']
+    Ti = phase[pore_temperature]
+    mu_i = phase[pore_viscosity]
     value = DABo*(Ti/To)*(mu_o/mu_i)
     return value
diff --git a/OpenPNM/Phases/models/electrical_conductivity.py b/OpenPNM/Phases/models/electrical_conductivity.py
index 251a52d67..9b8fd39be 100644
--- a/OpenPNM/Phases/models/electrical_conductivity.py
+++ b/OpenPNM/Phases/models/electrical_conductivity.py
@@ -1,11 +1,45 @@
 r"""
 ===============================================================================
-Submodule -- electrical_conductance
+Submodule -- electrical_conductivity
 ===============================================================================
 
 """
-import scipy as sp
+import scipy as _sp
 
 
-def na(**kwargs):
-    print('Not Implemented Yet')
+def percolating_continua(phase,
+                         phi_crit,
+                         tau,
+                         volume_fraction='pore.volume_fraction',
+                         bulk_property='pore.intrinsic_conductivity',
+                         **kwargs):
+    r'''
+    Calculates the effective property of a continua using percolation theory
+
+    Parameters
+    ----------
+    volume_fraction : string
+        The dictionary key in the Phase object containing the volume fraction
+        of the conducting component
+    bulk_property : string
+        The dictionary key in the Phase object containing the intrinsic
+        property of the conducting component
+    phi_crit : float
+        The volume fraction below which percolation does NOT occur
+    tau : float
+        The exponent of the percolation relationship
+
+    Notes
+    -----
+    This model uses the following standard percolation relationship:
+
+    .. math::
+
+        \sigma_{effective}=\sigma_{bulk}(\phi - \phi_{critical})^\lambda
+
+    '''
+    sigma = phase[bulk_property]
+    phi = phase[volume_fraction]
+    diff_phi = _sp.clip(phi - phi_crit, a_min=0, a_max=_sp.inf)
+    sigma_eff = sigma*(diff_phi)**tau
+    return sigma_eff
diff --git a/OpenPNM/Phases/models/misc.py b/OpenPNM/Phases/models/misc.py
index 9d63c0185..9887e7f5a 100644
--- a/OpenPNM/Phases/models/misc.py
+++ b/OpenPNM/Phases/models/misc.py
@@ -21,23 +21,23 @@ def random(phase, seed=None, **kwargs):
     r"""
     Assigns specified constant value
     """
-    _sp.random.seed(seed)
-    value = _sp.random.rand(_sp.shape(phase.pores())[0])
+    _sp.random.seed(seed=seed)
+    value = _sp.random.rand(phase.Np)
     return value
 
 
-def linear(phase, m, b, poreprop='pore.temperature', **kwargs):
+def linear(phase, m, b, poreprop, **kwargs):
     r"""
     Calculates a property as a linear function of a given property
 
     Parameters
     ----------
-    m, b: floats
+    m, b : floats
         Slope and intercept of the linear corelation
 
     poreprop : string
-        The property name of the independent variable or phase property.  The
-        default is 'pore.temperature'.
+        The dictionary key containing the independent variable or phase
+        property to be used in the correlation.
 
     """
     T = phase[poreprop]
@@ -45,20 +45,20 @@ def linear(phase, m, b, poreprop='pore.temperature', **kwargs):
     return value
 
 
-def polynomial(phase, a, poreprop='pore.temperature', **kwargs):
+def polynomial(phase, a, poreprop, **kwargs):
     r"""
     Calculates a property as a polynomial function of a given property
 
     Parameters
     ----------
-    a: list of floats
+    a : array_like
         A list containing the polynomial coefficients, where element 0 in the
         list corresponds to a0 and so on.  Note that no entries can be skipped
         so 0 coefficients must be sent as 0.
 
     poreprop : string
-        The property name of the independent variable or phase property.  The
-        default is 'pore.temperature'.
+        The dictionary key containing the independent variable or phase
+        property to be used in the correlation.
 
     """
     x = phase[poreprop]
@@ -68,16 +68,19 @@ def polynomial(phase, a, poreprop='pore.temperature', **kwargs):
     return value
 
 
-def ideal_mixture(phase, propname, pore_prop=None,
-                  composition='pore.mole_fraction', **kwargs):
+def ideal_mixture(phase, poreprop,
+                  composition='pore.mole_fraction',
+                  **kwargs):
     r"""
     Calcualtes a given mixture property as the composition weighted average
     of the pure compononent properties
 
     Parameters
     ----------
-    pore_prop : string
-        The name of the property on the pure component
+    poreprop : string
+        The dictionary key containing the independent variable or phase
+        property to be used in the correlation.
+
     composition : string, optional (default is 'pore.mole_fraction')
         The name of the pore property where the composition information
         is stored on each pure component
@@ -104,11 +107,9 @@ def ideal_mixture(phase, propname, pore_prop=None,
 
 
     """
-    if pore_prop is None:
-        pore_prop = propname
     value = _sp.zeros((phase.Np,))
     for comp in phase._phases:
-        value = value + comp[pore_prop]*comp[composition]
+        value = value + comp[poreprop]*comp[composition]
     return value
 
 
diff --git a/OpenPNM/Phases/models/molar_density.py b/OpenPNM/Phases/models/molar_density.py
index e0b827465..ca4285358 100644
--- a/OpenPNM/Phases/models/molar_density.py
+++ b/OpenPNM/Phases/models/molar_density.py
@@ -7,70 +7,119 @@ Submodule -- molar_density
 import scipy as sp
 
 
-def standard(phase, **kwargs):
+def standard(phase,
+             pore_MW='pore.molecular_weight',
+             pore_density='pore.density',
+             **kwargs):
     r"""
-    Calculates the molar density from the molecular with and mass density
+    Calculates the molar density from the molecular weight and mass density
+
+    Parameters
+    ----------
+    phase : OpenPNM Phase Object
+        The Phase object for which these values are being calculated.  This
+        controls the length of the calculated array, and also provides access
+        to other necessary thermofluid properties.
+
+    pore_MW : string
+        The dictionary key containing the molecular weight in kg/mol
+
+    pore_temperature : string
+        The dictionary key containing the density in kg/m3
     """
-    MW = phase['pore.molecular_weight']
-    rho = phase['pore.density']
+    MW = phase[pore_MW]
+    rho = phase[pore_density]
     value = rho/MW
     return value
 
 
-def ideal_gas(phase, **kwargs):
+def ideal_gas(phase,
+              pore_pressure='pore.pressure',
+              pore_temperature='pore.temperature',
+              **kwargs):
     r"""
     Uses ideal gas law to calculate the molar density of an ideal gas
 
     Parameters
     ----------
-    P, T, MW: float, array_like
-        P pressure of the gas in [Pa]
-        T temperature of the gas in [K]
-        MW molecular weight of the gas in [kg/kmole]
+    phase : OpenPNM Phase Object
+        The Phase object for which these values are being calculated.  This
+        controls the length of the calculated array, and also provides access
+        to other necessary thermofluid properties.
+
+    pore_temperature : string
+        The dictionary key containing the density in kg/m3
+
+    pore_pressure : string
+        The dictionary key containing the pressure values in Pascals (Pa)
 
     Returns
     -------
     rho, the density in [mol/m3]
 
-    """
+    Notes
+    -----
+    This method uses the SI value for the ideal gas constant, hence the need to
+    provide the temperature and pressure in SI.  In general, OpenPNM use SI
+    throughout for consistency.
 
-    P = phase['pore.pressure']
-    T = phase['pore.temperature']
+    """
     R = 8.31447
+    P = phase[pore_pressure]
+    T = phase[pore_temperature]
     value = P/(R*T)
     return value
 
 
-def vanderwaals(phase, P, T, Pc, Tc, MW, **kwargs):
+def vanderwaals(phase,
+                pore_P='pore.pressure',
+                pore_T='pore.temperature',
+                pore_Pc='pore.critical_pressure',
+                pore_Tc='pore.critical_temperature',
+                **kwargs):
     r"""
     Uses Van der Waals equation of state to calculate the density of a real gas
 
     Parameters
     ----------
-    P, T, Pc, Tc, MW: float, array_like
-        P pressure of the gas in [Pa]
-        T temperature of the gas in [K]
-        Pc critical pressure of the gas in [Pa]
-        T critical temperature of the gas in [K]
-        MW molecular weight of the gas in [kg/mol]
+    pore_P : string
+        The dictionary key containing the pressure values in Pascals (Pa)
+
+    pore_T : string
+        The dictionary key containing the temperature values in Kelvin (K)
+
+    pore_Pc : string
+        The dictionary key containing the critical pressure values in Pascals
+        (Pa)
+
+    pore_Tc : string
+        The dictionary key containing the critical temperature values in Kelvin
+        (K)
 
     Returns
     -------
     rho, the density in [mol/m3]
 
+    Notes
+    -----
+    This equation and its constant coefficients are taken [1]_ which uses the
+    cgs units system. All input parameters are expected in SI, then converted
+    in the method.
+
     """
 
-    P = phase['pore.pressure']
-    T = phase['pore.temperature']
-    Pc = phase['pore.criticalpressure']
-    Tc = phase['pore.criticaltemperature']
-    R = 8.314
+    P = phase[pore_P]/100000
+    T = phase[pore_T]
+    Pc = phase[pore_Pc]/100000  # convert to bars
+    Tc = phase[pore_Tc]
+    R = 83.1447
     a = 27*(R**2)*(Tc**2)/(64*Pc)
     b = R*Tc/(8*Pc)
-    a0 = 1
     a1 = -1/b
     a2 = (R*T+b*P)/(a*b)
     a3 = -P/(a*b)
-    density = sp.roots([a0, a1, a2, a3])
-    value = sp.real(density[2])
+    a0 = sp.ones(sp.shape(a1))
+    coeffs = sp.vstack((a0, a1, a2, a3)).T
+    density = sp.array([sp.roots(C) for C in coeffs])
+    value = sp.real(density[:, 2])*1e6  # Convert it to mol/m3
     return value
diff --git a/OpenPNM/Phases/models/molar_mass.py b/OpenPNM/Phases/models/molar_mass.py
index dd77ba529..027f25de5 100644
--- a/OpenPNM/Phases/models/molar_mass.py
+++ b/OpenPNM/Phases/models/molar_mass.py
@@ -7,8 +7,10 @@ Submodule -- molar_mass
 import scipy as sp
 
 
-def mixture(phase, molar_mass='pore.molar_mass',
-            mole_frac='pore.mole_fraction', **kwargs):
+def mixture(phase,
+            molar_mass='pore.molar_mass',
+            mole_frac='pore.mole_fraction',
+            **kwargs):
     r"""
     Calculates the average molecular weight of a mixture using mole fraction
     weighting
@@ -17,7 +19,7 @@ def mixture(phase, molar_mass='pore.molar_mass',
     ----------
     phase : OpenPNM Phase Object
         The phase for which the molar mass is to be calculated.  This phase
-        must have sub-phases
+        must have sub-phases.
 
     molar_mass : string, optional (default = 'pore.molar_mass')
         The property name for the molar masses of each sub-phase
diff --git a/OpenPNM/Phases/models/surface_tension.py b/OpenPNM/Phases/models/surface_tension.py
index 84ed73514..53797d7d2 100644
--- a/OpenPNM/Phases/models/surface_tension.py
+++ b/OpenPNM/Phases/models/surface_tension.py
@@ -8,16 +8,25 @@ Submodule -- surface_tension
 import scipy as sp
 
 
-def water(phase, **kwargs):
+def water(phase,
+          pore_T='pore.temperature',
+          pore_salinity='pore.salinity',
+          **kwargs):
     r"""
-    Calculates surface tension of pure water or seawater at atmospheric pressure
-    using Eq. (28) given by Sharqawy et. al [1]_. Values at temperature higher
-    than the normal boiling temperature are calculated at the saturation pressure.
+    Calculates surface tension of pure water or seawater at atmospheric
+    pressure using Eq. (28) given by Sharqawy et. al [1]_. Values at
+    temperature higher than the normal boiling temperature are calculated at
+    the saturation pressure.
 
     Parameters
     ----------
-    T, S: strings
-        Property names where phase temperature and salinity are located.
+    pore_temperature : string
+        The dictionary key containing the temperature values.  Temperature must
+        be in Kelvin for this emperical equation to work
+
+    pore_salinity : string
+        The dictionary key containing the salinity values.  Salinity must be
+        expressed in g of salt per kg of solution (ppt).
 
     Returns
     -------
@@ -32,8 +41,8 @@ def water(phase, **kwargs):
 
     References
     ----------
-    [1] Sharqawy M. H., Lienhard J. H., and Zubair, S. M., Desalination and Water
-    Treatment, 2010.
+    [1] Sharqawy M. H., Lienhard J. H., and Zubair, S. M., Desalination and
+    Water Treatment, 2010.
 
     """
     T = phase['pore.temperature']
@@ -51,48 +60,97 @@ def water(phase, **kwargs):
     return value
 
 
-def eotvos(phase, k, pore_molar_density='pore.molar_density', **kwargs):
+def eotvos(phase, k,
+           pore_T='pore.temperature',
+           pore_Tc='pore.critical_temperature',
+           pore_molar_density='pore.molar_density', **kwargs):
     r"""
-    Documentation for this method is being updated, we are sorry for the
-    inconvenience.
+    Missing description
+
+    Parameters
+    ----------
+    k : float
+        Constant parameter specific to fluid
+
+    pore_T : string
+        The dictionary key containing the temperature values (K)
+
+    pore_Tc : string
+        The dictionary key containing the critical temperature values (K)
+
+    pore_molar_density : string
+        The dictionary key containing the molar density values (K)
+
+
+    TODO: Needs description, and improve definition of k
+
     """
-    Tc = phase['pore.Tc']
-    T = phase['pore.temperature']
+    Tc = phase[pore_Tc]
+    T = phase[pore_T]
     Vm = 1/phase[pore_molar_density]
     value = k*(Tc-T)/(Vm**(2/3))
     return value
 
 
-def guggenheim_katayama(phase, K2, n, **kwargs):
+def guggenheim_katayama(phase, K2, n,
+                        pore_T='pore.temperature',
+                        pore_Tc='pore.critical_temperature',
+                        pore_Pc='pore.critical_pressure',
+                        **kwargs):
     r"""
-    Documentation for this method is being updated, we are sorry for the
-    inconvenience.
+    Missing description
+
+    Parameters
+    ----------
+    K2 : scalar
+        Fluid specific constant
+
+    n : scalar
+        Fluid specific constant
+
+    pore_T : string
+        The dictionary key containing the temperature values (K)
+
+    pore_Tc : string
+        The dictionary key containing the critical temperature values (K)
+
+    pore_Pc : string
+        The dictionary key containing the critical pressure values (K)
+
+    TODO: Needs description
     """
-    T = phase['pore.temperature']
-    Pc = phase['pore.Pc']
-    Tc = phase['pore.Tc']
+    T = phase[pore_T]
+    Pc = phase[pore_Pc]
+    Tc = phase[pore_Tc]
     sigma_o = K2*Tc**(1/3)*Pc**(2/3)
     value = sigma_o*(1-T/Tc)**n
     return value
 
 
-def brock_bird_scaling(phase, sigma_o, To, **params):
+def brock_bird_scaling(phase, sigma_o, To,
+                       pore_T='pore.temperature',
+                       pore_Tc='pore.critical_temperature',
+                       **params):
     r"""
     Uses Brock_Bird model to adjust surface tension from it's value at a given
     reference temperature to temperature of interest
 
     Parameters
     ----------
-    phase : OpenPNM Phase Object
+    To : float
+        Reference temperature (K)
 
     sigma_o : float
         Surface tension at reference temperature (N/m)
 
-    To : float
-        Temperature at reference conditions (K)
+    pore_T : string
+        The dictionary key containing the temperature values (K)
+
+    pore_Tc : string
+        The dictionary key containing the critical temperature values (K)
     """
-    Tc = phase['pore.Tc']
-    Ti = phase['pore.temperature']
+    Tc = phase[pore_Tc]
+    Ti = phase[pore_T]
     Tro = To/Tc
     Tri = Ti/Tc
     value = sigma_o*(1-Tri)**(11/9)/(1-Tro)**(11/9)
diff --git a/OpenPNM/Phases/models/thermal_conductivity.py b/OpenPNM/Phases/models/thermal_conductivity.py
index fb1bc8dbd..dbf75a841 100644
--- a/OpenPNM/Phases/models/thermal_conductivity.py
+++ b/OpenPNM/Phases/models/thermal_conductivity.py
@@ -8,16 +8,27 @@ Submodule -- thermal_conductance
 import scipy as sp
 
 
-def water(phase, **kwargs):
+def water(phase,
+          pore_T='pore.temperature',
+          pore_salinity='pore.salinity',
+          **kwargs):
     r"""
-    Calculates thermal conductivity of pure water or seawater at atmospheric pressure
-    using the correlation given in [1]_. Values at temperature higher
-    than the normal boiling temperature are calculated at the saturation pressure.
+    Calculates thermal conductivity of pure water or seawater at atmospheric
+    pressure using the correlation given in [1]_. Values at temperature higher
+    than the normal boiling temperature are calculated at the saturation
+    pressure.
 
     Parameters
     ----------
-    T, S: strings
-        Property names where phase temperature and salinity are located.
+    phase : OpenPNM Phase Object
+
+    pore_temperature : string
+        The dictionary key containing the temperature values.  Temperature must
+        be in Kelvin for this emperical equation to work
+
+    pore_salinity : string
+        The dictionary key containing the salinity values.  Salinity must be
+        expressed in g of salt per kg of solution (ppt).
 
     Returns
     -------
@@ -35,9 +46,9 @@ def water(phase, **kwargs):
     [1] D. T. Jamieson, and J. S. Tudhope, Desalination, 8, 393-401, 1970.
 
     """
-    T = phase['pore.temperature']
+    T = phase[pore_T]
     try:
-        S = phase['pore.salinity']
+        S = phase[pore_salinity]
     except:
         S = 0
     T68 = 1.00024*T  # convert from T_90 to T_68
@@ -49,25 +60,48 @@ def water(phase, **kwargs):
     return value
 
 
-def chung(phase, Cv, MW, acentric, pore_viscosity='pore.viscosity', **kwargs):
+def chung(phase,
+          pore_Cv='pore.heat_capacity',
+          pore_acentric='pore.acentric_factor',
+          pore_MW='pore.molecular_weight',
+          pore_viscosity='pore.viscosity',
+          pore_T='pore.temperature',
+          pore_Tc='pore.critical_temperature',
+          **kwargs):
     r"""
     Uses Chung et al. model to estimate thermal conductivity for gases with
     low pressure(<10 bar) from first principles at conditions of interest
 
     Parameters
     ----------
-    Cv :  float, array_like
-        Heat capacity at constant volume (J/(mol.K))
-    MW : float, array_like
-        Molecular weight of the component (kg/mol)
-    acentric : float, array_like
-        Acentric factor of the component
+    pore_acentric : string
+        Dictionary key containing the acentric factor of the component
+
+    pore_Cv :  string
+        Dictionary key containing the heat capacity at constant volume
+        (J/(mol.K))
+
+    pore_MW : string
+        Dictionary key containing the molecular weight of the component
+        (kg/mol)
+
+    pore_viscosity : string
+        The dictionary key containing the viscosity values (Pa.s)
+
+    pore_T : string
+        The dictionary key containing the temperature values (K)
+
+    pore_Tc: string
+        The dictionary key containing the critical temperature values (K)
 
     """
+    Cv = phase[pore_Cv]
+    acentric = phase[pore_acentric]
+    MW = phase[pore_MW]
     R = 8.314
-    T = phase['pore.temperature']
+    T = phase[pore_T]
     mu = phase[pore_viscosity]
-    Tc = phase['pore.Tc']
+    Tc = phase[pore_Tc]
     Tr = T/Tc
     z = 2.0 + 10.5*Tr**2
     beta = 0.7862 - 0.7109*acentric + 1.3168*acentric**2
@@ -78,22 +112,36 @@ def chung(phase, Cv, MW, acentric, pore_viscosity='pore.viscosity', **kwargs):
     return value
 
 
-def sato(phase, Tb, MW, **params):
+def sato(phase,
+         pore_MW='pore.molecular_weight',
+         pore_Tb='pore.boiling_point',
+         pore_T='pore.temperature',
+         pore_Tc='pore.critical_temperature',
+         **params):
     r"""
     Uses Sato et al. model to estimate thermal conductivity for pure liquids
     from first principles at conditions of interest
 
     Parameters
     ----------
-    Tb :  float, array_like
-        Boiling temperature of the component (K)
-    MW : float, array_like
-        Molecular weight of the component (kg/mol)
+    pore_Tb :  string
+        Dictionary key containing the toiling temperature of the component (K)
+
+    pore_MW : string
+        Dictionary key containing the molecular weight of the component
+        (kg/mol)
+
+    pore_T : string
+        The dictionary key containing the temperature values (K)
+
+    pore_Tc: string
+        The dictionary key containing the critical temperature values (K)
 
     """
-    T = phase['pore.temperature']
-    Tc = phase['pore.Tc']
-    Tbr = Tb/Tc
+    T = phase[pore_T]
+    Tc = phase[pore_Tc]
+    MW = phase[pore_MW]
+    Tbr = phase[pore_Tb]/Tc
     Tr = T/Tc
     value = (1.11/((MW*1e3)**0.5))*(3+20*(1-Tr)**(2/3))/(3+20*(1-Tbr)**(2/3))
     return value
diff --git a/OpenPNM/Phases/models/vapor_pressure.py b/OpenPNM/Phases/models/vapor_pressure.py
index 1e900f9d7..3198fbe04 100644
--- a/OpenPNM/Phases/models/vapor_pressure.py
+++ b/OpenPNM/Phases/models/vapor_pressure.py
@@ -10,40 +10,51 @@ Methods for predicing the vapor pressure of pure species
 import scipy as sp
 
 
-def antoine(phase, A, B, C, **kwargs):
+def antoine(phase, A, B, C, pore_temperature='pore.temperature', **kwargs):
     r"""
     Uses Antoine equation [1]_ to estimate vapor pressure of a pure component
 
     Parameters
     ----------
-    A, B, C :  float, array_like
-            Antoine vapor pressure coefficients for pure compounds. Note that
-            these coefficients should be converted such that the Temperature
-            is in [K] and the vapor pressure is in [Pa].
+    A, B, C :  scalars
+        Antoine vapor pressure coefficients for pure compounds. Since virtually
+        all Antoine coefficients are reported for units of mmHg and C for
+        historical reaons, this method assumes these A, B and C values are for
+        mmHg and C, but converts all properties internally to returrn Pascals.
+
+    pore_temperature : string
+        The dictionary key containing the phase temperature values in Kelvin
+        [K]
 
     [1] Antoine, C. (1888), Vapor Pressure: a new relationship between pressure
         and temperature, Comptes Rendus des Séances de l'Académie des Sciences
         (in French) 107: 681–684, 778–780, 836–837
 
     """
-    T = phase['pore.temperature']
-    value = (10**(A-B/(C+T)))
+    T = phase[pore_temperature] - 273.15
+    value = (10**(A-B/(C+T)))/760*101325
     return value
 
 
-def water(phase, **kwargs):
+def water(phase,
+          pore_temperature='pore.temperature',
+          pore_salinity='pore.salinity',
+          **kwargs):
     r"""
     Calculates vapor pressure of pure water or seawater given by [1]_ based on
     Raoult's law. The pure water vapor pressure is given by [2]_
 
     Parameters
     ----------
-    T, S: strings
-        Property names where phase temperature and salinity are located.
+    pore_temperature : strings
+        The dictionary key containing the phase temperature values
+
+    pore_salinity : strings
+        The dictionary key containing the phase salinity values
 
     Returns
     -------
-    Pv_sw, the vapor pressure of water/seawater in [Pa]
+    The vapor pressure of water/seawater in [Pa]
 
     Notes
     -----
@@ -54,14 +65,14 @@ def water(phase, **kwargs):
 
     References
     ----------
-    [1] Sharqawy M. H., Lienhard J. H., and Zubair, S. M., Desalination and Water
-        Treatment, 2010.
+    [1] Sharqawy M. H., Lienhard J. H., and Zubair, S. M., Desalination and
+    Water Treatment, 2010.
     [2] ASHRAE handbook: Fundamentals, ASHRAE; 2005.
 
     """
-    T = phase['pore.temperature']
+    T = phase[pore_temperature]
     try:
-        S = phase['pore.salinity']
+        S = phase[pore_salinity]
     except:
         S = 0
     a1 = -5.8002206E+03
diff --git a/OpenPNM/Phases/models/viscosity.py b/OpenPNM/Phases/models/viscosity.py
index 06dea4a5e..d8fe9ea0e 100644
--- a/OpenPNM/Phases/models/viscosity.py
+++ b/OpenPNM/Phases/models/viscosity.py
@@ -10,16 +10,27 @@ Models for predicting phase viscosity
 import scipy as sp
 
 
-def water(phase, **kwargs):
+def water(phase,
+          pore_T='pore.temperature',
+          pore_salinity='pore.salinity',
+          **kwargs):
     r"""
     Calculates viscosity of pure water or seawater at atmospheric pressure
     using Eq. (22) given by Sharqawy et. al [1]_. Values at temperature higher
-    than the normal boiling temperature are calculated at the saturation pressure.
+    than the normal boiling temperature are calculated at the saturation
+    pressure.
 
     Parameters
     ----------
-    T, S: strings
-        Property names where phase temperature and salinity are located.
+    phase : OpenPNM Phase Object
+
+    pore_temperature : string
+        The dictionary key containing the temperature values.  Temperature must
+        be in Kelvin for this emperical equation to work
+
+    pore_salinity : string
+        The dictionary key containing the salinity values.  Salinity must be
+        expressed in g of salt per kg of solution (ppt).
 
     Returns
     -------
@@ -34,13 +45,13 @@ def water(phase, **kwargs):
 
     References
     ----------
-    [1] Sharqawy M. H., Lienhard J. H., and Zubair, S. M., Desalination and Water
-        Treatment, 2010.
+    [1] Sharqawy M. H., Lienhard J. H., and Zubair, S. M., Desalination and
+        Water Treatment, 2010.
 
     """
-    T = phase['pore.temperature']
+    T = phase[pore_T]
     try:
-        S = phase['pore.salinity']
+        S = phase[pore_salinity]
     except:
         S = 0
     TC = T-273.15
@@ -63,7 +74,9 @@ def water(phase, **kwargs):
     return value
 
 
-def reynolds(phase, uo, b, **kwargs):
+def reynolds(phase, uo, b,
+             pore_T='pore.temperature',
+             **kwargs):
     r"""
     Uses exponential model by Reynolds [1]_ for the temperature dependance of
     shear viscosity
@@ -74,38 +87,50 @@ def reynolds(phase, uo, b, **kwargs):
             Coefficients of the viscosity exponential model (mu = u0*Exp(-b*T)
             where T is the temperature in Kelvin
 
+    pore_T : string
+        The dictionary key containing the temperature values (K)
+
     [1] Reynolds O. (1886). Phil Trans Royal Soc London, v. 177, p.157.
 
     """
-    T = phase['pore.temperature']
+    T = phase[pore_T]
     value = uo*sp.exp(b*T)
     return value
 
 
-def chung(phase, MW='molecular_weight', Tc='critical_temperature',
-          Vc='critical_volume', **kwargs):
+def chung(phase,
+          pore_T='pore.temperature',
+          pore_MW='pore.molecular_weight',
+          pore_Tc='pore.critical_temperature',
+          pore_Vc='pore.critical_volume',
+          **kwargs):
     r"""
-    Uses Chung et al. [2]_ model to estimate viscosity for gases with low pressure
-    (much less than the critical pressure) at conditions of interest
+    Uses Chung et al. [2]_ model to estimate viscosity for gases with low
+    pressure (much less than the critical pressure) at conditions of interest.
 
     Parameters
     ----------
-    Vc :  float, array_like
-        Critical volume of the gas (m3/kmol)
-    Tc :  float, array_like
-        Critical temperature of the gas (K)
-    MW : float, array_like
-        Molecular weight of the gas (kg/kmol)
+    pore_T : string
+        The dictionary key containing the temperature values (K)
+
+    pore_Tc : string
+        The dictionary key containing the temperature values (K)
+
+    pore_MW : string
+        The dictionary key containing the molecular weight values (kg/mol)
+
+    pore_Vc : string
+        The dictionary key containing the critical volume values (m3/kmol)
 
     [2] Chung, T.H., Lee, L.L., and Starling, K.E., Applications of Kinetic Gas
         Theories and Multiparameter Correlation for Prediction of Dilute Gas
         Viscosity and Thermal Conductivity”, Ind. Eng. Chem. Fundam.23:8, 1984.
 
     """
-    T = phase['pore.temperature']
-    MW = phase['pore.'+MW]
-    Tc = phase['pore.'+Tc]
-    Vc = phase['pore.'+Vc]
+    T = phase[pore_T]
+    MW = phase[pore_MW]
+    Tc = phase[pore_Tc]
+    Vc = phase[pore_Vc]
     Tr = T / Tc
     Tstar = 1.2593*Tr
     A = 1.161415
diff --git a/OpenPNM/Postprocessing/Plots.py b/OpenPNM/Postprocessing/Plots.py
index 51d0ce64c..eab2be8a0 100644
--- a/OpenPNM/Postprocessing/Plots.py
+++ b/OpenPNM/Postprocessing/Plots.py
@@ -280,29 +280,8 @@ def drainage_curves(inv_alg, fig=None, Pc='inv_Pc', sat='inv_sat',
     The invasion algorithm for which the graphs are desired
 
     timing : string
-    if algorithm keeps track of simulated time, insert string here.
-
-    Examples
-    --------
-    >>> import OpenPNM
-    >>> pn = OpenPNM.Network.TestNet()
-    >>> geo = OpenPNM.Geometry.TestGeometry(network=pn,
-    ...                                     pores=pn.pores(),
-    ...                                     throats=pn.throats())
-    >>> phase1 = OpenPNM.Phases.TestPhase(network=pn)
-    >>> phase2 = OpenPNM.Phases.TestPhase(network=pn)
-    >>> phys1 = OpenPNM.Physics.TestPhysics(network=pn, phase=phase1,
-    ...                                     pores=pn.pores(),throats=pn.throats())
-    >>> phys2 = OpenPNM.Physics.TestPhysics(network=pn, phase=phase2,
-    ...                                     pores=pn.pores(),throats=pn.throats())
-    >>> IP = OpenPNM.Algorithms.InvasionPercolationTimed(network=pn)
-    >>> IP.run(invading_phase=phase1, defending_phase=phase2,
-    ...        inlets=pn.pores('top'), outlets=pn.pores('bottom'))
-         IP algorithm at 0 % completion at 0.0 seconds
-         IP algorithm at 20 % completion at 0.0 seconds
-         IP algorithm at 40 % completion at 0.0 seconds
-         IP algorithm at 60 % completion at 0.0 seconds
-         IP algorithm at 100% completion at  0.0  seconds
+    if algorithm keeps track of simulated time, insert string here
+
     """
     inv_throats = inv_alg.toindices(inv_alg['throat.' + seq] > 0)
     sort_seq = _sp.argsort(inv_alg['throat.'+seq][inv_throats])
diff --git a/docs/index.rst b/docs/index.rst
index 066d0b981..dcfd9e4af 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -36,6 +36,7 @@ Documentation
 
    userguide/index.rst
    examples/index.rst
+   resources/index.rst
    docOpenPNM.rst
    docExtraSidebar.rst
 
diff --git a/docs/resources/index.rst b/docs/resources/index.rst
new file mode 100644
index 000000000..0cd4f39f3
--- /dev/null
+++ b/docs/resources/index.rst
@@ -0,0 +1,12 @@
+.. _resources:
+
+###############################################################################
+Resources
+###############################################################################
+
+This section contains recipes and various code snippets which may be useful but were considered too specific or arcane to include in the main package.
+
+.. toctree::
+   :maxdepth: 1
+   
+   pore_scale_models.rst
diff --git a/docs/resources/pore_scale_models.rst b/docs/resources/pore_scale_models.rst
new file mode 100644
index 000000000..123966956
--- /dev/null
+++ b/docs/resources/pore_scale_models.rst
@@ -0,0 +1,81 @@
+.. _pore_scale_models:
+
+###############################################################################
+Spatially correlated pore seeds using Perlin noise
+###############################################################################
+Perlin noise [ref?] is an algorithm for creating spatially correlated random
+noise.  
+
+.. code-block:: python
+
+    def perlin_noise(geometry, freq=1, octaves=4, mode='classic', **kwargs):
+        r"""
+        Generate pore seed values using the Perlin noise algorithm.  This approach
+        imparts some spatial clumpiness to the pore seeds.
+
+        Parameters
+        ----------
+        freq, octaves : int
+            Parameters that control the qualities of the noise.  Lower frequency
+            results in more smaller clumps.  Higher octaves gives a more textured
+            noise.
+        mode : {'classic','simplex'}
+            The algorithm to use when generating the noise.
+
+            * 'classic' : (default) The original algorithm developed by Perlin
+            * 'simplex' : A newer algorithm that is supposedly faster and results
+            in a more natural texture.
+
+        Returns
+        -------
+        A list of pore seeds values between 0:1 that are spatially correlated (i.e.
+        similar values are clumped together)
+
+        Notes
+        -----
+        - This method uses image analysis type tools, so only works on Cubic
+        networks
+        - This method requires the 'noise' module is installed
+
+        Examples
+        --------
+        >>> import OpenPNM
+        >>> pn = OpenPNM.Network.Cubic(shape=[50, 50, 50])
+        >>> geom = OpenPNM.Geometry.GenericGeometry(network=pn,pores=pn.Ps,throats=pn.Ts)
+        >>> geom.add_model(propname='pore.seed',
+        ...                model=OpenPNM.Geometry.models.pore_seed.perlin_noise)
+        >>> im = pn.asarray(geom['pore.seed'])
+
+        Visualizing the end result can be done with:
+
+        .. code-block:: python
+
+            matplotlib.pyplot.imshow(im[:, 25, :],interpolation='none')
+
+        """
+        from noise import pnoise3, snoise3
+        import scipy.stats as spst
+
+        net = geometry._net
+        if mode == 'classic':
+            model = pnoise3
+        elif mode == 'simplex':
+            model = snoise3
+        freq = freq * octaves
+        # The following will only work on Cubic networks
+        x = net._shape[0]
+        y = net._shape[1]
+        z = net._shape[2]
+        temp = _sp.ndarray((x, y, z))
+        for k in range(z):
+            for j in range(y):
+                for i in range(x):
+                    temp[i, j, k] = model(i / freq, j / freq, k / freq, octaves) + 0.5
+        # Assuming that the noise is normally distributed, find seeds of that dist
+        temp = _sp.reshape(temp, (temp.size,))
+        x_mean = _sp.mean(temp)
+        x_sigma = _sp.sqrt(1/(temp.size-1)*_sp.sum((temp - x_mean)**2))
+        fn1 = spst.norm(loc=x_mean, scale=x_sigma)
+        values = fn1.cdf(temp)
+        values = values[geometry.map_pores(target=net, pores=geometry.Ps)]
+        return values.flatten()
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 9da9e37e2..5daa826c5 100644
--- a/setup.py
+++ b/setup.py
@@ -40,9 +40,7 @@ setup(
     install_requires = [
         'numpy',
         'scipy>=0.14.0',
-        'matplotlib',
-        'noise',
-        'scikit-image'
+        'matplotlib'
     ],
     author = 'OpenPNM Team',
     author_email = 'jeff.gostick@mcgill.ca',


Test patch:
diff --git a/OpenPNM/Physics/__TestPhysics__.py b/OpenPNM/Physics/__TestPhysics__.py
index 8550c7988..f59a7886d 100644
--- a/OpenPNM/Physics/__TestPhysics__.py
+++ b/OpenPNM/Physics/__TestPhysics__.py
@@ -46,6 +46,3 @@ class TestPhysics(GenericPhysics):
                 self['throat.thermal_conductance'] = \
                     phase['throat.thermal_conductivity'] * \
                     self._net['throat.diameter'] / self._net['throat.length']
-
-if __name__ == '__main__':
-    print('none yet')
diff --git a/test/integration/script_test.py b/test/integration/script_test.py
index 0594eeeda..ec3434f39 100644
--- a/test/integration/script_test.py
+++ b/test/integration/script_test.py
@@ -3,13 +3,17 @@ import scipy as sp
 import OpenPNM
 import pytest
 
+
 def test_controller():
-    # The following tests check if the Controller's ghost, save, load and clear methods work
+    # The following tests check if the Controller's ghost, save, load and
+    # clear methods work
     ctrl = OpenPNM.Base.Controller()
     ctrl.clear()  # Clear controller to make sure it has no lingering objects
     ctrl.loglevel = 50
     pn = OpenPNM.Network.TestNet()
-    geom = OpenPNM.Geometry.TestGeometry(network=pn,pores=pn.Ps,throats=pn.Ts)
+    geom = OpenPNM.Geometry.TestGeometry(network=pn,
+                                         pores=pn.Ps,
+                                         throats=pn.Ts)
     loc1 = geom.__repr__()
     loc2 = ctrl[geom.name].__repr__()
     assert loc1 == loc2  # Ensure both locations hold same object
@@ -24,12 +28,12 @@ def test_controller():
     assert geom2._net is pn  # Ghost geometry has same ancestor Network
     assert geom2 not in pn._geometries  # Network is unaware of ghost Geometry
     # Test clear functionality
-    ctrl.save('test_net')  # Save current state
+    ctrl.save('test')  # Save current state
     ctrl.clear()  # Clear Controller
     assert ctrl.keys() == {}.keys()  # Empty dict
     assert pn.controller == {}  # Controller is now an empty dict
     assert geom.controller == pn.controller  # Both have empty dict's
-    ctrl.load('test_net')  # Load saved state
+    ctrl.load('test')  # Load saved state
     assert pn.name in ctrl.keys()  # Ensure loaded objects match originals
     assert geom.name in ctrl.keys()
     pn2 = ctrl[pn.name]  # Retrieve loaded Network from the Controller dict
@@ -41,217 +45,287 @@ def test_controller():
     assert geom2.name not in ctrl.keys()  # Geometry is purged from Controller
     assert 'pore.'+geom2.name not in pn2.keys()  # Geometry label is removed from Simulation objects too
 
+
 def test_cubic_standard_call():
-    pn = OpenPNM.Network.Cubic(shape=[3,4,5])
-    np.testing.assert_almost_equal(pn['pore.coords'][0], [0.5,0.5,0.5])
+    pn = OpenPNM.Network.Cubic(shape=[3, 4, 5])
+    np.testing.assert_almost_equal(pn['pore.coords'][0], [0.5, 0.5, 0.5])
+
 
 def test_trim_extend():
-    pn = OpenPNM.Network.Cubic(shape=[5,5,5])
-    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0),[ 1,  5, 25]))
-    assert [pn.Np,pn.Nt] == [125,300]
+    pn = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0), [1, 5, 25]))
+    assert [pn.Np, pn.Nt] == [125, 300]
     pn.trim(pores=[0])
-    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0),[ 1,  5, 25]))
-    assert [pn.Np,pn.Nt] == [124,297]
-    pn.extend(pore_coords=[0,0,0],throat_conns=[[124,0]])
-    assert [pn.Np,pn.Nt] == [125,298]
-    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0),[ 1,  5, 25, 124]))
+    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0), [1, 5, 25]))
+    assert [pn.Np, pn.Nt] == [124, 297]
+    pn.extend(pore_coords=[0, 0, 0], throat_conns=[[124, 0]])
+    assert [pn.Np, pn.Nt] == [125, 298]
+    assert sp.all(sp.in1d(pn.find_neighbor_pores(pores=0), [1, 5, 25, 124]))
+
 
 def test_stitch():
     ctrl = OpenPNM.Base.Controller()
-    [Nx,Ny,Nz]=[10,10,10]
-    pn = OpenPNM.Network.Cubic(shape=[Nx,Ny,Nz])
-    pn2 = OpenPNM.Network.Cubic(shape=[Nx,Ny,Nz])
-    pn2['pore.coords'][:,2] += Nz
-    pn.stitch(donor=pn2,P_network=pn.pores('top'),P_donor=pn2.pores('bottom'),len_max=1,method='nearest')
+    [Nx, Ny, Nz] = [10, 10, 10]
+    pn = OpenPNM.Network.Cubic(shape=[Nx, Ny, Nz])
+    pn2 = OpenPNM.Network.Cubic(shape=[Nx, Ny, Nz])
+    pn2['pore.coords'][:, 2] += Nz
+    pn.stitch(donor=pn2,
+              P_network=pn.pores('top'),
+              P_donor=pn2.pores('bottom'),
+              len_max=1,
+              method='nearest')
     assert pn.Np == 2*pn2.Np  # Ensure number of pores doubled
-    assert pn.Nt == 2*pn2.Nt + Nx*Ny  # Ensure correct number of new stitch throats
+    assert pn.Nt == (2*pn2.Nt + Nx*Ny)  # Ensure correct number of new throats
     assert pn2 not in ctrl.values()  # Donor Network is removed from Controller
     # Reuse the donor Network in another stitch
-    pn2['pore.coords'][:,2] -= 2*Nz
-    pn.stitch(donor=pn2,P_network=pn.pores('bottom'),P_donor=pn2.pores('top'),len_max=1,method='nearest')
+    pn2['pore.coords'][:, 2] -= 2*Nz
+    pn.stitch(donor=pn2,
+              P_network=pn.pores('bottom'),
+              P_donor=pn2.pores('top'),
+              len_max=1,
+              method='nearest')
     assert pn.Np == 3*pn2.Np  # Ensure number of pores increased again
-    assert pn.Nt == 3*pn2.Nt + 2*Nx*Ny  # Ensure correct number of new stitch throats
+    assert pn.Nt == (3*pn2.Nt + 2*Nx*Ny)  # Ensure correct number of new throats
     ctrl.clear()
 
+
 def test_linear_solvers():
-  pn = OpenPNM.Network.Cubic([1,40,30], spacing=0.0001)
-  geom = OpenPNM.Geometry.Toray090(network=pn,pores=pn.pores(),throats=pn.throats())
-  air = OpenPNM.Phases.Air(network=pn)
-  phys_air = OpenPNM.Physics.Standard(network=pn,phase=air,pores=pn.pores(),throats=pn.throats())
-
-  BC1_pores = pn.pores(labels=['left'])
-  BC2_pores = pn.pores(labels=['right'])
-
-  alg_1 = OpenPNM.Algorithms.FickianDiffusion(network=pn,phase=air)
-  alg_1.set_boundary_conditions(bctype='Dirichlet', bcvalue=1, pores=BC1_pores)
-  alg_1.set_boundary_conditions(bctype='Dirichlet', bcvalue=0, pores=BC2_pores)
-  alg_1.run()
-
-  alg_2 = OpenPNM.Algorithms.FickianDiffusion(network=pn,phase=air)
-  alg_2.set_boundary_conditions(bctype='Neumann', bcvalue = -1e-11, pores=BC1_pores)
-  alg_2.set_boundary_conditions(bctype='Dirichlet', bcvalue=0, pores=BC2_pores)
-  alg_2.run()
-
-  alg_3 = OpenPNM.Algorithms.FickianDiffusion(network=pn,phase=air)
-  alg_3.set_boundary_conditions(bctype='Neumann_group', bcvalue = -3e-10, pores=BC1_pores)
-  alg_3.set_boundary_conditions(bctype='Dirichlet', bcvalue=0, pores=BC2_pores)
-  alg_3.run()
-
-  alg_4 = OpenPNM.Algorithms.FickianDiffusion(network=pn,phase=air)
-  alg_4.set_boundary_conditions(bctype='Neumann_group', bcvalue = -3e-10, pores=BC1_pores)
-  alg_4.set_boundary_conditions(bctype='Dirichlet', bcvalue=0, pores=BC2_pores)
-  alg_4.setup()
-  alg_4.solve()
-
-  print( alg_1['pore.'+air.name+'_mole_fraction'][BC1_pores] )
-  print( alg_2['pore.'+air.name+'_mole_fraction'][BC1_pores] )
-  print( alg_3['pore.'+air.name+'_mole_fraction'][BC1_pores] )
-
-  assert round(sp.absolute(alg_1.rate(BC1_pores))[0],16) == round(sp.absolute(alg_1.rate(BC2_pores))[0],16)
-  assert round(sp.absolute(alg_2.rate(BC2_pores))[0],16) == round(sp.absolute(sp.unique(alg_2['pore.'+air.name+'_bcval_Neumann']))[0]*len(BC1_pores),16)
-  assert round(sp.absolute(alg_3.rate(BC2_pores))[0],16) == round(sp.absolute(sp.unique(alg_3['pore.'+air.name+'_bcval_Neumann_group']))[0],16)
-  assert round(sp.absolute(alg_4.rate(BC2_pores))[0],16) == round(sp.absolute(sp.unique(alg_4['pore.'+air.name+'_bcval_Neumann_group']))[0],16)
-
-  assert round(sp.absolute(sp.sum(alg_1.rate(BC1_pores,mode='single'))),16) == round(sp.absolute(alg_1.rate(BC1_pores))[0],16)
-  assert round(sp.absolute(sp.sum(alg_2.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_2.rate(BC2_pores))[0],16)
-  assert round(sp.absolute(sp.sum(alg_3.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_3.rate(BC2_pores))[0],16)
-  assert round(sp.absolute(sp.sum(alg_4.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_4.rate(BC2_pores))[0],16)
+    pn = OpenPNM.Network.Cubic([1, 40, 30], spacing=0.0001)
+    geom = OpenPNM.Geometry.Toray090(network=pn,
+                                     pores=pn.pores(),
+                                     throats=pn.throats())
+    air = OpenPNM.Phases.Air(network=pn)
+    phys_air = OpenPNM.Physics.Standard(network=pn,
+                                        phase=air,
+                                        pores=pn.pores(),
+                                        throats=pn.throats())
+
+    BC1_pores = pn.pores(labels=['left'])
+    BC2_pores = pn.pores(labels=['right'])
+
+    alg_1 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
+    alg_1.set_boundary_conditions(bctype='Dirichlet',
+                                  bcvalue=1,
+                                  pores=BC1_pores)
+    alg_1.set_boundary_conditions(bctype='Dirichlet',
+                                  bcvalue=0,
+                                  pores=BC2_pores)
+    alg_1.run()
+
+    alg_2 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
+    alg_2.set_boundary_conditions(bctype='Neumann',
+                                  bcvalue=-1e-11,
+                                  pores=BC1_pores)
+    alg_2.set_boundary_conditions(bctype='Dirichlet',
+                                  bcvalue=0,
+                                  pores=BC2_pores)
+    alg_2.run()
+
+    alg_3 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
+    alg_3.set_boundary_conditions(bctype='Neumann_group',
+                                  bcvalue=-3e-10,
+                                  pores=BC1_pores)
+    alg_3.set_boundary_conditions(bctype='Dirichlet',
+                                  bcvalue=0,
+                                  pores=BC2_pores)
+    alg_3.run()
+
+    alg_4 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
+    alg_4.set_boundary_conditions(bctype='Neumann_group',
+                                  bcvalue=-3e-10,
+                                  pores=BC1_pores)
+    alg_4.set_boundary_conditions(bctype='Dirichlet',
+                                  bcvalue=0,
+                                  pores=BC2_pores)
+    alg_4.setup()
+    alg_4.solve()
+
+    print(alg_1['pore.'+air.name+'_mole_fraction'][BC1_pores])
+    print(alg_2['pore.'+air.name+'_mole_fraction'][BC1_pores])
+    print(alg_3['pore.'+air.name+'_mole_fraction'][BC1_pores])
+
+    assert round(sp.absolute(alg_1.rate(BC1_pores))[0], 16) == round(sp.absolute(alg_1.rate(BC2_pores))[0], 16)
+    assert round(sp.absolute(alg_2.rate(BC2_pores))[0], 16) == round(sp.absolute(sp.unique(alg_2['pore.'+air.name+'_bcval_Neumann']))[0]*len(BC1_pores), 16)
+    assert round(sp.absolute(alg_3.rate(BC2_pores))[0], 16) == round(sp.absolute(sp.unique(alg_3['pore.'+air.name+'_bcval_Neumann_group']))[0], 16)
+    assert round(sp.absolute(alg_4.rate(BC2_pores))[0], 16) == round(sp.absolute(sp.unique(alg_4['pore.'+air.name+'_bcval_Neumann_group']))[0], 16)
+
+    assert round(sp.absolute(sp.sum(alg_1.rate(BC1_pores,mode='single'))),16) == round(sp.absolute(alg_1.rate(BC1_pores))[0],16)
+    assert round(sp.absolute(sp.sum(alg_2.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_2.rate(BC2_pores))[0],16)
+    assert round(sp.absolute(sp.sum(alg_3.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_3.rate(BC2_pores))[0],16)
+    assert round(sp.absolute(sp.sum(alg_4.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_4.rate(BC2_pores))[0],16)
 
 
 def test_add_boundary():
-  pn = OpenPNM.Network.Cubic(shape=[5,5,5])
-  pn.add_boundaries()
+    pn = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+    pn.add_boundaries()
 
-  keys_expected = {'pore.back', 'pore.bottom', 'pore.top_boundary',
-  'pore.right_boundary', 'throat.back_boundary', 'throat.all',
-  'throat.bottom_boundary', 'throat.front_boundary', 'pore.boundary',
-  'throat.left_boundary', 'throat.conns', 'throat.top_boundary',
-  'pore.back_boundary', 'pore.top', 'pore.front_boundary', 'pore.all',
-  'pore.front', 'pore.left_boundary', 'throat.boundary', 'pore.bottom_boundary',
-  'throat.right_boundary', 'pore.coords', 'pore.internal', 'pore.index',
-  'pore.left', 'pore.right'}
+    keys_expected = {'pore.back', 'pore.bottom', 'pore.top_boundary',
+                     'pore.right_boundary', 'throat.back_boundary',
+                     'throat.all', 'throat.bottom_boundary',
+                     'throat.front_boundary', 'pore.boundary',
+                     'throat.left_boundary', 'throat.conns',
+                     'throat.top_boundary', 'pore.back_boundary', 'pore.top',
+                     'pore.front_boundary', 'pore.all', 'pore.front',
+                     'pore.left_boundary', 'throat.boundary',
+                     'pore.bottom_boundary', 'throat.right_boundary',
+                     'pore.coords', 'pore.internal', 'pore.index', 'pore.left',
+                     'pore.right'}
 
-  keys_found = set(pn.keys())
+    keys_found = set(pn.keys())
+
+    symmetric_diff = keys_found ^ keys_expected
+    assert not symmetric_diff
 
-  symmetric_diff = keys_found ^ keys_expected
-  assert not symmetric_diff
 
 def test_open_air_diffusivity():
-    pn = OpenPNM.Network.Cubic([5,5,5], spacing=1)
+    pn = OpenPNM.Network.Cubic([5, 5, 5], spacing=1)
     pn.add_boundaries()
-    Ps = pn.pores('boundary',mode='not')
-    Ts = pn.find_neighbor_throats(pores=Ps,mode='intersection',flatten=True)
-    geom = OpenPNM.Geometry.Cube_and_Cuboid(network=pn,pores=Ps,throats=Ts)
+    Ps = pn.pores('boundary', mode='not')
+    Ts = pn.find_neighbor_throats(pores=Ps, mode='intersection', flatten=True)
+    geom = OpenPNM.Geometry.Cube_and_Cuboid(network=pn, pores=Ps, throats=Ts)
     geom['pore.diameter'] = 0.999999
     geom['throat.diameter'] = 0.999999
-    geom.regenerate(['pore.diameter','throat.diameter'],mode='exclude')
+    geom.regenerate(['pore.diameter', 'throat.diameter'], mode='exclude')
     Ps = pn.pores('boundary')
-    Ts = pn.find_neighbor_throats(pores=Ps,mode='not_intersection')
-    boun = OpenPNM.Geometry.Boundary(network=pn,pores=Ps,throats=Ts,shape='cubes')
+    Ts = pn.find_neighbor_throats(pores=Ps, mode='not_intersection')
+    boun = OpenPNM.Geometry.Boundary(network=pn,
+                                     pores=Ps,
+                                     throats=Ts,
+                                     shape='cubes')
     air = OpenPNM.Phases.Air(network=pn)
     Ps = pn.pores()
     Ts = pn.throats()
-    phys_air = OpenPNM.Physics.Standard(network=pn,phase=air,pores=Ps,throats=Ts)
+    phys_air = OpenPNM.Physics.Standard(network=pn,
+                                        phase=air,
+                                        pores=Ps,
+                                        throats=Ts)
     BC1_pores = pn.pores(labels=['top_boundary'])
     BC2_pores = pn.pores(labels=['bottom_boundary'])
-    Diff = OpenPNM.Algorithms.FickianDiffusion(network=pn,phase=air)
+    Diff = OpenPNM.Algorithms.FickianDiffusion(network=pn,
+                                               phase=air)
     # Assign Dirichlet boundary conditions to top and bottom surface pores
-    Diff.set_boundary_conditions(bctype='Dirichlet', bcvalue=0.6, pores=BC1_pores)
-    Diff.set_boundary_conditions(bctype='Dirichlet', bcvalue=0.4, pores=BC2_pores)
+    Diff.set_boundary_conditions(bctype='Dirichlet',
+                                 bcvalue=0.6,
+                                 pores=BC1_pores)
+    Diff.set_boundary_conditions(bctype='Dirichlet',
+                                 bcvalue=0.4,
+                                 pores=BC2_pores)
     Diff.run()
     Diff.return_results()
     Diff_deff = Diff.calc_eff_diffusivity()/np.mean(air['pore.diffusivity'])
-    assert np.round(Diff_deff,3) == 1
+    assert np.round(Diff_deff, 3) == 1
+
 
 def test_thermal_conduction():
-    #Generate Network and clean up boundaries (delete z-face pores)
-    divs = [10,50]
-    Lc   = 0.1  # cm
-    pn = OpenPNM.Network.Cubic(shape= divs, spacing = Lc)
+    # Generate Network and clean up boundaries (delete z-face pores)
+    divs = [10, 50]
+    Lc = 0.1  # cm
+    pn = OpenPNM.Network.Cubic(shape=divs, spacing=Lc)
     pn.add_boundaries()
-    pn.trim(pores=pn.pores(['top_boundary','bottom_boundary']))
-    #Generate Geometry objects for internal and boundary pores
+    pn.trim(pores=pn.pores(['top_boundary', 'bottom_boundary']))
+    # Generate Geometry objects for internal and boundary pores
     Ps = pn.pores('internal')
     Ts = pn.throats()
-    geom = OpenPNM.Geometry.GenericGeometry(network=pn,pores=Ps,throats=Ts)
-    geom['pore.area']     = Lc**2
+    geom = OpenPNM.Geometry.GenericGeometry(network=pn,
+                                            pores=Ps,
+                                            throats=Ts)
+    geom['pore.area'] = Lc**2
     geom['pore.diameter'] = Lc
     geom['throat.length'] = 1e-25
-    geom['throat.area']   = Lc**2
+    geom['throat.area'] = Lc**2
     Ps = pn.pores('boundary')
-    boun = OpenPNM.Geometry.GenericGeometry(network=pn,pores=Ps)
-    boun['pore.area']     = Lc**2
-    boun['pore.diameter'] =  1e-25
-    #Create Phase object and associate with a Physics object
+    boun = OpenPNM.Geometry.GenericGeometry(network=pn, pores=Ps)
+    boun['pore.area'] = Lc**2
+    boun['pore.diameter'] = 1e-25
+    # Create Phase object and associate with a Physics object
     Cu = OpenPNM.Phases.GenericPhase(network=pn)
     Cu['pore.thermal_conductivity'] = 1.0  # W/m.K
-    phys = OpenPNM.Physics.GenericPhysics(network=pn,phase=Cu,pores=pn.pores(),throats=pn.throats())
+    phys = OpenPNM.Physics.GenericPhysics(network=pn,
+                                          phase=Cu,
+                                          pores=pn.pores(),
+                                          throats=pn.throats())
     mod = OpenPNM.Physics.models.thermal_conductance.series_resistors
-    phys.add_model(propname='throat.thermal_conductance',model=mod)
+    phys.add_model(propname='throat.thermal_conductance', model=mod)
     phys.regenerate()  # Update the conductance values
-    #Setup Algorithm object
-    Fourier_alg = OpenPNM.Algorithms.FourierConduction(network=pn,phase=Cu)
+    # Setup Algorithm object
+    Fourier_alg = OpenPNM.Algorithms.FourierConduction(network=pn, phase=Cu)
     inlets = pn.pores('back_boundary')
-    outlets = pn.pores(['front_boundary','left_boundary','right_boundary'])
-    T_in = 30*sp.sin(sp.pi*pn['pore.coords'][inlets,1]/5)+50
-    Fourier_alg.set_boundary_conditions(bctype='Dirichlet',bcvalue=T_in,pores=inlets)
-    Fourier_alg.set_boundary_conditions(bctype='Dirichlet',bcvalue=50,pores=outlets)
+    outlets = pn.pores(['front_boundary', 'left_boundary', 'right_boundary'])
+    T_in = 30*sp.sin(sp.pi*pn['pore.coords'][inlets, 1]/5)+50
+    Fourier_alg.set_boundary_conditions(bctype='Dirichlet',
+                                        bcvalue=T_in,
+                                        pores=inlets)
+    Fourier_alg.set_boundary_conditions(bctype='Dirichlet',
+                                        bcvalue=50,
+                                        pores=outlets)
     Fourier_alg.run()
     Fourier_alg.return_results()
-    #Calculate analytical solution over the same domain spacing
-    Cu['pore.analytical_temp'] = 30*sp.sinh(sp.pi*pn['pore.coords'][:,0]/5)/sp.sinh(sp.pi/5)*sp.sin(sp.pi*pn['pore.coords'][:,1]/5) + 50
+    # Calculate analytical solution over the same domain spacing
+    Cu['pore.analytical_temp'] = 30*sp.sinh(sp.pi*pn['pore.coords'][:, 0]/5)/sp.sinh(sp.pi/5)*sp.sin(sp.pi*pn['pore.coords'][:, 1]/5) + 50
     b = Cu['pore.analytical_temp'][pn.pores(geom.name)]
     a = Cu['pore.temperature'][pn.pores(geom.name)]
-    a = sp.reshape(a,(divs[0],divs[1]))
-    b = sp.reshape(b,(divs[0],divs[1]))
+    a = sp.reshape(a, (divs[0], divs[1]))
+    b = sp.reshape(b, (divs[0], divs[1]))
     diff = a - b
     assert sp.amax(np.absolute(diff)) < 0.015
 
+
 def test_Darcy_alg():
-    #Generate Network and clean up some of boundaries
-    divs = [1,50,10]
+    # Generate Network and clean up some of boundaries
+    divs = [1, 50, 10]
     Lc = 0.00004
-    pn = OpenPNM.Network.Cubic(shape = divs, spacing = Lc)
+    pn = OpenPNM.Network.Cubic(shape=divs, spacing=Lc)
     pn.add_boundaries()
-    Ps = pn.pores(['front_boundary','back_boundary'])
+    Ps = pn.pores(['front_boundary', 'back_boundary'])
     pn.trim(pores=Ps)
-    #Generate Geometry objects for internal and boundary pores
-    Ps = pn.pores('boundary',mode='not')
-    Ts = pn.find_neighbor_throats(pores=Ps,mode='intersection',flatten=True)
-    geom = OpenPNM.Geometry.Toray090(network=pn,pores=Ps,throats=Ts)
+    # Generate Geometry objects for internal and boundary pores
+    Ps = pn.pores('boundary', mode='not')
+    Ts = pn.find_neighbor_throats(pores=Ps, mode='intersection', flatten=True)
+    geom = OpenPNM.Geometry.Toray090(network=pn, pores=Ps, throats=Ts)
     Ps = pn.pores('boundary')
-    Ts = pn.find_neighbor_throats(pores=Ps,mode='not_intersection')
-    boun = OpenPNM.Geometry.Boundary(network=pn,pores=Ps,throats=Ts)
-    #Create Phase object and associate with a Physics object
+    Ts = pn.find_neighbor_throats(pores=Ps, mode='not_intersection')
+    boun = OpenPNM.Geometry.Boundary(network=pn, pores=Ps, throats=Ts)
+    # Create Phase object and associate with a Physics object
     air = OpenPNM.Phases.Air(network=pn)
     Ps = pn.pores()
     Ts = pn.throats()
-    phys = OpenPNM.Physics.GenericPhysics(network=pn,phase=air,pores=Ps,throats=Ts)
+    phys = OpenPNM.Physics.GenericPhysics(network=pn,
+                                          phase=air,
+                                          pores=Ps,
+                                          throats=Ts)
     from OpenPNM.Physics import models as pm
     phys.add_model(propname='throat.hydraulic_conductance',
                    model=pm.hydraulic_conductance.hagen_poiseuille,
                    calc_pore_len=False)
     phys.regenerate()  # Update the conductance values
-    #Setup Algorithm objects
-    Darcy1 = OpenPNM.Algorithms.StokesFlow(network=pn,phase=air)
+    # Setup Algorithm objects
+    Darcy1 = OpenPNM.Algorithms.StokesFlow(network=pn, phase=air)
     inlets = pn.pores('bottom_boundary')
     Ps = pn.pores('top_boundary')
-    outlets = Ps[pn['pore.coords'][Ps,1]<(divs[1]*Lc/2)]
+    outlets = Ps[pn['pore.coords'][Ps, 1] < (divs[1]*Lc/2)]
     P_out = 0  # Pa
     Q_in = 0.6667*(Lc**2)*divs[1]*divs[0]  # m^3/s
-    Darcy1.set_boundary_conditions(bctype='Neumann_group',bcvalue=-Q_in,pores=inlets)
-    Darcy1.set_boundary_conditions(bctype='Dirichlet',bcvalue=P_out,pores=outlets)
+    Darcy1.set_boundary_conditions(bctype='Neumann_group',
+                                   bcvalue=-Q_in,
+                                   pores=inlets)
+    Darcy1.set_boundary_conditions(bctype='Dirichlet',
+                                   bcvalue=P_out,
+                                   pores=outlets)
     Darcy1.run()
     Darcy1.return_results()
     print('pore pressure for Darcy1 algorithm:')
     print(air['pore.pressure'])
-    Darcy2 = OpenPNM.Algorithms.StokesFlow(network=pn,phase=air)
+    Darcy2 = OpenPNM.Algorithms.StokesFlow(network=pn, phase=air)
     inlets = pn.pores('bottom_boundary')
     outlets = pn.pores('top_boundary')
     P_out = 10  # Pa
     P_in = 1000  # Pa
-    Darcy2.set_boundary_conditions(bctype='Dirichlet',bcvalue=P_in,pores=inlets)
-    Darcy2.set_boundary_conditions(bctype='Dirichlet',bcvalue=P_out,pores=outlets)
+    Darcy2.set_boundary_conditions(bctype='Dirichlet',
+                                   bcvalue=P_in,
+                                   pores=inlets)
+    Darcy2.set_boundary_conditions(bctype='Dirichlet',
+                                   bcvalue=P_out,
+                                   pores=outlets)
     Darcy2.run()
     print('pore pressure for Darcy2 algorithm:')
     print(Darcy2['pore.'+air.name+'_pressure'])
@@ -260,24 +334,31 @@ def test_Darcy_alg():
     Vp = sp.sum(pn['pore.volume']) + sp.sum(pn['throat.volume'])
     Vb = sp.prod(divs)*Lc**3
     e = Vp/Vb
-    print('Effective permeability: ',K,'- Porosity: ',e)
+    print('Effective permeability: ', K, '- Porosity: ', e)
+
+    a = round(sp.absolute(Darcy1.rate(outlets))[0], 16)
+    pore_prop = 'pore.'+air.name+'_bcval_Neumann_group'
+    b = round(sp.absolute(sp.unique(Darcy1[pore_prop]))[0], 16)
+    assert a == b
+
+    a = round(sp.absolute(Darcy2.rate(inlets))[0], 16)
+    b = round(sp.absolute(Darcy2.rate(outlets))[0], 16)
+    assert a == b
 
-    assert round(sp.absolute(Darcy1.rate(outlets))[0],16) == round(sp.absolute(sp.unique(Darcy1['pore.'+air.name+'_bcval_Neumann_group']))[0],16)
-    assert round(sp.absolute(Darcy2.rate(inlets))[0],16) == round(sp.absolute(Darcy2.rate(outlets))[0],16)
 
 def test_mapping():
     # Create small cubic network
-    pn = OpenPNM.Network.Cubic(shape=[3,3,3],spacing=0.0001)
+    pn = OpenPNM.Network.Cubic(shape=[3, 3, 3], spacing=0.0001)
     # Assign 3 different geometries to each layer in the z-direction
-    Pa = sp.arange(0,9)
+    Pa = sp.arange(0, 9)
     Ta = pn.find_neighbor_throats(Pa)
-    geom1 = OpenPNM.Geometry.GenericGeometry(network=pn,pores=Pa,throats=Ta)
-    Pc = sp.arange(18,27)
+    geom1 = OpenPNM.Geometry.GenericGeometry(network=pn, pores=Pa, throats=Ta)
+    Pc = sp.arange(18, 27)
     Tc = pn.find_neighbor_throats(Pc)
-    geom3 = OpenPNM.Geometry.GenericGeometry(network=pn,pores=Pc,throats=Tc)
-    Pb = sp.arange(9,18)
-    Tb = pn.find_neighbor_throats(pores=Pb,mode='intersection')
-    geom2 = OpenPNM.Geometry.GenericGeometry(network=pn,pores=Pb,throats=Tb)
+    geom3 = OpenPNM.Geometry.GenericGeometry(network=pn, pores=Pc, throats=Tc)
+    Pb = sp.arange(9, 18)
+    Tb = pn.find_neighbor_throats(pores=Pb, mode='intersection')
+    geom2 = OpenPNM.Geometry.GenericGeometry(network=pn, pores=Pb, throats=Tb)
     # Create an index in the Network
     pn['pore.num1'] = pn.Ps
     # Create the same index across each geom
@@ -288,62 +369,41 @@ def test_mapping():
     assert(sp.all(pn['pore.num1'] == pn['pore.num2']))
     # Send junk pores to ensure error is raised
     with pytest.raises(Exception):
-        pn.map_pores(pores=[0,pn.Np-1],target=geom1)
-        pn.map_pores(pores=[0,pn.Np+1],target=geom1)
-        pn.map_pores(pores=[pn.Np-1],target=geom1)
-        pn.map_pores(pores=[pn.Np+1],target=geom1)
-        geom1.map_pores(pores=[0,geom1.Np+1],target=pn)
-        geom1.map_pores(pores=[0,pn.Np+1],target=pn)
-        geom1.map_pores(pores=[geom1.Np+1],target=pn)
-        geom1.map_pores(pores=[pn.Np+1],target=pn)
-        geom2.map_pores(pores=[0],target=geom1)
-        geom2.map_pores(pores=[geom2.Np+1],target=geom1)
-        geom2.map_pores(pores=[0,geom2.Np-1],target=geom1)
-        geom2.map_pores(pores=[0,geom2.Np+1],target=geom1)
+        pn.map_pores(pores=[0, pn.Np-1], target=geom1)
+        pn.map_pores(pores=[0, pn.Np+1], target=geom1)
+        pn.map_pores(pores=[pn.Np-1], target=geom1)
+        pn.map_pores(pores=[pn.Np+1], target=geom1)
+        geom1.map_pores(pores=[0, geom1.Np+1], target=pn)
+        geom1.map_pores(pores=[0, pn.Np+1], target=pn)
+        geom1.map_pores(pores=[geom1.Np+1], target=pn)
+        geom1.map_pores(pores=[pn.Np+1], target=pn)
+        geom2.map_pores(pores=[0], target=geom1)
+        geom2.map_pores(pores=[geom2.Np+1], target=geom1)
+        geom2.map_pores(pores=[0, geom2.Np-1], target=geom1)
+        geom2.map_pores(pores=[0, geom2.Np+1], target=geom1)
     # Trim column from center of Network
-    pn.trim(pores=[4,13,22])
+    pn.trim(pores=[4, 13, 22])
     # Confirm index still match
     assert(sp.all(pn['pore.num1'] == pn['pore.num2']))
     # Check mapping between each Geometry object and in both directions
     # Check geom1
-    a = geom1.map_pores(pores=geom1.Ps,target=pn)
-    b = pn.map_pores(pores=a,target=geom1)
+    a = geom1.map_pores(pores=geom1.Ps, target=pn)
+    b = pn.map_pores(pores=a, target=geom1)
     assert(sp.all(b == geom1.Ps))
-    a = geom1.map_throats(throats=geom1.Ts,target=pn)
-    b = pn.map_throats(throats=a,target=geom1)
+    a = geom1.map_throats(throats=geom1.Ts, target=pn)
+    b = pn.map_throats(throats=a, target=geom1)
     assert(sp.all(b == geom1.Ts))
     # Check geom2
-    a = geom2.map_pores(pores=geom2.Ps,target=pn)
-    b = pn.map_pores(pores=a,target=geom2)
+    a = geom2.map_pores(pores=geom2.Ps, target=pn)
+    b = pn.map_pores(pores=a, target=geom2)
     assert(sp.all(b == geom2.Ps))
-    a = geom2.map_throats(throats=geom2.Ts,target=pn)
-    b = pn.map_throats(throats=a,target=geom2)
+    a = geom2.map_throats(throats=geom2.Ts, target=pn)
+    b = pn.map_throats(throats=a, target=geom2)
     assert(sp.all(b == geom2.Ts))
     # Check geom3
-    a = geom3.map_pores(pores=geom3.Ps,target=pn)
-    b = pn.map_pores(pores=a,target=geom3)
+    a = geom3.map_pores(pores=geom3.Ps, target=pn)
+    b = pn.map_pores(pores=a, target=geom3)
     assert(sp.all(b == geom3.Ps))
-    a = geom3.map_throats(throats=geom3.Ts,target=pn)
-    b = pn.map_throats(throats=a,target=geom3)
+    a = geom3.map_throats(throats=geom3.Ts, target=pn)
+    b = pn.map_throats(throats=a, target=geom3)
     assert(sp.all(b == geom3.Ts))
-
-def test_geometries():
-    ctrl = OpenPNM.Base.Controller()
-    pn = OpenPNM.Network.TestNet()
-    geom = OpenPNM.Geometry.Boundary(network=pn,pores=pn.Ps,throats=pn.Ts)
-    ctrl.clear()
-    pn = OpenPNM.Network.TestNet()
-    geom = OpenPNM.Geometry.Cube_and_Cuboid(network=pn,pores=pn.Ps,throats=pn.Ts)
-    ctrl.clear()
-    pn = OpenPNM.Network.TestNet()
-    geom = OpenPNM.Geometry.SGL10(network=pn,pores=pn.Ps,throats=pn.Ts)
-    ctrl.clear()
-    pn = OpenPNM.Network.TestNet()
-    geom = OpenPNM.Geometry.Stick_and_Ball(network=pn,pores=pn.Ps,throats=pn.Ts)
-    ctrl.clear()
-    pn = OpenPNM.Network.TestNet()
-    geom = OpenPNM.Geometry.Toray090(network=pn,pores=pn.Ps,throats=pn.Ts)
-    ctrl.clear()
-    pn = OpenPNM.Network.Delaunay(num_pores=20,domain_size=[1,1,1])
-    geom = OpenPNM.Geometry.Voronoi(network=pn,pores=pn.Ps,throats=pn.Ts)
-    ctrl.clear()
diff --git a/test/unit/Base/ControllerTest.py b/test/unit/Base/ControllerTest.py
index cf885b54f..23f93b1d5 100644
--- a/test/unit/Base/ControllerTest.py
+++ b/test/unit/Base/ControllerTest.py
@@ -1,27 +1,32 @@
+import OpenPNM
 from OpenPNM.Base import Controller
-from OpenPNM.Network import GenericNetwork
 
 
 class ControllerTest:
 
     def setup_class(self):
         self.controller = Controller()
-        self.pn = GenericNetwork(name='test_net')
-
-    def test_str(self):
-        actual_string = self.controller.__str__()
-        expected_string = \
-            '------------------------------------------------------------\n' + \
-            'Object:         Name                 (Class)\n' + \
-            '------------------------------------------------------------\n' + \
-            'Network:        test_net             (GenericNetwork)'
-
-        assert actual_string == expected_string
+        self.net = OpenPNM.Network.Cubic(shape=[10, 10, 10])
+        self.geo = OpenPNM.Geometry.TestGeometry(network=self.net,
+                                                 pores=self.net.Ps,
+                                                 throats=self.net.Ts)
 
     def test_get_log_level(self):
         self.controller.loglevel = 50
         assert self.controller.loglevel == 'Log level is currently set to: 50'
 
+    def test_save_and_load(self):
+        self.controller.save(self.net.name)
+        self.controller.clear()
+        assert self.controller == {}
+        self.controller.load(self.net.name)
+        assert self.net.name in self.controller.keys()
+        self.controller.save_simulation(self.net)
+        self.controller.purge_object(self.net, mode='complete')
+        self.controller.load_simulation(self.net.name)
+        self.controller.purge_object(self.net, mode='single')
+
     def teardown_class(self):
         del(self.controller)
-        del(self.pn)
+        del(self.net)
+        del(self.geo)
diff --git a/test/unit/Base/CoreTest.py b/test/unit/Base/CoreTest.py
index f4796f228..16ade506b 100644
--- a/test/unit/Base/CoreTest.py
+++ b/test/unit/Base/CoreTest.py
@@ -1,3 +1,312 @@
+import OpenPNM
+import scipy as sp
+
+
 class CoreTest:
-    def test_eq(self):
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[3, 3, 3])
+        self.geo = OpenPNM.Geometry.GenericGeometry(network=self.net,
+                                                    pores=self.net.Ps,
+                                                    throats=self.net.Ts)
+        self.geo['pore.diameter'] = sp.rand(self.net.Np)
+        self.geo.models.add(propname='pore.volume',
+                            model=OpenPNM.Geometry.models.pore_volume.sphere)
+        self.geo['throat.diameter'] = sp.rand(self.net.Nt)
+        self.geo.models.add(propname='throat.area',
+                            model=OpenPNM.Geometry.models.throat_area.cylinder)
+        self.geo['throat.label1'] = False
+        self.geo['throat.label2'] = False
+        self.geo['throat.label1'][0:6] = True
+        self.geo['throat.label2'][3:9] = True
+
+    def test_props_all(self):
+        a = self.geo.props()
+        assert sorted(a) == ['pore.diameter', 'pore.volume',
+                             'throat.area', 'throat.diameter']
+
+    def test_props_models(self):
+        a = self.geo.props(mode='models')
+        b = ['pore.volume', 'throat.area']
+        assert sorted(a) == sorted(b)
+
+    def test_props_constants(self):
+        a = self.geo.props(mode='constants')
+        b = ['pore.diameter', 'throat.diameter']
+        assert sorted(a) == sorted(b)
+
+    def test_props_pores_all(self):
+        a = self.geo.props(element='pores')
+        b = ['pore.diameter', 'pore.volume']
+        assert sorted(a) == sorted(b)
+
+    def test_props_pores_models(self):
+        a = self.geo.props(element='pores', mode='models')
+        b = ['pore.volume']
+        assert sorted(a) == sorted(b)
+
+    def test_props_pores_constants(self):
+        a = self.geo.props(element='pores', mode='constants')
+        b = ['pore.diameter']
+        assert sorted(a) == sorted(b)
+
+    def test_labels(self):
+        a = self.net.labels()
+        assert 'pore.top' in a
+
+    def test_labels_on_pores(self):
+        a = self.net.labels(element='pores')
+        b = ['pore.all', 'pore.back', 'pore.bottom', 'pore.front',
+             'pore.internal', 'pore.left', 'pore.right', 'pore.'+self.geo.name,
+             'pore.top']
+        assert sorted(a) == sorted(b)
+
+    def test_labels_on_throats(self):
+        a = self.net.labels(element='throats')
+        b = ['throat.all', 'throat.'+self.geo.name]
+        assert sorted(a) == sorted(b)
+
+    def test_labels_on_foo(self):
+        a = self.net.labels(element='foo')
+        assert a is None
+
+    def test_labels_on_all_pores(self):
+        a = self.net.labels(pores='all')
+        b = ['pore.all', 'pore.back', 'pore.bottom', 'pore.front',
+             'pore.internal', 'pore.left', 'pore.right', 'pore.'+self.geo.name,
+             'pore.top']
+        assert sorted(a) == sorted(b)
+
+    def test_labels_on_all_throats(self):
+        a = self.net.labels(throats='all')
+        b = ['throat.all', 'throat.'+self.geo.name]
+        assert sorted(a) == sorted(b)
+
+    def test_labels_on_one_pore(self):
+        a = self.net.labels(pores=0)
+        b = ['pore.all', 'pore.bottom', 'pore.front', 'pore.internal',
+             'pore.left', 'pore.'+self.geo.name]
+        assert sorted(a) == sorted(b)
+
+    def test_labels_on_list_of_pores(self):
+        a = self.net.labels(pores=[0, 1])
+        b = ['pore.all', 'pore.bottom', 'pore.front', 'pore.internal',
+             'pore.left', 'pore.'+self.geo.name]
+        assert sorted(a) == sorted(b)
+
+    def test_labels_locations_boolean(self):
+        ind = sp.zeros((self.net.Np), dtype=bool)
+        ind[[0, 1]] = True
+        a = self.net.labels(pores=ind)
+        b = ['pore.all', 'pore.bottom', 'pore.front', 'pore.internal',
+             'pore.left', 'pore.'+self.geo.name]
+        assert sorted(a) == sorted(b)
+
+    def test_labels_pores_mode_union(self):
+        a = self.net.labels(pores=[0, 1, 2], mode='union')
+        b = ['pore.all', 'pore.bottom', 'pore.front', 'pore.internal',
+             'pore.left', 'pore.'+self.geo.name, 'pore.top']
+        assert sorted(a) == sorted(b)
+
+    def test_labels_pores_mode_intersection(self):
+        a = self.net.labels(pores=[0, 1, 2], mode='intersection')
+        b = ['pore.all', 'pore.front', 'pore.internal', 'pore.left',
+             'pore.'+self.geo.name]
+        assert sorted(a) == sorted(b)
+
+    def test_labels_pores_mode_count(self):
+        a = self.net.labels(pores=[0, 1, 2], mode='count')
+        assert sp.all(a == [6, 5, 6])
+
+    def test_labels_pores_mode_mask(self):
+        a = self.net.labels(pores=[0, 1], mode='mask')
+        assert sp.sum(a) == 11
+
+    def test_labels_pores_mode_difference(self):
+        a = self.net.labels(pores=[0, 1, 2], mode='difference')
+        b = ['pore.back', 'pore.bottom', 'pore.right', 'pore.top']
+        assert sorted(a) == sorted(b)
+
+    def test_labels_pores_mode_none(self):
+        a = self.net.labels(pores=[0, 1], mode='none')
+        assert a[0] != a[1]
+
+    def test_labels_pores_mode_foo(self):
+        a = self.net.labels(pores=[0, 1], mode='foo')
+        assert a is None
+
+    def test_pores(self):
+        a = self.net.pores()
+        assert sp.all(a == sp.arange(0, self.net.Np))
+
+    def test_pores_one_label(self):
+        a = self.net.pores(labels='top')
+        assert sp.all(a == [2, 5, 8, 11, 14, 17, 20, 23, 26])
+
+    def test_pores_two_labels_union(self):
+        a = self.net.pores(labels=['top', 'front'], mode='union')
+        assert sp.all(a == [0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 14, 17, 20, 23, 26])
+
+    def test_pores_two_labels_intersection(self):
+        a = self.net.pores(labels=['top', 'front'], mode='intersection')
+        assert sp.all(a == [2, 5, 8])
+
+    def test_pores_two_labels_not_intersection(self):
+        a = self.net.pores(labels=['top', 'front'], mode='not_intersection')
+        assert sp.all(a == [0, 1, 3, 4, 6, 7, 11, 14, 17, 20, 23, 26])
+
+    def test_pores_two_labels_difference(self):
+        a = self.net.pores(labels=['top', 'front'], mode='difference')
+        assert sp.all(a == [9, 10, 12, 13, 15, 16, 18, 19, 21, 22, 24, 25])
+
+    def test_throats(self):
+        a = self.net.throats()
+        assert sp.all(a == sp.arange(0, self.net.Nt))
+
+    def test_throats_one_label(self):
+        a = self.net.throats(labels='label1')
+        assert sp.all(a == [0, 1, 2, 3, 4, 5])
+
+    def test_throats_two_labels_union(self):
+        a = self.net.throats(labels=['label1', 'label2'], mode='union')
+        assert sp.all(a == [0, 1, 2, 3, 4, 5, 6, 7, 8])
+
+    def test_throats_two_labels_intersection(self):
+        a = self.net.throats(labels=['label1', 'label2'], mode='intersection')
+        assert sp.all(a == [3, 4, 5])
+
+    def test_throats_two_labels_not_intersection(self):
+        a = self.net.throats(labels=['label1', 'label2'],
+                             mode='not_intersection')
+        assert sp.all(a == [0, 1, 2, 6, 7, 8])
+
+    def test_tomask_pores(self):
+        a = self.net.tomask(pores=self.net.pores('top'))
+        assert sp.sum(a) == 9
+
+    def test_tomask_throats(self):
+        a = self.net.tomask(throats=self.net.throats('label1'))
+        assert sp.sum(a) == 6
+
+    def test_toindices_pores(self):
+        mask = sp.zeros((self.net.Np), dtype=bool)
+        Ps = [0, 3, 6]
+        mask[Ps] = True
+        a = self.net.toindices(mask)
+        assert sp.all(a == Ps)
+
+    def test_toindices_throats(self):
+        mask = sp.zeros((self.net.Nt), dtype=bool)
+        Ts = [0, 3, 6]
+        mask[Ts] = True
+        a = self.net.toindices(mask)
+        assert sp.all(a == Ts)
+
+    def test_num_pores(self):
+        a = self.net.num_pores()
+        assert a == 27
+
+    def test_num_pores_one_label(self):
+        a = self.net.num_pores(labels='top')
+        assert a == 9
+
+    def test_num_pores_two_labels_union(self):
+        a = self.net.num_pores(labels=['top', 'front'], mode='union')
+        assert a == 15
+
+    def test_num_pores_two_labels_intersection(self):
+        a = self.net.num_pores(labels=['top', 'front'], mode='intersection')
+        assert a == 3
+
+    def test_num_pores_two_labels_notintersection(self):
+        a = self.net.num_pores(labels=['top', 'front'],
+                               mode='not_intersection')
+        assert a == 12
+
+    def test_num_pores_two_labels_difference(self):
+        a = self.net.num_pores(labels=['top', 'front'], mode='difference')
+        assert a == 12
+
+    def test_num_throats(self):
+        a = self.net.num_throats()
+        assert a == 54
+
+    def test_num_throats_one_label(self):
+        a = self.net.num_throats(labels='label1')
+        assert a == 6
+
+    def test_num_throats_two_labels_union(self):
+        a = self.net.num_throats(labels=['label1', 'label2'], mode='union')
+        assert a == 9
+
+    def test_num_throats_two_labels_intersection(self):
+        a = self.net.num_throats(labels=['label1', 'label2'],
+                                 mode='intersection')
+        assert a == 3
+
+    def test_num_throats_two_labels_notintersection(self):
+        a = self.net.num_throats(labels=['label1', 'label2'],
+                                 mode='not_intersection')
+        assert a == 6
+
+    def test_num_throats_two_labels_difference(self):
+        a = self.net.num_throats(labels=['label1', 'label2'],
+                                 mode='difference')
+        assert a == 45
+
+    def test_setitem_wrong_prefix(self):
+        self.geo['pore2.test'] = 0
+        assert 'pore2.test' not in self.geo.keys()
+
+    def test_setitem_wrong_length(self):
+        self.geo['pore.test'] = sp.ones((self.geo.Np+1))
+        assert 'pore.test' not in self.geo.keys()
+
+    def test_setitem_replace_all(self):
+        array_len = sp.size(self.geo['pore.all'])
+        self.geo['pore.all'] = sp.ones((self.geo.Np+1))
+        assert sp.size(self.geo['pore.all']) == array_len
+
+    def test_setitem_overwrite_into_all(self):
         pass
+        # This test will fail as there is currently no way to prevent this
+        # array_sum = sp.sum(self.geo['pore.all'])
+        # self.geo['pore.all'][0] = False
+        # assert sp.sum(self.geo['pore.all']) == array_sum
+
+    def test_object_name_name_conflict(self):
+        flag = False
+        try:
+            self.geo.name = self.net.name
+        except:
+            flag = True
+        assert flag
+
+    def test_object_name_array_conflict(self):
+        flag = False
+        try:
+            self.geo.name = 'coords'
+        except:
+            flag = True
+        assert flag
+
+        Np = self.geo.Np
+        Nt = self.geo.Nt
+        self.geo.clear()
+        assert self.geo.Np == Np
+        assert self.geo.Nt == Nt
+
+    def test_get_indices(self):
+        flag = False
+        temp = self.net.pop('pore.all')
+        try:
+            self.net._get_indices(element='pores', labels='blah')
+        except:
+            flag = True
+        assert flag
+        self.net.update({'pore.all': temp})
+
+    def test_get_indices_wildcard(self):
+        a = self.net._get_indices(element='pore', labels='ri*')
+        assert sp.all(a == [6, 7, 8, 15, 16, 17, 24, 25, 26])
+        b = self.net._get_indices(element='pore', labels='*ght')
+        assert sp.all(a == b)
diff --git a/test/unit/Base/ModelsDictTest.py b/test/unit/Base/ModelsDictTest.py
index 55b092491..bfb65e7b4 100644
--- a/test/unit/Base/ModelsDictTest.py
+++ b/test/unit/Base/ModelsDictTest.py
@@ -1,3 +1,4 @@
+import pytest
 from OpenPNM.Base import ModelsDict
 from OpenPNM.Base.__ModelsDict__ import ModelWrapper
 
@@ -28,3 +29,21 @@ class ModelsDictTest:
 
         def test_str_with_models(self):
             pass
+
+        def test_keys(self):
+            assert self.models_dict.keys() == []
+            self.models_dict['test'] = {'model': 'test content'}
+            assert self.models_dict.keys() == ['test']
+
+        def test_add_with_no_master(self):
+            with pytest.raises(Exception):
+                self.models_dict.add('test model', {})
+
+        def test_add_with_more_than_one_master(self):
+            pass
+
+        def test_remove(self):
+            pass
+
+        def test_find_master_with_more_than_one_master(self):
+            pass
diff --git a/test/unit/Geometry/GenericGeometryTest.py b/test/unit/Geometry/GenericGeometryTest.py
index 44f036fc8..da043584b 100644
--- a/test/unit/Geometry/GenericGeometryTest.py
+++ b/test/unit/Geometry/GenericGeometryTest.py
@@ -1,3 +1,48 @@
+import OpenPNM
+import scipy as sp
+
+
 class GenericGeometryTest:
-    def test_getitem(self):
-        pass
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+        P = self.net.pores('top', mode='not')
+        T = self.net.find_neighbor_throats(pores=P)
+        self.geo = OpenPNM.Geometry.GenericGeometry(network=self.net,
+                                                    pores=P,
+                                                    throats=T)
+        P = self.net.pores('top')
+        T = self.net.find_neighbor_throats(pores=P, mode='intersection')
+        self.geo2 = OpenPNM.Geometry.GenericGeometry(network=self.net,
+                                                     pores=P,
+                                                     throats=T)
+
+    def test_init_without_network(self):
+        standalone_geo = OpenPNM.Geometry.GenericGeometry()
+        assert standalone_geo.Np == 0
+        assert standalone_geo.Nt == 0
+
+    def test_getitem_coords_from_net(self):
+        a = self.geo['pore.coords']
+        b = self.net['pore.coords'][self.net.pores(self.geo.name)]
+        assert sp.all(a == b)
+
+    def test_getitem_geo_conns_includes_nans(self):
+        a = self.geo['throat.conns']
+        assert sp.all(len(a) == self.geo.Nt)
+        assert sp.amax(a) == (self.geo.Np - 1)
+        b = [pore for pore in a[:, 1] if pore is not sp.nan]
+        assert sp.size(b) > self.geo.Np - 25
+
+    def test_getitem_geo2_conns_has_no_nans(self):
+        # Because of the way geo2 was constructed it has not dangling throats
+        # to another Geometry, so should have no nans in its conns
+        a = self.geo2['throat.conns']
+        assert sp.all(len(a) == self.geo2.Nt)
+        assert sp.amax(a) == (self.geo2.Np - 1)
+        b = sp.where(~(a[:, 1] < sp.inf))[0]
+        assert sp.size(b) == 0
+
+    def test_getitem_own_pores_by_name(self):
+        a = self.geo.pores(self.geo.name)
+        b = self.geo.Ps
+        assert sp.all(a == b)
diff --git a/test/unit/Geometry/models/PoreSeedTest.py b/test/unit/Geometry/models/PoreSeedTest.py
index 02780f6d3..289cecd28 100644
--- a/test/unit/Geometry/models/PoreSeedTest.py
+++ b/test/unit/Geometry/models/PoreSeedTest.py
@@ -1,3 +1,27 @@
+import OpenPNM
+import scipy as sp
+
+
 class PoreSeedTest:
-    def test_perlin_noice(self):
-        pass
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+        self.geo = OpenPNM.Geometry.GenericGeometry(network=self.net,
+                                                    pores=self.net.Ps,
+                                                    throats=self.net.Ts)
+
+    def test_random(self):
+        f = OpenPNM.Geometry.models.pore_seed.random
+        self.geo.models.add(propname='pore.seed',
+                            model=f,
+                            seed=0,
+                            num_range=[0.1, 2])
+        assert sp.amax(self.geo['pore.seed']) > 1.9
+        assert sp.amin(self.geo['pore.seed']) > 0.1
+
+    def test_spatially_correlated(self):
+        f = OpenPNM.Geometry.models.pore_seed.spatially_correlated
+        self.geo.models.add(propname='pore.seed',
+                            model=f,
+                            weights=[2, 2, 2])
+        assert sp.amin(self.geo['pore.seed'] > 0)
+        assert sp.amax(self.geo['pore.seed'] < 1)
diff --git a/test/unit/Geometry/models/ThroatSeedTest.py b/test/unit/Geometry/models/ThroatSeedTest.py
new file mode 100644
index 000000000..6165c6415
--- /dev/null
+++ b/test/unit/Geometry/models/ThroatSeedTest.py
@@ -0,0 +1,33 @@
+import OpenPNM
+import scipy as sp
+
+
+class ThroatSeedTest:
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+        self.geo = OpenPNM.Geometry.GenericGeometry(network=self.net,
+                                                    pores=self.net.Ps,
+                                                    throats=self.net.Ts)
+        self.geo['pore.seed'] = sp.rand(self.net.Np)
+
+    def test_random(self):
+        self.geo.models.add(propname='throat.seed',
+                            model=OpenPNM.Geometry.models.throat_seed.random,
+                            seed=0,
+                            num_range=[0.1, 2])
+        assert sp.amax(self.geo['throat.seed']) > 1.9
+        assert sp.amin(self.geo['throat.seed']) > 0.1
+
+    def test_neighbor(self):
+        self.geo.models.add(propname='throat.seed_max',
+                            model=OpenPNM.Geometry.models.throat_seed.neighbor,
+                            mode='max')
+        self.geo.models.add(propname='throat.seed_min',
+                            model=OpenPNM.Geometry.models.throat_seed.neighbor,
+                            mode='min')
+        self.geo.models.add(propname='throat.seed_mean',
+                            model=OpenPNM.Geometry.models.throat_seed.neighbor,
+                            mode='mean')
+        assert sp.all(self.geo['throat.seed_min'] <= self.geo['throat.seed_max'])
+        assert sp.all(self.geo['throat.seed_min'] <= self.geo['throat.seed_mean'])
+        assert sp.all(self.geo['throat.seed_mean'] <= self.geo['throat.seed_max'])
diff --git a/test/unit/Network/CubicTest.py b/test/unit/Network/CubicTest.py
index 7b1b74618..2a8e6014b 100644
--- a/test/unit/Network/CubicTest.py
+++ b/test/unit/Network/CubicTest.py
@@ -1,3 +1,50 @@
+import OpenPNM
+import scipy as sp
+
+
 class CubicTest:
-    def test_label_surfaces(self):
-        pass
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[5, 5, 5], spacing=1)
+        self.net['pore.diameter'] = sp.rand(self.net.Np)
+
+    def test_connectivity_6(self):
+        net = OpenPNM.Network.Cubic(shape=[3, 3, 3], connectivity=6)
+        assert net.num_neighbors(pores=13) == 6
+
+    def test_connectivity_8(self):
+        net = OpenPNM.Network.Cubic(shape=[3, 3, 3], connectivity=8)
+        assert net.num_neighbors(pores=13) == 8
+
+    def test_connectivity_12(self):
+        net = OpenPNM.Network.Cubic(shape=[3, 3, 3], connectivity=12)
+        assert net.num_neighbors(pores=13) == 12
+
+    def test_connectivity_14(self):
+        net = OpenPNM.Network.Cubic(shape=[3, 3, 3], connectivity=14)
+        assert net.num_neighbors(pores=13) == 14
+
+    def test_connectivity_18(self):
+        net = OpenPNM.Network.Cubic(shape=[3, 3, 3], connectivity=18)
+        assert net.num_neighbors(pores=13) == 18
+
+    def test_connectivity_20(self):
+        net = OpenPNM.Network.Cubic(shape=[3, 3, 3], connectivity=20)
+        assert net.num_neighbors(pores=13) == 20
+
+    def test_connectivity_26(self):
+        net = OpenPNM.Network.Cubic(shape=[3, 3, 3], connectivity=26)
+        assert net.num_neighbors(pores=13) == 26
+
+    def test_domain_area(self):
+        A = self.net.domain_area(face=self.net.pores('top'))
+        assert sp.allclose(A, 25, rtol=0.1)
+
+    def test_domain_length(self):
+        L = self.net.domain_length(face_1=self.net.pores('top'),
+                                   face_2=self.net.pores('bottom'))
+        assert sp.allclose(L, 4, rtol=1e-02)
+
+    def test_from_array(self):
+        arr = sp.ones(self.net._shape, dtype=int)
+        self.net.fromarray(array=arr, propname='pore.test')
+        assert sp.sum(self.net['pore.test']) == self.net.Np
diff --git a/test/unit/Network/GenericNetworkTest.py b/test/unit/Network/GenericNetworkTest.py
index 04bd7ec06..f35c8f163 100644
--- a/test/unit/Network/GenericNetworkTest.py
+++ b/test/unit/Network/GenericNetworkTest.py
@@ -1,3 +1,178 @@
+import OpenPNM
+import scipy as sp
+
+
 class GenericNetworkTest:
-    def test_setitem(self):
-        pass
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[10, 10, 10])
+
+    def teardown_class(self):
+        ctrl = OpenPNM.Base.Controller()
+        ctrl.clear()
+
+    def test_find_connected_pores_numeric_not_flattend(self):
+        a = self.net.find_connected_pores(throats=[0, 1])
+        assert sp.all(a.flatten() == [0, 1, 1, 2])
+
+    def test_find_connected_pores_numeric_flattend(self):
+        a = self.net.find_connected_pores(throats=[0, 1], flatten=True)
+        assert sp.all(a == [0, 1, 2])
+
+    def test_find_connected_pores_boolean_flattend(self):
+        Tind = sp.zeros((self.net.Nt,), dtype=bool)
+        Tind[[0, 1]] = True
+        a = self.net.find_connected_pores(throats=Tind, flatten=True)
+        assert sp.all(a == [0, 1, 2])
+
+    def test_find_connected_pores_empty_flattend(self):
+        a = self.net.find_connected_pores(throats=[], flatten=True)
+        assert sp.shape(a) == (0, 2)
+
+    def test_find_neighbor_pores_numeric(self):
+        a = self.net.find_neighbor_pores(pores=[])
+        assert sp.size(a) == 0
+
+    def test_find_neighbor_pores_boolean(self):
+        Pind = sp.zeros((self.net.Np,), dtype=bool)
+        Pind[[0, 1]] = True
+        a = self.net.find_neighbor_pores(pores=Pind)
+        assert sp.all(a == [2, 10, 11, 100, 101])
+
+    def test_find_neighbor_pores_numeric_union(self):
+        a = self.net.find_neighbor_pores(pores=[0, 2],
+                                         mode='union')
+        assert sp.all(a == [1, 3, 10, 12, 100, 102])
+
+    def test_find_neighbor_pores_numeric_intersection(self):
+        a = self.net.find_neighbor_pores(pores=[0, 2],
+                                         mode='intersection')
+        assert sp.all(a == [1])
+
+    def test_find_neighbor_pores_numeric_notintersection(self):
+        a = self.net.find_neighbor_pores(pores=[0, 2],
+                                         mode='not_intersection')
+        assert sp.all(a == [3, 10, 12, 100, 102])
+
+    def test_find_neighbor_pores_numeric_union_incl_self(self):
+        a = self.net.find_neighbor_pores(pores=[0, 2],
+                                         mode='union',
+                                         excl_self=False)
+        assert sp.all(a == [0, 1, 2, 3, 10, 12, 100, 102])
+
+    def test_find_neighbor_pores_numeric_intersection_incl_self(self):
+        a = self.net.find_neighbor_pores(pores=[0, 2],
+                                         mode='intersection',
+                                         excl_self=False)
+        assert sp.all(a == [1])
+
+    def test_find_neighbor_pores_numeric_notintersection_incl_self(self):
+        a = self.net.find_neighbor_pores(pores=[0, 2],
+                                         mode='not_intersection',
+                                         excl_self=False)
+        assert sp.all(a == [0, 2, 3, 10, 12, 100, 102])
+
+    def test_find_neighbor_throats_empty(self):
+        a = self.net.find_neighbor_throats(pores=[])
+        assert sp.size(a) == 0
+
+    def test_find_neighbor_throats_boolean(self):
+        Pind = sp.zeros((self.net.Np,), dtype=bool)
+        Pind[[0, 1]] = True
+        a = self.net.find_neighbor_throats(pores=Pind)
+        assert sp.all(a == [0, 1, 900, 901, 1800, 1801])
+
+    def test_find_neighbor_throats_numeric_union(self):
+        a = self.net.find_neighbor_throats(pores=[0, 2], mode='union')
+        assert sp.all(a == [0, 1, 2, 900, 902, 1800, 1802])
+
+    def test_find_neighbor_throats_numeric_intersection(self):
+        a = self.net.find_neighbor_throats(pores=[0, 2], mode='intersection')
+        assert sp.size(a) == 0
+
+    def test_find_neighbor_throats_numeric_not_intersection(self):
+        a = self.net.find_neighbor_throats(pores=[0, 2],
+                                           mode='not_intersection')
+        assert sp.all(a == [0, 1, 2, 900, 902, 1800, 1802])
+
+    def test_num_neighbors_empty(self):
+        a = self.net.num_neighbors(pores=[])
+        assert sp.size(a) == 0
+
+    def test_num_neighbors_boolean(self):
+        Pind = sp.zeros((self.net.Np,), dtype=bool)
+        Pind[0] = True
+        a = self.net.num_neighbors(pores=Pind)
+        assert a == 3
+
+    def test_num_neighbors_numeric_flattened(self):
+        a = self.net.num_neighbors(pores=[0, 2], flatten=True)
+        assert a == 6
+        assert isinstance(a, int)
+
+    def test_num_neighbors_numeric_notflattened(self):
+        a = self.net.num_neighbors(pores=[0, 2], flatten=False)
+        assert sp.all(a == [3, 4])
+
+    def test_num_neighbors_single_pore_notflattened(self):
+        a = self.net.num_neighbors(pores=0, flatten=False)
+        assert sp.all(a == [3])
+        assert isinstance(a, sp.ndarray)
+
+    def test_find_interface_throats(self):
+        self.net['pore.domain1'] = False
+        self.net['pore.domain2'] = False
+        self.net['pore.domain3'] = False
+        self.net['pore.domain1'][[0, 1, 2]] = True
+        self.net['pore.domain2'][[100, 101, 102]] = True
+        self.net['pore.domain3'][900:] = True
+        a = self.net.find_interface_throats(labels=['domain1', 'domain2'])
+        assert sp.all(a == [1800, 1801, 1802])
+        a = self.net.find_interface_throats(labels=['domain1', 'domain3'])
+        assert sp.size(a) == 0
+
+    def test_check_network_health_healthy(self):
+        a = self.net.check_network_health()
+        assert len(list(a.values())) == 5
+        assert sp.size(list(a.values())) == 0
+
+    def test_check_network_isolated_pores(self):
+        net = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+        Ts = net.find_neighbor_throats(pores=0)
+        net.trim(throats=Ts)
+        a = net.check_network_health()
+        assert a['isolated_pores'] == 0
+        net.trim(a['trim_pores'])
+        a = net.check_network_health()
+        assert sp.size(list(a.values())) == 0
+
+    def test_check_network_health_duplicate_throat(self):
+        net = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+        P12 = net['throat.conns'][0]
+        net.extend(throat_conns=[P12])
+        a = net.check_network_health()
+        assert len(a['duplicate_throats'][0]) == 2
+
+    def test_check_network_health_triplicate_throats(self):
+        net = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+        P12 = net['throat.conns'][0]
+        net.extend(throat_conns=[P12])
+        net.extend(throat_conns=[P12])
+        a = net.check_network_health()
+        assert len(a['duplicate_throats'][0]) == 3
+
+    def test_check_network_health_multiple_duplicate_throats(self):
+        net = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+        P12 = net['throat.conns'][0]
+        net.extend(throat_conns=[P12])
+        P12 = net['throat.conns'][1]
+        net.extend(throat_conns=[P12])
+        a = net.check_network_health()
+        assert len(a['duplicate_throats'][1]) == 2
+
+    def test_check_network_health_bidirectional_throat(self):
+        net = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+        P12 = net['throat.conns'][0]
+        net['throat.conns'][0] = [P12[1], P12[0]]
+        a = net.check_network_health()
+        assert sp.size(a['bidirectional_throats']) == 1
+        assert sp.size(a['duplicate_throats']) == 0
diff --git a/test/unit/Phases/models/DensityTest.py b/test/unit/Phases/models/DensityTest.py
index 3059df672..69453a573 100644
--- a/test/unit/Phases/models/DensityTest.py
+++ b/test/unit/Phases/models/DensityTest.py
@@ -1,3 +1,35 @@
+import OpenPNM
+import scipy as sp
+
+
 class DensityTest:
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[3, 3, 3])
+        self.phase = OpenPNM.Phases.GenericPhase(network=self.net)
+        self.phase['pore.temperature'] = 298.0  # K
+        self.phase['pore.pressure'] = 101325.0  # Pa
+        self.phase['pore.molecular_weight'] = 0.018  # kg/mol
+        self.phase['pore.molar_density'] = 55539.0  # mol/m3
+        self.phase['pore.salinity'] = 0.0  # ppt
+
     def test_standard(self):
-        pass
+        # Liquid water
+        self.phase.models.add(propname='pore.density',
+                              model=OpenPNM.Phases.models.density.standard)
+        assert sp.allclose(self.phase['pore.density'], 999.702)
+
+    def test_ideal_gas(self):
+        # Water vapor
+        self.phase.models.add(propname='pore.density',
+                              model=OpenPNM.Phases.models.density.ideal_gas)
+        assert sp.allclose(self.phase['pore.density'], 0.73610248)
+
+    def test_water(self):
+        # Liquid water
+        self.phase.models.add(propname='pore.density',
+                              model=OpenPNM.Phases.models.density.water)
+        assert sp.allclose(self.phase['pore.density'], 996.9522)
+
+    def teardown_class(self):
+        del(self.phase)
+        del(self.net)
diff --git a/test/unit/Phases/models/DiffusivityTest.py b/test/unit/Phases/models/DiffusivityTest.py
index 502c0cd59..624a1b270 100644
--- a/test/unit/Phases/models/DiffusivityTest.py
+++ b/test/unit/Phases/models/DiffusivityTest.py
@@ -1,3 +1,49 @@
+import OpenPNM
+import scipy as sp
+
+
 class DiffusivityTest:
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[3, 3, 3])
+        self.phase = OpenPNM.Phases.GenericPhase(network=self.net)
+        self.phase['pore.temperature'] = 298.0  # K
+        self.phase['pore.pressure'] = 101325  # Pa
+        self.phase['pore.viscosity'] = 1.75e-5  # Pa.s
+
     def test_fuller(self):
-        pass
+        f = OpenPNM.Phases.models.diffusivity.fuller
+        self.phase.models.add(propname='pore.diffusivity',
+                              model=f,
+                              MA=0.032,
+                              MB=0.028,
+                              vA=16.6,
+                              vB=17.9)
+        assert sp.allclose(self.phase['pore.diffusivity'], 2.06754784e-05)
+
+    def test_fuller_scaling(self):
+        f = OpenPNM.Phases.models.diffusivity.fuller_scaling
+        self.phase.models.add(propname='pore.diffusivity',
+                              model=f,
+                              DABo=1.79712526e-05,
+                              Po=100000,
+                              To=273)
+        assert sp.allclose(self.phase['pore.diffusivity'], 2.06754784e-05)
+
+    def test_tyn_calus(self):
+        f = OpenPNM.Phases.models.diffusivity.tyn_calus
+        self.phase.models.add(propname='pore.diffusivity',
+                              model=f,
+                              VA=16.5,
+                              VB=17.9,
+                              sigma_A=1,
+                              sigma_B=1)
+        assert sp.allclose(self.phase['pore.diffusivity'], 9.84851806e-05)
+
+    def test_tyn_calus_scaling(self):
+        f = OpenPNM.Phases.models.diffusivity.tyn_calus_scaling
+        self.phase.models.add(propname='pore.diffusivity',
+                              model=f,
+                              DABo=5.26300839e-05,
+                              mu_o=3e-5,
+                              To=273)
+        assert sp.allclose(self.phase['pore.diffusivity'], 9.84851806e-05)
diff --git a/test/unit/Phases/models/ElectricalConductivityTest.py b/test/unit/Phases/models/ElectricalConductivityTest.py
new file mode 100644
index 000000000..ec756a4b9
--- /dev/null
+++ b/test/unit/Phases/models/ElectricalConductivityTest.py
@@ -0,0 +1,16 @@
+import OpenPNM
+
+
+class ElectricalConductivityTest:
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[3, 3, 3])
+        self.phase = OpenPNM.Phases.GenericPhase(network=self.net)
+        self.phase['pore.intrinsic_conductivity'] = 1
+        self.phase['pore.volume_fraction'] = 0.5
+
+    def test_percolating_continua(self):
+        f = OpenPNM.Phases.models.electrical_conductivity.percolating_continua
+        self.phase.models.add(propname='pore.effective_conductivity',
+                              model=f,
+                              phi_crit=0.25,
+                              tau=2)
diff --git a/test/unit/Phases/models/MolarDensityTest.py b/test/unit/Phases/models/MolarDensityTest.py
index 76767de24..bf7a3b5e4 100644
--- a/test/unit/Phases/models/MolarDensityTest.py
+++ b/test/unit/Phases/models/MolarDensityTest.py
@@ -1,3 +1,32 @@
+import OpenPNM
+import scipy as sp
+
+
 class MolarDensityTest:
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[3, 3, 3])
+        self.phase = OpenPNM.Phases.GenericPhase(network=self.net)
+        self.phase['pore.molecular_weight'] = 0.0291  # kg/mol
+        self.phase['pore.density'] = 1.19  # kg/m3
+        self.phase['pore.temperature'] = 298.0  # K
+        self.phase['pore.pressure'] = 101325  # Pa
+        self.phase['pore.critical_temperature'] = 132.65  # K
+        self.phase['pore.critical_pressure'] = 3771000.0  # Pa
+
     def test_standard(self):
-        pass
+        f = OpenPNM.Phases.models.molar_density.standard
+        self.phase.models.add(propname='pore.molar_density',
+                              model=f)
+        assert sp.allclose(self.phase['pore.molar_density'], 40.8934707)
+
+    def test_ideal_gas(self):
+        f = OpenPNM.Phases.models.molar_density.ideal_gas
+        self.phase.models.add(propname='pore.molar_density',
+                              model=f)
+        assert sp.allclose(self.phase['pore.molar_density'], 40.8945824)
+
+    def test_vanderwaals(self):
+        f = OpenPNM.Phases.models.molar_density.vanderwaals
+        self.phase.models.add(propname='pore.molar_density',
+                              model=f)
+        assert sp.allclose(self.phase['pore.molar_density'], 40.92524916)
diff --git a/test/unit/Phases/models/SurfaceTensionTest.py b/test/unit/Phases/models/SurfaceTensionTest.py
index e963a7b9c..57b5f6a79 100644
--- a/test/unit/Phases/models/SurfaceTensionTest.py
+++ b/test/unit/Phases/models/SurfaceTensionTest.py
@@ -1,3 +1,43 @@
+import OpenPNM
+import scipy as sp
+
+
 class SurfaceTensionTest:
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[3, 3, 3])
+        self.phase = OpenPNM.Phases.GenericPhase(network=self.net)
+        self.phase['pore.temperature'] = 298.0  # K
+        self.phase['pore.molecular_weight'] = 0.018  # kg/mol
+        self.phase['pore.critical_temperature'] = 647.15  # K
+        self.phase['pore.critical_pressure'] = 3771000.0  # Pa
+        self.phase['pore.salinity'] = 0  # g/kg
+        self.phase['pore.molar_density'] = 55.5  # mol/m3
+
     def test_water(self):
-        pass
+        f = OpenPNM.Phases.models.surface_tension.water
+        self.phase.models.add(propname='pore.surface_tension',
+                              model=f)
+        assert sp.allclose(self.phase['pore.surface_tension'], 0.07199533)
+
+    def test_eotvos(self):
+        f = OpenPNM.Phases.models.surface_tension.eotvos
+        self.phase.models.add(propname='pore.surface_tension',
+                              model=f,
+                              k=0.000014)
+        assert sp.allclose(self.phase['pore.surface_tension'], 0.07112169)
+
+    def test_guggenheim_katayama(self):
+        f = OpenPNM.Phases.models.surface_tension.guggenheim_katayama
+        self.phase.models.add(propname='pore.surface_tension',
+                              model=f,
+                              K2=0.0000014,
+                              n=0.1)
+        assert sp.allclose(self.phase['pore.surface_tension'], 0.27582571)
+
+    def test_brock_bird_scaling(self):
+        f = OpenPNM.Phases.models.surface_tension.brock_bird_scaling
+        self.phase.models.add(propname='pore.surface_tension',
+                              model=f,
+                              sigma_o=0.0608,
+                              To=363)
+        assert sp.allclose(self.phase['pore.surface_tension'], 0.07820761)
diff --git a/test/unit/Phases/models/ThermalConductivityTest.py b/test/unit/Phases/models/ThermalConductivityTest.py
index e8df5a57b..7f921da2e 100644
--- a/test/unit/Phases/models/ThermalConductivityTest.py
+++ b/test/unit/Phases/models/ThermalConductivityTest.py
@@ -1,3 +1,34 @@
+import OpenPNM
+import scipy as sp
+
+
 class ThermalConductivityTest:
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[3, 3, 3])
+        self.phase = OpenPNM.Phases.GenericPhase(network=self.net)
+        self.phase['pore.temperature'] = 298.0  # K
+        self.phase['pore.salinity'] = 0.0  # g/kg
+        self.phase['pore.viscosity'] = 0.001  # Pa.s
+        self.phase['pore.critical_temperature'] = 647.15  # K
+        self.phase['pore.molecular_weight'] = 0.018  # kg/mol
+        self.phase['pore.boiling_point'] = 373.15  # K
+        self.phase['pore.heat_capacity'] = 75.28  # J/mol K
+        self.phase['pore.acentric_factor'] = 11.5  # J/mol K
+
     def test_water(self):
-        pass
+        f = OpenPNM.Phases.models.thermal_conductivity.water
+        self.phase.models.add(propname='pore.thermal_conductivity',
+                              model=f)
+        assert sp.allclose(self.phase['pore.thermal_conductivity'], 0.61047611)
+
+    def test_chung(self):
+        f = OpenPNM.Phases.models.thermal_conductivity.chung
+        self.phase.models.add(propname='pore.thermal_conductivity',
+                              model=f)
+        assert sp.allclose(self.phase['pore.thermal_conductivity'], 0.62063913)
+
+    def test_sato(self):
+        f = OpenPNM.Phases.models.thermal_conductivity.sato
+        self.phase.models.add(propname='pore.thermal_conductivity',
+                              model=f)
+        assert sp.allclose(self.phase['pore.thermal_conductivity'], 0.29787023)
diff --git a/test/unit/Phases/models/VaporPressureTest.py b/test/unit/Phases/models/VaporPressureTest.py
index 4dfa71680..d56e6bde6 100644
--- a/test/unit/Phases/models/VaporPressureTest.py
+++ b/test/unit/Phases/models/VaporPressureTest.py
@@ -1,3 +1,28 @@
+import OpenPNM
+import scipy as sp
+
+
 class VaporPressureTest:
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[3, 3, 3])
+        self.phase = OpenPNM.Phases.GenericPhase(network=self.net)
+        self.phase['pore.temperature'] = 300*sp.ones(self.phase.Np,)
+        self.phase['pore.salinity'] = sp.zeros((self.phase.Np,))
+
     def test_antoine(self):
-        pass
+        f = OpenPNM.Phases.models.vapor_pressure.antoine
+        self.phase.models.add(propname='pore.test',
+                              model=f,
+                              pore_temperature='pore.temperature',
+                              A=8.07131,
+                              B=1730.63,
+                              C=233.426)
+        assert sp.allclose(self.phase['pore.test'], 3523.72641773)
+
+    def test_water(self):
+        f = OpenPNM.Phases.models.vapor_pressure.water
+        self.phase.models.add(propname='pore.test',
+                              model=f,
+                              pore_temperature='pore.temperature',
+                              pore_salinity='pore.salinity')
+        assert sp.allclose(self.phase['pore.test'], 3536.01)
diff --git a/test/unit/Phases/models/ViscosityTest.py b/test/unit/Phases/models/ViscosityTest.py
index 96414ff4c..72ac7b59a 100644
--- a/test/unit/Phases/models/ViscosityTest.py
+++ b/test/unit/Phases/models/ViscosityTest.py
@@ -1,3 +1,30 @@
+import OpenPNM
+import scipy as sp
+
+
 class ViscosityTest:
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[3, 3, 3])
+        self.phase = OpenPNM.Phases.GenericPhase(network=self.net)
+        self.phase['pore.temperature'] = 298.0  # K
+        self.phase['pore.molecular_weight'] = 0.018  # kg/mol
+        self.phase['pore.critical_temperature'] = 647.15  # K
+        self.phase['pore.critical_volume'] = 1.805e-5  # m3/kmol
+        self.phase['pore.salinity'] = 0  # g/kg
+
     def test_water(self):
-        pass
+        self.phase.models.add(propname='pore.viscosity',
+                              model=OpenPNM.Phases.models.viscosity.water)
+        assert sp.allclose(self.phase['pore.viscosity'], 0.00089319)
+
+    def test_reynolds(self):
+        self.phase.models.add(propname='pore.viscosity',
+                              model=OpenPNM.Phases.models.viscosity.reynolds,
+                              uo=0.001,
+                              b=0.001)
+        assert sp.allclose(self.phase['pore.viscosity'], 0.00134716)
+
+    def test_chung(self):
+        self.phase.models.add(propname='pore.viscosity',
+                              model=OpenPNM.Phases.models.viscosity.chung)
+        assert sp.allclose(self.phase['pore.viscosity'], 6.47289919e-05)
diff --git a/test/unit/Physics/GenericPhysicsTest.py b/test/unit/Physics/GenericPhysicsTest.py
index 4e04b1614..76e1523a0 100644
--- a/test/unit/Physics/GenericPhysicsTest.py
+++ b/test/unit/Physics/GenericPhysicsTest.py
@@ -1,3 +1,19 @@
+import OpenPNM
+
+
 class GenericPhysicsTest:
-    def test_setitem(self):
-        pass
+    def setup_class(self):
+        self.net = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+        self.net = OpenPNM.Network.Cubic(shape=[5, 5, 5])
+        P1 = self.net.pores('top', mode='not')
+        T1 = self.net.find_neighbor_throats(pores=P1)
+        self.geo1 = OpenPNM.Geometry.GenericGeometry(network=self.net,
+                                                     pores=P1,
+                                                     throats=T1)
+        P2 = self.net.pores('top')
+        T2 = self.net.find_neighbor_throats(pores=P2, mode='intersection')
+        self.geo2 = OpenPNM.Geometry.GenericGeometry(network=self.net,
+                                                     pores=P2,
+                                                     throats=T2)
+        self.phase1 = OpenPNM.Phases.GenericPhase(network=self.net)
+        self.phase2 = OpenPNM.Phases.GenericPhase(network=self.net)
diff --git a/test/unit/Physics/models/CapillaryPressureTest.py b/test/unit/Physics/models/CapillaryPressureTest.py
index 0458495f3..6ec0d157d 100644
--- a/test/unit/Physics/models/CapillaryPressureTest.py
+++ b/test/unit/Physics/models/CapillaryPressureTest.py
@@ -1,3 +1,12 @@
+import OpenPNM
+
+
 class CapillaryPressureTest:
+
+    def setup_class(self):
+        self.pn = OpenPNM.Network.TestNet()
+        self.phase = OpenPNM.Phases.TestPhase()
+        self.physics = OpenPNM.Physics.TestPhysics()
+
     def test_washburn(self):
-        self
+        pass