Repository: PMEAL/OpenPNM
Pull request #827 (2018-03-13 02:06:37)

Generalizing misc models
This PR re-arranges the models module to have a misc submodule, and moved all the pore_misc, throat_misc, etc, files there


Patch:
diff --git a/openpnm/core/ModelsMixin.py b/openpnm/core/ModelsMixin.py
index cb04a35e4..50a582d76 100644
--- a/openpnm/core/ModelsMixin.py
+++ b/openpnm/core/ModelsMixin.py
@@ -27,6 +27,50 @@ class ModelsDict(PrintableDict):
 class ModelsMixin():
 
     def add_model(self, propname, model, regen_mode='deferred', **kwargs):
+        r"""
+        Adds a new model to the models dictionary (``object.models``)
+
+        Parameters
+        ----------
+        propname : string
+            The name of the property to be calculated by the model.
+
+        model : function
+            A reference (handle) to the function to be used.
+
+        regen_mode : string
+            Controls how/when the model is run (See Notes for more details).
+            Options are:
+
+            *'normal'* or *'eager'* : The model is run directly upon being
+            assiged, and also run every time ``object.regenerate_models`` is
+            called.
+
+            *'constant'* : The model is run directly upon being assigned, but
+            is not called again, thus making it's data act like a constant.
+
+            *'deferred'* or *'lazy'*: (default) Is not run upon being assigned,
+            but is run the first time that it's data is requested.
+
+        Notes
+        -----
+        The difference between 'eager' and 'lazy' execution is useful to
+        understand.  In 'eager' mode the model is run as soon as it's assigned
+        to the object, which means that users must be careful to assign
+        dependent models first, or else a KeyError will be raised since needed
+        data is not present.  In 'lazy' mode models are not run until their
+        data is asked for, which creates a cascading call to all necessary
+        models that have not been run yet.  The latter behavior is more
+        convenient, but can be confusing since a lot happens behind the scenes.
+
+        For example: In 'eager' mode, if 'pore.volume' is assigned before
+        'pore.diameter' a KeyError will occur when the model attempts to
+        access the 'pore.diameter' data.  In 'lazy' mode the request for
+        `pore.volume` data will run the 'pore.volume' model, which attempts
+        to access `pore.diameter` data, and upon failing to find it will
+        attempt to run the `pore.diameter` model.
+
+        """
         # Add model and regen_mode to kwargs dictionary
         kwargs.update({'model': model, 'regen_mode': regen_mode})
         # Insepct model to extract arguments and default values
@@ -41,10 +85,55 @@ class ModelsMixin():
         # Store all keyword argumnents in model
         self.models[propname] = kwargs
         # Regenerate model values if necessary
-        if regen_mode != 'deferred':
+        if regen_mode not in ['deferred', 'lazy']:
             self._regen(propname)
 
+    def remove_model(self, propname=None, mode=['model', 'data']):
+        r"""
+        Removes model and data from object.
+
+        Parameters
+        ----------
+        propname : string or list of strings
+            The property or list of properties to remove
+
+        mode : list of strings
+            Controls what is removed.  Options are:
+
+            *'model'* : Removes the model but not any numerical data that may
+            already exist.
+
+            *'data'* : Removes the data but leaves the model.
+
+        The default is both.
+
+        """
+        if type(propname) is str:
+            propname = [propname]
+        for item in propname:
+            if 'model' in mode:
+                if item in self.models.keys():
+                    del self.models[item]
+            if 'data' in mode:
+                if item in self.keys():
+                    del self[item]
+
     def regenerate_models(self, propnames=None, exclude=[]):
+        r"""
+        Re-runs the specified models.
+
+        Parameters
+        ----------
+        propnames : string or list of strings
+            The list of property names to be regenerated.  If None are given
+            then ALL models are re-run (except for those whose ``regen_mode``
+            is 'constant').
+
+        exclude : list of strings
+            Since the default behavior is to run ALL models, this can be used
+            to exclude specific models.  It may be more convenient to supply
+            as list of 2 models to exclude than to specify 8 models include.
+        """
         # If no props given, then regenerate them all
         if propnames is None:
             propnames = list(self.models.keys())
diff --git a/openpnm/geometry/StickAndBall.py b/openpnm/geometry/StickAndBall.py
index bdfc17391..52ad1dacf 100644
--- a/openpnm/geometry/StickAndBall.py
+++ b/openpnm/geometry/StickAndBall.py
@@ -1,5 +1,4 @@
-# -*- coding: utf-8 -*-
-import openpnm.models.geometry as gm
+import openpnm.models as mods
 from openpnm.geometry import GenericGeometry
 
 
@@ -26,52 +25,52 @@ class StickAndBall(GenericGeometry):
         super().__init__(**kwargs)
 
         self.add_model(propname='pore.seed',
-                       model=gm.pore_misc.random,
+                       model=mods.misc.random,
+                       element='pore',
                        num_range=[0, 0.1],
                        seed=None)
 
         self.add_model(propname='pore.max_size',
-                       model=gm.pore_size.largest_sphere,
+                       model=mods.geometry.pore_size.largest_sphere,
                        iters=10)
 
         self.add_model(propname='pore.area',
-                       model=gm.pore_area.spherical,
+                       model=mods.geometry.pore_area.sphere,
                        pore_diameter='pore.diameter')
 
         self.add_model(propname='pore.volume',
-                       model=gm.pore_volume.sphere,
+                       model=mods.geometry.pore_volume.sphere,
                        pore_diameter='pore.diameter')
 
         self.add_model(propname='pore.diameter',
-                       model=gm.misc.product,
-                       arg1='pore.max_size',
-                       arg2='pore.seed')
+                       model=mods.misc.product,
+                       props=['pore.max_size', 'pore.seed'])
 
         self.add_model(propname='throat.max_size',
-                       model=gm.throat_misc.neighbor,
+                       model=mods.misc.from_neighbor_pores,
                        mode='min',
                        pore_prop='pore.diameter')
 
         self.add_model(propname='throat.diameter',
-                       model=gm.misc.scaled,
+                       model=mods.misc.scaled,
                        factor=0.5,
                        prop='throat.max_size')
 
         self.add_model(propname='throat.length',
-                       model=gm.throat_length.straight,
+                       model=mods.geometry.throat_length.straight,
                        L_negative=1e-09,
                        pore_diameter='pore.diameter')
 
         self.add_model(propname='throat.surface_area',
-                       model=gm.throat_surface_area.cylinder,
+                       model=mods.geometry.throat_surface_area.cylinder,
                        throat_diameter='throat.diameter',
                        throat_length='throat.length')
 
         self.add_model(propname='throat.volume',
-                       model=gm.throat_volume.cylinder,
+                       model=mods.geometry.throat_volume.cylinder,
                        throat_diameter='throat.diameter',
                        throat_length='throat.length')
 
         self.add_model(propname='throat.area',
-                       model=gm.throat_area.cylinder,
+                       model=mods.geometry.throat_area.cylinder,
                        throat_diameter='throat.diameter')
diff --git a/openpnm/io/VTK.py b/openpnm/io/VTK.py
index 2f96bda98..24b7d6136 100644
--- a/openpnm/io/VTK.py
+++ b/openpnm/io/VTK.py
@@ -101,13 +101,13 @@ class VTK(GenericIO):
             if element is not None:
                 cell_data_node.append(element)
 
+        tree = ET.ElementTree(root)
+        tree.write(filename)
+
         if filename == '':
             filename = project.name
         filename = cls._parse_filename(filename=filename, ext='vtp')
 
-        tree = ET.ElementTree(root)
-        tree.write(filename)
-
         with open(filename, 'r+') as f:
             string = f.read()
             string = string.replace('</DataArray>', '</DataArray>\n\t\t\t')
diff --git a/openpnm/models/__init__.py b/openpnm/models/__init__.py
index 361856d1f..9d3675702 100644
--- a/openpnm/models/__init__.py
+++ b/openpnm/models/__init__.py
@@ -1,3 +1,4 @@
+from . import misc
 from . import geometry
 from . import phase
 from . import physics
diff --git a/openpnm/models/geometry/__init__.py b/openpnm/models/geometry/__init__.py
index 2ef48a69d..04c3cd696 100644
--- a/openpnm/models/geometry/__init__.py
+++ b/openpnm/models/geometry/__init__.py
@@ -1,5 +1,4 @@
 from . import pore_size
-from . import pore_misc
 from . import pore_seed
 from . import pore_volume
 from . import pore_surface_area
@@ -7,9 +6,7 @@ from . import pore_area
 from . import throat_area
 from . import throat_size
 from . import throat_length
-from . import throat_misc
 from . import throat_perimeter
-from . import throat_seed
 from . import throat_surface_area
 from . import throat_volume
 from . import throat_shape_factor
diff --git a/openpnm/models/geometry/pore_area.py b/openpnm/models/geometry/pore_area.py
index f7e35b67f..68675ba9a 100644
--- a/openpnm/models/geometry/pore_area.py
+++ b/openpnm/models/geometry/pore_area.py
@@ -1,13 +1,7 @@
-r"""
-===============================================================================
-pore_area -- Models for cross-sectional area of a pore body
-===============================================================================
+from numpy import pi as _pi
 
-"""
-import scipy as _sp
 
-
-def spherical(target, pore_diameter='pore.diameter'):
+def sphere(target, pore_diameter='pore.diameter'):
     r"""
     Calculate cross-sectional area assuming the pore body is a sphere
 
@@ -24,11 +18,11 @@ def spherical(target, pore_diameter='pore.diameter'):
 
     """
     diams = target[pore_diameter]
-    value = _sp.pi/4*(diams)**2
+    value = _pi/4*(diams)**2
     return value
 
 
-def cubic(target, pore_diameter='pore.diameter'):
+def cube(target, pore_diameter='pore.diameter'):
     r"""
     Calculate cross-sectional area assuming the pore body is a cube
 
diff --git a/openpnm/models/geometry/pore_misc.py b/openpnm/models/geometry/pore_misc.py
deleted file mode 100644
index a71b6bb20..000000000
--- a/openpnm/models/geometry/pore_misc.py
+++ /dev/null
@@ -1,105 +0,0 @@
-r"""
-===============================================================================
-pore_misc -- miscillaneous and generic functions to apply to pores
-===============================================================================
-
-"""
-from . import misc as _misc
-import scipy as _sp
-
-
-def constant(target, value):
-    r"""
-    Assign specified constant value to all pores on the specifed Geometry
-    object.
-
-    Parameters
-    ----------
-    target : OpenPNM Object
-        The object which this model is associated with. This controls the
-        length of the calculated array, and also provides access to other
-        necessary properties.
-
-    value : scalar
-        The scalar value to be applied
-
-    Notes
-    -----
-    This function is redundant and could be accomplished with
-    geometry['pore.prop'] = value.
-    """
-    if _sp.shape(value):
-        raise Exception('Only scalar values can be assigned')
-    value = _sp.ones(target.num_pores(),)*value
-    return value
-
-
-def random(target, seed=None, num_range=[0, 1]):
-    r"""
-    Generate an array of random numbers using the Scipy random number
-    generator.
-
-    Parameters
-    ----------
-    target : OpenPNM Object
-        The object which this model is associated with. This controls the
-        length of the calculated array, and also provides access to other
-        necessary properties.
-
-    seed : int
-        The seed value to initialize the random number generator.  The default
-        is None which means different numbers will be produced each time this
-        function is run.
-
-    num_range : list
-        The min and max values of the returned random numbers.  The values will
-        be uniformly distributed over this range.
-
-    Notes
-    -----
-    This method seems trivial but it can be useful to add this as a pore-scale
-    model so that random number can be regenerated upon demand to create
-    new realization of the pore and throat size distributions (i.e. assuming
-    ``seed`` is not set to None).
-    """
-    values = _misc.random(target, element='pore', seed=seed,
-                          num_range=num_range)
-    return values
-
-
-def neighbor(target, throat_prop='throat.seed', mode='min'):
-    r"""
-    Adopt a value from the values found in neighboring throats
-
-    Parameters
-    ----------
-    target : OpenPNM Object
-        The object which this model is associated with. This controls the
-        length of the calculated array, and also provides access to other
-        necessary properties.
-
-    throat_prop : string
-        The dictionary key of the array containing the throat property to be
-        used in the calculation.  The default is 'throat.seed'.
-
-    mode : string
-        Controls how the pore property is calculated.  Options are 'min',
-        'max' and 'mean'.
-    """
-    network = target.project.network
-    Ps = target.pores()
-    data = target[throat_prop]
-    neighborTs = network.find_neighbor_throats(pores=Ps,
-                                               flatten=False,
-                                               mode='intersection')
-    values = _sp.ones((_sp.shape(Ps)[0],))*_sp.nan
-    if mode == 'min':
-        for pore in Ps:
-            values[pore] = _sp.amin(data[neighborTs[pore]])
-    if mode == 'max':
-        for pore in Ps:
-            values[pore] = _sp.amax(data[neighborTs[pore]])
-    if mode == 'mean':
-        for pore in Ps:
-            values[pore] = _sp.mean(data[neighborTs[pore]])
-    return values
diff --git a/openpnm/models/geometry/pore_seed.py b/openpnm/models/geometry/pore_seed.py
index 08ad7ca27..2c0ddd249 100644
--- a/openpnm/models/geometry/pore_seed.py
+++ b/openpnm/models/geometry/pore_seed.py
@@ -1,15 +1,7 @@
-# -*- coding: utf-8 -*-
-r"""
-===============================================================================
-pore_seed -- Methods for generating fields of values for use as seeds in
-statistical pore size distributions
-===============================================================================
-
-"""
-from . import misc as _misc
 import scipy as _sp
-from openpnm.core import logging
-_logger = logging.getLogger(__name__)
+from openpnm.models import misc as _misc
+from openpnm.core import logging as _logging
+_logger = _logging.getLogger(__name__)
 
 
 def random(target, seed=None, num_range=[0, 1]):
diff --git a/openpnm/models/geometry/pore_size.py b/openpnm/models/geometry/pore_size.py
index 653d9cf01..a25d2ea64 100644
--- a/openpnm/models/geometry/pore_size.py
+++ b/openpnm/models/geometry/pore_size.py
@@ -1,13 +1,7 @@
-r"""
-===============================================================================
-pore_diameter
-===============================================================================
-
-"""
-from openpnm.core import logging
-from . import misc as _misc
-import scipy as _sp
-_logger = logging.getLogger(__name__)
+from openpnm.core import logging as _logging
+from openpnm.models import misc as _misc
+import numpy as _np
+_logger = _logging.getLogger(__name__)
 
 
 def weibull(target, shape, scale, loc, seeds='pore.seed'):
@@ -34,34 +28,16 @@ generic.__doc__ = _misc.generic.__doc__
 
 
 def random(target, seed=None, num_range=[0, 1]):
-    r"""
-    Assign pore sizes from a random distribution
-
-    Parameters
-    ----------
-    target : OpenPNM Object
-        The object which this model is associated with. This controls the
-        length of the calculated array, and also provides access to other
-        necessary properties.
+    return _misc.random(target, element='pore', seed=seed, num_range=num_range)
 
-    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.  The default is [0, 1] but this can be adjusted to produce
-        pore sizes directly; for instance pores between 10 and 100 um can be
-        generated with ``num_range = [0.00001, 0.0001]``.
-    """
-    return _misc.random(target, element='pore', seed=seed, num_range=num_range)
+random.__doc__ = _misc.random.__doc__
 
 
 def largest_sphere(target, iters=10):
     r"""
-    Finds the maximum diameter pore that can be place in each location that
-    does not overlap with any neighbors.
+    Finds the maximum diameter pore that can be places in each location without
+    doverlapping any neighbors.
 
     Parameters
     ----------
@@ -89,17 +65,17 @@ def largest_sphere(target, iters=10):
 
     """
     network = target.project.network
-    D = _sp.zeros([network.Np, ], dtype=float)
+    D = _np.zeros([network.Np, ], dtype=float)
     Ps = network.pores(target.name)
     C1 = network['pore.coords'][network['throat.conns'][:, 0]]
     C2 = network['pore.coords'][network['throat.conns'][:, 1]]
-    L = _sp.sqrt(_sp.sum((C1 - C2)**2, axis=1))
+    L = _np.sqrt(_np.sum((C1 - C2)**2, axis=1))
     while iters >= 0:
         iters -= 1
-        Lt = L - _sp.sum(D[network['throat.conns']], axis=1)/2
+        Lt = L - _np.sum(D[network['throat.conns']], axis=1)/2
         am = network.create_adjacency_matrix(weights=Lt, fmt='lil')
-        D[Ps] = D[Ps] + _sp.array([_sp.amin(row) for row in am.data])[Ps]*0.95
-    if _sp.any(D < 0):
+        D[Ps] = D[Ps] + _np.array([_np.amin(row) for row in am.data])[Ps]*0.95
+    if _np.any(D < 0):
         _logger.warning('Negative pore diameters found!  Neighboring pores' +
                         ' must be larger than the pore spacing.')
     return D[network.pores(target.name)]
@@ -122,5 +98,5 @@ def equivalent_sphere(target, pore_volume='pore.volume'):
     """
     from scipy.special import cbrt
     pore_vols = target[pore_volume]
-    value = cbrt(6*pore_vols/_sp.pi)
+    value = cbrt(6*pore_vols/_np.pi)
     return value
diff --git a/openpnm/models/geometry/pore_surface_area.py b/openpnm/models/geometry/pore_surface_area.py
index 7d27c5864..d836d2de6 100644
--- a/openpnm/models/geometry/pore_surface_area.py
+++ b/openpnm/models/geometry/pore_surface_area.py
@@ -1,10 +1,4 @@
-r"""
-===============================================================================
-Submodule -- pore_surface_area
-===============================================================================
-
-"""
-import scipy as _sp
+import numpy as _np
 
 
 def sphere(target, pore_diameter='pore.diameter', throat_area='throat.area'):
@@ -31,19 +25,17 @@ def sphere(target, pore_diameter='pore.diameter', throat_area='throat.area'):
     """
     network = target.project.network
     R = target[pore_diameter]/2
-    Asurf = 4*_sp.constants.pi*R**2
+    Asurf = 4*_np.pi*R**2
     Tn = network.find_neighbor_throats(pores=target.Ps, flatten=False)
-    Tsurf = _sp.array([_sp.sum(network[throat_area][Ts]) for Ts in Tn])
+    Tsurf = _np.array([_np.sum(network[throat_area][Ts]) for Ts in Tn])
     value = Asurf - Tsurf
     return value
 
 
 def cube(target, pore_diameter='pore.diameter', throat_area='throat.area'):
     r"""
-    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.
+    Calculates internal surface area of pore bodies assuming they are cubes
+    then subtracts the area of the neighboring throats.
 
     Parameters
     ----------
@@ -62,6 +54,6 @@ def cube(target, pore_diameter='pore.diameter', throat_area='throat.area'):
     network = target.project.network
     D = target[pore_diameter]
     Tn = network.find_neighbor_throats(pores=target.Ps, flatten=False)
-    Tsurf = _sp.array([_sp.sum(network[throat_area][Ts]) for Ts in Tn])
+    Tsurf = _np.array([_np.sum(network[throat_area][Ts]) for Ts in Tn])
     value = 6*D**2 - Tsurf
     return value
diff --git a/openpnm/models/geometry/pore_volume.py b/openpnm/models/geometry/pore_volume.py
index d194f8223..9e1453216 100644
--- a/openpnm/models/geometry/pore_volume.py
+++ b/openpnm/models/geometry/pore_volume.py
@@ -1,10 +1,4 @@
-r"""
-===============================================================================
-pore_volume --
-===============================================================================
-
-"""
-import scipy as _sp
+from numpy import pi as _pi
 
 
 def sphere(target, pore_diameter='pore.diameter'):
@@ -19,11 +13,11 @@ def sphere(target, pore_diameter='pore.diameter'):
         necessary geometric properties.
 
     pore_diameter : string
-        The dictionary key to the pore diameter values
+        The dictionary key of the pore diameter values
 
     """
     diams = target[pore_diameter]
-    value = _sp.pi/6*diams**3
+    value = _pi/6*diams**3
     return value
 
 
@@ -39,7 +33,7 @@ def cube(target, pore_diameter='pore.diameter'):
         necessary geometric properties.
 
     pore_diameter : string
-        The dictionary key to the pore diameter values
+        The dictionary key of the pore diameter values
 
     """
     diams = target[pore_diameter]
diff --git a/openpnm/models/geometry/throat_area.py b/openpnm/models/geometry/throat_area.py
index 76502302a..109d2f9c2 100644
--- a/openpnm/models/geometry/throat_area.py
+++ b/openpnm/models/geometry/throat_area.py
@@ -1,10 +1,4 @@
-r"""
-===============================================================================
-Submodule -- throat_area
-===============================================================================
-
-"""
-import scipy as _sp
+from numpy import pi as _pi
 
 
 def cylinder(target, throat_diameter='throat.diameter'):
@@ -19,11 +13,11 @@ def cylinder(target, throat_diameter='throat.diameter'):
         necessary geometric properties.
 
     throat_diameter : string
-        Dictionary key to the throat diameter values
+        Dictionary key of the throat diameter values
 
     """
     diams = target[throat_diameter]
-    value = _sp.pi/4*(diams)**2
+    value = _pi/4*(diams)**2
     return value
 
 
@@ -39,7 +33,7 @@ def cuboid(target, throat_diameter='throat.diameter'):
         necessary geometric properties.
 
     throat_diameter : string
-        Dictionary key to the throat diameter values
+        Dictionary key of the throat diameter values
 
     """
     diams = target[throat_diameter]
diff --git a/openpnm/models/geometry/throat_length.py b/openpnm/models/geometry/throat_length.py
index 403b0b600..fb4f974a5 100644
--- a/openpnm/models/geometry/throat_length.py
+++ b/openpnm/models/geometry/throat_length.py
@@ -1,12 +1,6 @@
-r"""
-===============================================================================
-Submodule -- throat_length
-===============================================================================
-
-"""
 import scipy as _sp
-from openpnm.core import logging
-_logger = logging.getLogger(__name__)
+from openpnm.core import logging as _logging
+_logger = _logging.getLogger(__name__)
 
 
 def straight(target, pore_diameter='pore.diameter', L_negative=1e-9):
diff --git a/openpnm/models/geometry/throat_misc.py b/openpnm/models/geometry/throat_misc.py
deleted file mode 100644
index c5643ab4e..000000000
--- a/openpnm/models/geometry/throat_misc.py
+++ /dev/null
@@ -1,74 +0,0 @@
-r"""
-===============================================================================
-throat_misc --  Miscillaneous and generic functions to apply to throats
-===============================================================================
-
-"""
-import scipy as _sp
-from . import misc as _misc
-
-
-def random(target, seed=None, num_range=[0, 1]):
-    r"""
-    Generate an array of random numbers using the Scipy random number
-    generator.
-
-    Parameters
-    ----------
-    target : OpenPNM Object
-        The object which this model is associated with. This controls the
-        length of the calculated array, and also provides access to other
-        necessary properties.
-
-    seed : int
-        The seed value to initialize the random number generator.  The default
-        is None which means different numbers will be produced each time this
-        function is run.
-
-    num_range : list
-        The min and max values of the returned random numbers.  The values will
-        be uniformly distributed over this range.
-
-    Notes
-    -----
-    This method seems trivial but it can be useful to add this as a pore-scale
-    model so that random number can be regenerated upon demand to create
-    new realization of the pore and throat size distributions (i.e. assuming
-    ``seed`` is not set to None).
-    """
-    values = _misc.random(target, element='throat', seed=seed,
-                          num_range=num_range)
-    return values
-
-
-def neighbor(target, pore_prop='pore.seed', mode='min'):
-    r"""
-    Adopt a value based on the values in neighboring pores
-
-    Parameters
-    ----------
-    target : OpenPNM Object
-        The object which this model is associated with. This controls the
-        length of the calculated array, and also provides access to other
-        necessary properties.
-
-    pore_prop : string
-        The dictionary key to the array containing the pore property to be
-        used in the calculation.  Default is 'pore.seed'.
-
-    mode : string
-        Controls how the throat property is calculated.  Options are 'min',
-        'max' and 'mean'.
-
-    """
-    network = target.project.network
-    throats = network.throats(target.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/models/geometry/throat_perimeter.py b/openpnm/models/geometry/throat_perimeter.py
index 85513606c..30319af0f 100644
--- a/openpnm/models/geometry/throat_perimeter.py
+++ b/openpnm/models/geometry/throat_perimeter.py
@@ -1,9 +1,3 @@
-r"""
-===============================================================================
-Submodule -- throat_perimeter
-===============================================================================
-
-"""
 import scipy as _sp
 
 
diff --git a/openpnm/models/geometry/throat_seed.py b/openpnm/models/geometry/throat_seed.py
index 11a6d5ecd..9f3c67d7a 100644
--- a/openpnm/models/geometry/throat_seed.py
+++ b/openpnm/models/geometry/throat_seed.py
@@ -1,11 +1,4 @@
-r"""
-===============================================================================
-Submodule -- throat_seeds
-===============================================================================
-
-"""
-from . import misc as _misc
-import scipy as _sp
+from openpnm.models import misc as _misc
 
 
 def random(target, seed=None, num_range=[0, 1]):
@@ -16,36 +9,9 @@ def random(target, seed=None, num_range=[0, 1]):
 random.__doc__ = _misc.random.__doc__
 
 
-def neighbor(target, pore_prop='pore.seed', mode='min'):
-    r"""
-    Adopt a value based on the values in the neighboring pores
-
-    Parameters
-    ----------
-    target : OpenPNM Object
-        The object which this model is associated with. This controls the
-        length of the calculated array, and also provides access to other
-        necessary properties.
-
-    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.
-    """
-    network = target.project.network
-    throats = network.throats(target.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
+def from_neighbor_pores(target, pore_prop='pore.seed', mode='min'):
+    return _misc.from_neighbor_pores(target=target, pore_prop=pore_prop,
+                                     mode=mode)
+
+
+random.__doc__ = _misc.from_neighbor_pores.__doc__
diff --git a/openpnm/models/geometry/throat_shape_factor.py b/openpnm/models/geometry/throat_shape_factor.py
index 509982e78..f4ab3c4b7 100644
--- a/openpnm/models/geometry/throat_shape_factor.py
+++ b/openpnm/models/geometry/throat_shape_factor.py
@@ -1,9 +1,3 @@
-r"""
-===============================================================================
-Submodule -- throat_shape_factor
-===============================================================================
-
-"""
 import scipy as _sp
 
 
@@ -47,8 +41,10 @@ def mason_morrow(target, throat_perimeter='throat.perimeter',
                  throat_area='throat.area'):
     r"""
     Mason and Morrow relate the capillary pressure to the shaped factor in a
-    Similar way to Mortensen but for triangles.
-    Ref:
+    similar way to Mortensen but for triangles.
+
+    References
+    ----------
     Mason, G. and Morrow, N.R., 1991. Capillary behavior of a perfectly wetting
     liquid in irregular triangular tubes. Journal of Colloid and Interface
     Science, 141(1), pp.262-274.
@@ -68,7 +64,9 @@ def jenkins_rao(target, throat_perimeter='throat.perimeter',
     r"""
     Jenkins and Rao relate the capillary pressure in an eliptical throat to
     the aspect ratio
-    Ref:
+
+    References
+    ----------
     Jenkins, R.G. and Rao, M.B., 1984. The effect of elliptical pores on
     mercury porosimetry results. Powder technology, 38(2), pp.177-180.
     """
diff --git a/openpnm/models/geometry/throat_size.py b/openpnm/models/geometry/throat_size.py
index ac2f26ff1..0f829c93d 100644
--- a/openpnm/models/geometry/throat_size.py
+++ b/openpnm/models/geometry/throat_size.py
@@ -1,12 +1,5 @@
-r"""
-===============================================================================
-Submodule -- throat_diameter
-===============================================================================
-
-"""
-from . import misc as _misc
-from .throat_misc import neighbor as _neighbor
-import scipy as _sp
+from openpnm.models import misc as _misc
+import numpy as _np
 
 
 def weibull(target, shape, scale, loc, seeds='throat.seed'):
@@ -31,39 +24,22 @@ def generic(target, func, seeds='throat.seed'):
 generic.__doc__ = _misc.generic.__doc__
 
 
-def neighbor(target, pore_prop, mode='min'):
-    return _neighbor(target=target, pore_prop=pore_prop, mode=mode)
+def from_neighbor_pores(target, pore_prop, mode='min'):
+    return _misc.from_neighbor_pores(target=target, pore_prop=pore_prop,
+                                     mode=mode)
 
 
-neighbor.__doc__ = _neighbor.__doc__
+from_neighbor_pores.__doc__ = _misc.from_neighbor_pores.__doc__
 
 
 def random(target, seed=None, num_range=[0, 1]):
-    r"""
-    Assign throat sizes from a random
-
-    Parameters
-    ----------
-    target : OpenPNM Object
-        The object which this model is associated with. This controls the
-        length of the calculated array, and also provides access to other
-        necessary properties.
-
-    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.  The default is [0, 1] but this can be adjusted to produce
-        throat sizes directly; for instance pores between 10 and 100 um can be
-        generated with ``num_range = [0.00001, 0.0001]``.
-    """
     return _misc.random(target=target, element='throat', seed=seed,
                         num_range=num_range)
 
 
+random.__doc__ = _misc.random.__doc__
+
+
 def equivalent_circle(target, throat_area='throat.area'):
     r"""
     Calculates the diameter of a cirlce with same area as the throat.
@@ -79,5 +55,5 @@ def equivalent_circle(target, throat_area='throat.area'):
         The dictionary key to the throat area values
     """
     areas = target[throat_area]
-    value = 2*_sp.sqrt(areas/_sp.pi)
+    value = 2*_np.sqrt(areas/_np.pi)
     return value
diff --git a/openpnm/models/geometry/throat_surface_area.py b/openpnm/models/geometry/throat_surface_area.py
index 096f875e7..f0952c0e3 100644
--- a/openpnm/models/geometry/throat_surface_area.py
+++ b/openpnm/models/geometry/throat_surface_area.py
@@ -1,9 +1,3 @@
-r"""
-===============================================================================
-Submodule -- throat_surface_area
-===============================================================================
-
-"""
 import scipy as _sp
 
 
diff --git a/openpnm/models/geometry/throat_volume.py b/openpnm/models/geometry/throat_volume.py
index 440dfbb9a..49771aa86 100644
--- a/openpnm/models/geometry/throat_volume.py
+++ b/openpnm/models/geometry/throat_volume.py
@@ -1,9 +1,3 @@
-r"""
-===============================================================================
-Submodule -- throat_volume
-===============================================================================
-
-"""
 import scipy as _sp
 
 
diff --git a/openpnm/models/misc/__init__.py b/openpnm/models/misc/__init__.py
new file mode 100644
index 000000000..aea35ad45
--- /dev/null
+++ b/openpnm/models/misc/__init__.py
@@ -0,0 +1,11 @@
+from .misc import constant
+from .misc import product
+from .misc import random
+from .misc import scaled
+from .misc import linear
+from .misc import polynomial
+from .misc import weibull
+from .misc import normal
+from .misc import generic
+from .misc import from_neighbor_throats
+from .misc import from_neighbor_pores
diff --git a/openpnm/models/geometry/misc.py b/openpnm/models/misc/misc.py
similarity index 50%
rename from openpnm/models/geometry/misc.py
rename to openpnm/models/misc/misc.py
index 01224f3d5..b31ea69c8 100644
--- a/openpnm/models/geometry/misc.py
+++ b/openpnm/models/misc/misc.py
@@ -1,18 +1,46 @@
-r"""
-===============================================================================
-Submodule -- Miscillaneous functions
-===============================================================================
+import numpy as np
+import scipy.stats as spts
 
-"""
-import scipy as _sp
-import scipy.stats as _spts
 
+def constant(target, value):
+    r"""
+    Places a constant value into the target object
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The object which this model is associated with. This controls the
+        length of the calculated array, and also provides access to other
+        necessary properties.
+
+    value : scalar
+        The numerical value to apply
+
+    Notes
+    -----
+    This model is mostly useless and for testing purposes, but might be used
+    to 'reset' an array back to a default value.
+    """
+    return value
 
-def product(target, **props):
+
+def product(target, props):
     r"""
+    Calculates the product of multiple property values
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The object which this model is associated with. This controls the
+        length of the calculated array, and also provides access to other
+        necessary properties.
+
+    props : list of strings
+        The names of each property to be multiplied
+
     """
     value = 1.0
-    for item in props.values():
+    for item in props:
         value *= target[item]
     return value
 
@@ -40,8 +68,8 @@ def random(target, element, seed=None, num_range=[0, 1]):
     """
     range_size = num_range[1] - num_range[0]
     range_min = num_range[0]
-    _sp.random.seed(seed)
-    value = _sp.random.rand(target._count(element),)
+    np.random.seed(seed)
+    value = np.random.rand(target._count(element),)
     value = value*range_size + range_min
     return value
 
@@ -68,7 +96,59 @@ def scaled(target, prop, factor):
     return value
 
 
-def weibull(target, shape, scale, loc, seeds):
+def linear(target, m, b, prop):
+    r"""
+    Calculates a property as a linear function of a given property
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The 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.
+
+    m, b : floats
+        Slope and intercept of the linear corelation
+
+    prop : string
+        The dictionary key containing the independent variable or phase
+        property to be used in the correlation.
+
+    """
+    x = target[prop]
+    value = m*x + b
+    return value
+
+
+def polynomial(target, a, prop, **kwargs):
+    r"""
+    Calculates a property as a polynomial function of a given property
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The 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.
+
+    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 dictionary key containing the independent variable or phase
+        property to be used in the correlation.
+
+    """
+    x = target[prop]
+    value = 0.0
+    for i in range(0, len(a)):
+        value += a[i]*x**i
+    return value
+
+
+def weibull(target, seeds, shape, scale, loc):
     r"""
     Produces values from a Weibull distribution given a set of random numbers.
 
@@ -79,6 +159,10 @@ def weibull(target, shape, scale, loc, seeds):
         length of the calculated array, and also provides access to other
         necessary properties.
 
+    seeds : string, optional
+        The dictionary key on the Geometry object containing random seed values
+        (between 0 and 1) to use in the statistical distribution.
+
     shape : float
         This controls the skewness of the distribution, with 'shape' < 1 giving
         values clustered on the low end of the range with a long tail, and
@@ -92,12 +176,6 @@ def weibull(target, shape, scale, loc, seeds):
         Applies an offset to the distribution such that the smallest values are
         above this number.
 
-    seeds : string, optional
-        The dictionary key on the Geometry object containing random seed values
-        (between 0 and 1) to use in the statistical distribution.  If none is
-        specified, then an array of random numbers will be automatically
-        generated and stored on the Geometry object.
-
     Examples
     --------
     The following code illustrates the inner workings of this function,
@@ -112,11 +190,11 @@ def weibull(target, shape, scale, loc, seeds):
 
     """
     seeds = target[seeds]
-    value = _spts.weibull_min.ppf(q=seeds, c=shape, scale=scale, loc=loc)
+    value = spts.weibull_min.ppf(q=seeds, c=shape, scale=scale, loc=loc)
     return value
 
 
-def normal(target, scale, loc, seeds):
+def normal(target, seeds, scale, loc):
     r"""
     Produces values from a Weibull distribution given a set of random numbers.
 
@@ -128,18 +206,16 @@ def normal(target, scale, loc, seeds):
         geom.Nt) and (2) provide access to other necessary values
         (i.e. geom['pore.seed']).
 
+    seeds : string, optional
+        The dictionary key on the Geometry object containing random seed values
+        (between 0 and 1) to use in the statistical distribution.
+
     scale : float
         The standard deviation of the Normal distribution
 
     loc : float
         The mean of the Normal distribution
 
-    seeds : string, optional
-        The dictionary key on the Geometry object containing random seed values
-        (between 0 and 1) to use in the statistical distribution.  If none is
-        specified, then an array of random numbers will be automatically
-        generated and stored on teh Geometry object.
-
     Examples
     --------
     The following code illustrates the inner workings of this function,
@@ -153,11 +229,11 @@ def normal(target, scale, loc, seeds):
 
     """
     seeds = target[seeds]
-    value = _spts.norm.ppf(q=seeds, scale=scale, loc=loc)
+    value = spts.norm.ppf(q=seeds, scale=scale, loc=loc)
     return value
 
 
-def generic(target, func, seeds):
+def generic(target, seeds, func):
     r"""
     Accepts an 'rv_frozen' object from the Scipy.stats submodule and returns
     values from the distribution for the given seeds using the ``ppf`` method.
@@ -169,16 +245,14 @@ def generic(target, func, seeds):
         length of the calculated array, and also provides access to other
         necessary properties.
 
+    seeds : string, optional
+        The dictionary key on the Geometry object containing random seed values
+        (between 0 and 1) to use in the statistical distribution.
+
     func : object
         An 'rv_frozen' object from the Scipy.stats library with all of the
         parameters pre-specified.
 
-    seeds : string, optional
-        The dictionary key on the Geometry object containing random seed values
-        (between 0 and 1) to use in the statistical distribution.  If none is
-        specified, then an array of random numbers will be automatically
-        generated and stored on teh Geometry object.
-
     Examples
     --------
     The following code illustrates the process of obtaining a 'frozen' Scipy
@@ -194,3 +268,74 @@ def generic(target, func, seeds):
     seeds = target[seeds]
     value = func.ppf(seeds)
     return value
+
+
+def from_neighbor_throats(target, throat_prop='throat.seed', mode='min'):
+    r"""
+    Adopt a value from the values found in neighboring throats
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The object which this model is associated with. This controls the
+        length of the calculated array, and also provides access to other
+        necessary properties.
+
+    throat_prop : string
+        The dictionary key of the array containing the throat property to be
+        used in the calculation.  The default is 'throat.seed'.
+
+    mode : string
+        Controls how the pore property is calculated.  Options are 'min',
+        'max' and 'mean'.
+    """
+    network = target.project.network
+    Ps = target.pores()
+    data = target[throat_prop]
+    neighborTs = network.find_neighbor_throats(pores=Ps,
+                                               flatten=False,
+                                               mode='intersection')
+    values = np.ones((np.shape(Ps)[0],))*np.nan
+    if mode == 'min':
+        for pore in Ps:
+            values[pore] = np.amin(data[neighborTs[pore]])
+    if mode == 'max':
+        for pore in Ps:
+            values[pore] = np.amax(data[neighborTs[pore]])
+    if mode == 'mean':
+        for pore in Ps:
+            values[pore] = np.mean(data[neighborTs[pore]])
+    return values
+
+
+def from_neighbor_pores(target, pore_prop='pore.seed', mode='min'):
+    r"""
+    Adopt a value based on the values in neighboring pores
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The object which this model is associated with. This controls the
+        length of the calculated array, and also provides access to other
+        necessary properties.
+
+    pore_prop : string
+        The dictionary key to the array containing the pore property to be
+        used in the calculation.  Default is 'pore.seed'.
+
+    mode : string
+        Controls how the throat property is calculated.  Options are 'min',
+        'max' and 'mean'.
+
+    """
+    network = target.project.network
+    throats = network.throats(target.name)
+    P12 = network.find_connected_pores(throats)
+    pvalues = network[pore_prop][P12]
+    if mode == 'min':
+        value = np.amin(pvalues, axis=1)
+    if mode == 'max':
+        value = np.amax(pvalues, axis=1)
+    if mode == 'mean':
+        value = np.mean(pvalues, axis=1)
+    return value
diff --git a/openpnm/models/phase/__init__.py b/openpnm/models/phase/__init__.py
index 2aaa2728e..c6911b209 100644
--- a/openpnm/models/phase/__init__.py
+++ b/openpnm/models/phase/__init__.py
@@ -2,7 +2,6 @@ from . import contact_angle
 from . import density
 from . import diffusivity
 from . import electrical_conductivity
-from . import misc
 from . import molar_density
 from . import surface_tension
 from . import thermal_conductivity
diff --git a/openpnm/models/phase/density.py b/openpnm/models/phase/density.py
index 206dd10db..592b84f51 100644
--- a/openpnm/models/phase/density.py
+++ b/openpnm/models/phase/density.py
@@ -1,10 +1,3 @@
-r"""
-===============================================================================
-Submodule -- density
-===============================================================================
-
-"""
-import scipy as _sp
 
 
 def standard(target, mol_weight='pore.molecular_weight',
diff --git a/openpnm/models/phase/misc.py b/openpnm/models/phase/misc.py
deleted file mode 100644
index df4ab0574..000000000
--- a/openpnm/models/phase/misc.py
+++ /dev/null
@@ -1,61 +0,0 @@
-r"""
-===============================================================================
-Submodule -- miscillaneous
-===============================================================================
-
-Models for applying basic phase properties
-
-"""
-import scipy as _sp
-
-
-def linear(target, m, b, prop):
-    r"""
-    Calculates a property as a linear function of a given property
-
-    Parameters
-    ----------
-    target : OpenPNM Object
-        The 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.
-
-    m, b : floats
-        Slope and intercept of the linear corelation
-
-    prop : string
-        The dictionary key containing the independent variable or phase
-        property to be used in the correlation.
-
-    """
-    x = target[prop]
-    value = m*x + b
-    return value
-
-
-def polynomial(target, a, prop, **kwargs):
-    r"""
-    Calculates a property as a polynomial function of a given property
-
-    Parameters
-    ----------
-    target : OpenPNM Object
-        The 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.
-
-    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 dictionary key containing the independent variable or phase
-        property to be used in the correlation.
-
-    """
-    x = target[prop]
-    value = 0.0
-    for i in range(0, len(a)):
-        value += a[i]*x**i
-    return value
diff --git a/openpnm/models/phase/vapor_pressure.py b/openpnm/models/phase/vapor_pressure.py
index 777368129..c7bf71fd8 100644
--- a/openpnm/models/phase/vapor_pressure.py
+++ b/openpnm/models/phase/vapor_pressure.py
@@ -1,13 +1,4 @@
-# -*- coding: utf-8 -*-
-r"""
-===============================================================================
-Submodule -- vapor_pressure
-===============================================================================
-
-Methods for predicing the vapor pressure of pure species
-
-"""
-import scipy as sp
+import numpy as np
 
 
 def antoine(target, A, B, C, temperature='pore.temperature'):
@@ -88,7 +79,7 @@ def water(target, temperature='pore.temperature', salinity='pore.salinity'):
     a4 = 4.1764768E-05
     a5 = -1.4452093E-08
     a6 = 6.5459673E+00
-    Pv_w = sp.exp((a1/T) + a2 + a3*T + a4*T**2 + a5*T**3 + a6*sp.log(T))
+    Pv_w = np.exp((a1/T) + a2 + a3*T + a4*T**2 + a5*T**3 + a6*np.log(T))
     Pv_sw = Pv_w/(1+0.57357*(S/(1000-S)))
     value = Pv_sw
     return value
diff --git a/openpnm/models/physics/capillary_pressure.py b/openpnm/models/physics/capillary_pressure.py
index b2bacfb11..69a87f106 100644
--- a/openpnm/models/physics/capillary_pressure.py
+++ b/openpnm/models/physics/capillary_pressure.py
@@ -46,14 +46,24 @@ def washburn(target, surface_tension='pore.surface_tension',
     """
     network = target.project.network
     phase = target.project.find_phase(target)
-    if surface_tension.split('.')[0] == 'pore':
-        sigma = phase.interpolate_data(propname=surface_tension)
-    else:
-        sigma = phase[surface_tension]
-    if contact_angle.split('.')[0] == 'pore':
-        theta = phase.interpolate_data(propname=contact_angle)
-    else:
-        theta = phase[contact_angle]
+    if 'pore' in diameter:
+        if 'throat' in surface_tension:
+            sigma = phase.interpolate_data(propname=surface_tension)
+        else:
+            sigma = phase[surface_tension]
+        if 'throat' in contact_angle:
+            theta = phase.interpolate_data(propname=contact_angle)
+        else:
+            theta = phase[contact_angle]
+    if 'throat' in diameter:
+        if 'pore' in surface_tension:
+            sigma = phase.interpolate_data(propname=surface_tension)
+        else:
+            sigma = phase[surface_tension]
+        if 'pore' in contact_angle:
+            theta = phase.interpolate_data(propname=contact_angle)
+        else:
+            theta = phase[contact_angle]
     r = network[diameter]/2
     value = -2*sigma*_sp.cos(_sp.radians(theta))/r
     if diameter.split('.')[0] == 'throat':
@@ -114,16 +124,24 @@ def purcell(target, r_toroid, surface_tension='pore.surface_tension',
     """
     network = target.project.network
     phase = target.project.find_phase(target)
-    if surface_tension.split('.')[0] == 'pore':
-        sigma = phase[surface_tension]
-        sigma = phase.interpolate_data(propname=surface_tension)
-    else:
-        sigma = phase[surface_tension]
-    if contact_angle.split('.')[0] == 'pore':
-        theta = phase[contact_angle]
-        theta = phase.interpolate_data(propname=contact_angle)
-    else:
-        theta = phase[contact_angle]
+    if 'pore' in diameter:
+        if 'throat' in surface_tension:
+            sigma = phase.interpolate_data(propname=surface_tension)
+        else:
+            sigma = phase[surface_tension]
+        if 'throat' in contact_angle:
+            theta = phase.interpolate_data(propname=contact_angle)
+        else:
+            theta = phase[contact_angle]
+    if 'throat' in diameter:
+        if 'pore' in surface_tension:
+            sigma = phase.interpolate_data(propname=surface_tension)
+        else:
+            sigma = phase[surface_tension]
+        if 'pore' in contact_angle:
+            theta = phase.interpolate_data(propname=contact_angle)
+        else:
+            theta = phase[contact_angle]
     r = network[diameter]/2
     R = r_toroid
     alpha = theta - 180 + _sp.arcsin(_sp.sin(_sp.radians(theta)/(1+r/R)))
diff --git a/openpnm/models/physics/diffusive_conductance.py b/openpnm/models/physics/diffusive_conductance.py
index 815727139..df410853c 100644
--- a/openpnm/models/physics/diffusive_conductance.py
+++ b/openpnm/models/physics/diffusive_conductance.py
@@ -1,12 +1,5 @@
-r"""
-===============================================================================
-Submodule -- diffusive_conductance
-===============================================================================
-
-"""
-
+import openpnm as op
 import scipy as _sp
-import openpnm.utils.misc as misc
 
 
 def bulk_diffusion(target, molar_density='pore.molar_density',
@@ -55,7 +48,7 @@ def bulk_diffusion(target, molar_density='pore.molar_density',
     DABp = phase[diffusivity]
     DABt = phase.interpolate_data(propname=diffusivity)
     if calc_pore_len:
-        lengths = misc.conduit_lengths(network, mode='centroid')
+        lengths = op.utils.misc.conduit_lengths(network, mode='centroid')
         plen1 = lengths[:, 0]
         plen2 = lengths[:, 2]
     else:
diff --git a/openpnm/models/physics/hydraulic_conductance.py b/openpnm/models/physics/hydraulic_conductance.py
index fb03dc284..5d1605df1 100644
--- a/openpnm/models/physics/hydraulic_conductance.py
+++ b/openpnm/models/physics/hydraulic_conductance.py
@@ -1,12 +1,5 @@
-r"""
-===============================================================================
-Submodule -- hydraulic_conductance
-===============================================================================
-
-"""
-
+import openpnm as op
 import scipy as _sp
-import openpnm.utils.misc as misc
 
 
 def hagen_poiseuille(target, viscosity='pore.viscosity',
@@ -36,10 +29,15 @@ def hagen_poiseuille(target, viscosity='pore.viscosity',
     # Get Nt-by-2 list of pores connected to each throat
     Ps = network['throat.conns']
     # Get properties in every pore in the network
-    mut = phase[viscosity]
+    if 'pore' in viscosity:
+        mut = phase.interpolate_data(viscosity)
+        mup = phase[viscosity]
+    elif 'throat' in viscosity:
+        mup = phase.interpolate_data(viscosity)
+        mut = phase[viscosity]
     pdia = network[pore_diameter]
     if calc_pore_len:
-        lengths = misc.conduit_lengths(network, mode='centroid')
+        lengths = op.utils.misc.conduit_lengths(network, mode='centroid')
         plen1 = lengths[:, 0]
         plen2 = lengths[:, 2]
     else:
diff --git a/openpnm/models/physics/multiphase.py b/openpnm/models/physics/multiphase.py
index 1ee177735..46657faf2 100644
--- a/openpnm/models/physics/multiphase.py
+++ b/openpnm/models/physics/multiphase.py
@@ -8,10 +8,10 @@ Submodule -- diffusive_conductance
 import scipy as sp
 
 
-def conduit_conductance(physics, phase, network, throat_conductance,
+def conduit_conductance(target, throat_conductance,
                         throat_occupancy='throat.occupancy',
                         pore_occupancy='pore.occupancy',
-                        mode='strict', factor=1e-6, **kwargs):
+                        mode='strict', factor=1e-6):
     r"""
     Add a new multiphase conductance property to the conduits of network, where
     a conduit is ( 1/2 pore - full throat - 1/2 pore ) based on the areas.
@@ -52,7 +52,7 @@ def conduit_conductance(physics, phase, network, throat_conductance,
     calculated.
 
     """
-    throats = phase.Ts
+    throats = target.Ts
     if mode == 'loose':
         closed_conduits = ~sp.array(phase[throat_occupancy], dtype=bool)
     else:
diff --git a/openpnm/network/GenericNetwork.py b/openpnm/network/GenericNetwork.py
index d4c8a724a..da040a740 100644
--- a/openpnm/network/GenericNetwork.py
+++ b/openpnm/network/GenericNetwork.py
@@ -485,11 +485,8 @@ class GenericNetwork(Base, ModelsMixin):
         pores = self._parse_indices(pores)
         if sp.size(pores) == 0:
             return sp.array([], ndmin=1, dtype=int)
-        if 'lil' not in self._am.keys():
-            self.get_adjacency_matrix(fmt='lil')
         neighbors = topotools.find_neighbor_sites(sites=pores, logic=mode,
-                                                  am=self._am['lil'],
-                                                  flatten=flatten,
+                                                  am=self.am, flatten=flatten,
                                                   exclude_input=excl_self)
         return neighbors
 
@@ -542,11 +539,8 @@ class GenericNetwork(Base, ModelsMixin):
         pores = self._parse_indices(pores)
         if sp.size(pores) == 0:
             return sp.array([], ndmin=1, dtype=int)
-        if 'lil' not in self._im.keys():
-            self.get_incidence_matrix(fmt='lil')
         neighbors = topotools.find_neighbor_bonds(sites=pores, logic=mode,
-                                                  im=self._im['lil'],
-                                                  flatten=flatten)
+                                                  im=self.im, flatten=flatten)
         return neighbors
 
     def _find_neighbors(self, pores, element, **kwargs):
diff --git a/openpnm/phases/Air.py b/openpnm/phases/Air.py
index 6f80793ab..e4e00cc7f 100644
--- a/openpnm/phases/Air.py
+++ b/openpnm/phases/Air.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 from openpnm.phases import GenericPhase
-import openpnm.models.phase as fm
+import openpnm.models as mods
 
 
 class Air(GenericPhase):
@@ -33,8 +33,30 @@ class Air(GenericPhase):
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
 
-        self['pore.molecular_weight'] = 0.02896            # kg/mol
-        self['pore.critical_pressure'] = 3.786E6           # Pa
-        self['pore.critical_temperature'] = 132.5          # K
-        self['pore.critical_volume'] = 0.002917            # kg/m3
-        self['pore.contact_angle'] = 110.0                 # Degree
+        self['pore.molecular_weight'] = 0.02896
+        self['pore.critical_pressure'] = 3.786E6
+        self['pore.critical_temperature'] = 132.5
+        self['pore.critical_volume'] = 0.002917
+        self['pore.contact_angle'] = 110.0
+        self['pore.molecular_weight'] = 0.02896
+        self['pore.critical_pressure'] = 3.786E6
+        self['pore.critical_temperature'] = 132.5
+        self['pore.critical_volume'] = 0.002917
+        self['pore.contact_angle'] = 110.0
+        self.add_model(propname='pore.density',
+                       model=mods.phase.density.ideal_gas)
+        self.add_model(propname='pore.molar_density',
+                       model=mods.phase.molar_density.ideal_gas)
+        self.add_model(propname='pore.diffusivity',
+                       model=mods.phase.diffusivity.fuller,
+                       MA=0.032, MB=0.028,
+                       vA=16.6, vB=17.9)
+        self.add_model(propname='pore.thermal_conductivity',
+                       model=mods.misc.polynomial,
+                       prop='pore.temperature',
+                       a=[0.00422791, 0.0000789606, -1.56383E-08])
+        self.add_model(propname='pore.viscosity',
+                       model=mods.misc.polynomial,
+                       prop='pore.temperature',
+                       a=[0.00000182082, 6.51815E-08, -3.48553E-11,
+                          1.11409E-14])
diff --git a/openpnm/phases/Mercury.py b/openpnm/phases/Mercury.py
index 1e8e842de..480e4f841 100644
--- a/openpnm/phases/Mercury.py
+++ b/openpnm/phases/Mercury.py
@@ -1,6 +1,5 @@
-# -*- coding: utf-8 -*-
 from openpnm.phases import GenericPhase
-import openpnm.models.phase as fm
+import openpnm.models as mods
 
 
 class Mercury(GenericPhase):
@@ -32,30 +31,30 @@ class Mercury(GenericPhase):
     """
     def __init__(self, name=None, **kwargs):
         super().__init__(name=name, **kwargs)
-        self['pore.molecular_weight'] = 0.2006                         # kg/mol
-        self['pore.critical_pressure'] = 1.662E8                       # Pa
-        self['pore.critical_temperature'] = 1733                       # K
-        self['pore.critical_volume'] = 0.000189                        # kg/m3
-        self['pore.contact_angle'] = 140                               # Degree
-        self.add_model(propname='pore.vapor_pressure',                 # Pa
-                       model=fm.vapor_pressure.antoine,
+        self['pore.molecular_weight'] = 0.2006
+        self['pore.critical_pressure'] = 1.662E8
+        self['pore.critical_temperature'] = 1733
+        self['pore.critical_volume'] = 0.000189
+        self['pore.contact_angle'] = 140
+        self.add_model(propname='pore.vapor_pressure',
+                       model=mods.phase.vapor_pressure.antoine,
                        A=9.85767, B=3007.129, C=-10.001)
-        self.add_model(propname='pore.density',                        # kg/m3
-                       model=fm.misc.linear,
+        self.add_model(propname='pore.density',
+                       model=mods.misc.linear,
                        prop='pore.temperature',
                        b=14280.9, m=-2.47004)
         self.add_model(propname='pore.molar_density',
-                       model=fm.molar_density.standard)                # mol/m3
-        self.add_model(propname='pore.surface_tension',                # N/m
-                       model=fm.misc.linear,
+                       model=mods.phase.molar_density.standard)
+        self.add_model(propname='pore.surface_tension',
+                       model=mods.misc.linear,
                        prop='pore.temperature',
                        b=0.56254, m=-0.00028)
-        self.add_model(propname='pore.thermal_conductivity',           # W/m.K
-                       model=fm.misc.polynomial,
+        self.add_model(propname='pore.thermal_conductivity',
+                       model=mods.misc.polynomial,
                        prop='pore.temperature',
                        a=[3.98691, 0.0170967, -0.0000063862])
-        self.add_model(propname='pore.viscosity',                      # kg/m.s
-                       model=fm.misc.polynomial,
+        self.add_model(propname='pore.viscosity',
+                       model=mods.misc.polynomial,
                        prop='pore.temperature',
                        a=[0.00355837, -0.0000100131, 1.23684E-08, -5.1684E-12])
         self.regenerate_models()
diff --git a/openpnm/phases/Water.py b/openpnm/phases/Water.py
index a0b93f336..a3ef2ba50 100644
--- a/openpnm/phases/Water.py
+++ b/openpnm/phases/Water.py
@@ -1,5 +1,5 @@
 from openpnm.phases import GenericPhase
-import openpnm.models.phase as fm
+import openpnm.models as mods
 
 
 class Water(GenericPhase):
@@ -25,23 +25,23 @@ class Water(GenericPhase):
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
 
-        self['pore.molecular_weight'] = 0.01802               # kg/mol
-        self['pore.critical_pressure'] = 2.2064E7             # Pa
-        self['pore.critical_temperature'] = 647.1             # K
-        self['pore.critical_volume'] = 0.003106               # kg/m3
-        self['pore.contact_angle'] = 110.0                    # Degree
-        self['pore.diffusivity'] = 1e-9                       # m2/s
+        self['pore.molecular_weight'] = 0.01802
+        self['pore.critical_pressure'] = 2.2064E7
+        self['pore.critical_temperature'] = 647.1
+        self['pore.critical_volume'] = 0.003106
+        self['pore.contact_angle'] = 110.0
+        self['pore.diffusivity'] = 1e-9
         self.add_model(propname='pore.density',
-                       model=fm.density.water)               # kg/m3
+                       model=mods.phase.density.water)
         self.add_model(propname='pore.molar_density',
-                       model=fm.molar_density.standard)      # mol/m3
+                       model=mods.phase.molar_density.standard)
         self.add_model(propname='pore.surface_tension',
-                       model=fm.surface_tension.water)       # N/m
+                       model=mods.phase.surface_tension.water)
         self.add_model(propname='pore.thermal_conductivity',
-                       model=fm.thermal_conductivity.water)  # W/m.K
-        self.add_model(propname='pore.vapor_pressure',       # Pa
-                       model=fm.vapor_pressure.antoine,
+                       model=mods.phase.thermal_conductivity.water)
+        self.add_model(propname='pore.vapor_pressure',
+                       model=mods.phase.vapor_pressure.antoine,
                        A=8.088, B=1750.71, C=236.191)
         self.add_model(propname='pore.viscosity',
-                       model=fm.viscosity.water)             # kg/m.s
+                       model=mods.phase.viscosity.water)
         self.regenerate_models()


Test patch:
diff --git a/tests/temp/dummy b/tests/temp/dummy
deleted file mode 100644
index e69de29bb..000000000
diff --git a/tests/unit/models/geometry/PoreAreaTest.py b/tests/unit/models/geometry/PoreAreaTest.py
new file mode 100644
index 000000000..8db1db90e
--- /dev/null
+++ b/tests/unit/models/geometry/PoreAreaTest.py
@@ -0,0 +1,44 @@
+import openpnm as op
+import openpnm.models.geometry.pore_area as mods
+import numpy as np
+
+
+class PoreAreaTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5], spacing=1.0)
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.air,
+                                              geometry=self.geo)
+        self.geo['pore.diameter'] = 1.0
+        self.geo['throat.area'] = 0.1
+
+    def test_sphere(self):
+        self.geo.add_model(propname='pore.area',
+                           model=mods.sphere,
+                           regen_mode='normal')
+        a = np.array([0.78539816])
+        b = np.unique(self.geo['pore.area'])
+        assert np.allclose(a, b)
+
+    def test_cube(self):
+        self.geo.add_model(propname='pore.area',
+                           model=mods.cube,
+                           regen_mode='normal')
+        a = np.array([1.0])
+        b = np.unique(self.geo['pore.area'])
+        assert np.allclose(a, b)
+
+
+if __name__ == '__main__':
+
+    t = PoreAreaTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/geometry/PoreCentroidTest.py b/tests/unit/models/geometry/PoreCentroidTest.py
new file mode 100644
index 000000000..db8bd253f
--- /dev/null
+++ b/tests/unit/models/geometry/PoreCentroidTest.py
@@ -0,0 +1,3 @@
+class PoreCentroidTest:
+    def test_voronoi(self):
+        pass
diff --git a/tests/unit/models/geometry/PoreSeedTest.py b/tests/unit/models/geometry/PoreSeedTest.py
new file mode 100644
index 000000000..9466c8472
--- /dev/null
+++ b/tests/unit/models/geometry/PoreSeedTest.py
@@ -0,0 +1,71 @@
+import openpnm as op
+import scipy as sp
+import openpnm.models.geometry.pore_seed as mods
+
+
+class PoreSeedTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5])
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+
+    def test_random(self):
+        f = mods.random
+        self.geo.add_model(propname='pore.seed',
+                           model=f,
+                           seed=0,
+                           num_range=[0.1, 2])
+        assert sp.amax(self.geo['pore.seed']) > 1
+        assert sp.amin(self.geo['pore.seed']) < 1
+
+    def test_spatially_correlated(self):
+        f = mods.spatially_correlated
+        self.geo.add_model(propname='pore.seed',
+                           model=f,
+                           weights=[2, 2, 2],
+                           regen_mode='normal')
+        assert sp.amin(self.geo['pore.seed'] > 0)
+        assert sp.amax(self.geo['pore.seed'] < 1)
+
+    def test_spatially_correlated_zero_weights(self):
+        f = mods.spatially_correlated
+        self.geo.add_model(propname='pore.seed',
+                           model=f,
+                           weights=[0, 0, 0],
+                           regen_mode='normal')
+        assert sp.amin(self.geo['pore.seed'] > 0)
+        assert sp.amax(self.geo['pore.seed'] < 1)
+
+    def test_location_adjusted(self):
+        image1 = sp.ones([5, 5, 5])
+        image1[:, :, 0] = 0.5
+        f = mods.location_adjusted
+        self.geo.add_model(propname='pore.seed',
+                           model=f,
+                           image=image1)
+        Ps = self.net.pores('bottom')
+        a = sp.mean(self.geo["pore.seed"][Ps])
+        Other = self.net.pores('bottom', mode='not')
+        b = sp.mean(self.geo["pore.seed"][Other])
+        assert a < b
+        image2 = sp.ones([2, 2])
+        self.geo.add_model(propname='pore.seed',
+                           model=f,
+                           image=image2,
+                           regen_mode='normal')
+        self.geo.add_model(propname='pore.seed',
+                           model=f,
+                           image=-image1,
+                           regen_mode='normal')
+
+
+if __name__ == '__main__':
+
+    t = PoreSeedTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/geometry/PoreSizeTest.py b/tests/unit/models/geometry/PoreSizeTest.py
new file mode 100644
index 000000000..36e72677a
--- /dev/null
+++ b/tests/unit/models/geometry/PoreSizeTest.py
@@ -0,0 +1,84 @@
+import openpnm as op
+import scipy as sp
+import scipy.stats as spst
+import openpnm.models.geometry.pore_size as mods
+
+
+class PoreSizeTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5])
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.geo['pore.seed'] = sp.rand(self.geo.Np)
+
+    def test_normal(self):
+        self.geo.add_model(propname='pore.diameter',
+                           model=mods.normal,
+                           scale=0.01,
+                           loc=0.5,
+                           seeds='pore.seed')
+        assert 0.45 < sp.mean(self.geo['pore.diameter']) < 0.55
+        del self.geo['pore.diameter']
+
+    def test_weibull(self):
+        self.geo.add_model(propname='pore.diameter',
+                           model=mods.weibull,
+                           shape=1.5,
+                           scale=0.0001,
+                           loc=0.001,
+                           seeds='pore.seed')
+        assert sp.amin(self.geo['pore.diameter']) > 0.001
+        del self.geo['pore.diameter']
+
+    def test_generic(self):
+        func = spst.gamma(a=2, loc=0.001, scale=0.0001)
+        self.geo.add_model(propname='pore.diameter',
+                           model=mods.generic,
+                           func=func,
+                           seeds='pore.seed')
+        assert sp.amin(self.geo['pore.diameter']) > 0.001
+        del self.geo['pore.diameter']
+
+    def test_largest_sphere(self):
+        net = op.network.Cubic(shape=[5, 5, 5], spacing=[0.1, 0.2, 0.3])
+        geo = op.geometry.GenericGeometry(network=net, pores=net.Ps,
+                                          throats=net.Ts)
+        geo.add_model(propname='pore.diameter',
+                      model=mods.largest_sphere,
+                      iters=1)
+        dmin = sp.amin(geo['pore.diameter'])
+        assert dmin <= 0.1
+        geo.models['pore.diameter']['iters'] = 5
+        geo.regenerate_models()
+        assert dmin < sp.amin(geo['pore.diameter'])
+
+    def test_largest_sphere_multiple_geometries(self):
+        net = op.network.Cubic(shape=[10, 10, 10], spacing=[5, 5, 5])
+        net['pore.coords'][net.pores('top')] += [0, 0, -3]
+        geom2 = op.geometry.GenericGeometry(network=net,
+                                            pores=net.pores('top'))
+        geom2['pore.diameter'] = 1.0
+        Ps = net.pores('top', mode='not')
+        geom1 = op.geometry.GenericGeometry(network=net, pores=Ps,
+                                            throats=net.Ts)
+        mod = op.models.geometry.pore_size.largest_sphere
+        geom1.add_model(propname='pore.diameter',
+                        model=mod,
+                        iters=15)
+        assert sp.all(geom2['pore.diameter'] == 1.0)
+#        assert sp.all(sp.ceil(sp`.unique(geom1['pore.diameter'])) == [3.0, 5.0])
+
+
+if __name__ == '__main__':
+
+    t = PoreSizeTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            try:
+                t.__getattribute__(item)()
+            except TypeError:
+                t.__getattribute__(item)(tmpdir=py.path.local())
diff --git a/tests/unit/models/geometry/PoreSurfaceAreaTest.py b/tests/unit/models/geometry/PoreSurfaceAreaTest.py
new file mode 100644
index 000000000..1cabde6cb
--- /dev/null
+++ b/tests/unit/models/geometry/PoreSurfaceAreaTest.py
@@ -0,0 +1,44 @@
+import openpnm as op
+import openpnm.models.geometry.pore_surface_area as mods
+import numpy as np
+
+
+class PoreSurfaceAreaTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5], spacing=1.0)
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.air,
+                                              geometry=self.geo)
+        self.geo['pore.diameter'] = 1
+        self.geo['throat.area'] = 0.1
+
+    def test_sphere(self):
+        self.geo.add_model(propname='pore.area',
+                           model=mods.sphere,
+                           regen_mode='normal')
+        a = np.array([2.54159265, 2.64159265, 2.74159265, 2.84159265])
+        b = np.unique(self.geo['pore.area'])
+        assert np.allclose(a, b)
+
+    def test_cube(self):
+        self.geo.add_model(propname='pore.area',
+                           model=mods.cube,
+                           regen_mode='normal')
+        a = np.array([5.4, 5.5, 5.6, 5.7])
+        b = np.unique(self.geo['pore.area'])
+        assert np.allclose(a, b)
+
+
+if __name__ == '__main__':
+
+    t = PoreSurfaceAreaTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/geometry/PoreVolumeTest.py b/tests/unit/models/geometry/PoreVolumeTest.py
new file mode 100644
index 000000000..9722c9e36
--- /dev/null
+++ b/tests/unit/models/geometry/PoreVolumeTest.py
@@ -0,0 +1,44 @@
+import openpnm as op
+import openpnm.models.geometry.pore_volume as mods
+import numpy as np
+
+
+class PoreVolumeTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5], spacing=1.0)
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.air,
+                                              geometry=self.geo)
+        self.geo['pore.diameter'] = 1.0
+        self.geo['throat.area'] = 0.1
+
+    def test_sphere(self):
+        self.geo.add_model(propname='pore.volume',
+                           model=mods.sphere,
+                           regen_mode='normal')
+        a = np.array([0.52359878])
+        b = np.unique(self.geo['pore.volume'])
+        assert np.allclose(a, b)
+
+    def test_cube(self):
+        self.geo.add_model(propname='pore.volume',
+                           model=mods.cube,
+                           regen_mode='normal')
+        a = np.array([1.0])
+        b = np.unique(self.geo['pore.volume'])
+        assert np.allclose(a, b)
+
+
+if __name__ == '__main__':
+
+    t = PoreVolumeTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/geometry/ThroatAreaTest.py b/tests/unit/models/geometry/ThroatAreaTest.py
new file mode 100644
index 000000000..794891f38
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatAreaTest.py
@@ -0,0 +1,44 @@
+import openpnm as op
+import openpnm.models.geometry.throat_area as mods
+import numpy as np
+
+
+class ThroatSurfaceAreaTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5], spacing=1.0)
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.air,
+                                              geometry=self.geo)
+        self.geo['throat.length'] = 1
+        self.geo['throat.diameter'] = 0.1
+
+    def test_sphere(self):
+        self.geo.add_model(propname='throat.area',
+                           model=mods.cylinder,
+                           regen_mode='normal')
+        a = np.array([0.00785398])
+        b = np.unique(self.geo['throat.area'])
+        assert np.allclose(a, b)
+
+    def test_cube(self):
+        self.geo.add_model(propname='throat.area',
+                           model=mods.cuboid,
+                           regen_mode='normal')
+        a = np.array([0.01])
+        b = np.unique(self.geo['throat.area'])
+        assert np.allclose(a, b)
+
+
+if __name__ == '__main__':
+
+    t = ThroatSurfaceAreaTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/geometry/ThroatCentroidTest.py b/tests/unit/models/geometry/ThroatCentroidTest.py
new file mode 100644
index 000000000..c5008fef6
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatCentroidTest.py
@@ -0,0 +1,3 @@
+class ThroatCentroidTest:
+    def test_voronoi(self):
+        pass
diff --git a/tests/unit/models/geometry/ThroatDiameterTest.py b/tests/unit/models/geometry/ThroatDiameterTest.py
new file mode 100644
index 000000000..a0a8ea447
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatDiameterTest.py
@@ -0,0 +1,62 @@
+import openpnm as op
+import openpnm.models.geometry as gm
+import scipy as sp
+import scipy.stats as spst
+
+
+class ThroatSizeTest:
+
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[2, 2, 2])
+        self.net.add_boundary_pores()
+        Ps = self.net.pores('*boundary', mode='not')
+        Ts = self.net.throats('*boundary', mode='not')
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=Ps,
+                                               throats=Ts)
+        sp.random.RandomState(seed=0)
+        self.geo['throat.seed'] = sp.rand(self.geo.Nt)
+        BPs = self.net.pores('*boundary')
+        BTs = self.net.throats('*boundary')
+        self.boun = op.geometry.GenericGeometry(network=self.net,
+                                                pores=BPs,
+                                                throats=BTs)
+
+    def test_normal(self):
+        self.geo.add_model(propname='throat.diameter',
+                           model=gm.throat_size.normal,
+                           scale=0.01,
+                           loc=0.5,
+                           seeds='throat.seed')
+        assert 0.45 < sp.mean(self.geo['throat.diameter']) < 0.55
+        del self.geo['throat.diameter']
+
+    def test_weibull(self):
+        self.geo.add_model(propname='throat.diameter',
+                           model=gm.throat_size.weibull,
+                           shape=1.5,
+                           scale=0.0001,
+                           loc=0.001,
+                           seeds='throat.seed')
+        assert sp.amin(self.geo['throat.diameter']) > 0.001
+        del self.geo['throat.diameter']
+
+    def test_generic(self):
+        func = spst.gamma(a=2, loc=0.001, scale=0.0001)
+        self.geo.add_model(propname='throat.diameter',
+                           model=gm.throat_size.generic,
+                           func=func,
+                           seeds='throat.seed')
+        assert sp.amin(self.geo['throat.diameter']) > 0.001
+        del self.geo['throat.diameter']
+
+
+if __name__ == '__main__':
+
+    t = ThroatSizeTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/geometry/ThroatLengthTest.py b/tests/unit/models/geometry/ThroatLengthTest.py
new file mode 100644
index 000000000..0b214c856
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatLengthTest.py
@@ -0,0 +1,35 @@
+import openpnm as op
+import openpnm.models.geometry.throat_length as mods
+import numpy as np
+
+
+class ThroatStraightTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5], spacing=1.0)
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.air,
+                                              geometry=self.geo)
+        self.geo['pore.diameter'] = .05
+
+    def test_straight(self):
+        self.geo.add_model(propname='throat.length',
+                           model=mods.straight,
+                           regen_mode='normal')
+        a = np.array([0.95])
+        b = np.unique(self.geo['throat.length'])
+        assert np.allclose(a, b)
+
+
+if __name__ == '__main__':
+
+    t = ThroatStraightTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/geometry/ThroatNormalTest.py b/tests/unit/models/geometry/ThroatNormalTest.py
new file mode 100644
index 000000000..1204af316
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatNormalTest.py
@@ -0,0 +1,3 @@
+class ThroatNormalTest:
+    def test_voronoi(self):
+        pass
diff --git a/tests/unit/models/geometry/ThroatOffsetVerticesTest.py b/tests/unit/models/geometry/ThroatOffsetVerticesTest.py
new file mode 100644
index 000000000..6f82f0c10
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatOffsetVerticesTest.py
@@ -0,0 +1,3 @@
+class ThroatOffsetVerticesTest:
+    def test_voronoi(self):
+        pass
diff --git a/tests/unit/models/geometry/ThroatPerimeterTest.py b/tests/unit/models/geometry/ThroatPerimeterTest.py
new file mode 100644
index 000000000..c2234b570
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatPerimeterTest.py
@@ -0,0 +1,43 @@
+import openpnm as op
+import openpnm.models.geometry.throat_perimeter as mods
+import numpy as np
+
+
+class ThroatPerimeterTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5], spacing=1.0)
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.air,
+                                              geometry=self.geo)
+        self.geo['throat.diameter'] = 0.1
+
+    def test_cylinder(self):
+        self.geo.add_model(propname='throat.perimeter',
+                           model=mods.cylinder,
+                           regen_mode='normal')
+        a = np.array([0.31415927])
+        b = np.unique(self.geo['throat.perimeter'])
+        assert np.allclose(a, b)
+
+    def test_cube(self):
+        self.geo.add_model(propname='throat.perimeter',
+                           model=mods.cuboid,
+                           regen_mode='normal')
+        a = np.array([0.4])
+        b = np.unique(self.geo['throat.perimeter'])
+        assert np.allclose(a, b)
+
+
+if __name__ == '__main__':
+
+    t = ThroatPerimeterTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/geometry/ThroatSeedTest.py b/tests/unit/models/geometry/ThroatSeedTest.py
new file mode 100644
index 000000000..fa9b670b7
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatSeedTest.py
@@ -0,0 +1,47 @@
+import openpnm as op
+import scipy as sp
+import openpnm.models.geometry.throat_seed as mods
+
+
+class ThroatSeedTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5])
+        self.geo = op.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.add_model(propname='throat.seed',
+                           model=mods.random,
+                           seed=0,
+                           num_range=[0.1, 2])
+        self.geo.regenerate_models()
+        assert sp.amax(self.geo['throat.seed']) > 1.9
+        assert sp.amin(self.geo['throat.seed']) > 0.1
+
+    def test_neighbor(self):
+        self.geo.add_model(propname='throat.seed_max',
+                           model=mods.from_neighbor_pores,
+                           mode='max')
+        self.geo.add_model(propname='throat.seed_min',
+                           model=mods.from_neighbor_pores,
+                           mode='min')
+        self.geo.add_model(propname='throat.seed_mean',
+                           model=mods.from_neighbor_pores,
+                           mode='mean')
+        self.geo.regenerate_models()
+        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'])
+
+
+if __name__ == '__main__':
+
+    t = ThroatSeedTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/geometry/ThroatShapeFactorTest.py b/tests/unit/models/geometry/ThroatShapeFactorTest.py
new file mode 100644
index 000000000..52ab90eac
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatShapeFactorTest.py
@@ -0,0 +1,3 @@
+class ThroatShapeFactorTest:
+    def test_compactness(self):
+        pass
diff --git a/tests/unit/models/geometry/ThroatSurfaceAreaTest.py b/tests/unit/models/geometry/ThroatSurfaceAreaTest.py
new file mode 100644
index 000000000..fdbe327f8
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatSurfaceAreaTest.py
@@ -0,0 +1,45 @@
+import openpnm as op
+import scipy as sp
+
+
+class ThroatSurfaceAreaTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5])
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.geo['throat.diameter'] = sp.ones((self.geo.Nt, ))
+        self.geo['throat.length'] = sp.ones((self.geo.Nt, ))
+        self.geo['throat.perimeter'] = sp.ones((self.geo.Nt, ))
+
+    def test_cylinder(self):
+        f = op.models.geometry.throat_surface_area.cylinder
+        self.geo.add_model(propname='throat.surface_area',
+                           model=f)
+        self.geo.regenerate_models()
+        assert sp.all(self.geo['throat.surface_area'] == sp.pi)
+
+    def test_cuboid(self):
+        f = op.models.geometry.throat_surface_area.cuboid
+        self.geo.add_model(propname='throat.surface_area',
+                           model=f,)
+        self.geo.regenerate_models()
+        assert sp.all(self.geo['throat.surface_area'] == 4)
+
+    def test_extrusion(self):
+        f = op.models.geometry.throat_surface_area.extrusion
+        self.geo.add_model(propname='throat.surface_area',
+                           model=f)
+        self.geo.regenerate_models()
+        assert sp.all(self.geo['throat.surface_area'] == 1)
+
+
+if __name__ == '__main__':
+
+    t = ThroatSurfaceAreaTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/geometry/ThroatVectorTest.py b/tests/unit/models/geometry/ThroatVectorTest.py
new file mode 100644
index 000000000..d6878bf70
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatVectorTest.py
@@ -0,0 +1,3 @@
+class ThroatVectorTest:
+    def test_pore_to_pore(self):
+        pass
diff --git a/tests/unit/models/geometry/ThroatVolumeTest.py b/tests/unit/models/geometry/ThroatVolumeTest.py
new file mode 100644
index 000000000..773fdc087
--- /dev/null
+++ b/tests/unit/models/geometry/ThroatVolumeTest.py
@@ -0,0 +1,44 @@
+import openpnm as op
+import openpnm.models.geometry.throat_volume as mods
+import numpy as np
+
+
+class ThroatVolumeTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5], spacing=1.0)
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.air,
+                                              geometry=self.geo)
+        self.geo['throat.diameter'] = 0.1
+        self.geo['throat.length'] = 1.0
+
+    def test_cylinder(self):
+        self.geo.add_model(propname='throat.volume',
+                           model=mods.cylinder,
+                           regen_mode='normal')
+        a = np.array([0.00785398])
+        b = np.unique(self.geo['throat.volume'])
+        assert np.allclose(a, b)
+
+    def test_cube(self):
+        self.geo.add_model(propname='throat.volume',
+                           model=mods.cuboid,
+                           regen_mode='normal')
+        a = np.array([0.01])
+        b = np.unique(self.geo['throat.volume'])
+        assert np.allclose(a, b)
+
+
+if __name__ == '__main__':
+
+    t = ThroatVolumeTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/misc/MiscTest.py b/tests/unit/models/misc/MiscTest.py
new file mode 100644
index 000000000..c8cb0d767
--- /dev/null
+++ b/tests/unit/models/misc/MiscTest.py
@@ -0,0 +1,169 @@
+import openpnm as op
+import scipy as sp
+import openpnm.models.misc as mods
+
+
+class MiscTest:
+
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5])
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+
+    def test_constant(self):
+        self.geo.add_model(model=mods.constant,
+                           propname='pore.value',
+                           value=3.3)
+        assert sp.all(sp.unique(self.geo['pore.value']) == 3.3)
+
+    def test_product(self):
+        self.geo.add_model(model=mods.constant,
+                           propname='pore.value1',
+                           value=2)
+        self.geo.add_model(model=mods.constant,
+                           propname='pore.value2',
+                           value=4)
+        self.geo.add_model(model=mods.product,
+                           propname='pore.value3',
+                           props=['pore.value1', 'pore.value2'])
+        assert sp.all(sp.unique(self.geo['pore.value3']) == 8)
+
+    def test_scaled(self):
+        self.geo['pore.value4'] = 4
+        self.geo.add_model(model=mods.scaled,
+                           propname='pore.value5',
+                           prop='pore.value4',
+                           factor=2)
+        assert sp.all(sp.unique(self.geo['pore.value5']) == 8)
+
+    def test_linear(self):
+        self.geo['pore.value4'] = 4
+        self.geo.add_model(model=mods.linear,
+                           propname='pore.value6',
+                           prop='pore.value4',
+                           m=2, b=2)
+        assert sp.all(sp.unique(self.geo['pore.value6']) == 10)
+
+    def test_polynomial(self):
+        self.geo['pore.value4'] = 4
+        self.geo.add_model(model=mods.polynomial,
+                           propname='pore.value7',
+                           prop='pore.value4',
+                           a=[0, 2, 4, 6])
+        assert sp.all(sp.unique(self.geo['pore.value7']) == 456)
+
+    def test_random_no_seed(self):
+        self.geo.add_model(model=mods.random,
+                           propname='pore.seed',
+                           element='pore',
+                           seed=None)
+        temp1 = self.geo['pore.seed'].copy()
+        self.geo.regenerate_models()
+        temp2 = self.geo['pore.seed'].copy()
+        assert sp.all(~(temp1 == temp2))
+
+    def test_random_with_seed(self):
+        self.geo.add_model(model=mods.random,
+                           propname='pore.seed',
+                           element='pore',
+                           seed=0)
+        temp1 = self.geo['pore.seed'].copy()
+        self.geo.regenerate_models()
+        temp2 = self.geo['pore.seed'].copy()
+#        assert np.testing.assert_array_almost_equal_nulp(temp1, temp2)
+
+    def test_random_with_range(self):
+        self.geo.add_model(model=mods.random,
+                           propname='pore.seed',
+                           element='pore',
+                           num_range=[0.1, 0.9])
+        self.geo.regenerate_models()
+        assert sp.amax(self.geo['pore.seed']) <= 0.9
+        assert sp.amin(self.geo['pore.seed']) >= 0.1
+
+    def test_from_neighbor_throats_min(self):
+        catch = self.geo.pop('pore.seed', None)
+        catch = self.geo.models.pop('pore.seed', None)
+        catch = self.geo.models.pop('throat.seed', None)
+        self.geo['throat.seed'] = sp.rand(self.net.Nt,)
+        self.geo.add_model(model=mods.from_neighbor_throats,
+                           propname='pore.seed',
+                           throat_prop='throat.seed',
+                           mode='min')
+        assert sp.all(sp.in1d(self.geo['pore.seed'], self.geo['throat.seed']))
+        pmax = sp.amax(self.geo['pore.seed'])
+        tmax = sp.amax(self.geo['throat.seed'])
+        assert pmax <= tmax
+
+    def test_from_neighbor_throats_max(self):
+        catch = self.geo.pop('pore.seed', None)
+        catch = self.geo.models.pop('pore.seed', None)
+        catch = self.geo.models.pop('throat.seed', None)
+        self.geo['throat.seed'] = sp.rand(self.net.Nt,)
+        self.geo.add_model(model=mods.from_neighbor_throats,
+                           propname='pore.seed',
+                           throat_prop='throat.seed',
+                           mode='max')
+        assert sp.all(sp.in1d(self.geo['pore.seed'], self.geo['throat.seed']))
+        pmin = sp.amin(self.geo['pore.seed'])
+        tmin = sp.amin(self.geo['throat.seed'])
+        assert pmin >= tmin
+
+    def test_from_neighbor_throats_mean(self):
+        catch = self.geo.pop('pore.seed', None)
+        catch = self.geo.models.pop('pore.seed', None)
+        catch = self.geo.models.pop('throat.seed', None)
+        self.geo['throat.seed'] = sp.rand(self.net.Nt,)
+        self.geo.add_model(model=mods.from_neighbor_throats,
+                           propname='pore.seed',
+                           throat_prop='throat.seed',
+                           mode='mean')
+        tmax = sp.amax(self.geo['throat.seed'])
+        tmin = sp.amin(self.geo['throat.seed'])
+        assert sp.all(self.geo['pore.seed'] > tmin)
+        assert sp.all(self.geo['pore.seed'] < tmax)
+
+    def test_from_neighbor_pores_min(self):
+        self.geo.remove_model('throat.seed')
+        self.geo['pore.seed'] = sp.rand(self.net.Np,)
+        self.geo.add_model(model=mods.from_neighbor_pores,
+                           propname='throat.seed',
+                           pore_prop='pore.seed',
+                           mode='min')
+        P12 = self.net['throat.conns']
+        tseed = sp.amin(self.geo['pore.seed'][P12], axis=1)
+        assert sp.allclose(self.geo['throat.seed'], tseed)
+
+    def test_from_neighbor_pores_max(self):
+        self.geo.remove_model('throat.seed')
+        self.geo['pore.seed'] = sp.rand(self.net.Np,)
+        self.geo.add_model(model=mods.from_neighbor_pores,
+                           propname='throat.seed',
+                           pore_prop='pore.seed',
+                           mode='max')
+        P12 = self.net['throat.conns']
+        tseed = sp.amax(self.geo['pore.seed'][P12], axis=1)
+        assert sp.allclose(self.geo['throat.seed'], tseed)
+
+    def test_from_neighbor_pores_mean(self):
+        self.geo.remove_model('throat.seed')
+        self.geo['pore.seed'] = sp.rand(self.net.Np,)
+        self.geo.add_model(model=mods.from_neighbor_pores,
+                           propname='throat.seed',
+                           pore_prop='pore.seed',
+                           mode='mean')
+        P12 = self.net['throat.conns']
+        tseed = sp.mean(self.geo['pore.seed'][P12], axis=1)
+        assert sp.allclose(self.geo['throat.seed'], tseed)
+
+
+if __name__ == '__main__':
+
+    t = MiscTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/phase/ContactAngleTest.py b/tests/unit/models/phase/ContactAngleTest.py
new file mode 100644
index 000000000..c791c1cef
--- /dev/null
+++ b/tests/unit/models/phase/ContactAngleTest.py
@@ -0,0 +1,3 @@
+class ContactAngleTest:
+    def test_young(self):
+        pass
diff --git a/tests/unit/models/phase/DensityTest.py b/tests/unit/models/phase/DensityTest.py
new file mode 100644
index 000000000..0970bf517
--- /dev/null
+++ b/tests/unit/models/phase/DensityTest.py
@@ -0,0 +1,48 @@
+import openpnm as op
+import scipy as sp
+
+
+class DensityTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.phase = op.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):
+        # Liquid water
+        self.phase.add_model(propname='pore.density',
+                             model=op.models.phase.density.standard)
+        assert sp.allclose(self.phase['pore.density'], 999.702)
+
+    def test_ideal_gas(self):
+        # Water vapor
+        self.phase.add_model(propname='pore.density',
+                             model=op.models.phase.density.ideal_gas)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.density'], 0.73610248)
+
+    def test_water(self):
+        # Liquid water
+        self.phase.add_model(propname='pore.density',
+                             model=op.models.phase.density.water)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.density'], 996.9522)
+
+    def teardown_class(self):
+        del(self.phase)
+        del(self.net)
+
+
+if __name__ == '__main__':
+
+    t = DensityTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/phase/DiffusivityTest.py b/tests/unit/models/phase/DiffusivityTest.py
new file mode 100644
index 000000000..5eb0171a1
--- /dev/null
+++ b/tests/unit/models/phase/DiffusivityTest.py
@@ -0,0 +1,64 @@
+import openpnm as op
+import scipy as sp
+
+
+class DiffusivityTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.phase = op.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):
+        f = op.models.phase.diffusivity.fuller
+        self.phase.add_model(propname='pore.diffusivity',
+                             model=f,
+                             MA=0.032,
+                             MB=0.028,
+                             vA=16.6,
+                             vB=17.9)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.diffusivity'], 2.06754784e-05)
+
+    def test_fuller_scaling(self):
+        f = op.models.phase.diffusivity.fuller_scaling
+        self.phase.add_model(propname='pore.diffusivity',
+                             model=f,
+                             DABo=1.79712526e-05,
+                             Po=100000,
+                             To=273)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.diffusivity'], 2.06754784e-05)
+
+    def test_tyn_calus(self):
+        f = op.models.phase.diffusivity.tyn_calus
+        self.phase.add_model(propname='pore.diffusivity',
+                             model=f,
+                             VA=16.5,
+                             VB=17.9,
+                             sigma_A=1,
+                             sigma_B=1)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.diffusivity'], 9.84851806e-05)
+
+    def test_tyn_calus_scaling(self):
+        f = op.models.phase.diffusivity.tyn_calus_scaling
+        self.phase.add_model(propname='pore.diffusivity',
+                             model=f,
+                             DABo=5.26300839e-05,
+                             mu_o=3e-5,
+                             To=273)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.diffusivity'], 9.84851806e-05)
+
+
+if __name__ == '__main__':
+
+    t = DiffusivityTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/phase/ElectricalConductivityTest.py b/tests/unit/models/phase/ElectricalConductivityTest.py
new file mode 100644
index 000000000..cf9dc4b70
--- /dev/null
+++ b/tests/unit/models/phase/ElectricalConductivityTest.py
@@ -0,0 +1,27 @@
+import openpnm as op
+
+
+class ElectricalConductivityTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.phase = op.phases.GenericPhase(network=self.net)
+        self.phase['pore.intrinsic_conductivity'] = 1
+        self.phase['pore.volume_fraction'] = 0.5
+
+    def test_percolating_continua(self):
+        f = op.models.phase.electrical_conductivity.percolating_continua
+        self.phase.add_model(propname='pore.effective_conductivity',
+                             model=f,
+                             phi_crit=0.25,
+                             tau=2)
+
+
+if __name__ == '__main__':
+
+    t = ElectricalConductivityTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/phase/MolarDensityTest.py b/tests/unit/models/phase/MolarDensityTest.py
new file mode 100644
index 000000000..3254da562
--- /dev/null
+++ b/tests/unit/models/phase/MolarDensityTest.py
@@ -0,0 +1,43 @@
+import openpnm as op
+import scipy as sp
+
+
+class MolarDensityTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.phase = op.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):
+        f = op.models.phase.molar_density.standard
+        self.phase.add_model(propname='pore.molar_density', model=f)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.molar_density'], 40.8934707)
+
+    def test_ideal_gas(self):
+        f = op.models.phase.molar_density.ideal_gas
+        self.phase.add_model(propname='pore.molar_density', model=f)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.molar_density'], 40.8945824)
+
+    def test_vanderwaals(self):
+        f = op.models.phase.molar_density.vanderwaals
+        self.phase.add_model(propname='pore.molar_density', model=f)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.molar_density'], 40.92524916)
+
+
+if __name__ == '__main__':
+
+    t = MolarDensityTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/phase/MolarMassTest.py b/tests/unit/models/phase/MolarMassTest.py
new file mode 100644
index 000000000..6bc2a38fe
--- /dev/null
+++ b/tests/unit/models/phase/MolarMassTest.py
@@ -0,0 +1,3 @@
+class MolarMassTest:
+    def test_mixture(self):
+        pass
diff --git a/tests/unit/models/phase/SurfaceTensionTest.py b/tests/unit/models/phase/SurfaceTensionTest.py
new file mode 100644
index 000000000..abe3d6788
--- /dev/null
+++ b/tests/unit/models/phase/SurfaceTensionTest.py
@@ -0,0 +1,58 @@
+import openpnm as op
+import scipy as sp
+
+
+class SurfaceTensionTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.phase = op.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):
+        f = op.models.phase.surface_tension.water
+        self.phase.add_model(propname='pore.surface_tension',
+                             model=f)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.surface_tension'], 0.07199533)
+
+    def test_eotvos(self):
+        f = op.models.phase.surface_tension.eotvos
+        self.phase.add_model(propname='pore.surface_tension',
+                             model=f,
+                             k=0.000014)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.surface_tension'], 0.07112169)
+
+    def test_guggenheim_katayama(self):
+        f = op.models.phase.surface_tension.guggenheim_katayama
+        self.phase.add_model(propname='pore.surface_tension',
+                             model=f,
+                             K2=0.0000014,
+                             n=0.1)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.surface_tension'], 0.27582571)
+
+    def test_brock_bird_scaling(self):
+        f = op.models.phase.surface_tension.brock_bird_scaling
+        self.phase.add_model(propname='pore.surface_tension',
+                             model=f,
+                             sigma_o=0.0608,
+                             To=363)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.surface_tension'], 0.07820761)
+
+
+if __name__ == '__main__':
+
+    t = SurfaceTensionTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/phase/ThermalConductivityTest.py b/tests/unit/models/phase/ThermalConductivityTest.py
new file mode 100644
index 000000000..eb7a5df14
--- /dev/null
+++ b/tests/unit/models/phase/ThermalConductivityTest.py
@@ -0,0 +1,48 @@
+import openpnm as op
+import scipy as sp
+
+
+class ThermalConductivityTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.phase = op.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):
+        f = op.models.phase.thermal_conductivity.water
+        self.phase.add_model(propname='pore.thermal_conductivity',
+                             model=f)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.thermal_conductivity'], 0.61047611)
+
+    def test_chung(self):
+        f = op.models.phase.thermal_conductivity.chung
+        self.phase.add_model(propname='pore.thermal_conductivity',
+                             model=f)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.thermal_conductivity'], 0.62063913)
+
+    def test_sato(self):
+        f = op.models.phase.thermal_conductivity.sato
+        self.phase.add_model(propname='pore.thermal_conductivity',
+                             model=f)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.thermal_conductivity'], 0.29787023)
+
+
+if __name__ == '__main__':
+
+    t = ThermalConductivityTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/phase/VaporPressureTest.py b/tests/unit/models/phase/VaporPressureTest.py
new file mode 100644
index 000000000..73bc2b959
--- /dev/null
+++ b/tests/unit/models/phase/VaporPressureTest.py
@@ -0,0 +1,52 @@
+import openpnm as op
+import scipy as sp
+
+
+class VaporPressureTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.phase = op.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):
+        f = op.models.phase.vapor_pressure.antoine
+        self.phase.add_model(propname='pore.test',
+                             model=f,
+                             temperature='pore.temperature',
+                             A=8.088,
+                             B=1750.71,
+                             C=236.191)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.test'], 3607.850875)
+
+    def test_water(self):
+        f = op.models.phase.vapor_pressure.water
+        self.phase.add_model(propname='pore.test',
+                             model=f,
+                             temperature='pore.temperature',
+                             salinity='pore.salinity')
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.test'], 3536.01)
+
+    def test_water_no_salinity(self):
+        f = op.models.phase.vapor_pressure.water
+        del self.phase['pore.salinity']
+        self.phase.add_model(propname='pore.test',
+                             model=f,
+                             temperature='pore.temperature',
+                             salinity='pore.salinity')
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.test'], 3536.01)
+        self.phase['pore.salinity'] = sp.zeros((self.phase.Np,))
+
+
+if __name__ == '__main__':
+
+    t = VaporPressureTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/phase/ViscosityTest.py b/tests/unit/models/phase/ViscosityTest.py
new file mode 100644
index 000000000..5955ba673
--- /dev/null
+++ b/tests/unit/models/phase/ViscosityTest.py
@@ -0,0 +1,44 @@
+import openpnm as op
+import scipy as sp
+
+
+class ViscosityTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.phase = op.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):
+        self.phase.add_model(propname='pore.viscosity',
+                             model=op.models.phase.viscosity.water)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.viscosity'], 0.00089319)
+
+    def test_reynolds(self):
+        self.phase.add_model(propname='pore.viscosity',
+                             model=op.models.phase.viscosity.reynolds,
+                             u0=0.001,
+                             b=0.001)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.viscosity'], 0.00134716)
+
+    def test_chung(self):
+        self.phase.add_model(propname='pore.viscosity',
+                             model=op.models.phase.viscosity.chung)
+        self.phase.regenerate_models()
+        assert sp.allclose(self.phase['pore.viscosity'], 6.47289919e-05)
+
+
+if __name__ == '__main__':
+
+    t = ViscosityTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/physics/CapillaryPressureTest.py b/tests/unit/models/physics/CapillaryPressureTest.py
new file mode 100644
index 000000000..727bc608d
--- /dev/null
+++ b/tests/unit/models/physics/CapillaryPressureTest.py
@@ -0,0 +1,97 @@
+import openpnm as op
+import scipy as sp
+
+
+class CapillaryPressureTest:
+
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.geo['throat.diameter'] = 1
+        self.geo['pore.diameter'] = 1
+        self.water = op.phases.GenericPhase(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              geometry=self.geo,
+                                              phase=self.water)
+
+    def test_washburn_pore_values(self):
+        self.water['pore.surface_tension'] = 0.072
+        self.water['pore.contact_angle'] = 120
+        f = op.models.physics.capillary_pressure.washburn
+        self.phys.add_model(propname='pore.capillary_pressure',
+                            model=f,
+                            surface_tension='pore.surface_tension',
+                            contact_angle='pore.contact_angle',
+                            diameter='pore.diameter')
+        self.phys.add_model(propname='throat.capillary_pressure',
+                            model=f,
+                            surface_tension='pore.surface_tension',
+                            contact_angle='pore.contact_angle',
+                            diameter='throat.diameter')
+        self.phys.regenerate_models()
+        a = 0.14399999999999993
+        assert sp.allclose(self.water['pore.capillary_pressure'][0], a)
+        self.phys.remove_model('pore.capillary_pressure')
+
+    def test_washburn_throat_values(self):
+        self.water['throat.surface_tension'] = 0.072
+        self.water['throat.contact_angle'] = 120
+        f = op.models.physics.capillary_pressure.washburn
+        self.phys.add_model(propname='throat.capillary_pressure',
+                            model=f,
+                            surface_tension='throat.surface_tension',
+                            contact_angle='throat.contact_angle',
+                            diameter='throat.diameter')
+        self.phys.regenerate_models()
+        a = 0.14399999999999993
+        assert sp.allclose(self.water['throat.capillary_pressure'][0], a)
+        self.phys.remove_model('throat.capillary_pressure')
+
+    def test_purcell_pore_values(self):
+        self.water['pore.surface_tension'] = 0.072
+        self.water['pore.contact_angle'] = 120
+        f = op.models.physics.capillary_pressure.purcell
+        self.phys.add_model(propname='pore.capillary_pressure',
+                            model=f,
+                            r_toroid=0.1,
+                            surface_tension='pore.surface_tension',
+                            contact_angle='pore.contact_angle',
+                            diameter='pore.diameter')
+        self.phys.add_model(propname='throat.capillary_pressure',
+                            model=f,
+                            r_toroid=0.1,
+                            surface_tension='pore.surface_tension',
+                            contact_angle='pore.contact_angle',
+                            diameter='throat.diameter')
+        self.phys.regenerate_models()
+        a = 0.26206427646507374
+        assert sp.allclose(self.water['pore.capillary_pressure'][0], a)
+        self.phys.remove_model('pore.capillary_pressure')
+
+    def test_purcell_throat_values(self):
+        self.water['throat.surface_tension'] = 0.072
+        self.water['throat.contact_angle'] = 120
+        f = op.models.physics.capillary_pressure.purcell
+        self.phys.add_model(propname='throat.capillary_pressure',
+                            model=f,
+                            r_toroid=0.1,
+                            surface_tension='throat.surface_tension',
+                            contact_angle='throat.contact_angle',
+                            diameter='throat.diameter')
+        self.phys.regenerate_models()
+        a = 0.26206427646507374
+        assert sp.allclose(self.water['throat.capillary_pressure'][0], a)
+        self.phys.remove_model('throat.capillary_pressure')
+
+
+if __name__ == '__main__':
+
+    t = CapillaryPressureTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/physics/DiffusiveConductanceTest.py b/tests/unit/models/physics/DiffusiveConductanceTest.py
new file mode 100644
index 000000000..f9930f919
--- /dev/null
+++ b/tests/unit/models/physics/DiffusiveConductanceTest.py
@@ -0,0 +1,45 @@
+import openpnm as op
+import scipy as sp
+
+
+class DiffusiveConductanceTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5], spacing=1.0)
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.geo['pore.diameter'] = 1.0
+        self.geo['pore.area'] = 1.0
+        self.geo['throat.diameter'] = 1.0
+        self.geo['throat.length'] = 1e-9
+        self.geo['throat.area'] = 1
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.air,
+                                              geometry=self.geo)
+
+    def test_bulk_diffusion(self):
+        mod = op.models.physics.diffusive_conductance.bulk_diffusion
+        self.phys.add_model(propname='throat.conductance1',
+                            model=mod)
+        self.phys.regenerate_models()
+        assert sp.allclose(a=self.phys['throat.conductance1'][0],
+                           b=0.00084552)
+
+        self.phys.add_model(propname='throat.conductance2',
+                            model=mod,
+                            calc_pore_len=True)
+        self.phys.regenerate_models()
+        assert sp.allclose(a=self.phys['throat.conductance2'][0],
+                           b=0.00084552)
+
+
+if __name__ == '__main__':
+
+    t = DiffusiveConductanceTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/physics/ElectricalConductanceTest.py b/tests/unit/models/physics/ElectricalConductanceTest.py
new file mode 100644
index 000000000..6c3e465e2
--- /dev/null
+++ b/tests/unit/models/physics/ElectricalConductanceTest.py
@@ -0,0 +1,39 @@
+import openpnm as op
+import scipy as sp
+
+
+class ElectricalConductanceTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[4, 4, 4])
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.phase = op.phases.GenericPhase(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.phase,
+                                              geometry=self.geo)
+
+    def test_electrical_conductance(self):
+        self.phase['pore.conductivity'] = 1
+        self.geo['pore.area'] = 1
+        self.geo['pore.diameter'] = 1
+        self.geo['throat.area'] = 1
+        self.geo['throat.length'] = 0.0001
+        f = op.models.physics.electrical_conductance.series_resistors
+        self.phys.add_model(propname='throat.electrical_conductance',
+                            conductivity='pore.conductivity',
+                            model=f)
+        self.phys.regenerate_models()
+        a = 0.99990001
+        assert sp.allclose(self.phys['throat.electrical_conductance'][0], a)
+
+
+if __name__ == '__main__':
+
+    t = ElectricalConductanceTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/physics/GenericSourceTerm.py b/tests/unit/models/physics/GenericSourceTerm.py
new file mode 100644
index 000000000..782258e6d
--- /dev/null
+++ b/tests/unit/models/physics/GenericSourceTerm.py
@@ -0,0 +1,283 @@
+import openpnm as op
+import numpy as np
+import openpnm.models.physics as pm
+
+
+class GenericSourceTermTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5])
+        Ps = self.net.Ps
+        Ts = self.net.Ts
+        self.geo = op.geometry.GenericGeometry(network=self.net, pores=Ps,
+                                               throats=Ts)
+        self.phase = op.phases.GenericPhase(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.phase,
+                                              geometry=self.geo)
+        self.phys['throat.diffusive_conductance'] = 5e-8
+        self.phase['pore.mole_fraction'] = 0.
+        self.alg = op.algorithms.GenericTransport(network=self.net,
+                                                  phase=self.phase)
+        BC_pores = np.arange(20, 30)
+        self.S_pores = np.arange(55, 85)
+        self.alg.set_boundary_conditions(bctype='Dirichlet',
+                                         bcvalue=0.4,
+                                         pores=BC_pores)
+
+#    def test_linear(self):
+#        self.phys['pore.item1'] = 0.5e-11
+#        self.phys['pore.item2'] = 1.5e-12
+#        self.phys.add_model(propname='pore.source1',
+#                            model=pm.generic_source_term.linear,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            x='mole_fraction',
+#                            return_rate=False,
+#                            regen_mode='on_demand')
+#        self.phys.add_model(propname='pore.source2',
+#                            model=pm.generic_source_term.linear,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            x='mole_fraction',
+#                            return_rate=True,
+#                            regen_mode='on_demand')
+#        self.alg.set_source_term(source_name='pore.source1',
+#                                 pores=self.S_pores,
+#                                 mode='overwrite')
+#        self.alg.run(conductance='throat.diffusive_conductance',
+#                     quantity='pore.mole_fraction',
+#                     super_pore_conductance=None)
+#        self.alg.return_results()
+#        self.phys.regenerate(props='pore.source1')
+#        self.phys.regenerate(props='pore.source2')
+#        X = self.phase['pore.mole_fraction']
+#        r1 = np.round(np.sum(0.5e-11 * X[self.S_pores] + 1.5e-12), 20)
+#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+#        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
+#        assert r1 == r2
+#        assert r2 == -r3
+#
+#    def test_power_law(self):
+#        self.phys['pore.item1'] = 0.5e-12
+#        self.phys['pore.item2'] = 2.5
+#        self.phys['pore.item3'] = -1.4e-11
+#        self.phys.add_model(propname='pore.source1',
+#                            model=pm.generic_source_term.power_law,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            A3='pore.item3',
+#                            x='mole_fraction',
+#                            return_rate=False,
+#                            regen_mode='on_demand')
+#        self.phys.add_model(propname='pore.source2',
+#                            model=pm.generic_source_term.power_law,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            A3='pore.item3',
+#                            x='mole_fraction',
+#                            return_rate=True,
+#                            regen_mode='on_demand')
+#        self.alg.set_source_term(source_name='pore.source1',
+#                                 pores=self.S_pores,
+#                                 mode='overwrite')
+#        self.alg.run(conductance='throat.diffusive_conductance',
+#                     quantity='pore.mole_fraction',
+#                     super_pore_conductance=None)
+#        self.alg.return_results()
+#        self.phys.regenerate(props='pore.source1')
+#        self.phys.regenerate(props='pore.source2')
+#        X = self.phase['pore.mole_fraction']
+#        r1 = np.round(np.sum(0.5e-12 * X[self.S_pores] ** 2.5 - 1.4e-11), 20)
+#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+#        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
+#        assert r1 == r2
+#        assert r2 == -r3
+#
+#    def test_exponential(self):
+#        self.phys['pore.item1'] = 0.8e-11
+#        self.phys['pore.item2'] = 3
+#        self.phys['pore.item3'] = 0.5
+#        self.phys['pore.item4'] = 2
+#        self.phys['pore.item5'] = -0.34
+#        self.phys['pore.item6'] = 2e-14
+#        self.phys.add_model(propname='pore.source1',
+#                            model=pm.generic_source_term.exponential,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            A3='pore.item3',
+#                            A4='pore.item4',
+#                            A5='pore.item5',
+#                            A6='pore.item6',
+#                            x='mole_fraction',
+#                            return_rate=False,
+#                            regen_mode='on_demand')
+#        self.phys.add_model(propname='pore.source2',
+#                            model=pm.generic_source_term.exponential,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            A3='pore.item3',
+#                            A4='pore.item4',
+#                            A5='pore.item5',
+#                            A6='pore.item6',
+#                            x='mole_fraction',
+#                            return_rate=True,
+#                            regen_mode='on_demand')
+#        self.alg.set_source_term(source_name='pore.source1',
+#                                 pores=self.S_pores,
+#                                 mode='overwrite')
+#        self.alg.run(conductance='throat.diffusive_conductance',
+#                     quantity='pore.mole_fraction',
+#                     super_pore_conductance=None)
+#        self.alg.return_results()
+#        self.phys.regenerate(props='pore.source1')
+#        self.phys.regenerate(props='pore.source2')
+#        X = self.phase['pore.mole_fraction']
+#        r1 = np.round(np.sum(0.8e-11 * 3 ** (0.5 * X[self.S_pores] ** 2 -
+#                      0.34) + 2e-14), 20)
+#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+#        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
+#        assert r1 == r2
+#        assert r2 == -r3
+#
+#    def test_natural_exponential(self):
+#        self.phys['pore.item1'] = 0.8e-11
+#        self.phys['pore.item2'] = 0.5
+#        self.phys['pore.item3'] = 2
+#        self.phys['pore.item4'] = -0.34
+#        self.phys['pore.item5'] = 2e-14
+#        self.phys.add_model(propname='pore.source1',
+#                            model=pm.generic_source_term.natural_exponential,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            A3='pore.item3',
+#                            A4='pore.item4',
+#                            A5='pore.item5',
+#                            x='mole_fraction',
+#                            return_rate=False,
+#                            regen_mode='on_demand')
+#        self.phys.add_model(propname='pore.source2',
+#                            model=pm.generic_source_term.natural_exponential,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            A3='pore.item3',
+#                            A4='pore.item4',
+#                            A5='pore.item5',
+#                            x='mole_fraction',
+#                            return_rate=True,
+#                            regen_mode='on_demand')
+#        self.alg.set_source_term(source_name='pore.source1',
+#                                 pores=self.S_pores,
+#                                 mode='overwrite')
+#        self.alg.run(conductance='throat.diffusive_conductance',
+#                     quantity='pore.mole_fraction',
+#                     super_pore_conductance=None)
+#        self.alg.return_results()
+#        self.phys.regenerate(props='pore.source1')
+#        self.phys.regenerate(props='pore.source2')
+#        X = self.phase['pore.mole_fraction']
+#        r1 = np.round(np.sum(0.8e-11 * np.exp(0.5 * X[self.S_pores] ** 2 -
+#                      0.34) + 2e-14), 20)
+#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+#        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
+#        assert r1 == r2
+#        assert r2 == -r3
+#
+#    def test_logarithm(self):
+#        self.phys['pore.item1'] = 0.16e-13
+#        self.phys['pore.item2'] = 10
+#        self.phys['pore.item3'] = 4
+#        self.phys['pore.item4'] = 1.4
+#        self.phys['pore.item5'] = 0.133
+#        self.phys['pore.item6'] = -5.1e-13
+#        self.phys.add_model(propname='pore.source1',
+#                            model=pm.generic_source_term.logarithm,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            A3='pore.item3',
+#                            A4='pore.item4',
+#                            A5='pore.item5',
+#                            A6='pore.item6',
+#                            x='mole_fraction',
+#                            return_rate=False,
+#                            regen_mode='on_demand')
+#        self.phys.add_model(propname='pore.source2',
+#                            model=pm.generic_source_term.logarithm,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            A3='pore.item3',
+#                            A4='pore.item4',
+#                            A5='pore.item5',
+#                            A6='pore.item6',
+#                            x='mole_fraction',
+#                            return_rate=True,
+#                            regen_mode='on_demand')
+#        self.alg.set_source_term(source_name='pore.source1',
+#                                 pores=self.S_pores,
+#                                 mode='overwrite')
+#        self.alg.run(conductance='throat.diffusive_conductance',
+#                     quantity='pore.mole_fraction',
+#                     super_pore_conductance=None)
+#        self.alg.return_results()
+#        self.phys.regenerate(props='pore.source1')
+#        self.phys.regenerate(props='pore.source2')
+#        X = self.phase['pore.mole_fraction']
+#        r1 = np.round(np.sum(0.16e-13 * np.log(4 * X[self.S_pores] ** (1.4) +
+#                             0.133) / np.log(10) - 5.1e-13), 20)
+#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+#        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
+#        assert r1 == r2
+#        assert r2 == -r3
+#
+#    def test_natural_logarithm(self):
+#        self.phys['pore.item1'] = 0.16e-14
+#        self.phys['pore.item2'] = 4
+#        self.phys['pore.item3'] = 1.4
+#        self.phys['pore.item4'] = 0.133
+#        self.phys['pore.item5'] = -5.1e-14
+#        self.phys.add_model(propname='pore.source1',
+#                            model=pm.generic_source_term.natural_logarithm,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            A3='pore.item3',
+#                            A4='pore.item4',
+#                            A5='pore.item5',
+#                            x='mole_fraction',
+#                            return_rate=False,
+#                            regen_mode='on_demand')
+#        self.phys.add_model(propname='pore.source2',
+#                            model=pm.generic_source_term.natural_logarithm,
+#                            A1='pore.item1',
+#                            A2='pore.item2',
+#                            A3='pore.item3',
+#                            A4='pore.item4',
+#                            A5='pore.item5',
+#                            x='mole_fraction',
+#                            return_rate=True,
+#                            regen_mode='on_demand')
+#        self.alg.set_source_term(source_name='pore.source1',
+#                                 pores=self.S_pores,
+#                                 mode='overwrite')
+#        self.alg.run(conductance='throat.diffusive_conductance',
+#                     quantity='pore.mole_fraction',
+#                     super_pore_conductance=None)
+#        self.alg.return_results()
+#        self.phys.regenerate(props='pore.source1')
+#        self.phys.regenerate(props='pore.source2')
+#        X = self.phase['pore.mole_fraction']
+#        r1 = np.round(np.sum(0.16e-14 * np.log(4 * X[self.S_pores] ** (1.4) +
+#                             0.133) - 5.1e-14), 20)
+#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+#        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
+#        assert r1 == r2
+#        assert r2 == -r3
+
+
+if __name__ == '__main__':
+
+    t = GenericSourceTermTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/physics/HydraulicConductanceTest.py b/tests/unit/models/physics/HydraulicConductanceTest.py
new file mode 100644
index 000000000..f964a3886
--- /dev/null
+++ b/tests/unit/models/physics/HydraulicConductanceTest.py
@@ -0,0 +1,45 @@
+import openpnm as op
+import scipy as sp
+
+
+class HydraulicConductanceTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5], spacing=1.0)
+        self.geo = op.geometry.GenericGeometry(network=self.net,
+                                               pores=self.net.Ps,
+                                               throats=self.net.Ts)
+        self.geo['pore.diameter'] = 1.0
+        self.geo['throat.diameter'] = 1.0
+        self.geo['throat.length'] = 1.0e-9
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.air,
+                                              geometry=self.geo)
+
+    def teardown_class(self):
+        mgr = op.core.Workspace()
+        mgr.clear()
+
+    def test_hagen_poiseuille(self):
+        mod = op.models.physics.hydraulic_conductance.hagen_poiseuille
+        self.phys.add_model(propname='throat.conductance1',
+                            model=mod)
+        assert sp.allclose(a=self.phys['throat.conductance1'][0],
+                           b=1330.68207684)
+
+        self.phys.add_model(propname='throat.conductance2',
+                            model=mod,
+                            calc_pore_len=True)
+        assert sp.allclose(a=self.phys['throat.conductance2'][0],
+                           b=1330.68207684)
+
+
+if __name__ == '__main__':
+
+    t = HydraulicConductanceTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/physics/MultiPhaseTest.py b/tests/unit/models/physics/MultiPhaseTest.py
new file mode 100644
index 000000000..579aa572b
--- /dev/null
+++ b/tests/unit/models/physics/MultiPhaseTest.py
@@ -0,0 +1,80 @@
+import openpnm as op
+import scipy as sp
+import openpnm.models.physics as pm
+
+
+class MultiPhaseTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        Ps = self.net.Ps
+        Ts = self.net.Ts
+        self.geo = op.geometry.GenericGeometry(network=self.net, pores=Ps,
+                                               throats=Ts)
+        self.phase = op.phases.GenericPhase(network=self.net)
+        self.phase['pore.occupancy'] = sp.ones(self.net.Np)
+        self.phase['throat.occupancy'] = sp.ones(self.net.Nt)
+        self.phase['pore.occupancy'] = 0.
+        self.phase['throat.occupancy'] = 0.
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.phase,
+                                              geometry=self.geo)
+        self.phys['throat.capillary_pressure'] = 5000
+        self.phys['throat.diffusive_conductance'] = 2
+
+#    def test_conduit_conductance_strict(self):
+#        self.phase['pore.occupancy'][[19, 25]] = 0
+#        t1 = self.net.Ts[self.phase['throat.occupancy'] == 0]
+#        t2 = self.net.Ts[~sp.in1d(self.net.Ts, t1)]
+#        self.phys.add_model(propname='throat.cond_conductance',
+#                            throat_conductance='throat.diffusive_conductance',
+#                            model=pm.multiphase.conduit_conductance,
+#                            mode='strict', factor=0)
+#        assert ~sp.all(self.phase['throat.cond_conductance'] == 0)
+#        assert sp.all(self.phase['throat.cond_conductance'][t1] == 0)
+#        assert ~sp.all(self.phase['throat.cond_conductance'][t2] != 0)
+#
+#    def test_conduit_conductance_medium(self):
+#        self.phase['pore.occupancy'][[19, 25]] = 0
+#        self.phase['pore.occupancy'][6] = 0
+#        t1 = self.net.Ts[self.phase['throat.occupancy'] == 0.]
+#        t2 = self.net.find_connecting_throat(6, 7)[0]
+#        t3 = self.net.Ts[~sp.in1d(self.net.Ts, sp.concatenate((t1, t2)))]
+#        self.phys.add_model(propname='throat.cond_conductance',
+#                            throat_conductance='throat.diffusive_conductance',
+#                            model=pm.multiphase.conduit_conductance,
+#                            mode='medium', factor=0)
+#        assert sp.all(self.phase['throat.cond_conductance'][t1] == 0)
+#        assert self.phase['throat.cond_conductance'][t2] == 0
+#        assert sp.all(self.phase['throat.cond_conductance'][t3] != 0)
+#
+#    def test_conduit_conductance_loose(self):
+#        self.phase['pore.occupancy'][[19, 20]] = 0
+#        t1 = self.net.Ts[self.phase['throat.occupancy'] == 0]
+#        t2 = self.net.Ts[~sp.in1d(self.net.Ts, t1)]
+#        self.phys.add_model(propname='throat.cond_conductance',
+#                            throat_conductance='throat.diffusive_conductance',
+#                            model=pm.multiphase.conduit_conductance,
+#                            mode='loose', factor=0)
+#        assert sp.all(self.phase['throat.cond_conductance'][t1] == 0)
+#        assert sp.all(self.phase['throat.cond_conductance'][t2] != 0)
+#
+#    def test_late_throat_filling(self):
+#        pass
+#
+#    def test_late_pore_filling(self):
+#        self.phys.add_model(propname='pore.late_p',
+#                            model=pm.multiphase.late_pore_filling,
+#                            Pc=5500)
+#        p = self.phase['pore.occupancy'] > 0
+#        assert sp.allclose(self.phase['pore.late_p'][p], 0.84973704)
+
+
+if __name__ == '__main__':
+
+    t = MultiPhaseTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/phases/AirTest.py b/tests/unit/phases/AirTest.py
new file mode 100644
index 000000000..301fd1456
--- /dev/null
+++ b/tests/unit/phases/AirTest.py
@@ -0,0 +1,3 @@
+class AirTest:
+    def test_generate(self):
+        pass
diff --git a/tests/unit/phase/GenericPhaseTest.py b/tests/unit/phases/GenericPhaseTest.py
similarity index 100%
rename from tests/unit/phase/GenericPhaseTest.py
rename to tests/unit/phases/GenericPhaseTest.py
diff --git a/tests/unit/phases/MercuryTest.py b/tests/unit/phases/MercuryTest.py
new file mode 100644
index 000000000..13bd694cf
--- /dev/null
+++ b/tests/unit/phases/MercuryTest.py
@@ -0,0 +1,3 @@
+class MercuryTest:
+    def test_generate(self):
+        pass
diff --git a/tests/unit/phases/WaterTest.py b/tests/unit/phases/WaterTest.py
new file mode 100644
index 000000000..d13ece022
--- /dev/null
+++ b/tests/unit/phases/WaterTest.py
@@ -0,0 +1,30 @@
+import openpnm as op
+import numpy.testing as nptest
+
+
+class WaterTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5])
+        self.water = op.phases.Water(network=self.net)
+
+    def test_change_temperature(self):
+        assert self.water['pore.temperature'][0] == 298.0
+        self.water.regenerate_models()
+        mu = self.water['pore.viscosity'][0]
+        nptest.assert_array_almost_equal_nulp(mu, 0.000893190947459112)
+        self.water['pore.temperature'][0] = 333.0
+        assert self.water['pore.temperature'][0] == 333.0
+        self.water.regenerate_models()
+        mu = self.water['pore.viscosity'][0]
+        nptest.assert_array_almost_equal_nulp(mu, 0.00046734380929754017)
+
+
+if __name__ == '__main__':
+
+    t = WaterTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()