Repository: PMEAL/OpenPNM
Pull request #1166 (2018-12-19 03:31:30)

Improving docs
I just merged this into master by accident :-(  I've reversed it, and now am doing the PR on the right branch.


Patch:
diff --git a/.gitignore b/.gitignore
index fafcb1de9..2f988de0e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,5 @@
 .spyproject/
 .pytest_cache/
 .vscode/
+
+docs/_build/
diff --git a/docs/getting_started/index.rst b/docs/getting_started/index.rst
index 450ca53ee..9f3f38dda 100644
--- a/docs/getting_started/index.rst
+++ b/docs/getting_started/index.rst
@@ -11,3 +11,5 @@ Contents:
     installation.rst
     quick_start.rst
     examples.rst
+    visualize_in_paraview.rst
+    migrating_to_V2.rst
diff --git a/docs/getting_started/migrating_to_V2.rst b/docs/getting_started/migrating_to_V2.rst
new file mode 100644
index 000000000..63d56bbf2
--- /dev/null
+++ b/docs/getting_started/migrating_to_V2.rst
@@ -0,0 +1,204 @@
+.. _migrating_to_V2:
+
+================================================================================
+Migrating from V1 to V2
+================================================================================
+
+--------------------------------------------------------------------------------
+Renaming Modules, and Reorganizing
+--------------------------------------------------------------------------------
+
+One of the main objectives of V2 was to 'reorganize' the package.  When we started building OpenPNM in 2012 we were new to Python and made some errors in the layout of the code that we fixed in V2.  Specifically, V2 fully adopts the pep8 naming conventions, which suggests that modules be named with lowercase while classes be uppercase.
+
+.. code-block:: python
+
+    # Old way
+    import OpenPNM as op
+    pn = op.Network.Cubic(shape=[10, 10, 10])
+
+    # New way
+    import openpnm as op
+    pn = op.network.Cubic(shape=[10, 10, 10])
+
+Note that the package itself was renamed to lowercase, as well as all the module names.  The tables below give an overview of the package comparing V1 and V2:
+
+The package itself was renamed using all lowercase letters:
+
++--------------------------------------+---------------------------------------+
+| Version 1                            | Version 2                             |
++======================================+=======================================+
+| OpenPNM                              | openpnm                               |
++--------------------------------------+---------------------------------------+
+
+All the submodules were renamed using lowercase letters:
+
++--------------------------------------+---------------------------------------+
+| Version 1                            | Version 2                             |
++======================================+=======================================+
+| Network                              | network                               |
++--------------------------------------+---------------------------------------+
+| Phases                               | phases                                |
++--------------------------------------+---------------------------------------+
+| Geometry                             | geometry                              |
++--------------------------------------+---------------------------------------+
+| Physics                              | physics                               |
++--------------------------------------+---------------------------------------+
+| Algorithms                           | algorithms                            |
++--------------------------------------+---------------------------------------+
+| Utilities                            | utils                                 |
++--------------------------------------+---------------------------------------+
+
+The terms Base and Core were switched; the Base module is now called core and the Core class is now Base.  This is more consistent with other packages, where all the 'core' classes are stored in a core module, and the 'Base' class is the main class from which all other classes inherit (e.g. Cubic <-- GenericNetwork <-- Base):
+
++--------------------------------------+---------------------------------------+
+| Version 1                            | Version 2                             |
++======================================+=======================================+
+| Base                                 | core                                  |
++--------------------------------------+---------------------------------------+
+| Base.Core                            | core.Base                             |
++--------------------------------------+---------------------------------------+
+
+Pore-scale models are now stored in a top-level module to improve findability:
+
++--------------------------------------+---------------------------------------+
+| Version 1                            | Version 2                             |
++======================================+=======================================+
+| Phases.models                        | models.phases                         |
++--------------------------------------+---------------------------------------+
+| Geometry.models                      | models.geometry                       |
++--------------------------------------+---------------------------------------+
+| Physics.models                       | models.physics                        |
++--------------------------------------+---------------------------------------+
+
+A new module was added called topotools to house all functions pertaining to manipulating and mangling the topology of the networks.  All functions formerly found under models and/or tools have been moved:
+
++--------------------------------------+---------------------------------------+
+| Version 1                            | Version 2                             |
++======================================+=======================================+
+| Network.models                       | *removed*                             |
++--------------------------------------+---------------------------------------+
+| Network.tools                        | topotools                             |
++--------------------------------------+---------------------------------------+
+
+Finally, the Postprocessing module was removed and the few functions it had worth keeping have been relocated:
+
++--------------------------------------+---------------------------------------+
+| Version 1                            | Version 2                             |
++======================================+=======================================+
+| Postprocessing                       | *removed*                             |
++--------------------------------------+---------------------------------------+
+
+--------------------------------------------------------------------------------
+Cleaning up Classes
+--------------------------------------------------------------------------------
+
+In V2, we have made an effort to reduce the number of methods on each class, particularly GenericNetwork. Many of the functions that were previously found on the various networks (trim, extend) are now stored in topotools, and instead of being class methods are simple functions that accept a network as an argument an operate on the network 'in-place' meaning they do not return a network, but the received network is changed.
+
+For instance:
+
+.. code-block:: python
+
+    # Old way:
+    pn = op.Network.Cubic(shape=[10, 10, 10])
+    pn.trim(pores=[1, 4, 7])
+
+    # New way:
+    pn = op.network.Cubic(shape=[10, 10, 10])
+    op.topotools.trim(network=pn, pores=[1, 4, 7])
+
+--------------------------------------------------------------------------------
+Instantiating Physics Objects
+--------------------------------------------------------------------------------
+
+In version 2 we have moved to a more strict assignment of Physics objects to different regions, so instead of assigning them to specific pores and throats (as is still done for Geometry objects), each Physics is associated 1-to-1 with a Geometry.  For example:
+
+.. code-block:: python
+
+    # Old way
+    import OpenPNM as op
+    pn = op.Network.Cubic(shape=[10, 10, 10])
+    geom = op.Geometry.GenericGeometry(network=pn, pores=pn.Ps, throats=pn.Ts)
+    water = op.Phases.Water(network=pn)
+    phys = op.Physics.GenericPhysics(network=pn, phase=water, pores=pn.Ps, throats=pn.Ts)
+
+    # New way:
+    import openpnm as op
+    pn = op.network.Cubic(shape=[10, 10, 10])
+    geom = op.geometry.GenericGeometry(network=pn, pores=pn.Ps, throats=pn.Ts)
+    water = op.phases.Water(network=pn)
+    phys = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom)
+
+The changes to the Physics instantiation were motivated by a new way organize the simulations called "The Grid".  The grid is a 2D table of objects where each column header is a Phase, each row label is a Geometry, and each row-column intersection contains the corresponding Physics.  By forcing a 1-to-1 associated with Geometry and Physics, we maintain the grid structure which is helpful for object look-ups.  The grid can be seen using the new Project object as follows:
+
+.. code-block:: python
+
+    import openpnm as op
+    pn = op.network.Cubic(shape=[10, 10, 10])
+    geom = op.geometry.GenericGeometry(network=pn, pores=pn.Ps, throats=pn.Ts)
+    water = op.phases.Water(network=pn)
+    phys = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom)
+    proj = pn.project  # (or geom.project, or water.project or phys.project)
+    print(proj.grid)
+
+Which will produce the following:
+
+::
+
+    ――――――――――――――――――――――――――――――――
+    |    net_01    |   phase_01    |
+    ――――――――――――――――――――――――――――――――
+    |    geo_01    |    phys_01    |
+    ――――――――――――――――――――――――――――――――
+
+--------------------------------------------------------------------------------
+Object Lookups
+--------------------------------------------------------------------------------
+
+Object look-ups have changed substantially.  In V1 every object possessed a handle to it's associated objects, such that ``water.network`` would actually contain a handle directly to the network object. This association was causing problems such as memory leaks, large objects with circular references, and complicated object deletion.  In V2 we have moved away from this completely, and now all object look-ups are done with a "top down" approach.  This means that you must ask the Project object to find the associated object.  For example:
+
+.. code-block:: python
+
+    # Old way
+    phys.phases()
+
+    # New way
+    proj = phys.project
+    proj.find_phase(physics=phys)
+
+The actual looking up is done by checking the labels on each object.  When an object such as ``phys`` was is created, it is given a name, such as ``phys_01``, and the Phase with which it is associated is updated to have the labels ``'pore.phys_01'`` and ``'throat.phys_01'`` with ``True`` values where the ``phys`` object was assigned.  This 'top-down' object look up will search each Phase object until it finds the label corresponding to ``phys``, and once it does, the associated Phase is found.  As similar approach is done to find which geometry is associated with which physics and so forth.  
+
+--------------------------------------------------------------------------------
+Dealing with Pore-Scale Models
+--------------------------------------------------------------------------------
+
+Pore-scale models are one of the more confusing aspects for new users of OpenPNM, so we have attempted to simplify them in several ways.  Each object still has a ``models`` attribute which is a dictionary containing all the model definitions, but this dictionary, called the ModelsDict, is much simpler to work with.
+
+1.  Models are now run in order of dependencies, so if a user adds a model for 'pore.volume' before adding 'pore.diameter', OpenPNM will detect this and run the 'pore.diameter' model first so the values are available to subsequent models.  In the event that OpenPNM cannot run models in the correct order (e.g. a needed model is missing), it will report a warning, and can be remedied by adding the missing model and calling ``regnerate_models``.
+
+2.  ``add_model`` and ``regnerate_models`` are now methods of the Base class rather than being methods of the models dict.  This was done to increase their visibility and make them easier to call.  For example:
+
+.. code-block:: python
+
+    # Old way
+    geom.models.add(<args go here>)
+    geom.models.regnerate()
+
+    # New way
+    geom.add_model(<args go here>)
+    geom.regnerate_models()
+
+3.  Models themselves are now much simpler and more flexible.  Any given model can be assigned to any object. The aim where was to allow, for instance, geometrical models like 'pore.diameter' to be assigned directly to a Network object.  The point of Geometry objects is to allow complex, heterogeneous materials (e.g. layers with different pore sizes), but this is not always needed, so users can bypass this feature.  Having no Geometry objects precludes the use of Physics objects (see previous section on instantiating Physics objects), which means that pore-scale physics models (e.g. capillary pressure) must be assigned to Phase objects, which is also possible under this new flexible scheme.  In all pore scale models, the object to which the model is attached is referred to as the ``target``.
+
+--------------------------------------------------------------------------------
+Other Breaking Changes
+--------------------------------------------------------------------------------
+
+1. Neumann and Dirichlet boundary conditions have been renamed to Rate and Value to be more descriptive.  Neumann in particular was a problem since we were not actually specifying the gradient in the BC, but the gradient multiplied by the conductance, hence the move to Rate.  Dirichlet was renamed to Value was for consistency with Rate.
+
+2. The behavior of the various lookup methods (pores, find_neighbor_pores) now all use the same ``mode`` keywords. The new keyswords are taken for logic theory and include 'and', 'or', 'xor', 'xnor'.  The functions also accept synonyms from set theory ('intersection', 'union', 'exclusive_or', while xnor has no equivalent).
+
+.. note:: The meaning of intersection has changed
+
+  **One very important change** is the meaning of ``mode='intersection'`` in the ``find_neighbor_pores`` method.  In version 1 this was incorrectly being used to find neighbors that shared *more than one* input pore, while in version 2 this means pores that share *all the* input pores.  The mode 'xnor' replaces the old 'intersection'.  This change needed to be made, but is problematic since 'intersection' is still an accepted mode but returns a different result
+
+3. Every object now has a ``settings`` attribute, which is a dictionary of key:value pairs.  This was added so that we could stop storing various flags and values as attributes (sometimes hidden) on objects, and keep all such information in one place.  All Algorithm objects in particular use this settings dictionary, and the ``setup`` method essentially just passes any received arguments onto the ``settings`` dict.
diff --git a/docs/getting_started/quick_start.rst b/docs/getting_started/quick_start.rst
index ca6669347..0658fed25 100644
--- a/docs/getting_started/quick_start.rst
+++ b/docs/getting_started/quick_start.rst
@@ -27,7 +27,7 @@ The resulting network can be quickly visualized using ``plot_coordinates`` and
     :width: 800px
     :align: center
 
-The ``network`` module has a number of network types to chose from, and they can be found on the :ref:`network_index` page.
+The ``network`` module has a number of network types to chose from, and they can be found on the :ref:`network_index` page. You can also import networks from various outside sources, include networks that have been extracted from tomographic images using our `PoreSpy Package <http:\\porespy.com>`_.
 
 --------------------------------------------------------------------------------
 Adding Geometrical Properties
@@ -47,13 +47,13 @@ In this case the `StickAndBall` class was used, which has preset pore-scale mode
 
     (1) The Network (``pn``) was passed to the Geometry class as an argument, so that ``geo`` knows which network it's associated with.
 
-    (2) the Geometry was assigned to specified pores (``Ps``) and throats (``Ts``), in this case it was 'all' of them but it's possible to use several different Geometry objects for different subsets of the domain.  This explained in more detail here XXX.
+    (2) the Geometry was assigned to specified pores (``Ps``) and throats (``Ts``), in this case it was 'all' of them but it's possible to use several different Geometry objects for different subsets of the domain.
 
 --------------------------------------------------------------------------------
 Creating Phases
 --------------------------------------------------------------------------------
 
-Phases must created to calculate the thermophysical properties of the fluids (and solids) used in the simulations:
+Phases must be created to calculate the thermophysical properties of the fluids (and solids) used in the simulations:
 
 .. code-block:: python
 
@@ -145,11 +145,17 @@ Similarly for the permeability calculation:
     >>> perm.set_value_BC(pores=pn.pores('right'), values=0)
     >>> # perm.run()
 
-The above code solves for the pressure in each pore and stores the result as ``perm['pore.pressure']``.  To find the permeability of the network, there is a ``calc_eff_permeability`` method on the StokeFlow class:
+The above code solves for the pressure in each pore and stores the result as ``perm['pore.pressure']``.  To find the permeability of the network, there is a ``calc_eff_permeability`` method on the StokeFlow class.  Start by telling the algorithm the area and length of the domain (unfortunately there is no sure way to get these accurately, though values will be guesstimated if not provided)
 
     >>> perm.domain_area = (10*0.0001)**2
     >>> perm.domain_length = (10*0.0001)
-    >>> # K = perm.calc_eff_permeability()
+
+Finally to run the algorithm and calculate network permeability use:
+
+.. code-block:: python
+
+    perm.run()
+    K = perm.calc_effective_permeability()
 
 .. note::
 
diff --git a/docs/getting_started/visualize_in_paraview.rst b/docs/getting_started/visualize_in_paraview.rst
new file mode 100644
index 000000000..f27d9a994
--- /dev/null
+++ b/docs/getting_started/visualize_in_paraview.rst
@@ -0,0 +1,77 @@
+.. _visualize:
+
+================================================================================
+Using Paraview to Visualize Output
+================================================================================
+
+This tutorial will show step by step how to plot the networks into Paraview. This post-processing software is open-source and can be downloaded from the `Paraview Website <http://www.paraview.org/>`_.
+
+
+--------------------------------------------------------------------------------
+Export a VTK File from OpenPNM
+--------------------------------------------------------------------------------
+
+Assume the following project has been created and you wish to visualize it:
+
+.. code-block:: python
+
+    >>> import scipy as sp
+    >>> import openpnm as op
+    >>> pn = op.network.Cubic(shape=[5, 5, 5])
+    >>> geo = op.geometry.StickAndBall(network=pn, pores=pn.Ps, throats=pn.Ts)
+    >>> water = op.phases.Water(network=pn)
+    >>> phys = op.physics.Standard(network=pn, phase=water, geometry=geo)
+
+Now we can export it into a variety of formats using the ``io`` module, but the
+VTK format is easiest for Paraview.
+
+.. note::
+
+    The XDMF format is actually must better for using in Paraview if the network
+    is of appreciable size since it's performance is far better, but it is
+    slightly more complicated since it splits the data into 2 files.  The
+    information in this tutorial applies equally well to both formats.
+
+Now that the project is setup we can export it:
+
+.. code-block:: python
+
+    >>> op.io.VTK.save(network=pn, phases=water, filename='test_file')
+
+This will create a file in your current working directory called 'test_file.vtp'.  There are several things to note about this line.  Firstly,  the 'vtp' file extension is used rather than 'vtk'. The last letter indicates the type of data in the file, where 'p' indicates PolyData.  The VTK accronym stands for Visualization Tool Kit, which is a general reference to all similar files.  Secondly, current working directory can be found from the interactive python terminal by importing the ``os`` package, and typing ``os.get_cwd()``.  Finally, sending Phase data to the VTK file is optional, and you must specify which phase data to write.
+
+--------------------------------------------------------------------------------
+Open the VTK file created by OpenPNM
+--------------------------------------------------------------------------------
+
+Once the OpenPNM script (run_script.py) have been run, a output file named net.vtp is created. This file can be opened in Paraview and contains all the parameters computed by OpenPNM. A summary of them is seen under the Properties bar. Once the file is open, click on Apply, this action should be done after each modification.
+
+1. Open the 'test_file.vtp' file using the normal GUI functions in Paraview:
+
+.. image:: http://i.imgur.com/gmPVmRM.png
+
+2. Click and apply, the network should be seen in the view port looking somthing like this:
+
+.. image:: http://i.imgur.com/wFinDmX.png
+
+--------------------------------------------------------------------------------
+Plot the Dore Data
+--------------------------------------------------------------------------------
+
+To visualize the pore data, we need to add some glyphs. First click on the Glyph button in the tool bar. Then, you can plot the pore data as spheres, where their size and/or their color can be mapped to some variables. In the images below spherical glyphs are assigned to the pores, where the diameter is linked to the pore diameters and the color to the concentration. Clicking on the Apply button renders these settings.
+
+.. image:: http://i.imgur.com/bnFuH3r.png
+
+--------------------------------------------------------------------------------
+Plot the Throat Data
+--------------------------------------------------------------------------------
+
+.. note::
+
+  Plotting throats as nicely oriented cylinders with controlled diameter and color remains an annoyance, although we are working on getting the necessary information into the VTK file for throats be added as simply as adding cylindrical glyphs.
+
+To visualize the throat data (like diameters or molar flow rates) we need to set up the following filters in Paraview. First, use the Shrink filter and set it up to 1. Then, the cell data needs to be transposed to the data point with CellDatatoPointdata filter. Then extract the surface with the filter ExtractSurface. Finally, the data may be plotted as tube by using the Tube filter. As previously for the pore data, either the Tube radius or color can be linked to throat data.
+
+Throat data plotted with tubes radius proportional to throat diameter and tubes color linked to the throat mole rate:
+
+.. image:: http://i.imgur.com/SX5YeVj.png
diff --git a/docs/modules/algorithms/advection_diffusion.rst b/docs/modules/algorithms/advection_diffusion.rst
new file mode 100644
index 000000000..0f0ef57e3
--- /dev/null
+++ b/docs/modules/algorithms/advection_diffusion.rst
@@ -0,0 +1,9 @@
+.. _advection_diffusion_api:
+
+--------------------------------------------------------------------------------
+AdvectionDiffusion
+--------------------------------------------------------------------------------
+
+.. autoclass:: openpnm.algorithms.AdvectionDiffusion
+   :members:
+   :show-inheritance:
diff --git a/docs/modules/algorithms/fickian_diffusion.rst b/docs/modules/algorithms/fickian_diffusion.rst
new file mode 100644
index 000000000..1d2caaf6d
--- /dev/null
+++ b/docs/modules/algorithms/fickian_diffusion.rst
@@ -0,0 +1,9 @@
+.. _fickian_diffusion_api:
+
+--------------------------------------------------------------------------------
+FickianDiffusion
+--------------------------------------------------------------------------------
+
+.. autoclass:: openpnm.algorithms.FickianDiffusion
+   :members:
+   :show-inheritance:
diff --git a/docs/modules/algorithms/fourier_conduction.rst b/docs/modules/algorithms/fourier_conduction.rst
new file mode 100644
index 000000000..9a6414de6
--- /dev/null
+++ b/docs/modules/algorithms/fourier_conduction.rst
@@ -0,0 +1,9 @@
+.. _fourier_conduction_api:
+
+--------------------------------------------------------------------------------
+FourierConduction
+--------------------------------------------------------------------------------
+
+.. autoclass:: openpnm.algorithms.FourierConduction
+   :members:
+   :show-inheritance:
diff --git a/docs/modules/algorithms/index.rst b/docs/modules/algorithms/index.rst
index 8e13ae785..e4f5fa3f2 100644
--- a/docs/modules/algorithms/index.rst
+++ b/docs/modules/algorithms/index.rst
@@ -14,7 +14,7 @@ Overview of Module
 .. automodule:: openpnm.algorithms
 
 --------------------------------------------------------------------------------
-Links to Detailed Class Documentation
+Generic Classes
 --------------------------------------------------------------------------------
 
 .. toctree::
@@ -24,7 +24,28 @@ Links to Detailed Class Documentation
    generic_transport.rst
    reactive_transport.rst
    transient_reactive_transport.rst
-   ordinary_percolation.rst
-   porosimetry.rst
-   invasion_percolation.rst
-   mixed_percolation.rst
+
+--------------------------------------------------------------------------------
+Percolation Classes
+--------------------------------------------------------------------------------
+
+.. toctree::
+  :maxdepth: 3
+
+  ordinary_percolation.rst
+  porosimetry.rst
+  invasion_percolation.rst
+  mixed_percolation.rst
+
+--------------------------------------------------------------------------------
+Transport Classes
+--------------------------------------------------------------------------------
+
+.. toctree::
+  :maxdepth: 3
+
+  fickian_diffusion.rst
+  fourier_conduction.rst
+  ohmic_conduction.rst
+  stokes_flow.rst
+  advection_diffusion.rst
diff --git a/docs/modules/algorithms/ohmic_conduction.rst b/docs/modules/algorithms/ohmic_conduction.rst
new file mode 100644
index 000000000..5d86efc79
--- /dev/null
+++ b/docs/modules/algorithms/ohmic_conduction.rst
@@ -0,0 +1,9 @@
+.. _ohmic_conduction_api:
+
+--------------------------------------------------------------------------------
+OhmicConduction
+--------------------------------------------------------------------------------
+
+.. autoclass:: openpnm.algorithms.OhmicConduction
+   :members:
+   :show-inheritance:
diff --git a/docs/modules/materials/cubic_berea.rst b/docs/modules/algorithms/stokes_flow.rst
similarity index 64%
rename from docs/modules/materials/cubic_berea.rst
rename to docs/modules/algorithms/stokes_flow.rst
index 656129c6f..b7e9cb76e 100644
--- a/docs/modules/materials/cubic_berea.rst
+++ b/docs/modules/algorithms/stokes_flow.rst
@@ -1,8 +1,9 @@
-.. _berea_cubic_api:
+.. _stokes_flow_api:
 
 --------------------------------------------------------------------------------
-BereaCubic
+StokesFlow
 --------------------------------------------------------------------------------
 
-.. autoclass:: openpnm.materials.BereaCubic
+.. autoclass:: openpnm.algorithms.StokesFlow
    :members:
+   :show-inheritance:
diff --git a/docs/modules/materials/index.rst b/docs/modules/materials/index.rst
index 8f3644ac0..de1b676fc 100644
--- a/docs/modules/materials/index.rst
+++ b/docs/modules/materials/index.rst
@@ -21,4 +21,3 @@ Links to Detailed Class Documentation
    :maxdepth: 3
 
    voronoi_fibers.rst
-   cubic_berea.rst
diff --git a/docs/modules/models/geometry.rst b/docs/modules/models/geometry.rst
index d9e4d2ac7..e1fed99a6 100644
--- a/docs/modules/models/geometry.rst
+++ b/docs/modules/models/geometry.rst
@@ -21,11 +21,7 @@ Detailed Model Documentation
 Pore Size
 ................................................................................
 
-.. autofunction:: openpnm.models.geometry.pore_size.weibull
-.. autofunction:: openpnm.models.geometry.pore_size.normal
-.. autofunction:: openpnm.models.geometry.pore_size.generic
-.. autofunction:: openpnm.models.geometry.pore_size.largest_sphere
-.. autofunction:: openpnm.models.geometry.pore_size.equivalent_diameter
+.. automodule:: openpnm.models.geometry.pore_size
 
 ----
 
@@ -33,8 +29,7 @@ Pore Size
 Pore Volume
 ................................................................................
 
-.. autofunction:: openpnm.models.geometry.pore_volume.sphere
-.. autofunction:: openpnm.models.geometry.pore_volume.cube
+.. automodule:: openpnm.models.geometry.pore_volume
 
 ----
 
@@ -42,8 +37,7 @@ Pore Volume
 Pore Seed
 ................................................................................
 
-.. autofunction:: openpnm.models.geometry.pore_seed.random
-.. autofunction:: openpnm.models.geometry.pore_seed.spatially_correlated
+.. automodule:: openpnm.models.geometry.pore_seed
 
 ----
 
@@ -51,8 +45,7 @@ Pore Seed
 Pore Surface Area
 ................................................................................
 
-.. autofunction:: openpnm.models.geometry.pore_surface_area.sphere
-.. autofunction:: openpnm.models.geometry.pore_surface_area.cube
+.. automodule:: openpnm.models.geometry.pore_surface_area
 
 ----
 
@@ -60,8 +53,7 @@ Pore Surface Area
 Throat Area
 ................................................................................
 
-.. autofunction:: openpnm.models.geometry.throat_area.cylinder
-.. autofunction:: openpnm.models.geometry.throat_area.cuboid
+.. automodule:: openpnm.models.geometry.throat_area
 
 ----
 
@@ -69,11 +61,7 @@ Throat Area
 Throat Size
 ................................................................................
 
-.. autofunction:: openpnm.models.geometry.throat_size.weibull
-.. autofunction:: openpnm.models.geometry.throat_size.normal
-.. autofunction:: openpnm.models.geometry.throat_size.generic
-.. autofunction:: openpnm.models.geometry.throat_size.from_neighbor_pores
-.. autofunction:: openpnm.models.geometry.throat_size.equivalent_diameter
+.. automodule:: openpnm.models.geometry.throat_size
 
 ----
 
@@ -81,7 +69,7 @@ Throat Size
 Throat Length
 ................................................................................
 
-.. autofunction:: openpnm.models.geometry.throat_length.piecewise
+.. automodule:: openpnm.models.geometry.throat_length
 
 ----
 
@@ -89,8 +77,7 @@ Throat Length
 Throat Perimeter
 ................................................................................
 
-.. autofunction:: openpnm.models.geometry.throat_perimeter.cylinder
-.. autofunction:: openpnm.models.geometry.throat_perimeter.cuboid
+.. automodule:: openpnm.models.geometry.throat_perimeter
 
 ----
 
@@ -98,9 +85,15 @@ Throat Perimeter
 Throat Surface Area
 ................................................................................
 
-.. autofunction:: openpnm.models.geometry.throat_surface_area.cylinder
-.. autofunction:: openpnm.models.geometry.throat_surface_area.cuboid
-.. autofunction:: openpnm.models.geometry.throat_surface_area.extrusion
+.. automodule:: openpnm.models.geometry.throat_surface_area
+
+----
+
+................................................................................
+Throat Vector
+................................................................................
+
+.. automodule:: openpnm.models.geometry.throat_vector
 
 ----
 
@@ -108,9 +101,7 @@ Throat Surface Area
 Throat Volume
 ................................................................................
 
-.. autofunction:: openpnm.models.geometry.throat_volume.cylinder
-.. autofunction:: openpnm.models.geometry.throat_volume.cuboid
-.. autofunction:: openpnm.models.geometry.throat_volume.extrusion
+.. automodule:: openpnm.models.geometry.throat_volume
 
 ----
 
@@ -118,6 +109,4 @@ Throat Volume
 Throat Shape Factor
 ................................................................................
 
-.. autofunction:: openpnm.models.geometry.throat_shape_factor.compactness
-.. autofunction:: openpnm.models.geometry.throat_shape_factor.mason_morrow
-.. autofunction:: openpnm.models.geometry.throat_shape_factor.jenkins_rao
+.. automodule:: openpnm.models.geometry.throat_shape_factor
diff --git a/docs/modules/models/misc.rst b/docs/modules/models/misc.rst
index f24e54377..b4b43d171 100644
--- a/docs/modules/models/misc.rst
+++ b/docs/modules/models/misc.rst
@@ -17,25 +17,19 @@ Overview of Submodule
 Detailed Model Documentation
 --------------------------------------------------------------------------------
 
-.. currentmodule:: openpnm.models.misc
-
 ................................................................................
-Simple Calculations
+Basic Math and Calculations
 ................................................................................
 
-.. autofunction:: constant
-.. autofunction:: product
-.. autofunction:: scaled
-.. autofunction:: random
+.. automodule:: openpnm.models.misc.basic_math
 
 ----
 
 ................................................................................
-Generic Equations
+Simple Equations
 ................................................................................
 
-.. autofunction:: linear
-.. autofunction:: polynomial
+.. automodule:: openpnm.models.misc.simple_equations
 
 ----
 
@@ -43,9 +37,7 @@ Generic Equations
 Statistical Distributions
 ................................................................................
 
-.. autofunction:: weibull
-.. autofunction:: normal
-.. autofunction:: generic
+.. automodule:: openpnm.models.misc.statistical_distributions
 
 ----
 
@@ -53,5 +45,4 @@ Statistical Distributions
 Neighbor Lookups
 ................................................................................
 
-.. autofunction:: from_neighbor_throats
-.. autofunction:: from_neighbor_pores
+.. automodule:: openpnm.models.misc.neighbor_lookups
diff --git a/docs/modules/models/phase.rst b/docs/modules/models/phase.rst
index b472a749a..638d1a5ad 100644
--- a/docs/modules/models/phase.rst
+++ b/docs/modules/models/phase.rst
@@ -21,9 +21,7 @@ Detailed Model Documentation
 Density
 ................................................................................
 
-.. autofunction:: openpnm.models.phases.density.water
-.. autofunction:: openpnm.models.phases.density.standard
-.. autofunction:: openpnm.models.phases.density.ideal_gas
+.. automodule:: openpnm.models.phases.density
 
 ----
 
@@ -31,10 +29,7 @@ Density
 Diffusivity
 ................................................................................
 
-.. autofunction:: openpnm.models.phases.diffusivity.fuller
-.. autofunction:: openpnm.models.phases.diffusivity.fuller_scaling
-.. autofunction:: openpnm.models.phases.diffusivity.tyn_calus
-.. autofunction:: openpnm.models.phases.diffusivity.tyn_calus_scaling
+.. automodule:: openpnm.models.phases.diffusivity
 
 ----
 
@@ -42,7 +37,7 @@ Diffusivity
 Electrical Conductivity
 ................................................................................
 
-.. autofunction:: openpnm.models.phases.electrical_conductivity.percolating_continua
+.. automodule:: openpnm.models.phases.electrical_conductivity
 
 ----
 
@@ -50,9 +45,7 @@ Electrical Conductivity
 Molar Density
 ................................................................................
 
-.. autofunction:: openpnm.models.phases.molar_density.standard
-.. autofunction:: openpnm.models.phases.molar_density.ideal_gas
-.. autofunction:: openpnm.models.phases.molar_density.vanderwaals
+.. automodule:: openpnm.models.phases.molar_density
 
 ----
 
@@ -60,10 +53,7 @@ Molar Density
 Surface Tension
 ................................................................................
 
-.. autofunction:: openpnm.models.phases.surface_tension.water
-.. autofunction:: openpnm.models.phases.surface_tension.eotvos
-.. autofunction:: openpnm.models.phases.surface_tension.guggenheim_katayama
-.. autofunction:: openpnm.models.phases.surface_tension.brock_bird_scaling
+.. automodule:: openpnm.models.phases.surface_tension
 
 ----
 
@@ -71,9 +61,7 @@ Surface Tension
 Thermal Conductivity
 ................................................................................
 
-.. autofunction:: openpnm.models.phases.thermal_conductivity.water
-.. autofunction:: openpnm.models.phases.thermal_conductivity.chung
-.. autofunction:: openpnm.models.phases.thermal_conductivity.sato
+.. automodule:: openpnm.models.phases.thermal_conductivity
 
 ----
 
@@ -81,8 +69,7 @@ Thermal Conductivity
 Vapor Pressure
 ................................................................................
 
-.. autofunction:: openpnm.models.phases.vapor_pressure.water
-.. autofunction:: openpnm.models.phases.vapor_pressure.antoine
+.. automodule:: openpnm.models.phases.vapor_pressure
 
 ----
 
@@ -90,6 +77,4 @@ Vapor Pressure
 Viscosity
 ................................................................................
 
-.. autofunction:: openpnm.models.phases.viscosity.water
-.. autofunction:: openpnm.models.phases.viscosity.reynolds
-.. autofunction:: openpnm.models.phases.viscosity.chung
+.. automodule:: openpnm.models.phases.viscosity
diff --git a/docs/modules/models/physics.rst b/docs/modules/models/physics.rst
index 12e9d4b3d..7000e74c0 100644
--- a/docs/modules/models/physics.rst
+++ b/docs/modules/models/physics.rst
@@ -21,8 +21,7 @@ Detailed Model Documentation
 Capillary Pressure
 ................................................................................
 
-.. autofunction:: openpnm.models.physics.capillary_pressure.washburn
-.. autofunction:: openpnm.models.physics.capillary_pressure.purcell
+.. automodule:: openpnm.models.physics.capillary_pressure
 
 ----
 
@@ -30,7 +29,7 @@ Capillary Pressure
 Diffusive Conductance
 ................................................................................
 
-.. autofunction:: openpnm.models.physics.diffusive_conductance.ordinary_diffusion
+.. automodule:: openpnm.models.physics.diffusive_conductance
 
 ----
 
@@ -38,7 +37,7 @@ Diffusive Conductance
 Electrical Conductance
 ................................................................................
 
-.. autofunction:: openpnm.models.physics.electrical_conductance.series_resistors
+.. automodule:: openpnm.models.physics.electrical_conductance
 
 ----
 
@@ -46,7 +45,7 @@ Electrical Conductance
 Hydraulic Conductance
 ................................................................................
 
-.. autofunction:: openpnm.models.physics.hydraulic_conductance.hagen_poiseuille
+.. automodule:: openpnm.models.physics.hydraulic_conductance
 
 ----
 
@@ -54,7 +53,7 @@ Hydraulic Conductance
 Thermal Conductance
 ................................................................................
 
-.. autofunction:: openpnm.models.physics.thermal_conductance.series_resistors
+.. automodule:: openpnm.models.physics.thermal_conductance
 
 ----
 
@@ -62,5 +61,28 @@ Thermal Conductance
 Multiphase
 ................................................................................
 
-.. autofunction:: openpnm.models.physics.multiphase.conduit_conductance
-.. autofunction:: openpnm.models.physics.multiphase.late_filling
+.. automodule:: openpnm.models.physics.multiphase
+
+---
+
+................................................................................
+Poisson Shape Factors
+................................................................................
+
+.. automodule:: openpnm.models.physics.poisson_shape_factors
+
+---
+
+................................................................................
+Flow Shape Factors
+................................................................................
+
+.. automodule:: openpnm.models.physics.flow_shape_factors
+
+---
+
+................................................................................
+Source Terms
+................................................................................
+
+.. automodule:: openpnm.models.physics.generic_source_term
diff --git a/docs/userguide/customizing.rst b/docs/userguide/customizing.rst
index 17e08e750..d23446970 100644
--- a/docs/userguide/customizing.rst
+++ b/docs/userguide/customizing.rst
@@ -107,7 +107,7 @@ For example, let's create a custom Phase object for an oil with temperature depe
 
   class Oil(GenericPhase):
       def __init__(self, **kwargs):
-          super().__init__(**kwargs)
+          super().__init__(**kwargs)  # This is a python thing, and calls the init of the parent class
 
     		  self.add_model(propname='pore.viscosity',
              						 model=op.models.misc.polynomial,
@@ -129,11 +129,13 @@ Creating Customized Networks
 
 Unlike Geometry, Phase and Physics objects, a Network object requires more than a collection models.  In fact, Networks typically have no models.  Creating a custom Network type is all about defining the locations of the pores, and the connections between the throats.  Consider the following basic graph:
 
-4 ----- 3
-| \   / |
-5 - 0   |
-|     \ |
-2 - 6 - 1
+::
+
+    4 ----- 3
+    | \   / |
+    5 - 0   |
+    |     \ |
+    2 - 6 - 1
 
 Implementing this as a 'custom' Network can be done as by noting that the pore coordinates are:
 
@@ -178,9 +180,9 @@ To create an actual Network subclass, you can add the following to your ``'my_ne
 
 
   class MyNetwork(GenericNetwork):
-    def __init__(self, **kwargs):
-        # Place pore locations and throat connecting generation code here
-        super().__init__(conns=conns, coords=coords, **kwargs)
+      def __init__(self, **kwargs):
+          # Place pore locations and throat connecting generation code here
+          super().__init__(conns=conns, coords=coords, **kwargs)
 
 
 --------------------------------------------------------------------------------
diff --git a/openpnm/algorithms/InvasionPercolation.py b/openpnm/algorithms/InvasionPercolation.py
index 13d06df82..ccfda5b74 100644
--- a/openpnm/algorithms/InvasionPercolation.py
+++ b/openpnm/algorithms/InvasionPercolation.py
@@ -261,20 +261,23 @@ class InvasionPercolation(GenericAlgorithm):
         reverse assessing the occupancy of pore neighbors. Consider the
         following scenario when running standard IP without trapping,
         3 situations can happen after each invasion step:
-            The number of defending clusters stays the same and clusters can
+
+            * The number of defending clusters stays the same and clusters can
             shrink
-            A cluster of size one is suppressed
-            A cluster is split into multiple clusters
+            * A cluster of size one is suppressed
+            * A cluster is split into multiple clusters
+
         In reverse the following opposite situations can happen:
-            The number of defending clusters stays the same and clusters can
+
+            * The number of defending clusters stays the same and clusters can
             grow
-            A cluster of size one is created
-            Mutliple clusters merge into one cluster
-        With trapping the reversed rules are adjusted so that:
-            Only clusters that do not connect to a sink can grow and merge.
-            At the point that a neighbor connected to a sink is touched the
-            trapped cluster stops growing as this is the point of trapping in
-            forward invasion time.
+            * A cluster of size one is created
+            * Mutliple clusters merge into one cluster
+
+        With trapping the reversed rules are adjusted so that only clusters
+        that do not connect to a sink can grow and merge. At the point that a
+        neighbor connected to a sink is touched the trapped cluster stops
+        growing as this is the point of trapping in forward invasion time.
 
         Logger info displays the invasion sequence and pore index and a message
         with condition number based on the modified trapping rules and the
diff --git a/openpnm/algorithms/MixedInvasionPercolation.py b/openpnm/algorithms/MixedInvasionPercolation.py
index a5a94dac1..8c37acbd4 100644
--- a/openpnm/algorithms/MixedInvasionPercolation.py
+++ b/openpnm/algorithms/MixedInvasionPercolation.py
@@ -557,7 +557,6 @@ class MixedInvasionPercolation(GenericAlgorithm):
         partial : boolean
             Indicating whether partially filled network
 
-
         Notes
         -----
         It is applied as a post-process and runs the percolation algorithm in
diff --git a/openpnm/algorithms/__init__.py b/openpnm/algorithms/__init__.py
index 34478bc8e..e08de1a0d 100644
--- a/openpnm/algorithms/__init__.py
+++ b/openpnm/algorithms/__init__.py
@@ -1,8 +1,12 @@
 r"""
+
 **openpnm.algorithms**
+
 ----
+
 The ``algorithms`` module contains classes for conducting transport simulations
 on pore networks.
+
 """
 
 from .GenericAlgorithm import GenericAlgorithm
diff --git a/openpnm/io/JSONGraphFormat.py b/openpnm/io/JSONGraphFormat.py
index 70fee9522..dfcc3f30a 100644
--- a/openpnm/io/JSONGraphFormat.py
+++ b/openpnm/io/JSONGraphFormat.py
@@ -25,7 +25,8 @@ class JSONGraphFormat(GenericIO):
     The JGF standard must contain data formatted according to
     http://jsongraphformat.info and enforced by JSON schema validation.
 
-    Users must transfer any phase data to the network manually if they wish to keep it.
+    Users must transfer any phase data to the network manually if they wish to
+    keep it.
     """
 
     @classmethod
@@ -141,7 +142,7 @@ class JSONGraphFormat(GenericIO):
         with open(filename, 'r') as file:
             json_file = json.load(file)
             if not self.__validate_json__(json_file):
-                raise Exception(f'Error - {filename} is not in the JSON Graph Format.')
+                raise Exception('FIle is not in the JSON Graph Format')
 
         # Extract graph metadata from JSON
         number_of_nodes = json_file['graph']['metadata']['number_of_nodes']
diff --git a/openpnm/models/geometry/pore_seed.py b/openpnm/models/geometry/pore_seed.py
index 862838489..ba3d9a71c 100644
--- a/openpnm/models/geometry/pore_seed.py
+++ b/openpnm/models/geometry/pore_seed.py
@@ -1,3 +1,14 @@
+r"""
+
+Pore seed models are use to calculate random numbers for each pore, which can
+subsequently be used in statistical distributions to calculate actual pore
+sizes.
+
+.. autofunction:: openpnm.models.geometry.pore_seed.random
+.. autofunction:: openpnm.models.geometry.pore_seed.spatially_correlated
+
+"""
+
 import scipy as _sp
 from openpnm.models import misc as _misc
 from openpnm.utils import logging as _logging
@@ -35,8 +46,10 @@ def spatially_correlated(target, weights=None, strel=None):
         and extent of correlations.  The following would achieve a basic
         correlation in the z-direction:
 
-        strel = sp.array([[[0, 0, 0], [0, 0, 0], [0, 0, 0]], \
-                          [[0, 0, 0], [1, 1, 1], [0, 0, 0]], \
+    ::
+
+        strel = sp.array([[[0, 0, 0], [0, 0, 0], [0, 0, 0]],
+                          [[0, 0, 0], [1, 1, 1], [0, 0, 0]],
                           [[0, 0, 0], [0, 0, 0], [0, 0, 0]]])
 
     Notes
diff --git a/openpnm/models/geometry/pore_size.py b/openpnm/models/geometry/pore_size.py
index 8a8470219..c347ee165 100644
--- a/openpnm/models/geometry/pore_size.py
+++ b/openpnm/models/geometry/pore_size.py
@@ -1,3 +1,19 @@
+r"""
+The pore-size models in this sub-module are used to apply desired pore-size
+distributions to a pore network.  Most of the models accept pore seeds, and
+then look-up pore-sizes from cumuative distribtions functions.  There is also
+a model for finding the largest possible sphere that can be placed on each
+site.
+
+.. autofunction:: openpnm.models.geometry.pore_size.weibull
+.. autofunction:: openpnm.models.geometry.pore_size.normal
+.. autofunction:: openpnm.models.geometry.pore_size.random
+.. autofunction:: openpnm.models.geometry.pore_size.generic_distribution
+.. autofunction:: openpnm.models.geometry.pore_size.largest_sphere
+.. autofunction:: openpnm.models.geometry.pore_size.equivalent_diameter
+
+"""
+
 from openpnm.utils import logging as _logging
 from openpnm.models import misc as _misc
 import numpy as _np
@@ -25,7 +41,7 @@ def random(target, seed=None, num_range=[0, 1]):
                         num_range=num_range)
 
 
-normal.__doc__ = _misc.normal.__doc__
+random.__doc__ = _misc.random.__doc__
 
 
 def generic_distribution(target, func, seeds='pore.seed'):
diff --git a/openpnm/models/geometry/pore_surface_area.py b/openpnm/models/geometry/pore_surface_area.py
index d836d2de6..0c1c3c425 100644
--- a/openpnm/models/geometry/pore_surface_area.py
+++ b/openpnm/models/geometry/pore_surface_area.py
@@ -1,3 +1,9 @@
+r"""
+
+.. autofunction:: openpnm.models.geometry.pore_surface_area.sphere
+.. autofunction:: openpnm.models.geometry.pore_surface_area.cube
+
+"""
 import numpy as _np
 
 
diff --git a/openpnm/models/geometry/pore_volume.py b/openpnm/models/geometry/pore_volume.py
index 9e1453216..b2dd19bf9 100644
--- a/openpnm/models/geometry/pore_volume.py
+++ b/openpnm/models/geometry/pore_volume.py
@@ -1,3 +1,11 @@
+r"""
+These models calculate pore volumes depending on the specified shape
+
+.. autofunction:: openpnm.models.geometry.pore_volume.sphere
+.. autofunction:: openpnm.models.geometry.pore_volume.cube
+
+"""
+
 from numpy import pi as _pi
 
 
diff --git a/openpnm/models/geometry/throat_area.py b/openpnm/models/geometry/throat_area.py
index 404218cde..9389227eb 100644
--- a/openpnm/models/geometry/throat_area.py
+++ b/openpnm/models/geometry/throat_area.py
@@ -1,3 +1,9 @@
+r"""
+
+.. autofunction:: openpnm.models.geometry.throat_area.cylinder
+.. autofunction:: openpnm.models.geometry.throat_area.cuboid
+
+"""
 from numpy import pi as _pi
 
 
diff --git a/openpnm/models/geometry/throat_length.py b/openpnm/models/geometry/throat_length.py
index 2c602a170..6710329fa 100644
--- a/openpnm/models/geometry/throat_length.py
+++ b/openpnm/models/geometry/throat_length.py
@@ -1,3 +1,10 @@
+r"""
+
+.. autofunction:: openpnm.models.geometry.throat_length.piecewise
+.. autofunction:: openpnm.models.geometry.throat_length.ctc
+.. autofunction:: openpnm.models.geometry.throat_length.conduit_lengths
+
+"""
 import numpy as _np
 from scipy import sqrt as _sqrt
 from openpnm.utils import logging as _logging
diff --git a/openpnm/models/geometry/throat_perimeter.py b/openpnm/models/geometry/throat_perimeter.py
index 30319af0f..5d16b0a49 100644
--- a/openpnm/models/geometry/throat_perimeter.py
+++ b/openpnm/models/geometry/throat_perimeter.py
@@ -1,3 +1,9 @@
+r"""
+
+.. autofunction:: openpnm.models.geometry.throat_perimeter.cylinder
+.. autofunction:: openpnm.models.geometry.throat_perimeter.cuboid
+
+"""
 import scipy as _sp
 
 
diff --git a/openpnm/models/geometry/throat_shape_factor.py b/openpnm/models/geometry/throat_shape_factor.py
index 67d468c29..3556bfa1c 100644
--- a/openpnm/models/geometry/throat_shape_factor.py
+++ b/openpnm/models/geometry/throat_shape_factor.py
@@ -1,3 +1,11 @@
+r"""
+
+.. autofunction:: openpnm.models.geometry.throat_shape_factor.compactness
+.. autofunction:: openpnm.models.geometry.throat_shape_factor.mason_morrow
+.. autofunction:: openpnm.models.geometry.throat_shape_factor.jenkins_rao
+
+"""
+
 import scipy as _sp
 
 
diff --git a/openpnm/models/geometry/throat_size.py b/openpnm/models/geometry/throat_size.py
index 1e0f25fb3..a1232d484 100644
--- a/openpnm/models/geometry/throat_size.py
+++ b/openpnm/models/geometry/throat_size.py
@@ -1,3 +1,13 @@
+r"""
+
+.. autofunction:: openpnm.models.geometry.throat_size.weibull
+.. autofunction:: openpnm.models.geometry.throat_size.normal
+.. autofunction:: openpnm.models.geometry.throat_size.random
+.. autofunction:: openpnm.models.geometry.throat_size.generic_distribution
+.. autofunction:: openpnm.models.geometry.throat_size.from_neighbor_pores
+.. autofunction:: openpnm.models.geometry.throat_size.equivalent_diameter
+
+"""
 from openpnm.models import misc as _misc
 import numpy as _np
 
@@ -24,6 +34,14 @@ def generic_distribution(target, func, seeds='throat.seed'):
 generic_distribution.__doc__ = _misc.generic_distribution.__doc__
 
 
+def random(target, seed=None, num_range=[0, 1]):
+    return _misc.random(target=target, element='throat', seed=seed,
+                        num_range=num_range)
+
+
+random.__doc__ = _misc.random.__doc__
+
+
 def from_neighbor_pores(target, pore_prop='pore.diameter', mode='min'):
     return _misc.from_neighbor_pores(target=target, pore_prop=pore_prop,
                                      mode=mode)
diff --git a/openpnm/models/geometry/throat_surface_area.py b/openpnm/models/geometry/throat_surface_area.py
index f0952c0e3..18161a428 100644
--- a/openpnm/models/geometry/throat_surface_area.py
+++ b/openpnm/models/geometry/throat_surface_area.py
@@ -1,3 +1,10 @@
+r"""
+
+.. autofunction:: openpnm.models.geometry.throat_surface_area.cylinder
+.. autofunction:: openpnm.models.geometry.throat_surface_area.cuboid
+.. autofunction:: openpnm.models.geometry.throat_surface_area.extrusion
+
+"""
 import scipy as _sp
 
 
diff --git a/openpnm/models/geometry/throat_vector.py b/openpnm/models/geometry/throat_vector.py
index ecbd49768..02f51fd76 100644
--- a/openpnm/models/geometry/throat_vector.py
+++ b/openpnm/models/geometry/throat_vector.py
@@ -1,7 +1,6 @@
 r"""
-===============================================================================
-Submodule -- throat_vector
-===============================================================================
+
+.. autofunction:: openpnm.models.geometry.throat_vector.pore_to_pore
 
 """
 from transforms3d import _gohlketransforms as tr
@@ -13,7 +12,7 @@ def pore_to_pore(target):
 
     Parameters
     ----------
-    geometry : OpenPNM Geometry object
+    target: OpenPNM object
         The object containing the geometrical properties of the throats
 
     Notes
diff --git a/openpnm/models/geometry/throat_volume.py b/openpnm/models/geometry/throat_volume.py
index 49771aa86..d87bcfd12 100644
--- a/openpnm/models/geometry/throat_volume.py
+++ b/openpnm/models/geometry/throat_volume.py
@@ -1,3 +1,10 @@
+r"""
+
+.. autofunction:: openpnm.models.geometry.throat_volume.cylinder
+.. autofunction:: openpnm.models.geometry.throat_volume.cuboid
+.. autofunction:: openpnm.models.geometry.throat_volume.extrusion
+
+"""
 import scipy as _sp
 
 
diff --git a/openpnm/models/misc/__init__.py b/openpnm/models/misc/__init__.py
index dce5b1b8b..1799224b1 100644
--- a/openpnm/models/misc/__init__.py
+++ b/openpnm/models/misc/__init__.py
@@ -10,15 +10,7 @@ generating an array of random numbers.
 
 """
 
-from .misc import generic_function
-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_distribution
-from .misc import from_neighbor_throats
-from .misc import from_neighbor_pores
+from .statistical_distributions import generic_distribution, normal, weibull, random
+from .simple_equations import linear, polynomial, generic_function
+from .basic_math import constant, product, scaled
+from .neighbor_lookups import from_neighbor_throats, from_neighbor_pores
diff --git a/openpnm/models/misc/basic_math.py b/openpnm/models/misc/basic_math.py
new file mode 100644
index 000000000..651a57220
--- /dev/null
+++ b/openpnm/models/misc/basic_math.py
@@ -0,0 +1,85 @@
+r"""
+
+.. autofunction:: openpnm.models.misc.basic_math.constant
+.. autofunction:: openpnm.models.misc.basic_math.product
+.. autofunction:: openpnm.models.misc.basic_math.scaled
+
+"""
+
+import numpy as np
+import scipy.stats as spts
+from openpnm.utils import logging
+logger = logging.getLogger(__name__)
+
+
+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, prop1, prop2, **kwargs):
+    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.
+
+    prop1 : string
+        The name of the first argument
+
+    prop2 : string
+        The name of the second argument
+
+    Notes
+    -----
+    Additional properties can be specified beyond just ``prop1`` and ``prop2``
+    by including additional arguments in the function call (i.e. ``prop3 =
+    'pore.foo'``).
+    """
+    value = target[prop1]*target[prop2]
+    for item in kwargs.values():
+        value *= target[item]
+    return value
+
+
+def scaled(target, prop, factor):
+    r"""
+    Scales an existing value by a factor.  Useful for constricting some throat
+    property.
+
+    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.
+
+    prop : string
+        The dictionary key of the array containing the values to be scaled.
+
+    factor : scalar
+        The factor by which the values should be scaled.
+    """
+    value = target[prop]*factor
+    return value
diff --git a/openpnm/models/misc/misc.py b/openpnm/models/misc/neighbor_lookups.py
similarity index 100%
rename from openpnm/models/misc/misc.py
rename to openpnm/models/misc/neighbor_lookups.py
diff --git a/openpnm/models/misc/simple_equations.py b/openpnm/models/misc/simple_equations.py
new file mode 100644
index 000000000..b2b54a3a7
--- /dev/null
+++ b/openpnm/models/misc/simple_equations.py
@@ -0,0 +1,108 @@
+r""""
+
+.. autofunction:: openpnm.models.misc.simple_equations.generic_function
+.. autofunction:: openpnm.models.misc.simple_equations.polynomial
+.. autofunction:: openpnm.models.misc.simple_equations.linear
+
+"""
+
+import numpy as np
+from openpnm.utils import logging
+logger = logging.getLogger(__name__)
+
+
+def generic_function(target, prop, func, **kwargs):
+    r"""
+    Runs an arbitrary function on the given data
+
+    This allows users to place a customized calculation into the automatated
+    model regeneration pipeline.
+
+    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.
+
+    prop : string
+        The dictionary key containing the array to be operated on
+
+    func : Numpy function
+        A handle to the function to apply
+
+    kwargs : keyward arguments
+        All arguments required by the specific Numpy function
+
+    Examples
+    --------
+    The following example shows how to use a Numpy function, but any function
+    can be used, as long as it returns an array object:
+
+    >>> import openpnm as op
+    >>> import numpy as np
+    >>> pn = op.network.Cubic(shape=[5, 5, 5])
+    >>> geo = op.geometry.GenericGeometry(network=pn, pores=pn.Ps, throats=pn.Ts)
+    >>> geo['pore.rand'] = np.random.rand(geo.Np)
+    >>> geo.add_model(propname='pore.cos',
+    ...               model=op.models.misc.generic_function,
+    ...               func=np.cos,
+    ...               prop='pore.rand')
+    """
+    values = target[prop]
+    result = func(values, **kwargs)
+    if not isinstance(result, np.ndarray):
+        logger.warning('Given function must return a Numpy array')
+    return result
+
+
+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.
+
+    prop : string
+        The dictionary key containing the independent variable or phase
+        property to be used in the polynomial.
+
+    """
+    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/misc/statistical_distributions.py b/openpnm/models/misc/statistical_distributions.py
new file mode 100644
index 000000000..4379b895f
--- /dev/null
+++ b/openpnm/models/misc/statistical_distributions.py
@@ -0,0 +1,180 @@
+r"""
+
+.. autofunction:: openpnm.models.misc.statistical_distributions.random
+.. autofunction:: openpnm.models.misc.statistical_distributions.weibull
+.. autofunction:: openpnm.models.misc.statistical_distributions.normal
+.. autofunction:: openpnm.models.misc.statistical_distributions.generic_distribution
+
+"""
+
+import numpy as np
+import scipy.stats as spts
+from openpnm.utils import logging
+logger = logging.getLogger(__name__)
+
+
+def random(target, element, seed=None, num_range=[0, 1]):
+    r"""
+    Create an array of random numbers of a specified size.
+
+    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.
+
+    """
+    range_size = num_range[1] - num_range[0]
+    range_min = num_range[0]
+    if seed is not None:
+        np.random.seed(seed)
+    value = np.random.rand(target._count(element),)
+    value = value*range_size + range_min
+    return value
+
+
+def weibull(target, seeds, shape, scale, loc):
+    r"""
+    Produces values from a Weibull distribution given a set of random numbers.
+
+    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.
+
+    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
+        'shape' > 1 giving a more symmetrical distribution.
+
+    scale : float
+        This controls the width of the distribution with most of values falling
+        below this number.
+
+    loc : float
+        Applies an offset to the distribution such that the smallest values are
+        above this number.
+
+    Examples
+    --------
+    The following code illustrates the inner workings of this function,
+    which uses the 'weibull_min' method of the scipy.stats module.  This can
+    be used to find suitable values of 'shape', 'scale'` and 'loc'.  Note that
+    'shape' is represented by 'c' in the actual function call.
+
+    >>> import scipy
+    >>> func = scipy.stats.weibull_min(c=1.5, scale=0.0001, loc=0)
+    >>> import matplotlib.pyplot as plt
+    >>> fig = plt.hist(func.ppf(q=scipy.rand(10000)), bins=50)
+
+    """
+    seeds = target[seeds]
+    value = spts.weibull_min.ppf(q=seeds, c=shape, scale=scale, loc=loc)
+    return value
+
+
+def normal(target, seeds, scale, loc):
+    r"""
+    Produces values from a Weibull distribution given a set of random numbers.
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The object with which this function as associated.  This argument
+        is required to (1) set number of values to generate (geom.Np or
+        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
+
+    Examples
+    --------
+    The following code illustrates the inner workings of this function,
+    which uses the 'norm' method of the scipy.stats module.  This can
+    be used to find suitable values of 'scale' and 'loc'.
+
+    >>> import scipy
+    >>> func = scipy.stats.norm(scale=.0001, loc=0.001)
+    >>> import matplotlib.pyplot as plt
+    >>> fig = plt.hist(func.ppf(q=scipy.rand(10000)), bins=50)
+
+    """
+    seeds = target[seeds]
+    value = spts.norm.ppf(q=seeds, scale=scale, loc=loc)
+    return value
+
+
+def generic_distribution(target, seeds, func):
+    r"""
+    Accepts an 'rv_frozen' object from the Scipy.stats submodule and returns
+    values from the distribution for the given seeds
+
+    This uses the ``ppf`` method of the stats 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.
+
+    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.
+
+    Examples
+    --------
+    The following code illustrates the process of obtaining a 'frozen' Scipy
+    stats object and adding it as a model:
+
+    >>> import scipy
+    >>> import openpnm as op
+    >>> pn = op.network.Cubic(shape=[3, 3, 3])
+    >>> geo = op.geometry.GenericGeometry(network=pn, pores=pn.Ps, throats=pn.Ts)
+    >>> geo.add_model(propname='pore.seed',
+    ...               model=op.models.geometry.pore_seed.random)
+
+    Now retrieve the stats distribution and add to ``geo`` as a model:
+
+    >>> stats_obj = scipy.stats.weibull_min(c=2, scale=.0001, loc=0)
+    >>> geo.add_model(propname='pore.size',
+    ...               model=op.models.geometry.pore_size.generic_distribution,
+    ...               seeds='pore.seed',
+    ...               func=stats_obj)
+
+
+    >>> import matplotlib.pyplot as plt
+    >>> fig = plt.hist(stats_obj.ppf(q=scipy.rand(1000)), bins=50)
+
+    """
+    seeds = target[seeds]
+    value = func.ppf(seeds)
+    return value
diff --git a/openpnm/models/phases/density.py b/openpnm/models/phases/density.py
index ddd6017b9..fe1585c27 100644
--- a/openpnm/models/phases/density.py
+++ b/openpnm/models/phases/density.py
@@ -1,3 +1,12 @@
+r"""
+
+.. autofunction:: openpnm.models.phases.density.water
+.. autofunction:: openpnm.models.phases.density.standard
+.. autofunction:: openpnm.models.phases.density.ideal_gas
+
+"""
+
+
 def standard(target, mol_weight='pore.molecular_weight',
              molar_density='pore.molar_density'):
     r"""
diff --git a/openpnm/models/phases/diffusivity.py b/openpnm/models/phases/diffusivity.py
index 25fe6e09e..6bfa884c0 100644
--- a/openpnm/models/phases/diffusivity.py
+++ b/openpnm/models/phases/diffusivity.py
@@ -1,3 +1,12 @@
+r"""
+
+.. autofunction:: openpnm.models.phases.diffusivity.fuller
+.. autofunction:: openpnm.models.phases.diffusivity.fuller_scaling
+.. autofunction:: openpnm.models.phases.diffusivity.tyn_calus
+.. autofunction:: openpnm.models.phases.diffusivity.tyn_calus_scaling
+
+"""
+
 import scipy as _sp
 
 
diff --git a/openpnm/models/phases/electrical_conductivity.py b/openpnm/models/phases/electrical_conductivity.py
index b24adbe60..48771262c 100644
--- a/openpnm/models/phases/electrical_conductivity.py
+++ b/openpnm/models/phases/electrical_conductivity.py
@@ -1,3 +1,8 @@
+r"""
+
+.. autofunction:: openpnm.models.phases.electrical_conductivity.percolating_continua
+
+"""
 import scipy as _sp
 
 
diff --git a/openpnm/models/phases/molar_density.py b/openpnm/models/phases/molar_density.py
index 11ac388d5..99d79bc3c 100644
--- a/openpnm/models/phases/molar_density.py
+++ b/openpnm/models/phases/molar_density.py
@@ -1,3 +1,11 @@
+r"""
+
+.. autofunction:: openpnm.models.phases.molar_density.standard
+.. autofunction:: openpnm.models.phases.molar_density.ideal_gas
+.. autofunction:: openpnm.models.phases.molar_density.vanderwaals
+
+"""
+
 import scipy as sp
 
 
diff --git a/openpnm/models/phases/surface_tension.py b/openpnm/models/phases/surface_tension.py
index 2aeb8dae6..2ea82e60f 100644
--- a/openpnm/models/phases/surface_tension.py
+++ b/openpnm/models/phases/surface_tension.py
@@ -1,3 +1,12 @@
+r"""
+
+.. autofunction:: openpnm.models.phases.surface_tension.water
+.. autofunction:: openpnm.models.phases.surface_tension.eotvos
+.. autofunction:: openpnm.models.phases.surface_tension.guggenheim_katayama
+.. autofunction:: openpnm.models.phases.surface_tension.brock_bird_scaling
+
+"""
+
 import scipy as sp
 
 
diff --git a/openpnm/models/phases/thermal_conductivity.py b/openpnm/models/phases/thermal_conductivity.py
index 617ca9ac1..13d86818b 100644
--- a/openpnm/models/phases/thermal_conductivity.py
+++ b/openpnm/models/phases/thermal_conductivity.py
@@ -1,3 +1,11 @@
+r"""
+
+.. autofunction:: openpnm.models.phases.thermal_conductivity.water
+.. autofunction:: openpnm.models.phases.thermal_conductivity.chung
+.. autofunction:: openpnm.models.phases.thermal_conductivity.sato
+
+"""
+
 import scipy as sp
 
 
diff --git a/openpnm/models/phases/vapor_pressure.py b/openpnm/models/phases/vapor_pressure.py
index 751017280..6d89336f7 100644
--- a/openpnm/models/phases/vapor_pressure.py
+++ b/openpnm/models/phases/vapor_pressure.py
@@ -1,3 +1,10 @@
+r"""
+
+.. autofunction:: openpnm.models.phases.vapor_pressure.water
+.. autofunction:: openpnm.models.phases.vapor_pressure.antoine
+
+"""
+
 import numpy as np
 
 
diff --git a/openpnm/models/phases/viscosity.py b/openpnm/models/phases/viscosity.py
index 117b362c9..1acf665f1 100644
--- a/openpnm/models/phases/viscosity.py
+++ b/openpnm/models/phases/viscosity.py
@@ -1,3 +1,11 @@
+r"""
+
+.. autofunction:: openpnm.models.phases.viscosity.water
+.. autofunction:: openpnm.models.phases.viscosity.reynolds
+.. autofunction:: openpnm.models.phases.viscosity.chung
+
+"""
+
 import scipy as sp
 
 
diff --git a/openpnm/models/physics/capillary_pressure.py b/openpnm/models/physics/capillary_pressure.py
index 630a2c25e..e9d0b386c 100644
--- a/openpnm/models/physics/capillary_pressure.py
+++ b/openpnm/models/physics/capillary_pressure.py
@@ -1,3 +1,13 @@
+r"""
+
+.. autofunction:: openpnm.models.physics.capillary_pressure.washburn
+.. autofunction:: openpnm.models.physics.capillary_pressure.purcell
+.. autofunction:: openpnm.models.physics.capillary_pressure.purcell_bidirectional
+.. autofunction:: openpnm.models.physics.capillary_pressure.ransohoff_snap_off
+.. autofunction:: openpnm.models.physics.capillary_pressure.sinusoidal_bidirectional
+
+"""
+
 import scipy as _sp
 import sympy as syp
 import numpy as np
diff --git a/openpnm/models/physics/diffusive_conductance.py b/openpnm/models/physics/diffusive_conductance.py
index c08fb99f5..a05bf58de 100644
--- a/openpnm/models/physics/diffusive_conductance.py
+++ b/openpnm/models/physics/diffusive_conductance.py
@@ -1,3 +1,9 @@
+r"""
+
+.. autofunction:: openpnm.models.physics.diffusive_conductance.ordinary_diffusion
+
+"""
+
 from .misc import generic_conductance
 
 
diff --git a/openpnm/models/physics/electrical_conductance.py b/openpnm/models/physics/electrical_conductance.py
index 5d33d9d5b..92b762450 100644
--- a/openpnm/models/physics/electrical_conductance.py
+++ b/openpnm/models/physics/electrical_conductance.py
@@ -1,3 +1,9 @@
+r"""
+
+.. autofunction:: openpnm.models.physics.electrical_conductance.series_resistors
+
+"""
+
 from .misc import generic_conductance
 
 
diff --git a/openpnm/models/physics/flow_shape_factors.py b/openpnm/models/physics/flow_shape_factors.py
index cba6ac36c..211406dcf 100644
--- a/openpnm/models/physics/flow_shape_factors.py
+++ b/openpnm/models/physics/flow_shape_factors.py
@@ -1,3 +1,10 @@
+r"""
+
+.. autofunction:: openpnm.models.physics.flow_shape_factors.ball_and_stick
+.. autofunction:: openpnm.models.physics.flow_shape_factors.conical_frustum_and_stick
+
+"""
+
 from scipy import pi as _pi
 from scipy import arctanh as _atanh
 import scipy as _sp
diff --git a/openpnm/models/physics/generic_source_term.py b/openpnm/models/physics/generic_source_term.py
index dc038be45..0aec5410f 100644
--- a/openpnm/models/physics/generic_source_term.py
+++ b/openpnm/models/physics/generic_source_term.py
@@ -1,3 +1,16 @@
+r"""
+
+.. autofunction:: openpnm.models.physics.generic_source_term.standard_kinetics
+.. autofunction:: openpnm.models.physics.generic_source_term.linear
+.. autofunction:: openpnm.models.physics.generic_source_term.power_law
+.. autofunction:: openpnm.models.physics.generic_source_term.exponential
+.. autofunction:: openpnm.models.physics.generic_source_term.natural_exponential
+.. autofunction:: openpnm.models.physics.generic_source_term.logarithm
+.. autofunction:: openpnm.models.physics.generic_source_term.natural_logarithm
+.. autofunction:: openpnm.models.physics.generic_source_term.general_symbolic
+
+"""
+
 import scipy as _sp
 import sympy as _syp
 
diff --git a/openpnm/models/physics/hydraulic_conductance.py b/openpnm/models/physics/hydraulic_conductance.py
index 820f887a4..2d2f29340 100644
--- a/openpnm/models/physics/hydraulic_conductance.py
+++ b/openpnm/models/physics/hydraulic_conductance.py
@@ -1,3 +1,9 @@
+r"""
+
+.. autofunction:: openpnm.models.physics.hydraulic_conductance.hagen_poiseuille
+
+"""
+
 from .misc import generic_conductance
 
 
diff --git a/openpnm/models/physics/meniscus.py b/openpnm/models/physics/meniscus.py
index c38a1b816..edcc29bcc 100644
--- a/openpnm/models/physics/meniscus.py
+++ b/openpnm/models/physics/meniscus.py
@@ -1,3 +1,11 @@
+r"""
+
+.. autofunction:: openpnm.models.physics.meniscus.sinusoidal
+.. autofunction:: openpnm.models.physics.meniscus.purcell
+.. autofunction:: openpnm.models.physics.meniscus.general_toroidal
+
+"""
+
 import numpy as np
 import logging
 import sympy as syp
diff --git a/openpnm/models/physics/multiphase.py b/openpnm/models/physics/multiphase.py
index fa3100510..e9af2ead5 100644
--- a/openpnm/models/physics/multiphase.py
+++ b/openpnm/models/physics/multiphase.py
@@ -1,3 +1,10 @@
+r"""
+
+.. autofunction:: openpnm.models.physics.multiphase.conduit_conductance
+.. autofunction:: openpnm.models.physics.multiphase.late_filling
+
+"""
+
 import scipy as sp
 
 
diff --git a/openpnm/models/physics/poisson_shape_factors.py b/openpnm/models/physics/poisson_shape_factors.py
index ea528f62e..921e35262 100644
--- a/openpnm/models/physics/poisson_shape_factors.py
+++ b/openpnm/models/physics/poisson_shape_factors.py
@@ -1,3 +1,11 @@
+r"""
+
+.. autofunction:: openpnm.models.physics.poisson_shape_factors.ball_and_stick
+.. autofunction:: openpnm.models.physics.poisson_shape_factors.conical_frustum_and_stick
+
+"""
+
+
 from scipy import pi as _pi
 from scipy import arctanh as _atanh
 import scipy as _sp
diff --git a/openpnm/models/physics/thermal_conductance.py b/openpnm/models/physics/thermal_conductance.py
index 68a4afd94..c4e2a37dd 100644
--- a/openpnm/models/physics/thermal_conductance.py
+++ b/openpnm/models/physics/thermal_conductance.py
@@ -1,3 +1,9 @@
+r"""
+
+.. autofunction:: openpnm.models.physics.thermal_conductance.series_resistors
+
+"""
+
 from .misc import generic_conductance
 
 
diff --git a/openpnm/network/Bravais.py b/openpnm/network/Bravais.py
index 5af26ee7f..e984db8ca 100644
--- a/openpnm/network/Bravais.py
+++ b/openpnm/network/Bravais.py
@@ -33,19 +33,19 @@ class Bravais(GenericNetwork):
     mode : string
         The type of lattice to create.  Options are:
 
-        - **'sc'** : Simple cubic (Same as ``Cubic``)
-        - **'bcc'** : Body-centered cubic lattice
-        - **'fcc'** : Face-centered cubic lattice
-        - **'hcp'** : Hexagonal close packed (Note Implemented Yet)
+        - 'sc' : Simple cubic (Same as ``Cubic``)
+        - 'bcc' : Body-centered cubic lattice
+        - 'fcc' : Face-centered cubic lattice
+        - 'hcp' : Hexagonal close packed (Note Implemented Yet)
 
     name : string
         An optional name for the object to help identify it.  If not given,
         one will be generated.
 
     project : OpenPNM Project object, optional
-        Each OpenPNM object must be part of a *Project*.  If none is supplied
+        Each OpenPNM object must be part of a Project.  If none is supplied
         then one will be created and this Network will be automatically
-        assigned to it.  To create a *Project* use ``openpnm.Project()``.
+        assigned to it.  To create a Project use ``openpnm.Project()``.
 
     See Also
     --------
diff --git a/openpnm/network/Cubic.py b/openpnm/network/Cubic.py
index 2265abec7..272c4dfba 100644
--- a/openpnm/network/Cubic.py
+++ b/openpnm/network/Cubic.py
@@ -149,13 +149,27 @@ class Cubic(GenericNetwork):
 
         self['pore.coords'] = points
         self['throat.conns'] = pairs
-        self['pore.internal'] = True
-        self['throat.internal'] = True
-        topotools.label_faces(self)
-        Ps = self.pores('surface')
-        Ts = self.find_neighbor_throats(pores=Ps, mode='xnor')
+        # Label faces
+        x, y, z = self['pore.coords'].T
+        labels = ['front', 'back', 'left', 'right', 'bottom', 'top']
+        for label in labels:
+            if 'pore.'+label not in self.keys():
+                self['pore.'+label] = False
+        self['pore.front'][x <= x.min()] = True
+        self['pore.back'][x >= x.max()] = True
+        self['pore.left'][y <= y.min()] = True
+        self['pore.right'][y >= y.max()] = True
+        self['pore.bottom'][z <= z.min()] = True
+        self['pore.top'][z >= z.max()] = True
+        # Label surface pores
+        self['pore.surface'] = False
+        Ps = self.pores(labels)
+        self['pore.surface'][Ps] = True
         self['throat.surface'] = False
+        Ts = self.find_neighbor_throats(pores=Ps, mode='xnor')
         self['throat.surface'][Ts] = True
+        self['pore.internal'] = True
+        self['throat.internal'] = True
         # Scale network to requested spacing
         self['pore.coords'] *= spacing
 
diff --git a/openpnm/network/CubicTemplate.py b/openpnm/network/CubicTemplate.py
index 1f5d53d72..11a139e1c 100644
--- a/openpnm/network/CubicTemplate.py
+++ b/openpnm/network/CubicTemplate.py
@@ -76,8 +76,3 @@ class CubicTemplate(Cubic):
         self['pore.drop'] = template.flatten() == 0
         topotools.trim(network=self, pores=self.pores('drop'))
         del self['pore.drop']
-        # remove labels pertaining to surface pores, then redo post-trim
-        self.clear(mode='labels')
-        self['pore.internal'] = True
-        self['throat.internal'] = True
-        topotools.find_surface_pores(self)
diff --git a/openpnm/topotools/__init__.py b/openpnm/topotools/__init__.py
index d81cfe440..a34a2fad2 100644
--- a/openpnm/topotools/__init__.py
+++ b/openpnm/topotools/__init__.py
@@ -14,7 +14,6 @@ from .topotools import add_boundary_pores
 from .topotools import bond_percolation
 from .topotools import clone_pores
 from .topotools import connect_pores
-from .topotools import dimensionality
 from .topotools import extend
 from .topotools import find_path
 from .topotools import find_surface_pores
diff --git a/openpnm/topotools/topotools.py b/openpnm/topotools/topotools.py
index 01237b405..8597460de 100644
--- a/openpnm/topotools/topotools.py
+++ b/openpnm/topotools/topotools.py
@@ -1030,12 +1030,10 @@ def label_faces(network, tol=0.0, label='surface'):
 
     label : string
         An identifying label to isolate the pores on the faces of the network.
-        The default is 'surface'.  Surface pores can be found using
-        ``find_surface_pores``.
+        default is 'surface'.
 
     """
-    label = label.split('.', 1)[-1]
-    if 'pore.'+label not in network.labels():
+    if label not in network.labels():
         find_surface_pores(network, label=label)
     Psurf = network['pore.'+label]
     crds = network['pore.coords']
@@ -1045,16 +1043,12 @@ def label_faces(network, tol=0.0, label='surface'):
     yspan = ymax - ymin
     zmin, zmax = sp.amin(crds[:, 2]), sp.amax(crds[:, 2])
     zspan = zmax - zmin
-    dims = dimensionality(network)
-    if dims[0]:
-        network['pore.front'] = (crds[:, 0] <= (xmin + tol*xspan)) * Psurf
-        network['pore.back'] = (crds[:, 0] >= (xmax - tol*xspan)) * Psurf
-    if dims[1]:
-        network['pore.left'] = (crds[:, 1] <= (ymin + tol*yspan)) * Psurf
-        network['pore.right'] = (crds[:, 1] >= (ymax - tol*yspan)) * Psurf
-    if dims[2]:
-        network['pore.top'] = (crds[:, 2] >= (zmax - tol*zspan)) * Psurf
-        network['pore.bottom'] = (crds[:, 2] <= (zmin + tol*zspan)) * Psurf
+    network['pore.back'] = (crds[:, 0] >= (xmax - tol*xspan)) * Psurf
+    network['pore.right'] = (crds[:, 1] >= (ymax - tol*yspan)) * Psurf
+    network['pore.top'] = (crds[:, 2] >= (zmax - tol*zspan)) * Psurf
+    network['pore.front'] = (crds[:, 0] <= (xmin + tol*xspan)) * Psurf
+    network['pore.left'] = (crds[:, 1] <= (ymin + tol*yspan)) * Psurf
+    network['pore.bottom'] = (crds[:, 2] <= (zmin + tol*zspan)) * Psurf
 
 
 def find_surface_pores(network, markers=None, label='surface'):
@@ -1105,39 +1099,26 @@ def find_surface_pores(network, markers=None, label='surface'):
     >>> net.num_pores(['top','bottom', 'left', 'right', 'front','back'])
     98
 
+    This function is mostly useful for unique networks such as spheres, random
+    topology, or networks that have been subdivied.
+
     """
     import scipy.spatial as sptl
     if markers is None:
-        dims = dimensionality(network)
-        coords = network['pore.coords'][:, dims]
-        # normalize coords to a 1 unit cube centered on origin
-        coords -= sp.amin(coords, axis=0)
-        coords /= sp.amax(coords, axis=0)
-        coords -= 0.5
-        npts = max((network.Np/10, 100))
-        if sum(dims) == 1:
-            network['pore.'+label] = True
-            return
-        if sum(dims) == 2:
-            r = 0.75
-            theta = sp.linspace(0, 2*sp.pi, npts, dtype=float)
-            x = r*sp.cos(theta)
-            y = r*sp.sin(theta)
-            markers = sp.vstack((x, y)).T
-        if sum(dims) == 3:
-            r = 1.00
-            indices = sp.arange(0, npts, dtype=float) + 0.5
-            phi = sp.arccos(1 - 2*indices/npts)
-            theta = sp.pi * (1 + 5**0.5) * indices
-            x = r*sp.cos(theta) * sp.sin(phi)
-            y = r*sp.sin(theta) * sp.sin(phi)
-            z = r*sp.cos(phi)
-            markers = sp.vstack((x, y, z)).T
-    else:
-        coords = network['pore.coords']
-        markers = sp.atleast_2d(markers)
-    pts = sp.vstack((coords, markers))
-    tri = sptl.Delaunay(pts, incremental=False)
+        (xmax, ymax, zmax) = sp.amax(network['pore.coords'], axis=0)
+        (xmin, ymin, zmin) = sp.amin(network['pore.coords'], axis=0)
+        xave = (xmin+xmax)/2
+        yave = (ymin+ymax)/2
+        zave = (zmin+zmax)/2
+        markers = [[xmax + xave, yave, zave],
+                   [xmin - xave, yave, zave],
+                   [xave, ymax + yave, zave],
+                   [xave, ymin - yave, zave],
+                   [xave, yave, zmax + zave],
+                   [xave, yave, zmin - zave]]
+    markers = sp.atleast_2d(markers)
+    tri = sptl.Delaunay(network['pore.coords'], incremental=True)
+    tri.add_points(markers)
     (indices, indptr) = tri.vertex_neighbor_vertices
     for k in range(network.Np, tri.npoints):
         neighbors = indptr[indices[k]:indices[k+1]]
@@ -1148,30 +1129,8 @@ def find_surface_pores(network, markers=None, label='surface'):
         network['pore.'+label][neighbors] = True
 
 
-def dimensionality(network):
-    r"""
-    Checks the dimensionality of the network
-
-    Parameters
-    ----------
-    network : OpenPNM Network object
-        The network whose dimensionality is to be checked
-
-    Returns
-    -------
-    Returns an 3-by-1 array containing ``True`` for each axis that contains
-    multiple values, indicating that the pores are spatially distributed
-    in that direction.
-    """
-    v = [True, True, True]
-    for i in [0, 1, 2]:
-        v[i] = not sp.allclose(network['pore.coords'][:, i],
-                               network['pore.coords'][1, i])
-    return v
-
-
 def clone_pores(network, pores, labels=['clone'], mode='parents'):
-    r"""
+    r'''
     Clones the specified pores and adds them to the network
 
     Parameters
@@ -1192,7 +1151,7 @@ def clone_pores(network, pores, labels=['clone'], mode='parents'):
         - 'siblings': Clones are only connected to each other in the same
                       manner as parents were connected
         - 'isolated': No connections between parents or siblings
-    """
+    '''
     if len(network.project.geometries()) > 0:
         logger.warning('Network has active Geometries, new pores must be \
                         assigned a Geometry')
@@ -1453,6 +1412,9 @@ def connect_pores(network, pores1, pores2, labels=[], add_conns=True):
     in which case it consecutively connects corresponding members of the two
     lists in a 1-to-1 fashion. Example: pores1 = [[0, 1], [2, 3]] and
     pores2 = [[5], [7, 9]] leads to creation of the following connections:
+
+    ::
+
         0 --> 5     2 --> 7     3 --> 7
         1 --> 5     2 --> 9     3 --> 9
 
@@ -1460,7 +1422,7 @@ def connect_pores(network, pores1, pores2, labels=[], add_conns=True):
     within ``pores1`` and ``pores2`` are of type list or ndarray.
 
     (3) It creates the connections in a format which is acceptable by
-    the default OpenPNM connection ('throat.conns') and either adds them to
+    the default OpenPNM connections ('throat.conns') and either adds them to
     the network or returns them.
 
     Examples
@@ -1511,7 +1473,7 @@ def connect_pores(network, pores1, pores2, labels=[], add_conns=True):
 
 def find_pore_to_pore_distance(network, pores1=None, pores2=None):
     r'''
-    Find the distance between all pores on set 1 to each pore in set 2
+    Find the distance between all pores on set one to each pore in set 2
 
     Parameters
     ----------
@@ -1531,13 +1493,6 @@ def find_pore_to_pore_distance(network, pores1=None, pores2=None):
     The distance between pore *i* in ``pores1`` and *j* in ``pores2`` is
     located at *(i, j)* and *(j, i)* in the distance matrix.
 
-    Notes
-    -----
-    This function computes and returns a distance matrix, which is a dense
-    matrix of size Np_1 by Np_2, so can get large.  For distances between
-    larger sets a KD-tree approach would be better, which is available in
-    ``scipy.spatial``.
-
     '''
     from scipy.spatial.distance import cdist
     p1 = sp.array(pores1, ndmin=1)
@@ -1836,7 +1791,7 @@ def _template_sphere_disc(dim, outer_radius, inner_radius):
     return img
 
 
-def template_sphere_shell(outer_radius, inner_radius=0, dim=3):
+def template_sphere_shell(outer_radius, inner_radius=0):
     r"""
     This method generates an image array of a sphere-shell. It is useful for
     passing to Cubic networks as a ``template`` to make spherical shaped
@@ -1851,17 +1806,13 @@ def template_sphere_shell(outer_radius, inner_radius=0, dim=3):
         Number of nodes in the inner radius of the shell.  a value of 0 will
         result in a solid sphere.
 
-    dim : scalar
-        Controls the number of dimensions of the result.  3 returns a sphere,
-        while 2 returns a disk.
-
     Returns
     -------
     A Numpy array containing 1's to demarcate the sphere-shell, and 0's
     elsewhere.
 
     """
-    img = _template_sphere_disc(dim=dim, outer_radius=outer_radius,
+    img = _template_sphere_disc(dim=3, outer_radius=outer_radius,
                                 inner_radius=inner_radius)
     return img
 


Test patch:
diff --git a/tests/unit/io/JSONGraphFormatTest.py b/tests/unit/io/JSONGraphFormatTest.py
index b539f4ffa..45fa900f7 100644
--- a/tests/unit/io/JSONGraphFormatTest.py
+++ b/tests/unit/io/JSONGraphFormatTest.py
@@ -143,10 +143,8 @@ class JSONGraphFormatTest:
         filename = Path(path.resolve(), 'invalid.json')
 
         # Ensure an exception was thrown
-        with pytest.raises(Exception) as e_info:
+        with pytest.raises(Exception):
             op.io.JSONGraphFormat.load(filename)
-        expected_error = f'Error - {filename} is not in the JSON Graph Format.'
-        assert expected_error in str(e_info.value)
 
     def test_load_success(self):
         # Load JSON file and ensure project integrity
diff --git a/tests/unit/network/CubicDualTest.py b/tests/unit/network/CubicDualTest.py
index d20f766f6..6e78da3b9 100644
--- a/tests/unit/network/CubicDualTest.py
+++ b/tests/unit/network/CubicDualTest.py
@@ -46,6 +46,8 @@ class CubicDualTest:
         assert net.num_pores('right') == 9
         assert net.num_pores('front') == 9
         assert net.num_pores('back') == 9
+        assert net.num_pores('top') == 41
+        assert net.num_pores('bottom') == 41
         assert net.num_throats('interconnect') == 96
 
     def test_generation_2D_XZ(self):
@@ -53,6 +55,8 @@ class CubicDualTest:
                                    label_2='secondary')
         assert net.Np == 57
         assert net.Nt == 176
+        assert net.num_pores('left') == 41
+        assert net.num_pores('right') == 41
         assert net.num_pores('front') == 9
         assert net.num_pores('back') == 9
         assert net.num_pores('top') == 9
@@ -66,6 +70,8 @@ class CubicDualTest:
         assert net.Nt == 176
         assert net.num_pores('left') == 9
         assert net.num_pores('right') == 9
+        assert net.num_pores('front') == 41
+        assert net.num_pores('back') == 41
         assert net.num_pores('top') == 9
         assert net.num_pores('bottom') == 9
         assert net.num_throats('interconnect') == 96
@@ -79,6 +85,8 @@ class CubicDualTest:
         assert net.num_pores('right') == 9
         assert net.num_pores('front') == 9
         assert net.num_pores('back') == 9
+        assert net.num_pores('top') == 41
+        assert net.num_pores('bottom') == 41
         assert net.num_throats('interconnect') == 96
 
 
diff --git a/tests/unit/network/LabelTest.py b/tests/unit/network/LabelTest.py
index bf7d5ad87..309834ef4 100644
--- a/tests/unit/network/LabelTest.py
+++ b/tests/unit/network/LabelTest.py
@@ -30,19 +30,6 @@ class LabelTest:
         assert net.num_pores('left') == 16
         assert net.num_pores(['left', 'top'], mode='xnor') == 4
 
-    def test_cubic_2D(self):
-        net = op.network.Cubic(shape=[5, 5, 1])
-        assert 'pore.top' not in net.labels()
-        assert 'pore.bottom' not in net.labels()
-        assert net.num_pores('surface') == 16
-        net = op.network.Cubic(shape=[5, 1, 5])
-        assert 'pore.left' not in net.labels()
-        assert 'pore.right' not in net.labels()
-        assert net.num_pores('surface') == 16
-        net = op.network.Cubic(shape=[1, 5, 5])
-        assert 'pore.front' not in net.labels()
-        assert 'pore.back' not in net.labels()
-        assert net.num_pores('surface') == 16
 
 if __name__ == '__main__':
 
diff --git a/tests/unit/topotools/TopotoolsTest.py b/tests/unit/topotools/TopotoolsTest.py
index 4c9c12ca0..f0479890f 100644
--- a/tests/unit/topotools/TopotoolsTest.py
+++ b/tests/unit/topotools/TopotoolsTest.py
@@ -56,7 +56,7 @@ class TopotoolsTest:
         net.clear(mode='labels')
         assert net.labels() == ['pore.all', 'throat.all']
         topotools.find_surface_pores(network=net)
-        assert net.num_pores('surface') == 66
+        assert net.num_pores('surface') == 42
 
     def test_find_pore_to_pore_distance(self):
         net = op.network.Cubic(shape=[3, 3, 3], connectivity=6)