Repository: PMEAL/OpenPNM
Pull request #838 (2018-04-02 05:00:08)

Reworking algs
After having ruminated on the changes made to algs, it's time to finalize them.


Patch:
diff --git a/example_script.py b/example_script.py
index 9e6ba55fe..cac66dbf9 100644
--- a/example_script.py
+++ b/example_script.py
@@ -39,23 +39,26 @@ water['pore.k'] = 2
 water.add_model(propname='pore.reaction',
                 model=op.models.physics.generic_source_term.standard_kinetics,
                 prefactor='pore.A', exponent='pore.k',
-                quantity='pore.pressure', regen_mode='deferred')
+                quantity='pore.pressure', regen_mode='normal')
 
 s = {'conductance': 'throat.conductance',
      'quantity': 'pore.pressure'}
-alg1 = op.algorithms.GenericTransport(network=pn, phase=water, settings=s)
+alg1 = op.algorithms.GenericTransport(network=pn)
+# You can specify phase and settings using the setup method
+alg1.setup(phase=water, **s)
 alg1.set_dirichlet_BC(pores=pn.pores('inlets'), values=1)
 alg1.set_dirichlet_BC(pores=pn.pores('outlets'), values=0)
 alg1.run()
 
+# You can also specify phase and settings during initialization
 alg2 = op.algorithms.ReactiveTransport(network=pn, phase=water, settings=s)
 alg2.set_dirichlet_BC(pores=pn.pores('inlets'), values=1)
-alg2.set_source_term(propname='pore.reaction', pores=pn.pores('outlets'))
+alg2.set_source(propname='pore.reaction', pores=pn.pores('outlets'))
 alg2.run()
 water.update(alg2.results())
 
 alg3 = op.algorithms.TransientTransport(network=pn, phase=water)
-# You can also set the settings afterwards.  Note that some of these
+# You can also set the settings directly.  Note that some of these
 # will have defaults when finally subclassed (i.e. quantity = pressure)
 alg3.settings.update({'t_initial': 0,
                       't_final': 1,
@@ -66,10 +69,9 @@ alg3.set_dirichlet_BC(pores=pn.pores('inlets'), values=1)
 alg3.set_IC(values=0)
 alg3.run()
 
-
 alg4 = op.algorithms.TransientReactiveTransport(network=pn, phase=water)
 alg4.settings.update(alg3.settings)  # Just copy settings from another alg
 alg4.set_dirichlet_BC(pores=pn.pores('inlets'), values=1)
 alg4.set_IC(values=0)
-alg4.set_source_term(propname='pore.reaction', pores=pn.pores('bottom'))
+alg4.set_source(propname='pore.reaction', pores=pn.pores('bottom'))
 alg4.run()
diff --git a/openpnm/algorithms/CrankNicholson.py b/openpnm/algorithms/CrankNicholson.py
deleted file mode 100644
index b630c0eaf..000000000
--- a/openpnm/algorithms/CrankNicholson.py
+++ /dev/null
@@ -1,19 +0,0 @@
-from openpnm.core import logging
-from openpnm.algorithms import TransientTransport
-logger = logging.getLogger(__name__)
-
-
-class CrankNicholson(TransientTransport):
-    r"""
-    """
-
-    def build_B(self):
-        r"""
-        """
-        print('CrankNicholsonMixin: build_B')
-
-    def update_A(self):
-        r"""
-        """
-        print('CrankNicholsonMixin: update_A')
-        self.build_B()
diff --git a/openpnm/algorithms/FickianDiffusion.py b/openpnm/algorithms/FickianDiffusion.py
index e2782da35..1dc729229 100644
--- a/openpnm/algorithms/FickianDiffusion.py
+++ b/openpnm/algorithms/FickianDiffusion.py
@@ -1,10 +1,10 @@
 import scipy as sp
-from openpnm.algorithms import GenericTransport
+from openpnm.algorithms import ReactiveTransport
 from openpnm.core import logging
 logger = logging.getLogger(__name__)
 
 
-class FickianDiffusion(GenericTransport):
+class FickianDiffusion(ReactiveTransport):
     r"""
     A subclass of GenericLinearTransport to simulate binary diffusion. The 2
     main roles of this subclass are to set the default property names and to
@@ -19,12 +19,49 @@ class FickianDiffusion(GenericTransport):
                               'conductance': 'throat.diffusive_conductance'})
         self.settings.update(settings)
 
+    def setup(self, phase=None, quantity='', conductance='', **kwargs):
+        r"""
+        This method takes several arguments that are essential to running the
+        algorithm and adds them to the settings.
+
+        Parameters
+        ----------
+        phase : OpenPNM Phase object
+            The phase on which the algorithm is to be run.  If no value is
+            given, the existing value is kept.
+
+        quantity : string
+            The name of the physical quantity to be calcualted.  If no value is
+            given, the existing value is kept.  The default value is
+            ``'pore.mole_fraction'``.
+
+        conductance : string
+            The name of the pore-scale transport conductance values.  These
+            are typically calculate by a model attached to a *Physics* object
+            associated with the given *Phase*.  If no value is given, the
+            existing value is kept.  The default value is
+            ``'throat.diffusive_conductance'``.
+
+        Notes
+        -----
+        Any additional arguments are added to the ``settings`` dictionary of
+        the object.
+
+        """
+        if phase:
+            self.settings['phase'] = phase.name
+        if quantity:
+            self.settings['quantity': quantity]
+        if conductance:
+            self.settings['conductance': conductance]
+        super().setup(**kwargs)
+
     def calc_eff_diffusivity(self):
         r"""
         This calculates the effective diffusivity in this linear transport
         algorithm.
         """
-        phase = self.project.phases[self.settings['phase']]
+        phase = self.project.phases()[self.settings['phase']]
         d_normal = self._calc_eff_prop()
         self._eff_property = d_normal / sp.mean(phase['pore.molar_density'])
         return self._eff_property
diff --git a/openpnm/algorithms/FourierConduction.py b/openpnm/algorithms/FourierConduction.py
index 1608e6eb5..e00396a60 100644
--- a/openpnm/algorithms/FourierConduction.py
+++ b/openpnm/algorithms/FourierConduction.py
@@ -1,9 +1,9 @@
-from openpnm.algorithms import GenericTransport
+from openpnm.algorithms import ReactiveTransport
 from openpnm.core import logging
 logger = logging.getLogger(__name__)
 
 
-class FourierConduction(GenericTransport):
+class FourierConduction(ReactiveTransport):
     r"""
     A subclass of GenericLinearTransport to simulate heat conduction.  The 2
     main roles of this subclass are to set the default property names and to
@@ -18,6 +18,43 @@ class FourierConduction(GenericTransport):
                               'conductance': 'throat.thermal_conductance'})
         self.settings.update(settings)
 
+    def setup(self, phase=None, quantity='', conductance='', **kwargs):
+        r"""
+        This method takes several arguments that are essential to running the
+        algorithm and adds them to the settings.
+
+        Parameters
+        ----------
+        phase : OpenPNM Phase object
+            The phase on which the algorithm is to be run.  If no value is
+            given, the existing value is kept.
+
+        quantity : string
+            The name of the physical quantity to be calcualted.  If no value is
+            given, the existing value is kept.  The default value is
+            ``'pore.temperature'``.
+
+        conductance : string
+            The name of the pore-scale transport conductance values.  These
+            are typically calculate by a model attached to a *Physics* object
+            associated with the given *Phase*.  If no value is given, the
+            existing value is kept.  The default value is
+            ``'throat.thermal_conductance'``.
+
+        Notes
+        -----
+        Any additional arguments are added to the ``settings`` dictionary of
+        the object.
+
+        """
+        if phase:
+            self.settings['phase'] = phase.name
+        if quantity:
+            self.settings['quantity': quantity]
+        if conductance:
+            self.settings['conductance': conductance]
+        super().setup(**kwargs)
+
     def calc_effective_conductivity(self):
         r"""
         This calculates the effective thermal conductivity.
diff --git a/openpnm/algorithms/GenericLinearTransport.py b/openpnm/algorithms/GenericLinearTransport.py
deleted file mode 100644
index 732d5b1da..000000000
--- a/openpnm/algorithms/GenericLinearTransport.py
+++ /dev/null
@@ -1,326 +0,0 @@
-"""
-===============================================================================
-module GenericLinearTransport: Class for solving linear transport processes
-===============================================================================
-
-"""
-import scipy as sp
-import scipy.sparse as sprs
-import scipy.sparse.csgraph as spgr
-from openpnm.algorithms import GenericAlgorithm
-from openpnm.core import logging
-logger = logging.getLogger(__name__)
-
-
-class GenericLinearTransport(GenericAlgorithm):
-    r"""
-    """
-
-    def __init__(self, network, phase, settings={}, **kwargs):
-        self.settings.update({'phase': phase.name,
-                              'conductance': None,
-                              'quantity': None,
-                              'solver': 'spsolve',
-                              'sources': [],
-                              'tolerance': 0.001})
-        self.settings.update(settings)
-        super().__init__(Np=network.Np, Nt=network.Nt, network=network, **kwargs)
-        # Initialize label 'all' in the object's own info dictionaries
-        self['pore._id'] = network['pore._id']
-        self['throat._id'] = network['throat._id']
-
-    def set_dirchlet_BC(self, pores, values):
-        r"""
-        """
-        self.set_BC(pores=pores, bctype='dirichlet', bcvalues=values,
-                    mode='merge')
-
-    def set_neumann_BC(self, pores, values):
-        r"""
-        """
-        self.set_BC(pores=pores, bctype='neumann', bcvalues=values,
-                    mode='merge')
-
-    def set_source(self, source):
-        self.settings['sources'].append(source.name)
-        source.setup(algorithm=self)
-
-    def set_BC(self, pores, bctype, bcvalues=None, mode='merge'):
-        r"""
-        Apply boundary conditions to specified pores
-
-        Parameters
-        ----------
-        pores : array_like
-            The pores where the boundary conditions should be applied
-
-        bctype : string
-            Specifies the type or the name of boundary condition to apply. The
-            types can be one one of the following:
-
-            - *'dirichlet'* : Specify the quantity in each location
-            - *'neumann'* : Specify the flow rate into each location
-
-        bcvalues : int or array_like
-            The boundary value to apply, such as concentration or rate.  If
-            a single value is given, it's assumed to apply to all locations.
-            Different values can be applied to all pores in the form of an
-            array of the same length as ``pores``.
-
-        mode : string, optional
-            Controls how the conditions are applied.  Options are:
-
-            - *'merge'*: (Default) Adds supplied boundary conditions to already
-            existing conditions.
-
-        Notes
-        -----
-        It is not possible to have multiple boundary conditions for a
-        specified location in just one algorithm. Use ``mode='remove'`` to
-        clear existing BCs before applying new ones or ``mode='overwrite'``
-        which removes all existing BC's before applying the new ones.
-
-        Instead of using ``mode='remove'`` you can also set certain locations
-        to NaN using ``mode='merge'``, which is equivalent to removing the BCs
-        from those locations.
-
-        """
-        # Hijack the parse_mode function to verify bctype argument
-        bctype = self._parse_mode(bctype, allowed=['dirichlet', 'neumann',
-                                                   'neumann_group'],
-                                  single=True)
-        mode = self._parse_mode(mode, allowed=['merge', 'overwrite', 'remove'],
-                                single=True)
-        pores = self._parse_indices(pores)
-
-        values = sp.array(bcvalues)
-        if values.size > 1 and values.size != pores.size:
-            raise Exception('The number of boundary values must match the ' +
-                            'number of locations')
-
-        # Label pores where a boundary condition will be applied
-        if ('pore.'+bctype not in self.keys()) or (mode == 'overwrite'):
-            self['pore.'+bctype] = False
-        self['pore.'+bctype][pores] = True
-
-        # Store boundary values
-        if ('pore.'+bctype+'_value' not in self.keys()) or \
-           (mode == 'overwrite'):
-            self['pore.'+bctype+'_value'] = sp.nan
-        self['pore.'+bctype+'_value'][pores] = values
-
-    def remove_BC(self, pores=None):
-        r"""
-        Removes all boundary conditions assigned to the specified pores
-
-        Parameters
-        ----------
-        pores : array_like
-            The pores from which boundary conditions are to be removed.  If no
-            pores are specified, then BCs are removed from all pores. No error
-            is thrown if the provided pores do not have any BCs assigned.
-        """
-        if pores is None:
-            pores = self.Ps
-        if 'pore.dirichlet' in self.keys():
-            self['pore.dirichlet'][pores] = False
-            self['pore.dirichlet_value'][pores] = sp.nan
-        if 'pore.neumann' in self.keys():
-            self['pore.neumann'][pores] = False
-            self['pore.neumann_value'][pores] = sp.nan
-
-    def setup(self, conductance, quantity):
-        r"""
-        Specifies the necessary parameters
-
-        Parameters
-        ----------
-        conductance : string
-            The dictionary key containing the calculated pore-scale
-            conductances.  For example, for StokesFlow this is
-            'throat.hydraulic_conductance' by default.
-
-        quantity : string
-            The dictionary key where the values computed by this algorithm are
-            stored.  For exaple, for StokesFLow this is 'pore.pressure' by
-            default.
-
-        """
-        # Providing conductance values for the algorithm from the Physics name
-        self.settings['conductance'] = self._parse_prop(conductance, 'throat')
-        self.settings['quantity'] = self._parse_prop(quantity, 'pore')
-
-        # Check health of conductance vector
-        phase = self.project.phases()[self.settings['phase']]
-        if sp.any(sp.isnan(phase[self.settings['conductance']])):
-            raise Exception('The provided throat conductance contains NaNs')
-
-    def build_A(self):
-        r"""
-        """
-        network = self.project.network
-        phase = self.project.phases()[self.settings['phase']]
-        g = phase[self.settings['conductance']]
-        am = network.create_adjacency_matrix(weights=-g, fmt='coo')
-        A = spgr.laplacian(am)
-        if 'pore.neumann' in self.keys():
-            pass  # Do nothing to A, only b changes
-        if 'pore.dirichlet' in self.keys():
-            # Find all entries on rows associated with dirichlet pores
-            P_bc = self.toindices(self['pore.dirichlet'])
-            indrow = sp.in1d(A.row, P_bc)
-            A.data[indrow] = 0  # Remove entries from A for all BC rows
-            datadiag = A.diagonal()  # Add diagonal entries back into A
-            datadiag[P_bc] = sp.ones_like(P_bc, dtype=float)
-            A.setdiag(datadiag)
-            A.eliminate_zeros()  # Remove 0 entries
-        self.A = A
-        return A
-
-    def build_b(self):
-        r"""
-        """
-        b = sp.zeros(shape=(self.Np, ), dtype=float)  # Create b matrix of 0's
-        if 'pore.dirichlet' in self.keys():
-            ind = self['pore.dirichlet']
-            b[ind] = self['pore.dirichlet_value'][ind]
-        if 'pore.neumann' in self.keys():
-            ind = self['pore.neumann']
-            b[ind] = self['pore.neumann_value'][ind]
-        self.b = b
-        return b
-
-    def run(self, x=None):
-        r"""
-        Builds the A and b matrices, and calls the solver specified in the
-        ``settings`` attribute.
-
-        Parameters
-        ----------
-        x : ND-array
-            Initial guess of unknown variable
-
-        """
-        self.build_A()
-        self.build_b()
-        if x is None:
-            x = sp.zeros(shape=[self.Np, ], dtype=float)
-        # Scan through all source objects registered on algorithm
-        for item in self.settings['sources']:
-            source = self.project[item]  # Obtain handle to source object
-            source.apply()  # Apply source object to update A and b
-        x_new = self.solve()
-        self[self.settings['quantity']] = x_new
-        res = sp.sum(sp.absolute(x**2 - x_new**2))
-        if res < self.settings['tolerance']:
-            print('Tolerance met, solution converged')
-            return
-        else:
-            print('Tolerance not met: ', res)
-            self.run(x=x_new)
-
-    def solve(self, A=None, b=None):
-        r"""
-        Sends the A and b matrices to the specified solver.
-
-        Parameters
-        ----------
-        A : sparse matrix
-            The coefficient matrix in sparse format
-
-        b : ND-array
-            The RHS matrix in any format
-
-        Notes
-        -----
-        The solver used here is specified in the ``settings`` attribute of the
-        algorithm.
-
-        """
-        if A is None:
-            A = self.A
-        if b is None:
-            b = self.b
-        solver = getattr(sprs.linalg, self.settings['solver'])
-        x = solver(A=A.tocsr(), b=b)
-        return x
-
-    def rate(self, pores=None, mode='group'):
-        r"""
-        Calculates the net rate of material moving into a given set of pores.
-
-        Parameters
-        ----------
-        pores : array_like
-            The pores for which the rate should be calculated
-
-        mode : string, optional
-            Controls how to return the rate.  Options are:
-
-            *'group'*: (default) Teturns the cumulative rate of material
-            moving into the given set of pores
-
-            *'single'* : Calculates the rate for each pore individually
-
-        Notes
-        -----
-        A negative rate indicates material moving into the pore or pores, such
-        as material being consumed.
-        """
-        network = self.project.network
-        phase = self.project.phases[self.settings['phase']]
-        conductance = phase[self.settings['conductance']]
-        quantity = self[self.settings['quantity']]
-        pores = self._parse_indices(pores)
-        R = []
-        if mode == 'group':
-            t = network.find_neighbor_throats(pores, flatten=True,
-                                              mode='not_intersection')
-            throat_group_num = 1
-        elif mode == 'single':
-            t = network.find_neighbor_throats(pores, flatten=False,
-                                              mode='not_intersection')
-            throat_group_num = sp.shape(t)[0]
-        for i in sp.r_[0: throat_group_num]:
-            if mode == 'group':
-                throats = t
-                P = pores
-            elif mode == 'single':
-                throats = t[i]
-                P = pores[i]
-            p1 = network.find_connected_pores(throats)[:, 0]
-            p2 = network.find_connected_pores(throats)[:, 1]
-            pores1 = sp.copy(p1)
-            pores2 = sp.copy(p2)
-            # Changes to pores1 and pores2 to make them as inner/outer pores
-            pores1[~sp.in1d(p1, P)] = p2[~sp.in1d(p1, P)]
-            pores2[~sp.in1d(p1, P)] = p1[~sp.in1d(p1, P)]
-            X1 = quantity[pores1]
-            X2 = quantity[pores2]
-            g = conductance[throats]
-            R.append(sp.sum(sp.multiply(g, (X2 - X1))))
-        return sp.array(R, ndmin=1)
-
-    def _calc_eff_prop(self):
-        r"""
-        Returns the main parameters for calculating the effective
-        property in a linear transport equation.  It also checks for the
-        proper boundary conditions, inlets and outlets.
-        """
-        network = self.project.network
-        if self.settings['quantity'] not in self.keys():
-            raise Exception('The algorit hm has not been run yet. Cannot ' +
-                            'calculate effective property.')
-
-        # Determine boundary conditions by analyzing algorithm object
-        Ps = self.pores('pore.dirichlet')
-        BCs = sp.unique(self['pore.dirichlet_value'][Ps])
-        inlets = sp.where(self['pore.dirichlet_value'] == sp.amax(BCs))[0]
-        outlets = sp.where(self['pore.dirichlet_value'] == sp.amin(BCs))[0]
-
-        # Fetch area and length of domain
-        A = network.domain_area(face=inlets)
-        L = network.domain_length(face_1=inlets, face_2=outlets)
-        flow = self.rate(pores=inlets)
-        D = sp.sum(flow)*L/A/(BCs[0] - BCs[1])
-        return D
diff --git a/openpnm/algorithms/GenericTransport.py b/openpnm/algorithms/GenericTransport.py
index 1df64b24c..109074425 100644
--- a/openpnm/algorithms/GenericTransport.py
+++ b/openpnm/algorithms/GenericTransport.py
@@ -1,6 +1,9 @@
 import numpy as np
 import scipy.sparse as sprs
 import scipy.sparse.csgraph as spgr
+from scipy.spatial import ConvexHull
+from scipy.spatial import cKDTree
+from openpnm.topotools import iscoplanar
 from openpnm.algorithms import GenericAlgorithm
 from openpnm.core import logging
 logger = logging.getLogger(__name__)
@@ -17,28 +20,81 @@ class GenericTransport(GenericAlgorithm):
                               'conductance': None,
                               'quantity': None,
                               'solver': 'spsolve'})
-        self.settings.update(settings)
-        if phase is not None:
-            self.settings['phase'] = phase.name
-
+        self.setup(phase=phase, **settings)
+        # If network given, get project, otherwise let parent class create it
         if network is not None:
             project = network.project
-
         super().__init__(project=project, **kwargs)
+        # Create some instance attributes
+        self.A = None
+        self._pure_A = None
+        self.b = None
+        self._pure_b = None
+
+    def setup(self, phase=None, **kwargs):
+        r"""
+        This method takes several arguments that are essential to running the
+        algorithm and adds them to the settings.
+
+        Notes
+        -----
+        This generic version should be subclassed, and the arguments given
+        suitable default names.
+        """
+        if phase:
+            self.settings['phase'] = phase.name
+        self.settings.update(kwargs)
 
     def set_dirichlet_BC(self, pores, values):
         r"""
+        Apply *Dirichlet* boundary conditons to the specified pore locations.
+        Dirichlet conditions refer to constant *quantity* (e.g. pressure).
+
+        Parameters
+        ----------
+        pores : array_like
+            The pore indices where the condition should be applied
+
+        values : scalar or array_like
+            The value to of the boundary condition.  If a scalar is supplied
+            it is assigne to all locations, and if a vector is applied it
+            corresponds directy to the locations given in ``pores``.
+
+        Notes
+        -----
+        The definition of ``quantity`` is specified in the algorithm's
+        ``settings``, e.g. ``alg.settings['quentity'] = 'pore.pressure'``.
         """
-        self.set_BC(pores=pores, bctype='dirichlet', bcvalues=values,
-                    mode='merge')
+        self._set_BC(pores=pores, bctype='dirichlet', bcvalues=values,
+                     mode='merge')
 
     def set_neumann_BC(self, pores, values):
         r"""
+        Apply *Neumann*-type boundary conditons to the specified pore
+        locations. Neumann conditions technnically refer to the *gradient* of
+        the *quantity* (e.g. dP/dz), while in OpenPNM this means the flow-rate
+        of the quantity, which is the gradient multiplied by the conductance
+        (e.g n = g dP/dz).
+
+        Parameters
+        ----------
+        pores : array_like
+            The pore indices where the condition should be applied
+
+        values : scalar or array_like
+            The value to of the boundary condition.  If a scalar is supplied
+            it is assigne to all locations, and if a vector is applied it
+            corresponds directy to the locations given in ``pores``.
+
+        Notes
+        -----
+        The definition of ``quantity`` is specified in the algorithm's
+        ``settings``, e.g. ``alg.settings['quentity'] = 'pore.pressure'``.
         """
-        self.set_BC(pores=pores, bctype='neumann', bcvalues=values,
-                    mode='merge')
+        self._set_BC(pores=pores, bctype='neumann', bcvalues=values,
+                     mode='merge')
 
-    def set_BC(self, pores, bctype, bcvalues=None, mode='merge'):
+    def _set_BC(self, pores, bctype, bcvalues=None, mode='merge'):
         r"""
         Apply boundary conditions to specified pores
 
@@ -69,7 +125,7 @@ class GenericTransport(GenericAlgorithm):
         Notes
         -----
         It is not possible to have multiple boundary conditions for a
-        specified location in just one algorithm. Use ``mode='remove'`` to
+        specified location in one algorithm. Use ``mode='remove'`` to
         clear existing BCs before applying new ones or ``mode='overwrite'``
         which removes all existing BC's before applying the new ones.
 
@@ -79,8 +135,7 @@ class GenericTransport(GenericAlgorithm):
 
         """
         # Hijack the parse_mode function to verify bctype argument
-        bctype = self._parse_mode(bctype, allowed=['dirichlet', 'neumann',
-                                                   'neumann_group'],
+        bctype = self._parse_mode(bctype, allowed=['dirichlet', 'neumann'],
                                   single=True)
         mode = self._parse_mode(mode, allowed=['merge', 'overwrite', 'remove'],
                                 single=True)
@@ -122,30 +177,64 @@ class GenericTransport(GenericAlgorithm):
             self['pore.neumann'][pores] = False
             self['pore.neumann_value'][pores] = np.nan
 
-    def build_A(self):
+    def _build_A(self, force=False):
         r"""
-        """
-        network = self.project.network
-        phase = self.project.phases()[self.settings['phase']]
-        g = phase[self.settings['conductance']]
-        am = network.create_adjacency_matrix(weights=g, fmt='coo')
-        A = spgr.laplacian(am)
-        self.A = A
-        return A
+        Builds the coefficient matrix based on conductances between pores.
+        The conductance to use is specified in the algorithm's ``settings``
+        under ``conductance``.  In subclasses (e.g. ``FickianDiffusion``)
+        this is set by default, though it can be overwritten.
 
-    def build_b(self):
+        Parameters
+        ----------
+        force : Boolean (default is ``False)
+            If set to ``True`` then the A matrix is built from new.  If
+            ``False`` (the default), a cached version of A is returned.  The
+            cached version is *clean* in the sense that no boundary conditions
+            or sources terms have been added to it.
+        """
+        if force:
+            self._pure_A = None
+        if self._pure_A is None:
+            network = self.project.network
+            phase = self.project.phases()[self.settings['phase']]
+            g = phase[self.settings['conductance']]
+            am = network.create_adjacency_matrix(weights=g, fmt='coo')
+            self._pure_A = spgr.laplacian(am)
+        self.A = self._pure_A.copy()
+
+    def _build_b(self, force=False):
         r"""
+        Builds the RHS matrix, without applying any boundary conditions or
+        source terms. This method is trivial an basically creates a column
+        vector of 0's.
+
+        Parameters
+        ----------
+        force : Boolean (default is ``False``)
+            If set to ``True`` then the b matrix is built from new.  If
+            ``False`` (the default), a cached version of b is returned.  The
+            cached version is *clean* in the sense that no boundary conditions
+            or sources terms have been added to it.
         """
-        b = np.zeros(shape=(self.Np, ), dtype=float)  # Create b matrix of 0's
-        self.b = b
-        return b
+        if force:
+            self._pure_b = None
+        if self._pure_b is None:
+            b = np.zeros(shape=(self.Np, ), dtype=float)  # Create vector of 0s
+            self._pure_b = b
+        self.b = self._pure_b.copy()
+
+    def _apply_BCs(self):
+        r"""
+        Applies all the boundary conditions that have been specified, by
+        adding values to the *A* and *b* matrices.
 
-    def apply_BCs(self):
+        """
+        self._build_A()
+        self._build_b()
         if 'pore.neumann' in self.keys():
             # Update b
             ind = self['pore.neumann']
             self.b[ind] = self['pore.neumann_value'][ind]
-
         if 'pore.dirichlet' in self.keys():
             # Update b
             ind = self['pore.dirichlet']
@@ -160,13 +249,6 @@ class GenericTransport(GenericAlgorithm):
             self.A.setdiag(datadiag)
             self.A.eliminate_zeros()  # Remove 0 entries
 
-    def setup(self):
-        r"""
-        """
-        self.build_A()
-        self.build_b()
-        self.apply_BCs()
-
     def run(self):
         r"""
         Builds the A and b matrices, and calls the solver specified in the
@@ -177,27 +259,38 @@ class GenericTransport(GenericAlgorithm):
         x : ND-array
             Initial guess of unknown variable
 
+        Returns
+        -------
+        Nothing is returned...the solution is stored on the objecxt under
+        ``pore.quantity`` where *quantity* is specified in the ``settings``
+        attribute.
+
         """
         print('―'*80)
         print('Running GenericTransport')
-        self.setup()
         self._run_generic()
 
     def _run_generic(self):
+        self._apply_BCs()
         x_new = self._solve()
         self[self.settings['quantity']] = x_new
 
     def _solve(self, A=None, b=None):
         r"""
-        Sends the A and b matrices to the specified solver.
+        Sends the A and b matrices to the specified solver, and solves for *x*
+        given the boundary conditions, and source terms based on the present
+        value of *x*.  This method does NOT iterate to solve for non-linear
+        source terms or march time steps.
 
         Parameters
         ----------
         A : sparse matrix
-            The coefficient matrix in sparse format
+            The coefficient matrix in sparse format. If not specified, then
+            it uses  the ``A`` matrix attached to the object.
 
         b : ND-array
-            The RHS matrix in any format
+            The RHS matrix in any format.  If not specified, then it uses
+            the ``b`` matrix attached to the object.
 
         Notes
         -----
@@ -207,14 +300,20 @@ class GenericTransport(GenericAlgorithm):
         """
         if A is None:
             A = self.A
+            if A is None:
+                raise Exception('The A matrix has not been built yet')
         if b is None:
             b = self.b
+            if b is None:
+                raise Exception('The b matrix has not been built yet')
         solver = getattr(sprs.linalg, self.settings['solver'])
         x = solver(A=A.tocsr(), b=b)
         return x
 
     def results(self):
         r"""
+        Fetches the calculated quantity from the algorithm and returns it as
+        an array.
         """
         quantity = self.settings['quantity']
         d = {quantity: self[quantity]}
@@ -232,7 +331,7 @@ class GenericTransport(GenericAlgorithm):
         mode : string, optional
             Controls how to return the rate.  Options are:
 
-            *'group'*: (default) Teturns the cumulative rate of material
+            *'group'*: (default) Returns the cumulative rate of material
             moving into the given set of pores
 
             *'single'* : Calculates the rate for each pore individually
@@ -242,60 +341,113 @@ class GenericTransport(GenericAlgorithm):
         A negative rate indicates material moving into the pore or pores, such
         as material being consumed.
         """
+        pores = self._parse_indices(pores)
+
         network = self.project.network
-        phase = self.project.phases[self.settings['phase']]
+        phase = self.project.phases()[self.settings['phase']]
         conductance = phase[self.settings['conductance']]
         quantity = self[self.settings['quantity']]
-        pores = self._parse_indices(pores)
-        R = []
-        if mode == 'group':
-            t = network.find_neighbor_throats(pores, flatten=True,
-                                              mode='not_intersection')
-            throat_group_num = 1
-        elif mode == 'single':
-            t = network.find_neighbor_throats(pores, flatten=False,
-                                              mode='not_intersection')
-            throat_group_num = np.shape(t)[0]
-        for i in np.r_[0: throat_group_num]:
-            if mode == 'group':
-                throats = t
-                P = pores
-            elif mode == 'single':
-                throats = t[i]
-                P = pores[i]
-            p1 = network.find_connected_pores(throats)[:, 0]
-            p2 = network.find_connected_pores(throats)[:, 1]
-            pores1 = np.copy(p1)
-            pores2 = np.copy(p2)
-            # Changes to pores1 and pores2 to make them as inner/outer pores
-            pores1[~np.in1d(p1, P)] = p2[~np.in1d(p1, P)]
-            pores2[~np.in1d(p1, P)] = p1[~np.in1d(p1, P)]
-            X1 = quantity[pores1]
-            X2 = quantity[pores2]
-            g = conductance[throats]
-            R.append(np.sum(np.multiply(g, (X2 - X1))))
+
+        P12 = network['throat.conns']
+        X12 = quantity[P12]
+        f = (-1)**np.argsort(X12, axis=1)[:, 1]
+        g = conductance
+        Dx = np.abs(np.diff(X12, axis=1).squeeze())
+        Qt = f*g*Dx
+
+        if mode == 'single':
+            Qp = np.zeros((self.Np, ))
+            np.add.at(Qp, P12[:, 0], Qt)
+            np.add.at(Qp, P12[:, 1], -Qt)
+            R = Qp[pores]
+        elif mode == 'group':
+            Ts = network.find_neighbor_throats(pores=pores,
+                                               mode='exclusive_or')
+            R = np.sum(Qt[Ts])
         return np.array(R, ndmin=1)
 
     def _calc_eff_prop(self):
         r"""
-        Returns the main parameters for calculating the effective
-        property in a linear transport equation.  It also checks for the
-        proper boundary conditions, inlets and outlets.
+        Returns the main parameters for calculating the effective property
+        in a linear transport equation.  It also checks for the proper
+        boundary conditions, inlets and outlets.
         """
         network = self.project.network
         if self.settings['quantity'] not in self.keys():
-            raise Exception('The algorit hm has not been run yet. Cannot ' +
+            raise Exception('The algorithm has not been run yet. Cannot ' +
                             'calculate effective property.')
 
         # Determine boundary conditions by analyzing algorithm object
+        inlets, outlets = self._get_inlets_and_outlets()
         Ps = self.pores('pore.dirichlet')
         BCs = np.unique(self['pore.dirichlet_value'][Ps])
-        inlets = np.where(self['pore.dirichlet_value'] == np.amax(BCs))[0]
-        outlets = np.where(self['pore.dirichlet_value'] == np.amin(BCs))[0]
 
         # Fetch area and length of domain
-        A = network.domain_area(face=inlets)
-        L = network.domain_length(face_1=inlets, face_2=outlets)
+        A = self.domain_area
+        L = self.domain_length
         flow = self.rate(pores=inlets)
         D = np.sum(flow)*L/A/(BCs[0] - BCs[1])
         return D
+
+    def _get_inlets_and_outlets(self):
+        # Determine boundary conditions by analyzing algorithm object
+        Ps = self.pores('pore.dirichlet')
+        BCs = np.unique(self['pore.dirichlet_value'][Ps])
+        inlets = np.where(self['pore.dirichlet_value'] == np.amax(BCs))[0]
+        outlets = np.where(self['pore.dirichlet_value'] == np.amin(BCs))[0]
+        return (inlets, outlets)
+
+    def _get_domain_area(self):
+        if not hasattr(self, '_area'):
+            logger.warning('Attempting to estimate inlet area...will be low')
+            network = self.project.network
+            Pin, Pout = self._get_inlets_and_outlets()
+            inlets = network['pore.coords'][Pin]
+            outlets = network['pore.coords'][Pout]
+            if not iscoplanar(inlets):
+                raise Exception('Detected inlet pores are not coplanar')
+            if not iscoplanar(outlets):
+                raise Exception('Detected outlet pores are not coplanar')
+            Nin = np.ptp(inlets, axis=0) > 0
+            if Nin.all():
+                raise Exception('Detected inlets are not oriented along a ' +
+                                'principle axis')
+            Nout = np.ptp(outlets, axis=0) > 0
+            if Nout.all():
+                raise Exception('Detected outlets are not oriented along a ' +
+                                'principle axis')
+            hull_in = ConvexHull(points=inlets[:, Nin])
+            hull_out = ConvexHull(points=outlets[:, Nout])
+            if hull_in.volume != hull_out.volume:
+                raise Exception('Inlet and outlet faces are different area')
+            self._area = hull_in.volume  # In 2D volume=area, area=perimeter
+        return self._area
+
+    def _set_domain_area(self, area):
+        self._area = area
+
+    domain_area = property(fget=_get_domain_area, fset=_set_domain_area)
+
+    def _get_domain_length(self):
+        if not hasattr(self, '_length'):
+            logger.warning('Attempting to estimate domain length... ' +
+                           'could be low if boundary pores were not added')
+            network = self.project.network
+            Pin, Pout = self._get_inlets_and_outlets()
+            inlets = network['pore.coords'][Pin]
+            outlets = network['pore.coords'][Pout]
+            if not iscoplanar(inlets):
+                raise Exception('Detected inlet pores are not coplanar')
+            if not iscoplanar(outlets):
+                raise Exception('Detected inlet pores are not coplanar')
+            tree = cKDTree(data=inlets)
+            Ls = np.unique(np.around(tree.query(x=outlets)[0], decimals=5))
+            if np.size(Ls) != 1:
+                raise Exception('A unique value of length could not be found')
+            self._length = Ls[0]
+        return self._length
+
+    def _set_domain_length(self, length):
+        self._length = length
+
+    domain_length = property(fget=_get_domain_length, fset=_set_domain_length)
diff --git a/openpnm/algorithms/OhmicConduction.py b/openpnm/algorithms/OhmicConduction.py
index 23f17a9e0..f9c8d83cc 100644
--- a/openpnm/algorithms/OhmicConduction.py
+++ b/openpnm/algorithms/OhmicConduction.py
@@ -1,9 +1,9 @@
-from openpnm.algorithms import GenericTransport
+from openpnm.algorithms import ReactiveTransport
 from openpnm.core import logging
 logger = logging.getLogger(__name__)
 
 
-class OhmicConduction(GenericTransport):
+class OhmicConduction(ReactiveTransport):
     r"""
     A subclass of GenericLinearTransport to simulate electron and ionic
     conduction.  The 2 main roles of this subclass are to set the default
@@ -18,6 +18,43 @@ class OhmicConduction(GenericTransport):
                               'conductance': 'throat.electrical_conductance'})
         self.settings.update(settings)
 
+    def setup(self, phase=None, quantity='', conductance='', **kwargs):
+        r"""
+        This method takes several arguments that are essential to running the
+        algorithm and adds them to the settings.
+
+        Parameters
+        ----------
+        phase : OpenPNM Phase object
+            The phase on which the algorithm is to be run.  If no value is
+            given, the existing value is kept.
+
+        quantity : string
+            The name of the physical quantity to be calcualted.  If no value is
+            given, the existing value is kept.  The default value is
+            ``'pore.voltage'``.
+
+        conductance : string
+            The name of the pore-scale transport conductance values.  These
+            are typically calculate by a model attached to a *Physics* object
+            associated with the given *Phase*.  If no value is given, the
+            existing value is kept.  The default value is
+            ``'throat.electrical_conductance'``.
+
+        Notes
+        -----
+        Any additional arguments are added to the ``settings`` dictionary of
+        the object.
+
+        """
+        if phase:
+            self.settings['phase'] = phase.name
+        if quantity:
+            self.settings['quantity': quantity]
+        if conductance:
+            self.settings['conductance': conductance]
+        super().setup(**kwargs)
+
     def calc_effective_conductivity(self):
         r"""
         This calculates the effective electrical conductivity.
diff --git a/openpnm/algorithms/Reaction.py b/openpnm/algorithms/Reaction.py
deleted file mode 100644
index 834b33d75..000000000
--- a/openpnm/algorithms/Reaction.py
+++ /dev/null
@@ -1,37 +0,0 @@
-import scipy as sp
-from openpnm.core import ModelsMixin
-from openpnm.algorithms import GenericAlgorithm
-
-
-class GenericReaction(GenericAlgorithm, ModelsMixin):
-
-    def __init__(self, network, pores, settings={}, **kwargs):
-        self.settings.update({'rate_model': 'pore.rate',
-                              'prefix': 'rxn'})
-        self.settings.update(settings)
-        super().__init__(Np=sp.size(pores), network=network, **kwargs)
-        self['pore._id'] = network['pore._id'][pores]
-
-    def setup(self, algorithm):
-        self.settings['algorithm'] = algorithm.name
-        self.settings['quantity'] = algorithm.settings['quantity']
-
-    def apply(self, A=None, b=None):
-        net = self.project.network
-        Ps = net.map_pores(self['pore._id'])
-        # Fetch algorithm object from project
-        alg = self.project[self.settings['algorithm']]
-        if A is None:
-            A = alg.A
-        if b is None:
-            b = alg.b
-        quantity = alg.settings['quantity']
-        x = alg[quantity].copy()
-        self[quantity] = x[Ps]
-        # Regenerate models with new guess
-        self.regenerate_models()
-        # Add S1 to diagonal of A
-        datadiag = A.diagonal()
-        datadiag[Ps] = datadiag[Ps] + self[self.settings['rate_model']][:, 1]
-        A.setdiag(datadiag)
-        b[Ps] = b[Ps] - self[self.settings['rate_model']][:, 2]
diff --git a/openpnm/algorithms/ReactiveTransport.py b/openpnm/algorithms/ReactiveTransport.py
index 2aa1086b2..931e4f35d 100644
--- a/openpnm/algorithms/ReactiveTransport.py
+++ b/openpnm/algorithms/ReactiveTransport.py
@@ -13,29 +13,43 @@ class ReactiveTransport(GenericTransport):
                               'tolerance': 0.001})
         super().__init__(**kwargs)
 
-    def set_source_term(self, propname, pores):
+    def set_source(self, propname, pores):
         r"""
+        Applies a given source term to the specified pores
+
+        Parameters
+        ----------
+        propname : string
+            The property name of the source term model to be applied
+
+        pores : array_like
+            The pore indices where the source term should be applied
+
         """
         self.settings['sources'].append(propname)
         self[propname] = self.tomask(pores=pores)
 
-    def apply_sources(self):
+    def _apply_sources(self):
         phase = self.project.phases()[self.settings['phase']]
+        physics = self.project.find_physics(phase=phase)
         for item in self.settings['sources']:
             Ps = self.pores(item)
             # Regenerate models with new guess
             quantity = self.settings['quantity']
             # Put quantity on phase so physics finds it when regenerating
             phase[quantity] = self[quantity]
+            # Regenerate models, on either phase or physics
             phase.regenerate_models(propnames=item)
+            for phys in physics:
+                phys.regenerate_models(propnames=item)
             # Add S1 to diagonal of A
             # TODO: We need this to NOT overwrite the A and b, but create
             # copy, otherwise we have to regenerate A and b on each loop
             datadiag = self.A.diagonal()
-            datadiag[Ps] = datadiag[Ps] + phase[item+'_S1'][Ps]
+            datadiag[Ps] = datadiag[Ps] + phase[item+'.'+'S1'][Ps]
             self.A.setdiag(datadiag)
             # Add S2 to b
-            self.b[Ps] = self.b[Ps] - phase[item+'_S2'][Ps]
+            self.b[Ps] = self.b[Ps] - phase[item+'.'+'S2'][Ps]
 
     def run(self, x=None):
         r"""
@@ -50,14 +64,14 @@ class ReactiveTransport(GenericTransport):
         """
         print('―'*80)
         print('Running ReactiveTransport')
-        self.setup()
-        self._run_reactive(x=x)
+        x = self._run_reactive(x=x)
+        return x
 
     def _run_reactive(self, x):
         if self.settings['quantity'] not in self.keys():
             self[self.settings['quantity']] = 0
-        self.setup()
-        self.apply_sources()
+        self._apply_BCs()
+        self._apply_sources()
         if x is None:
             x = np.zeros(shape=[self.Np, ], dtype=float)
         x_new = self._solve()
@@ -65,7 +79,7 @@ class ReactiveTransport(GenericTransport):
         res = np.sum(np.absolute(x**2 - x_new**2))
         if res < self.settings['tolerance']:
             print('Solution converged: ' + str(res))
-            return
+            return x_new
         else:
             print('Tolerance not met: ' + str(res))
             self._run_reactive(x=x_new)
diff --git a/openpnm/algorithms/StokesFlow.py b/openpnm/algorithms/StokesFlow.py
index 875244a50..63ca7f868 100644
--- a/openpnm/algorithms/StokesFlow.py
+++ b/openpnm/algorithms/StokesFlow.py
@@ -1,10 +1,10 @@
 import scipy as sp
-from openpnm.algorithms import GenericTransport
+from openpnm.algorithms import ReactiveTransport
 from openpnm.core import logging
 logger = logging.getLogger(__name__)
 
 
-class StokesFlow(GenericTransport):
+class StokesFlow(ReactiveTransport):
     r"""
     A subclass of GenericLinearTransport to simulate viscous flow.  The 2
     main roles of this subclass are to set the default property names and to
@@ -19,12 +19,49 @@ class StokesFlow(GenericTransport):
                               'conductance': 'throat.hydraulic_conductance'})
         self.settings.update(settings)
 
+    def setup(self, phase=None, quantity='', conductance='', **kwargs):
+        r"""
+        This method takes several arguments that are essential to running the
+        algorithm and adds them to the settings.
+
+        Parameters
+        ----------
+        phase : OpenPNM Phase object
+            The phase on which the algorithm is to be run.  If no value is
+            given, the existing value is kept.
+
+        quantity : string
+            The name of the physical quantity to be calcualted.  If no value is
+            given, the existing value is kept.  The default value is
+            ``'pore.pressure'``.
+
+        conductance : string
+            The name of the pore-scale transport conductance values.  These
+            are typically calculate by a model attached to a *Physics* object
+            associated with the given *Phase*.  If no value is given, the
+            existing value is kept.  The default value is
+            ``'throat.hydraulic_conductance'``.
+
+        Notes
+        -----
+        Any additional arguments are added to the ``settings`` dictionary of
+        the object.
+
+        """
+        if phase:
+            self.settings['phase'] = phase.name
+        if quantity:
+            self.settings['quantity': quantity]
+        if conductance:
+            self.settings['conductance': conductance]
+        super().setup(**kwargs)
+
     def calc_eff_permeability(self):
         r"""
-        This calculates the effective permeability in this linear
-        transport algorithm.
+        This calculates the effective permeability in this linear transport
+        algorithm.
         """
-        phase = self.project.phases[self.settings['phase']]
+        phase = self.project.phases()[self.settings['phase']]
         d_normal = self._calc_eff_prop()
-        self._eff_property = d_normal / sp.mean(phase['pore.viscosity'])
+        self._eff_property = d_normal * sp.mean(phase['pore.viscosity'])
         return self._eff_property
diff --git a/openpnm/algorithms/TransientReactiveTransport.py b/openpnm/algorithms/TransientReactiveTransport.py
index da8a7d6d6..43cd17cde 100644
--- a/openpnm/algorithms/TransientReactiveTransport.py
+++ b/openpnm/algorithms/TransientReactiveTransport.py
@@ -17,12 +17,12 @@ class TransientReactiveTransport(ReactiveTransport, TransientTransport):
         """
         print('―'*80)
         print('Running TransientReactiveTransport')
-        self.setup()
         x = self._run_transient_reactive(x=x, t=t)
 
     def _run_transient_reactive(self, x, t):
-        self.update_A()
-        self.update_b()
+        self._apply_BCs()
+        self._update_A()
+        self._update_b()
         x = self._run_reactive(x=x)
         if t < self.settings['t_final']:
             print('Current time step: '+str(t))
diff --git a/openpnm/algorithms/TransientTransport.py b/openpnm/algorithms/TransientTransport.py
index 2f51c66cf..f2b60920a 100644
--- a/openpnm/algorithms/TransientTransport.py
+++ b/openpnm/algorithms/TransientTransport.py
@@ -18,12 +18,12 @@ class TransientTransport(GenericTransport):
         """
         self[self.settings['quantity']] = values
 
-    def update_A(self):
+    def _update_A(self):
         r"""
         """
         pass
 
-    def update_b(self):
+    def _update_b(self):
         r"""
         """
         pass
@@ -33,12 +33,12 @@ class TransientTransport(GenericTransport):
         """
         print('―'*80)
         print('Running TransientTransport')
-        self.setup()
         self._run_transient(t=t)
 
     def _run_transient(self, t):
-        self.update_A()
-        self.update_b()
+        self._apply_BCs()
+        self._update_A()
+        self._update_b()
         x_new = self._solve()
         self[self.settings['quantity']] = x_new
         if t < self.settings['t_final']:
diff --git a/openpnm/algorithms/__init__.py b/openpnm/algorithms/__init__.py
index 0be732895..153486f7b 100644
--- a/openpnm/algorithms/__init__.py
+++ b/openpnm/algorithms/__init__.py
@@ -3,7 +3,6 @@ from .GenericTransport import GenericTransport
 from .ReactiveTransport import ReactiveTransport
 from .TransientTransport import TransientTransport
 from .TransientReactiveTransport import TransientReactiveTransport
-from .CrankNicholson import CrankNicholson
 from .StokesFlow import StokesFlow
 from .FickianDiffusion import FickianDiffusion
 from .FourierConduction import FourierConduction
diff --git a/openpnm/algorithms/models.py b/openpnm/algorithms/models.py
deleted file mode 100644
index 46ce69ecd..000000000
--- a/openpnm/algorithms/models.py
+++ /dev/null
@@ -1,16 +0,0 @@
-import scipy as sp
-
-
-def standard_kinetics(target, quantity, prefactor, exponent):
-    r"""
-
-    """
-    X = target[quantity]
-    A = target[prefactor]
-    b = target[exponent]
-
-    r = A*(X**b)
-    S1 = A*b*(X**(b - 1))
-    S2 = A*(1 - b)*(X**b)
-    values = sp.vstack((r, S1, S2)).T
-    return values
diff --git a/openpnm/core/Base.py b/openpnm/core/Base.py
index 2958f3ed3..d87a35e7a 100644
--- a/openpnm/core/Base.py
+++ b/openpnm/core/Base.py
@@ -30,9 +30,7 @@ class Base(dict):
         self.update({'throat.all': sp.ones(shape=(Nt, ), dtype=bool)})
 
     def __repr__(self):
-        return '<%s object at %s>' % (
-            self.__class__.__module__,
-            hex(id(self)))
+        return '<%s object at %s>' % (self.__class__.__module__, hex(id(self)))
 
     def __eq__(self, other):
         if hex(id(self)) == hex(id(other)):
@@ -53,7 +51,7 @@ class Base(dict):
         if hasattr(value, 'keys'):
             for item in value.keys():
                 prop = item.replace('pore.', '').replace('throat.', '')
-                self.__setitem__(key+'_'+prop, value[item])
+                self.__setitem__(key+'.'+prop, value[item])
             return
 
         value = sp.array(value, ndmin=1)  # Convert value to an ndarray
@@ -222,9 +220,10 @@ class Base(dict):
         if v is None:
             if hasattr(self, 'models'):
                 if key in self.models.keys():
-                    logger.info(key+' not found, regenerating model')
-                    self.regenerate_models(propnames=key)
-                    v = self.get(key=key, value=value)
+                    if self.models[key]['regen_mode'] in ['deferred', 'lazy']:
+                        logger.info(key+' not found, regenerating model')
+                        self.regenerate_models(propnames=key)
+                        v = self.get(key=key, value=value)
         return v
 
     # -------------------------------------------------------------------------
diff --git a/openpnm/core/ModelsMixin.py b/openpnm/core/ModelsMixin.py
index f38218532..3e56b5e18 100644
--- a/openpnm/core/ModelsMixin.py
+++ b/openpnm/core/ModelsMixin.py
@@ -160,6 +160,8 @@ class ModelsMixin():
         else:
             # Re-create propnames to ensure it's in correct order
             propnames = [i for i in all_props if i in propnames]
+        if len(propnames) == 0:
+            logger.info('List of propnames to regenerate is empty')
         # Scan through list of propnames and regenerate each one
         for item in propnames:
             logger.info('Regenerating model: '+item)
diff --git a/openpnm/core/Subdomain.py b/openpnm/core/Subdomain.py
index 671da5a5a..bc73b243f 100644
--- a/openpnm/core/Subdomain.py
+++ b/openpnm/core/Subdomain.py
@@ -11,16 +11,9 @@ class Subdomain(Base):
             boss = self.project.find_phase(self)
         else:
             boss = self.project.network
-        # Deal with a few special key items
-        if key.split('.')[-1] == '_id':
-            inds = boss._get_indices(element=element, labels=self.name)
-            return boss[element+'._id'][inds]
-        # Convert self.name into 'all'
-        elif key.split('.')[-1] in [self.name]:
-            return self[element+'.all']
-        # Now get values if present, or regenerate them
+        # Get values if present, or regenerate them
         vals = self.get(key)
-        # If still not found, check with boss object (interleave data)
+        # If still not found, check with boss object
         if vals is None:
             inds = boss._get_indices(element=element, labels=self.name)
             vals = boss[key][inds]
diff --git a/openpnm/models/physics/generic_source_term.py b/openpnm/models/physics/generic_source_term.py
index 7663f88f3..2ec020d3a 100644
--- a/openpnm/models/physics/generic_source_term.py
+++ b/openpnm/models/physics/generic_source_term.py
@@ -1,10 +1,3 @@
-r"""
-===============================================================================
-Submodule -- generic_source_term
-===============================================================================
-
-"""
-
 import scipy as _sp
 
 
@@ -19,19 +12,23 @@ def standard_kinetics(target, quantity, prefactor, exponent):
     r = A*(X**b)
     S1 = A*b*(X**(b - 1))
     S2 = A*(1 - b)*(X**b)
-    values = {'pore.S1': S1, 'pore.S2': S2, 'pore.rate': r}
+    values = {'S1': S1, 'S2': S2, 'rate': r}
     return values
 
 
-def linear(physics, phase, A1='', A2='', x='', return_rate=True, **kwargs):
+def linear(target, A1='', A2='', quantity=''):
     r"""
     For the following source term:
+
         .. math::
             r = A_{1}   x  +  A_{2}
+
     If return_rate is True, it returns the value of source term for the
     provided x in each pore.
+
     If return_rate is False, it calculates the slope and intercept for the
     following linear form :
+
         .. math::
             r = S_{1}   x  +  S_{2}
 
@@ -40,75 +37,46 @@ def linear(physics, phase, A1='', A2='', x='', return_rate=True, **kwargs):
     A1 , A2 : string
         The property name of the coefficients in the source term model.
         With A2 set to zero this equation takes on the familiar for of r=kx.
-    x : string or float/int or array/list
-        The property name or numerical value or array for the main quantity
-    Notes
-    -----
-    Because this source term is linear in concentration (x) is it not necessary
-    to iterate during the solver step.  Thus, when using the
-    ``set_source_term`` method for an algorithm, it is recommended to set the
-    ``maxiter``
-    argument to 0.  This will save 1 unncessary solution of the system, since
-    the solution would coverge after the first pass anyway.
+
+    quantity : string
+        The property name for the main quantity being produced or consumed.
+
+    Returns
+    -------
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
+
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
+
+    Starting in V2, OpenPNM allows the assigning of such ``dicts`` to a Core
+    object ``dict``, but it breaks the received ``dict`` into separate arrays
+    and prepends the assigned propname.  For instance:
+
+        obj['pore.source_term']  = d
+
+    is equivalent to:
+
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    if x is '':
-        X = _sp.ones(physics.Np) * _sp.nan
-    else:
-        if type(x) == str:
-            x = 'pore.' + x.split('.')[-1]
-            try:
-                X = physics[x]
-            except KeyError:
-                raise Exception(physics.name +
-                                ' does not have the pore property :' + x + '!')
-        else:
-            X = _sp.array(x)
-
-    length_X = _sp.size(X)
-    if length_X != physics.Np:
-        if length_X == 1:
-            X = X * _sp.ones(physics.Np)
-        elif length_X >= phase.Np:
-            X = X[physics.map_pores()]
-        else:
-            raise Exception('Wrong size for the numerical array of x!')
-
-    a = {}
-    source_params = [A1, A2]
-    for ind in _sp.arange(_sp.size(source_params)):
-        A = source_params[ind]
-        if A is '':
-            a[str(ind+1)] = 0
-        else:
-            if type(A) == str:
-                A = 'pore.' + A.split('.')[-1]
-                try:
-                    a[str(ind+1)] = physics[A]
-                except KeyError:
-                    raise Exception(physics.name + '/' + phase.name +
-                                    ' does not have the pore property :' +
-                                    A + '!')
-            else:
-                raise Exception('source_term parameters can only be string '
-                                'type!')
-
-    if return_rate:
-        return(a['1'] * X + a['2'])
-    else:
-        S1 = a['1']
-        S2 = a['2']
-        return(_sp.vstack((S1, S2)).T)
-
-
-def power_law(physics, phase, A1='', A2='', A3='', x='',
-              return_rate=True, **kwargs):
+    X = target[quantity]
+    r = target[A1] * X + target[A2]
+    S1 = target[A1]
+    S2 = target[A2]
+    values = {'S1': S1, 'S2': S2, 'rate': r}
+    return values
+
+
+def power_law(target, A1='', A2='', A3='', quantity=''):
     r"""
     For the following source term:
         .. math::
             r = A_{1}   x^{A_{2}}  +  A_{3}
+
     If return_rate is True, it returns the value of source term for the
     provided x in each pore.
+
     If return_rate is False, it calculates the slope and intercept for the
     following linear form :
         .. math::
@@ -118,69 +86,48 @@ def power_law(physics, phase, A1='', A2='', A3='', x='',
     ----------
     A1 -> A3 : string
         The property name of the coefficients in the source term model
-    x : string or float/int or array/list
-        The property name or numerical value or array for the main quantity
-    Notes
-    -----
+
+    quantity : string
+        The property name for the main quantity being produced or consumed.
+
+    Returns
+    -------
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
+
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
+
+    Starting in V2, OpenPNM allows the assigning of such ``dicts`` to a Core
+    object ``dict``, but it breaks the received ``dict`` into separate arrays
+    and prepends the assigned propname.  For instance:
+
+        obj['pore.source_term']  = d
+
+    is equivalent to:
+
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    if x is '':
-        X = _sp.ones(physics.Np) * _sp.nan
-    else:
-        if type(x) == str:
-            x = 'pore.' + x.split('.')[-1]
-            try:
-                X = physics[x]
-            except KeyError:
-                raise Exception(physics.name +
-                                ' does not have the pore property :' + x + '!')
-        else:
-            X = _sp.array(x)
-
-    length_X = _sp.size(X)
-    if length_X != physics.Np:
-        if length_X == 1:
-            X = X * _sp.ones(physics.Np)
-        elif length_X >= phase.Np:
-            X = X[physics.map_pores()]
-        else:
-            raise Exception('Wrong size for the numerical array of x!')
-
-    a = {}
-    source_params = [A1, A2, A3]
-    for ind in _sp.arange(_sp.size(source_params)):
-        A = source_params[ind]
-        if A is '':
-            a[str(ind+1)] = 0
-        else:
-            if type(A) == str:
-                A = 'pore.' + A.split('.')[-1]
-                try:
-                    a[str(ind+1)] = physics[A]
-                except KeyError:
-                    raise Exception(physics.name + '/' + phase.name +
-                                    ' does not have the pore property :' +
-                                    A + '!')
-            else:
-                raise Exception('source_term parameters can only be string '
-                                'type!')
-
-    if return_rate:
-        return(a['1'] * X ** a['2'] + a['3'])
-    else:
-        S1 = a['1'] * a['2'] * X ** (a['2'] - 1)
-        S2 = a['1'] * X ** a['2'] * (1 - a['2']) + a['3']
-        return(_sp.vstack((S1, S2)).T)
-
-
-def exponential(physics, phase, A1='', A2='', A3='', A4='', A5='', A6='',
-                x='', return_rate=True, **kwargs):
+    X = target[quantity]
+
+    r = target[A1] * X ** target[A2] + target[A3]
+    S1 = target[A1] * target[A2] * X ** (target[A2] - 1)
+    S2 = target[A1] * X ** target[A2] * (1 - target[A2]) + target[A3]
+    values = {'S1': S1, 'S2': S2, 'rate': r}
+    return values
+
+
+def exponential(target, A1='', A2='', A3='', A4='', A5='', A6='',
+                quantity='', return_rate=True):
     r"""
     For the following source term:
         .. math::
             r =  A_{1} A_{2}^{( A_{3} x^{ A_{4} } + A_{5})} + A_{6}
+
     If return_rate is True, it returns the value of source term for the
     provided x in each pore.
+
     If return_rate is False, it calculates the slope and intercept for the
     following linear form :
         .. math::
@@ -190,75 +137,56 @@ def exponential(physics, phase, A1='', A2='', A3='', A4='', A5='', A6='',
     ----------
     A1 -> A6 : string
         The property name of the coefficients in the source term model
-    x : string or float/int or array/list
-        The property name or numerical value or array for the main quantity
-    Notes
-    -----
+
+    quantity : string
+        The property name for the main quantity being produced or consumed.
+
+    Returns
+    -------
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
+
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
+
+    Starting in V2, OpenPNM allows the assigning of such ``dicts`` to a Core
+    object ``dict``, but it breaks the received ``dict`` into separate arrays
+    and prepends the assigned propname.  For instance:
+
+        obj['pore.source_term']  = d
+
+    is equivalent to:
+
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    if x is '':
-        X = _sp.ones(physics.Np) * _sp.nan
-    else:
-        if type(x) == str:
-            x = 'pore.'+x.split('.')[-1]
-            try:
-                X = physics[x]
-            except KeyError:
-                raise Exception(physics.name +
-                                ' does not have the pore property :' + x + '!')
-        else:
-            X = _sp.array(x)
-
-    length_X = _sp.size(X)
-    if length_X != physics.Np:
-        if length_X == 1:
-            X = X * _sp.ones(physics.Np)
-        elif length_X >= phase.Np:
-            X = X[physics.map_pores()]
-        else:
-            raise Exception('Wrong size for the numerical array of x!')
-
-    a = {}
-    source_params = [A1, A2, A3, A4, A5, A6]
-    for ind in _sp.arange(_sp.size(source_params)):
-        A = source_params[ind]
-        if A is '':
-            if ind == 0:
-                a[str(ind+1)] = 1
-            else:
-                a[str(ind+1)] = 0
+    X = target[quantity]
+    a = {'1': A1, '2': A2, '3': A3, '4': A4, '5': A5, '6': A6}
+    for item in a.keys():
+        if a[item] == '':
+            a[item] = 0
         else:
-            if type(A) == str:
-                A = 'pore.' + A.split('.')[-1]
-                try:
-                    a[str(ind+1)] = physics[A]
-                except KeyError:
-                    raise Exception(physics.name + '/' + phase.name +
-                                    ' does not have the pore property :' +
-                                    A + '!')
-            else:
-                raise Exception('source_term parameters can only be string '
-                                'type!')
-
-    if return_rate:
-        return a['1'] * a['2'] ** (a['3'] * X ** a['4'] + a['5']) + a['6']
-    else:
-        S1 = a['1'] * a['3'] * a['4'] * \
-            X ** (a['4'] - 1) * _sp.log(a['2']) * \
-            a['2'] ** (a['3'] * X ** a['4'] + a['5'])
-        S2 = a['1'] * a['2'] ** (a['3'] * X ** a['4'] + a['5']) * \
-            (1 - a['3'] * a['4'] * _sp.log(a['2']) * X ** a['4']) + a['6']
-        return(_sp.vstack((S1, S2)).T)
-
-
-def natural_exponential(physics, phase, A1='', A2='', A3='', A4='', A5='',
-                        x='', return_rate=True, **kwargs):
+            a[item] = target[a[item]]
+
+    r = a['1'] * a['2'] ** (a['3'] * X ** a['4'] + a['5']) + a['6']
+    S1 = a['1'] * a['3'] * a['4'] * X ** (a['4'] - 1) * _sp.log(a['2']) *\
+        a['2'] ** (a['3'] * X ** a['4'] + a['5'])
+    S2 = a['1'] * a['2'] ** (a['3'] * X ** a['4'] + a['5']) * \
+        (1 - a['3'] * a['4'] * _sp.log(a['2']) * X ** a['4']) + a['6']
+    values = {'S1': S1, 'S2': S2, 'rate': r}
+    return values
+
+
+def natural_exponential(target, A1='', A2='', A3='', A4='', A5='',
+                        quantity=''):
     r"""
     For the following source term:
         .. math::
             r =   A_{1} exp( A_{2}  x^{ A_{3} } + A_{4} )+ A_{5}
+
     If return_rate is True, it returns the value of source term for the
     provided x in each pore.
+
     If return_rate is False, it calculates the slope and intercept for the
     following linear form :
         .. math::
@@ -268,75 +196,58 @@ def natural_exponential(physics, phase, A1='', A2='', A3='', A4='', A5='',
     ----------
     A1 -> A5 : string
         The property name of the coefficients in the source term model
-    x : string or float/int or array/list
-        The property name or numerical value or array for the main quantity
-    Notes
-    -----
+
+    quantity : string
+        The property name for the main quantity being produced or consumed.
+
+    Returns
+    -------
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
+
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
+
+    Starting in V2, OpenPNM allows the assigning of such ``dicts`` to a Core
+    object ``dict``, but it breaks the received ``dict`` into separate arrays
+    and prepends the assigned propname.  For instance:
+
+        obj['pore.source_term']  = d
+
+    is equivalent to:
+
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    if x is '':
-        X = _sp.ones(physics.Np)*_sp.nan
-    else:
-        if type(x) == str:
-            x = 'pore.'+x.split('.')[-1]
-            try:
-                X = physics[x]
-            except KeyError:
-                raise Exception(physics.name +
-                                ' does not have the pore property :' + x + '!')
-        else:
-            X = _sp.array(x)
-
-    length_X = _sp.size(X)
-    if length_X != physics.Np:
-        if length_X == 1:
-            X = X * _sp.ones(physics.Np)
-        elif length_X >= phase.Np:
-            X = X[physics.map_pores()]
-        else:
-            raise Exception('Wrong size for the numerical array of x!')
-
-    a = {}
-    source_params = [A1, A2, A3, A4, A5]
-    for ind in _sp.arange(_sp.size(source_params)):
-        A = source_params[ind]
-        if A is '':
-            if ind == 0:
-                a[str(ind+1)] = 1
-            else:
-                a[str(ind+1)] = 0
+    X = target[quantity]
+    a = {'1': A1, '2': A2, '3': A3, '4': A4, '5': A5}
+    for item in a.keys():
+        if a[item] == '':
+            a[item] = 0
         else:
-            if type(A) == str:
-                A = 'pore.' + A.split('.')[-1]
-                try:
-                    a[str(ind+1)] = physics[A]
-                except KeyError:
-                    raise Exception(physics.name + '/' + phase.name +
-                                    ' does not have the pore property :' +
-                                    A + '!')
-            else:
-                raise Exception('source_term parameters can only be string '
-                                'type!')
-
-    if return_rate:
-        return(a['1'] * _sp.exp(a['2'] * X ** a['3'] + a['4']) + a['5'])
-    else:
-        S1 = a['1'] * a['2'] * \
-            a['3'] * X ** (a['3'] - 1) * \
-            _sp.exp(a['2'] * X ** a['3'] + a['4'])
-        S2 = a['1'] * (1 - a['2'] * a['3'] * X ** a['3']) * \
-            _sp.exp(a['2'] * X ** a['3'] + a['4']) + a['5']
-        return(_sp.vstack((S1, S2)).T)
-
-
-def logarithm(physics, phase, A1='', A2='', A3='', A4='', A5='', A6='',
-              x='', return_rate=True, **kwargs):
+            a[item] = target[a[item]]
+
+    r = a['1'] * _sp.exp(a['2'] * X ** a['3'] + a['4']) + a['5']
+
+    S1 = a['1'] * a['2'] * a['3'] * X ** (a['3'] - 1) * \
+        _sp.exp(a['2'] * X ** a['3'] + a['4'])
+
+    S2 = a['1'] * (1 - a['2'] * a['3'] * X ** a['3']) * \
+        _sp.exp(a['2'] * X ** a['3'] + a['4']) + a['5']
+
+    values = {'S1': S1, 'S2': S2, 'rate': r}
+    return values
+
+
+def logarithm(target, A1='', A2='', A3='', A4='', A5='', A6='', quantity=''):
     r"""
     For the following source term:
         .. math::
             r =  A_{1}   Log_{ A_{2} }( A_{3} x^{ A_{4} }+ A_{5})+ A_{6}
+
     If return_rate is True, it returns the value of source term for the
     provided x in each pore.
+
     If return_rate is False, it calculates the slope and intercept for the
     following linear form :
         .. math::
@@ -346,145 +257,112 @@ def logarithm(physics, phase, A1='', A2='', A3='', A4='', A5='', A6='',
     ----------
     A1 -> A6 : string
         The property name of the coefficients in the source term model
-    x : string or float/int or array/list
-        The property name or numerical value or array for the main quantity
-    Notes
-    -----
+
+    quantity : string
+        The property name for the main quantity being produced or consumed.
+
+    Returns
+    -------
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
+
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
+
+    Starting in V2, OpenPNM allows the assigning of such ``dicts`` to a Core
+    object ``dict``, but it breaks the received ``dict`` into separate arrays
+    and prepends the assigned propname.  For instance:
+
+        obj['pore.source_term']  = d
+
+    is equivalent to:
+
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    if x is '':
-        X = _sp.ones(physics.Np)*_sp.nan
-    else:
-        if type(x) == str:
-            x = 'pore.' + x.split('.')[-1]
-            try:
-                X = physics[x]
-            except KeyError:
-                raise Exception(physics.name +
-                                ' does not have the pore property :' + x + '!')
-        else:
-            X = _sp.array(x)
-
-    length_X = _sp.size(X)
-    if length_X != physics.Np:
-        if length_X == 1:
-            X = X * _sp.ones(physics.Np)
-        elif length_X >= phase.Np:
-            X = X[physics.map_pores()]
-        else:
-            raise Exception('Wrong size for the numerical array of x!')
-
-    a = {}
-    source_params = [A1, A2, A3, A4, A5, A6]
-    for ind in _sp.arange(_sp.size(source_params)):
-        A = source_params[ind]
-        if A is '':
-            if ind == 0:
-                a[str(ind+1)] = 1
-            else:
-                a[str(ind+1)] = 0
+    X = target[quantity]
+    a = {'1': A1, '2': A2, '3': A3, '4': A4, '5': A5, '6': A6}
+    for item in a.keys():
+        if a[item] == '':
+            a[item] = 0
         else:
-            if type(A) == str:
-                A = 'pore.' + A.split('.')[-1]
-                try:
-                    a[str(ind+1)] = physics[A]
-                except KeyError:
-                    raise Exception(physics.name + '/' + phase.name +
-                                    ' does not have the pore property :' +
-                                    A + '!')
-            else:
-                raise Exception('source_term parameters can only be string '
-                                'type!')
-
-    if return_rate:
-        return(a['1'] * _sp.log(a['3'] * X ** a['4'] + a['5']) /
-               _sp.log(a['2']) + a['6'])
-    else:
-        S1 = a['1'] * a['3'] * a['4'] * \
-            X ** (a['4'] - 1) / \
-            (_sp.log(a['2']) * (a['3'] * X ** a['4'] + a['5']))
-        S2 = a['1'] * _sp.log(a['3'] * X ** a['4'] + a['5']) / \
-            _sp.log(a['2']) + a['6'] - a['1'] * a['3'] * \
-            a['4'] * X ** a['4'] / \
-            (_sp.log(a['2']) * (a['3'] * X ** a['4'] + a['5']))
-        return(_sp.vstack((S1, S2)).T)
-
-
-def natural_logarithm(physics, phase, A1='', A2='', A3='', A4='', A5='',
-                      x='', return_rate=True, **kwargs):
+            a[item] = target[a[item]]
+
+    r = a['1'] * _sp.log(a['3'] * X ** a['4'] + a['5']) / _sp.log(a['2']) + \
+        a['6']
+
+    S1 = a['1'] * a['3'] * a['4'] * X ** (a['4'] - 1) / \
+        (_sp.log(a['2']) * (a['3'] * X ** a['4'] + a['5']))
+
+    S2 = a['1'] * _sp.log(a['3'] * X ** a['4'] + a['5']) / _sp.log(a['2']) + \
+        a['6'] - a['1'] * a['3'] * a['4'] * X ** a['4'] / \
+        (_sp.log(a['2']) * (a['3'] * X ** a['4'] + a['5']))
+
+    values = {'S1': S1, 'S2': S2, 'rate': r}
+    return values
+
+
+def natural_logarithm(target, A1='', A2='', A3='', A4='', A5='', quantity=''):
     r"""
     For the following source term:
+
         .. math::
             r =   A_{1}  Ln( A_{2} x^{ A_{3} }+ A_{4})+ A_{5}
-    If return_rate is True, it returns the value of source term for the
-    provided x in each pore.
-    If return_rate is False, it calculates the slope and intercept for the
-    following linear form :
+
+    This model calculates the rate of generation or consumption based on the
+    current value of ``quantity``, as well as the slope and interecpt of the
+    above equation at th current value of ``quantity`` giving the following
+    linear form:
+
         .. math::
             r = S_{1}   x  +  S_{2}
 
+    which can be used in linear solvers to produce a better guess of
+    ``quantity``, which can
+
     Parameters
     ----------
     A1 -> A5 : string
         The property name of the coefficients in the source term model
-    x : string or float/int or array/list
-        The property name or numerical value or array for the main quantity
-    Notes
-    -----
+
+    quantity : string
+        The property name for the main quantity being produced or consumed.
+
+    Returns
+    -------
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
+
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
+
+    Starting in V2, OpenPNM allows the assigning of such ``dicts`` to a Core
+    object ``dict``, but it breaks the received ``dict`` into separate arrays
+    and prepends the assigned propname.  For instance:
+
+        obj['pore.source_term']  = d
+
+    is equivalent to:
+
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    if x is '':
-        X = _sp.ones(physics.Np)*_sp.nan
-    else:
-        if type(x) == str:
-            x = 'pore.' + x.split('.')[-1]
-            try:
-                X = physics[x]
-            except KeyError:
-                raise Exception(physics.name +
-                                ' does not have the pore property :' + x + '!')
+    X = target[quantity]
+    a = {'1': A1, '2': A2, '3': A3, '4': A4, '5': A5}
+    for item in a.keys():
+        if a[item] == '':
+            a[item] = 0
         else:
-            X = _sp.array(x)
+            a[item] = target[a[item]]
 
-    length_X = _sp.size(X)
+    r = a['1'] * _sp.log(a['2'] * X ** a['3'] + a['4']) + a['5']
 
-    if length_X != physics.Np:
-        if length_X == 1:
-            X = X * _sp.ones(physics.Np)
-        elif length_X >= phase.Np:
-            X = X[physics.map_pores()]
-        else:
-            raise Exception('Wrong size for the numerical array of x!')
-
-    a = {}
-    source_params = [A1, A2, A3, A4, A5]
-    for ind in _sp.arange(_sp.size(source_params)):
-        A = source_params[ind]
-        if A is '':
-            if ind == 0:
-                a[str(ind+1)] = 1
-            else:
-                a[str(ind+1)] = 0
-        else:
-            if type(A) == str:
-                A = 'pore.' + A.split('.')[-1]
-                try:
-                    a[str(ind+1)] = physics[A]
-                except KeyError:
-                    raise Exception(physics.name + '/' + phase.name +
-                                    ' does not have the pore property :' +
-                                    A + '!')
-            else:
-                raise Exception('source_term parameters can only be string '
-                                'type!')
-
-    if return_rate:
-        return(a['1'] * _sp.log(a['2'] * X ** a['3'] + a['4']) + a['5'])
-    else:
-        S1 = a['1'] * a['2'] * a['3'] * \
-            X ** (a['3'] - 1) / \
-            (a['2'] * X ** a['3'] + a['4'])
-        S2 = a['1'] * _sp.log(a['2'] * X ** a['3'] + a['4']) + \
-            a['5'] - a['1'] * a['2'] * a['3'] * \
-            X ** a['3'] / (a['2'] * X ** a['3'] + a['4'])
-        return(_sp.vstack((S1, S2)).T)
+    S1 = a['1'] * a['2'] * a['3'] * X ** (a['3'] - 1) / \
+        (a['2'] * X ** a['3'] + a['4'])
+
+    S2 = a['1'] * _sp.log(a['2'] * X ** a['3'] + a['4']) + \
+        a['5'] - a['1'] * a['2'] * a['3'] * X ** a['3'] / \
+        (a['2'] * X ** a['3'] + a['4'])
+
+    values = {'S1': S1, 'S2': S2, 'rate': r}
+    return values
diff --git a/openpnm/models/physics/thermal_conductance.py b/openpnm/models/physics/thermal_conductance.py
index 977dc4569..58c1b0b49 100644
--- a/openpnm/models/physics/thermal_conductance.py
+++ b/openpnm/models/physics/thermal_conductance.py
@@ -8,11 +8,9 @@ Submodule -- thermal_conductance
 import scipy as _sp
 
 
-def series_resistors(physics, phase, network,
-                     thermal_conductivity='pore.thermal_conductivity',
+def series_resistors(target, thermal_conductivity='pore.thermal_conductivity',
                      pore_diameter='pore.diameter', pore_area='pore.area',
-                     throat_area='throat.area', throat_length='throat.length',
-                     **kwargs):
+                     throat_area='throat.area', throat_length='throat.length'):
     r"""
     Calculate the thermal conductance of void conduits in network ( 1/2 pore - full
     throat - 1/2 pore ) based on size (assuming cylindrical geometry)
@@ -33,11 +31,12 @@ def series_resistors(physics, phase, network,
     network then extracts the values for the appropriate throats at the end.
 
     """
+    network = target.project.network
+    phase = target.project.find_phase(target)
     # Get Nt-by-2 list of pores connected to each throat
     Ps = network['throat.conns']
     # Get properties in every pore in the network
-    kp = phase[thermal_conductivity]
-    kt = phase.interpolate_data(kp)
+    kt = phase.interpolate_data(thermal_conductivity)
     # Find g for half of pore 1
     pdia = network[pore_diameter]
     parea = network[pore_area]
@@ -58,5 +57,5 @@ def series_resistors(physics, phase, network,
     tlen[tlen <= 0] = 1e-12
     gt = kt*tarea/tlen
     value = (1/gt + 1/gp1 + 1/gp2)**(-1)
-    value = value[phase.throats(physics.name)]
+    value = value[phase.throats(target.name)]
     return value
diff --git a/openpnm/network/GenericNetwork.py b/openpnm/network/GenericNetwork.py
index 3f21e2259..1bb4c4591 100644
--- a/openpnm/network/GenericNetwork.py
+++ b/openpnm/network/GenericNetwork.py
@@ -48,7 +48,7 @@ class GenericNetwork(Base, ModelsMixin):
             self._gen_ids()
         # Now get values if present, or regenerate them
         vals = self.get(key)
-        if vals is None:
+        if vals is None:  # Invoke interleave data
             logger.debug(key + ' not on Network, check on Geometries')
             geoms = self.project.geometries().values()
             vals = self._interleave_data(key, geoms)
@@ -56,6 +56,7 @@ class GenericNetwork(Base, ModelsMixin):
 
     def _gen_ids(self):
         if ('pore._id' not in self.keys()):
+            logger.info('Generating pore IDs, please wait')
             self['pore._id'] = [str(uuid.uuid4()) for i in self.Ps]
         else:
             IDs = super().__getitem__('pore._id')
@@ -66,6 +67,7 @@ class GenericNetwork(Base, ModelsMixin):
                 IDs[inds] = temp
                 self['pore._id'] = IDs
         if ('throat._id' not in self.keys()):
+            logger.info('Generating throat IDs, please wait')
             self['throat._id'] = [str(uuid.uuid4()) for i in self.Ts]
         else:
             IDs = super().__getitem__('throat._id')
diff --git a/openpnm/topotools/__init__.py b/openpnm/topotools/__init__.py
index ed5259bf6..69bbfd4e1 100644
--- a/openpnm/topotools/__init__.py
+++ b/openpnm/topotools/__init__.py
@@ -9,6 +9,7 @@ from .topotools import find_connected_sites
 from .topotools import find_connecting_bonds
 from .topotools import find_clusters
 from .topotools import generate_base_points
+from .topotools import iscoplanar
 from .topotools import label_faces
 from .topotools import merge_pores
 from .topotools import plot_connections


Test patch:
diff --git a/tests/integration/AlgorithmTest.py b/tests/integration/AlgorithmTest.py
new file mode 100644
index 000000000..8e684bba1
--- /dev/null
+++ b/tests/integration/AlgorithmTest.py
@@ -0,0 +1,120 @@
+import openpnm as op
+import numpy as np
+
+
+def test_thermal_conduction():
+    # Generate Network and clean up boundaries (delete z-face pores)
+    divs = [10, 50]
+    Lc = 1/divs[0]
+    pn = op.network.Cubic(shape=divs, spacing=Lc)
+    pn.add_boundary_pores()
+    op.topotools.trim(network=pn,
+                      pores=pn.pores(['top_boundary', 'bottom_boundary']))
+    # Generate Geometry objects for internal and boundary pores
+    Ps = pn.pores('internal')
+    Ts = pn.throats()
+    geom = op.geometry.GenericGeometry(network=pn, pores=Ps, throats=Ts)
+    geom['pore.area'] = Lc**2
+    geom['pore.diameter'] = Lc
+    geom['throat.length'] = 1e-25
+    geom['throat.area'] = Lc**2
+    Ps = pn.pores('*boundary')
+    boun = op.geometry.GenericGeometry(network=pn, pores=Ps)
+    boun['pore.area'] = Lc**2
+    boun['pore.diameter'] = 1e-25
+    # Create Phase object and associate with a Physics object
+    Cu = op.phases.GenericPhase(network=pn)
+    Cu['pore.thermal_conductivity'] = 1.0  # W/m.K
+    phys1 = op.physics.GenericPhysics(network=pn, phase=Cu, geometry=geom)
+    phys2 = op.physics.GenericPhysics(network=pn, phase=Cu, geometry=boun)
+    mod = op.models.physics.thermal_conductance.series_resistors
+    phys1.add_model(propname='throat.thermal_conductance', model=mod)
+    phys2.models = phys1.models.copy()
+    phys1.regenerate_models()  # Update the conductance values
+    phys2.regenerate_models()  # Update the conductance values
+    # Setup Algorithm object
+    alg = op.algorithms.FourierConduction(network=pn, phase=Cu)
+    inlets = pn.pores('back_boundary')
+    outlets = pn.pores(['front_boundary', 'left_boundary', 'right_boundary'])
+    T_in = 30*np.sin(np.pi*pn['pore.coords'][inlets, 1]/5)+50
+    alg.set_dirichlet_BC(values=T_in, pores=inlets)
+    alg.set_dirichlet_BC(values=50, pores=outlets)
+    alg.run()
+    Cu.update(alg.results())
+    # Calculate analytical solution over the same domain spacing
+    T = 30*np.sinh(np.pi*pn['pore.coords'][:, 0]/5)/np.sinh(np.pi/5) * \
+        np.sin(np.pi*pn['pore.coords'][:, 1]/5) + 50
+    Cu['pore.analytical_temp'] = T
+    b = Cu['pore.analytical_temp'][pn.pores(geom.name)]
+    a = Cu['pore.temperature'][pn.pores(geom.name)]
+    a = np.reshape(a, (divs[0], divs[1]))
+    b = np.reshape(b, (divs[0], divs[1]))
+    diff = a - b
+    assert np.amax(np.absolute(diff)) < 0.015
+
+
+def test_open_air_diffusivity():
+    pass
+#    pn = op.network.Cubic([5, 5, 5], spacing=1)
+#    pn.add_boundary_pores()
+#    air = op.phases.Air(network=pn)
+#    Dab = np.mean(air['pore.diffusivity'])
+#    c = np.mean(air['pore.molar_density'])
+#    air['throat.diffusive_conductance'] = Dab * c
+#    BC1_pores = pn.pores(labels=['top_boundary'])
+#    BC2_pores = pn.pores(labels=['bottom_boundary'])
+#    Diff = op.algorithms.FickianDiffusion(network=pn, phase=air)
+#    # Assign Dirichlet boundary conditions to top and bottom surface pores
+#    Diff.set_dirichlet_BC(values=0.6, pores=BC1_pores)
+#    Diff.set_dirichlet_BC(values=0.4, pores=BC2_pores)
+#    Diff.run()
+#    Diff.domain_area = 25
+#    Diff.domain_length = 5
+#    Diff_deff = Diff.calc_eff_diffusivity()/Dab
+#    assert np.round(Diff_deff, 3) == 1
+
+
+def test_Darcy_alg():
+    pass
+#    # Generate Network and clean up some of boundaries
+#    divs = [1, 50, 10]
+#    Lc = 0.00004
+#    pn = op.network.Cubic(shape=divs, spacing=Lc)
+#    # Generate Geometry objects for internal and boundary pores
+#    geom = op.geometry.StickAndBall(network=pn, pores=pn.Ps, throats=pn.Ts)
+#    # Create Phase object and associate with a Physics object
+#    air = op.phases.Air(network=pn)
+#    phys = op.physics.GenericPhysics(network=pn, phase=air, geometry=geom)
+#    mod = op.models.physics.hydraulic_conductance.hagen_poiseuille
+#    phys.add_model(propname='throat.hydraulic_conductance',
+#                   model=mod,
+#                   calc_pore_len=False)
+#    phys.regenerate_models()  # Update the conductance values
+#    # Setup Algorithm objects
+#    alg1 = op.algorithms.StokesFlow(network=pn, phase=air)
+#    inlets = pn.pores('bottom')
+#    outlets = pn.pores('top')
+#    P_out = 0  # Pa
+#    Q_in = 0.6667*(Lc**2)*divs[1]*divs[0]  # m^3/s
+#    alg1.set_neumann_BC(values=-Q_in, pores=inlets)
+#    alg1.set_dirichlet_BC(values=P_out, pores=outlets)
+#    alg1.run()
+#    air.update(alg1.results())
+#    a = round(np.absolute(alg1.rate(outlets))[0], 16)
+#    b = round(np.absolute(np.sum(alg1['pore.neumann_value'][inlets])), 16)
+#    assert a == b
+#
+#    alg2 = op.algorithms.StokesFlow(network=pn, phase=air)
+#    inlets = pn.pores('bottom')
+#    outlets = pn.pores('top')
+#    P_out = 0  # Pa
+#    P_in = 1000  # Pa
+#    alg2.set_dirichlet_BC(values=P_in, pores=inlets)
+#    alg2.set_dirichlet_BC(values=P_out, pores=outlets)
+#    alg2.run()
+#    a = round(np.absolute(alg2.rate(inlets))[0], 16)
+#    b = round(np.absolute(alg2.rate(outlets))[0], 16)
+#    assert a == b
+#    Q = -alg2.rate(inlets)
+#    K = Q*air['pore.viscosity'][0]*divs[2]*Lc/(divs[0]*divs[1]*Lc**2*(P_in-P_out))
+#    K_alg = alg2.calc_eff_permeability()
diff --git a/tests/unit/algorithms/DrainageTest.py b/tests/unit/algorithms/DrainageTest.py
new file mode 100644
index 000000000..9eb86e30b
--- /dev/null
+++ b/tests/unit/algorithms/DrainageTest.py
@@ -0,0 +1,214 @@
+import openpnm as op
+import scipy as sp
+import pytest
+mgr = op.Workspace()
+
+
+class DrainageTest:
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[5, 5, 5], spacing=0.0005)
+        self.geo = op.geometry.StickAndBall(network=self.net,
+                                            pores=self.net.Ps,
+                                            throats=self.net.Ts)
+        self.water = op.phases.Water(network=self.net)
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.water,
+                                              geometry=self.geo)
+        mod = op.models.physics.capillary_pressure.washburn
+        self.phys.add_model(propname='throat.capillary_pressure',
+                            model=mod)
+        self.alg = op.algorithms.Drainage(network=self.net)
+
+    def test_set_inlets_modes(self):
+        self.alg.setup(invading_phase=self.water, defending_phase=self.air)
+
+        self.alg.set_inlets(pores=self.net.pores('top'), mode='add')
+        assert sp.sum(self.alg['pore.inlets']) == 25
+
+        self.alg.set_inlets(pores=self.net.pores('bottom'), mode='add')
+        assert sp.sum(self.alg['pore.inlets']) == 50
+
+        self.alg.set_inlets(pores=self.net.pores('top'), mode='overwrite')
+        assert sp.sum(self.alg['pore.inlets']) == 25
+
+        self.alg.set_inlets(pores=self.net.pores('top'), mode='remove')
+        assert sp.sum(self.alg['pore.inlets']) == 0
+
+        self.alg.set_inlets(pores=self.net.pores('top'), mode='add')
+        self.alg.set_inlets(mode='clear')
+        assert sp.sum(self.alg['pore.inlets']) == 0
+
+    def test_set_inlets_conflicting_with_outlets(self):
+        self.alg.setup(invading_phase=self.water, defending_phase=self.air)
+        self.alg['pore.outlets'][self.net.pores('top')] = True
+        with pytest.raises(Exception):
+            self.alg.set_inlets(pores=self.net.pores('top'), mode='add')
+
+    def test_set_outlets_conflicting_with_inlets(self):
+        self.alg.setup(invading_phase=self.water,
+                       defending_phase=self.air,
+                       trapping=True)
+        self.alg['pore.inlets'][self.net.pores('top')] = True
+        try:
+            self.alg.set_outlets(pores=self.net.pores('top'))
+        except:
+            flag = True
+        assert flag
+
+    def test_set_outlets_without_trapping(self):
+        self.alg.setup(invading_phase=self.water,
+                       defending_phase=self.air)
+        with pytest.raises(Exception):
+            self.alg.set_outlets(pores=self.net.pores('top'))
+
+    def test_set_outlets_modes(self):
+        self.alg.setup(invading_phase=self.water,
+                       defending_phase=self.air,
+                       trapping=True)
+
+        self.alg.set_outlets(pores=self.net.pores('top'), mode='add')
+        assert sp.sum(self.alg['pore.outlets']) == 25
+
+        self.alg.set_outlets(pores=self.net.pores('bottom'), mode='add')
+        assert sp.sum(self.alg['pore.outlets']) == 50
+
+        self.alg.set_outlets(pores=self.net.pores('top'), mode='overwrite')
+        assert sp.sum(self.alg['pore.outlets']) == 25
+
+        self.alg.set_outlets(pores=self.net.pores('top'), mode='remove')
+        assert sp.sum(self.alg['pore.outlets']) == 0
+
+        self.alg.set_outlets(pores=self.net.pores('top'), mode='add')
+        self.alg.set_outlets(mode='clear')
+        assert sp.sum(self.alg['pore.outlets']) == 0
+
+    def test_set_residual_modes(self):
+        self.alg.setup(invading_phase=self.water, defending_phase=self.air)
+
+        Ps = sp.random.randint(0, self.net.Np, 10)
+        Ts = self.net.find_neighbor_pores(pores=Ps)
+        self.alg.set_residual(pores=Ps, throats=Ts, mode='add')
+        assert sp.sum(self.alg['pore.residual']) == sp.size(sp.unique(Ps))
+        assert sp.sum(self.alg['throat.residual']) == sp.size(sp.unique(Ts))
+
+        Ps = sp.random.randint(0, self.net.Np, 10)
+        Ts = self.net.find_neighbor_pores(pores=Ps)
+        self.alg.set_residual(pores=Ps, throats=Ts, mode='add')
+        assert sp.sum(self.alg['pore.residual']) > sp.size(sp.unique(Ps))
+        assert sp.sum(self.alg['throat.residual']) > sp.size(sp.unique(Ts))
+
+        Ps = sp.random.randint(0, self.net.Np, 10)
+        Ts = self.net.find_neighbor_pores(pores=Ps)
+        self.alg.set_residual(pores=Ps, throats=Ts, mode='overwrite')
+        assert sp.sum(self.alg['pore.residual']) == sp.size(sp.unique(Ps))
+        assert sp.sum(self.alg['throat.residual']) == sp.size(sp.unique(Ts))
+
+        self.alg.set_residual(pores=Ps, throats=Ts, mode='remove')
+        assert sp.sum(self.alg['pore.residual']) == 0
+
+        self.alg.set_residual(pores=Ps, throats=Ts, mode='add')
+        self.alg.set_residual(mode='clear')
+        assert sp.sum(self.alg['pore.residual']) == 0
+
+    def test_set_boundary_conditions_clear(self):
+        self.alg.setup(invading_phase=self.water,
+                       defending_phase=self.air,
+                       trapping=True)
+
+        self.alg['pore.inlets'] = True
+        self.alg['pore.outlets'] = True
+        self.alg['pore.residual'] = True
+        self.alg['throat.residual'] = True
+
+        self.alg.set_boundary_conditions(mode='clear')
+        assert sp.sum(self.alg['pore.inlets']) == 0
+        assert sp.sum(self.alg['pore.outlets']) == 0
+        assert sp.sum(self.alg['pore.residual']) == 0
+        assert sp.sum(self.alg['throat.residual']) == 0
+
+    def test_set_boundary_conditions_bctypes(self):
+        self.alg.setup(invading_phase=self.water,
+                       defending_phase=self.air,
+                       trapping=True)
+        Ps = sp.random.randint(0, self.net.Np, 10)
+
+        self.alg.set_boundary_conditions(pores=Ps, bc_type='inlets')
+        assert sp.sum(self.alg['pore.inlets']) == sp.size(sp.unique(Ps))
+        self.alg['pore.inlets'] = False
+
+        self.alg.set_boundary_conditions(pores=Ps, bc_type='outlets')
+        assert sp.sum(self.alg['pore.outlets']) == sp.size(sp.unique(Ps))
+        self.alg['pore.outlets'] = False
+
+        self.alg.set_boundary_conditions(pores=Ps, bc_type='residual')
+        assert sp.sum(self.alg['pore.residual']) == sp.size(sp.unique(Ps))
+        self.alg['pore.residual'] = False
+
+        with pytest.raises(Exception):
+            self.alg.set_boundary_conditions(pores=Ps, bc_type='bad_type')
+
+        with pytest.raises(Exception):
+            self.alg.set_boundary_conditions(bc_type=None, mode='bad_type')
+
+    def test_run_npts(self):
+        self.alg.setup(invading_phase=self.water, defending_phase=self.air)
+        Ps = sp.random.randint(0, self.net.Np, 10)
+        self.alg.set_boundary_conditions(pores=Ps, bc_type='inlets')
+        self.alg.run(npts=20)
+
+    def test_run_inv_pressures(self):
+        self.alg.setup(invading_phase=self.water, defending_phase=self.air)
+        Ps = sp.random.randint(0, self.net.Np, 10)
+        self.alg.set_boundary_conditions(pores=Ps, bc_type='inlets')
+        self.alg.run(inv_pressures=range(0, 20000, 1000))
+        assert sp.all(self.alg._inv_points == range(0, 20000, 1000))
+
+    def test_run_no_inlets(self):
+        self.alg.setup(invading_phase=self.water, defending_phase=self.air)
+        with pytest.raises(Exception):
+            alg.run()
+
+    def test_run_w_trapping_but_no_outlets(self):
+        self.alg.setup(invading_phase=self.water,
+                       defending_phase=self.air,
+                       trapping=True)
+        Ps = sp.random.randint(0, self.net.Np, 10)
+        self.alg.set_boundary_conditions(pores=Ps, bc_type='inlets')
+        with pytest.raises(Exception):
+            alg.run()
+
+    def test_run_w_trapping(self):
+        self.alg.setup(invading_phase=self.water,
+                       defending_phase=self.air,
+                       trapping=True)
+        self.alg.set_boundary_conditions(pores=self.net.pores('top'),
+                                         bc_type='inlets')
+        self.alg.set_boundary_conditions(pores=self.net.pores('bottom'),
+                                         bc_type='outlets')
+#        self.alg.run()
+#        data = self.alg.get_drainage_data()
+#        assert 'capillary_pressure' in data.keys()
+#        assert 'invading_phase_saturation' in data.keys()
+
+    def test_run_w_residual_pores_and_throats(self):
+        self.alg.setup(invading_phase=self.water, defending_phase=self.air)
+        self.alg.set_boundary_conditions(pores=self.net.pores('top'),
+                                         bc_type='inlets')
+        self.alg.set_boundary_conditions(pores=self.net.pores('bottom'),
+                                         bc_type='residual')
+        self.alg.run()
+        data = self.alg.get_drainage_data()
+        assert 'capillary_pressure' in data.keys()
+        assert 'invading_phase_saturation' in data.keys()
+
+
+if __name__ == '__main__':
+
+    t = DrainageTest()
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
+    self = t
diff --git a/tests/unit/algorithms/GenericTransportTest.py b/tests/unit/algorithms/GenericTransportTest.py
new file mode 100644
index 000000000..50605ccb8
--- /dev/null
+++ b/tests/unit/algorithms/GenericTransportTest.py
@@ -0,0 +1,83 @@
+import openpnm as op
+import scipy as sp
+import pytest
+
+
+class GenericTransportTest:
+
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[9, 9, 9])
+        self.geo = op.geometry.StickAndBall(network=self.net,
+                                            pores=self.net.Ps,
+                                            throats=self.net.Ts)
+        self.phase = op.phases.Air(network=self.net)
+        self.phase['pore.mole_fraction'] = 0
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.phase,
+                                              geometry=self.geo)
+        self.phys['throat.diffusive_conductance'] = 1.0
+
+    def test_results(self):
+        alg = op.algorithms.GenericTransport(network=self.net,
+                                             phase=self.phase)
+        with pytest.raises(Exception):
+            alg.results()
+
+    def test_remove_boundary_conditions(self):
+        alg = op.algorithms.GenericTransport(network=self.net,
+                                             phase=self.phase)
+        alg.set_dirichlet_BC(pores=self.net.pores('top'), values=1)
+        alg.set_dirichlet_BC(pores=self.net.pores('bottom'), values=0)
+        assert sp.sum(alg['pore.dirichlet']) > 0
+        alg.remove_BC(pores=self.net.pores('top'))
+        assert sp.sum(alg['pore.dirichlet']) > 0
+        alg.remove_BC(pores=self.net.pores('bottom'))
+        assert sp.sum(alg['pore.dirichlet']) == 0
+
+    def test_generic_transport(self):
+        alg = op.algorithms.GenericTransport(network=self.net,
+                                             phase=self.phase)
+        alg.settings['conductance'] = 'throat.diffusive_conductance'
+        alg.settings['quantity'] = 'pore.mole_fraction'
+        alg.set_dirichlet_BC(pores=self.net.pores('top'), values=1)
+        alg.set_dirichlet_BC(pores=self.net.pores('bottom'), values=0)
+        alg.run()
+
+    def test_two_dirichlet_conditions(self):
+        alg = op.algorithms.GenericTransport(network=self.net,
+                                             phase=self.phase)
+        alg.settings['conductance'] = 'throat.diffusive_conductance'
+        alg.settings['quantity'] = 'pore.mole_fraction'
+        alg.set_dirichlet_BC(pores=self.net.pores('top'), values=1)
+        alg.set_dirichlet_BC(pores=self.net.pores('bottom'), values=0)
+        alg.run()
+        x = [0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0]
+        y = sp.unique(sp.around(alg['pore.mole_fraction'], decimals=3))
+        assert sp.all(x == y)
+
+    def test_one_dirichlet_one_neumann(self):
+        alg = op.algorithms.GenericTransport(network=self.net,
+                                             phase=self.phase)
+        alg.settings['conductance'] = 'throat.diffusive_conductance'
+        alg.settings['quantity'] = 'pore.mole_fraction'
+        alg.set_neumann_BC(pores=self.net.pores('bottom'), values=1)
+        alg.set_dirichlet_BC(pores=self.net.pores('top'), values=0)
+        alg.run()
+        x = [0., 1., 2., 3., 4., 5., 6., 7., 8.]
+        y = sp.unique(sp.around(alg['pore.mole_fraction'], decimals=3))
+        assert sp.all(x == y)
+
+    def teardown_class(self):
+        ws = op.core.Workspace()
+        ws.clear()
+
+
+if __name__ == '__main__':
+
+    t = GenericTransportTest()
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
+    self = t
diff --git a/tests/unit/algorithms/ReactiveTransportTest.py b/tests/unit/algorithms/ReactiveTransportTest.py
new file mode 100644
index 000000000..69373d0b3
--- /dev/null
+++ b/tests/unit/algorithms/ReactiveTransportTest.py
@@ -0,0 +1,54 @@
+import openpnm as op
+import scipy as sp
+import pytest
+
+
+class GenericTransportTest:
+
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[9, 9, 9])
+        self.geo = op.geometry.StickAndBall(network=self.net,
+                                            pores=self.net.Ps,
+                                            throats=self.net.Ts)
+        self.phase = op.phases.Air(network=self.net)
+        self.phase['pore.mole_fraction'] = 0
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.phase,
+                                              geometry=self.geo)
+        self.phys['throat.diffusive_conductance'] = 1.0
+        self.phys['pore.A'] = 1e6
+        self.phys['pore.k'] = 2
+        mod = op.models.physics.generic_source_term.standard_kinetics
+        self.phys.add_model(propname='pore.reaction',
+                            model=mod,
+                            prefactor='pore.A',
+                            exponent='pore.k',
+                            quantity='pore.mole_fraction',
+                            regen_mode='normal')
+
+    def test_one_dirichlet_one_source(self):
+        alg = op.algorithms.ReactiveTransport(network=self.net,
+                                              phase=self.phase)
+        alg.settings['conductance'] = 'throat.diffusive_conductance'
+        alg.settings['quantity'] = 'pore.mole_fraction'
+        alg.set_source(pores=self.net.pores('bottom'), propname='pore.reaction')
+        alg.set_dirichlet_BC(pores=self.net.pores('top'), values=1)
+        alg.run()
+        x = [0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0]
+        y = sp.unique(sp.around(alg['pore.mole_fraction'], decimals=3))
+        assert sp.all(x == y)
+
+    def teardown_class(self):
+        ws = op.core.Workspace()
+        ws.clear()
+
+
+if __name__ == '__main__':
+
+    t = GenericTransportTest()
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
+    self = t
diff --git a/tests/unit/algorithms/SubclassedTransportTest.py b/tests/unit/algorithms/SubclassedTransportTest.py
new file mode 100644
index 000000000..d05d6a987
--- /dev/null
+++ b/tests/unit/algorithms/SubclassedTransportTest.py
@@ -0,0 +1,77 @@
+import openpnm as op
+import scipy as sp
+import pytest
+
+
+class SubclassedTransportTest:
+
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[9, 9, 9])
+        self.geo = op.geometry.StickAndBall(network=self.net,
+                                            pores=self.net.Ps,
+                                            throats=self.net.Ts)
+        self.phase = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.phase,
+                                              geometry=self.geo)
+
+    def test_fickian_diffusion(self):
+        alg = op.algorithms.FickianDiffusion(network=self.net,
+                                             phase=self.phase)
+        self.phys['throat.diffusive_conductance'] = 1
+        alg.set_dirichlet_BC(pores=self.net.pores('top'), values=1)
+        alg.set_dirichlet_BC(pores=self.net.pores('bottom'), values=0)
+        alg.run()
+        alg.domain_area = 81
+        alg.domain_length = 9
+        Deff = alg.calc_eff_diffusivity()
+        assert sp.around(Deff, decimals=10) == 0.0275097564
+
+    def test_stokes_flow(self):
+        alg = op.algorithms.StokesFlow(network=self.net, phase=self.phase)
+        self.phys['throat.hydraulic_conductance'] = 1
+        alg.set_dirichlet_BC(pores=self.net.pores('top'), values=101325)
+        alg.set_dirichlet_BC(pores=self.net.pores('bottom'), values=0)
+        alg.run()
+        alg.domain_area = 81
+        alg.domain_length = 9
+        Keff = alg.calc_eff_permeability()
+        assert sp.around(Keff, decimals=10) == 2.075e-05
+
+    def test_forurier_conduction(self):
+        alg = op.algorithms.FourierConduction(network=self.net,
+                                              phase=self.phase)
+        self.phys['throat.thermal_conductance'] = 1
+        alg.set_dirichlet_BC(pores=self.net.pores('top'), values=101325)
+        alg.set_dirichlet_BC(pores=self.net.pores('bottom'), values=0)
+        alg.run()
+        alg.domain_area = 81
+        alg.domain_length = 9
+        Keff = alg.calc_effective_conductivity()
+        assert sp.around(Keff, decimals=10) == 1.125
+
+    def test_ohmic_conduction(self):
+        alg = op.algorithms.OhmicConduction(network=self.net, phase=self.phase)
+        self.phys['throat.electrical_conductance'] = 1
+        alg.set_dirichlet_BC(pores=self.net.pores('top'), values=101325)
+        alg.set_dirichlet_BC(pores=self.net.pores('bottom'), values=0)
+        alg.run()
+        alg.domain_area = 81
+        alg.domain_length = 9
+        Keff = alg.calc_effective_conductivity()
+        assert sp.around(Keff, decimals=10) == 1.125
+
+    def teardown_class(self):
+        ws = op.core.Workspace()
+        ws.clear()
+
+
+if __name__ == '__main__':
+
+    t = SubclassedTransportTest()
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
+    self = t
diff --git a/tests/unit/core/BaseTest.py b/tests/unit/core/BaseTest.py
index 61fb2797c..c3e55413b 100644
--- a/tests/unit/core/BaseTest.py
+++ b/tests/unit/core/BaseTest.py
@@ -563,9 +563,9 @@ class BaseTest:
 
     def test_write_dict(self):
         self.net['pore.test_dict'] = {'test1': 1, 'test2': self.net.Ps}
-        assert 'pore.test_dict_test1' in self.net.keys()
-        assert 'pore.test_dict_test2' in self.net.keys()
-        self.net['pore.test_dict'] = {'test1': 1, 'test2': self.net.Ps}
+        assert 'pore.test_dict.test1' in self.net.keys()
+        assert self.net['pore.test_dict.test1'].shape == (self.net.Np, )
+        assert 'pore.test_dict.test2' in self.net.keys()
 
     def test_pore_mapping(self):
         a = self.geo21['pore._id']
diff --git a/tests/unit/core/ExtrasTest.py b/tests/unit/core/ExtrasTest.py
new file mode 100644
index 000000000..53c5dd22a
--- /dev/null
+++ b/tests/unit/core/ExtrasTest.py
@@ -0,0 +1,66 @@
+import openpnm as op
+import scipy as sp
+import pytest
+
+
+class ExtrasTest:
+
+    def setup_class(self):
+        ws = op.Workspace()
+
+    def teardown_class(self):
+        ws = op.core.Workspace()
+        ws.clear()
+
+    def test_initialize_GenericNetwork_without_args(self):
+        net = op.network.GenericNetwork()
+        assert set(net.keys()) == set(['pore.all', 'throat.all'])
+        assert net.Np == 0
+        assert net.Nt == 0
+
+    def test_initialize_GenericGeometry_without_args(self):
+        obj = op.geometry.GenericGeometry()
+        assert set(obj.keys()) == set(['pore.all', 'throat.all'])
+        assert obj.Np == 0
+        assert obj.Nt == 0
+        assert len(obj.project) == 1
+
+    def test_initialize_StickAndBall_without_args(self):
+        obj = op.geometry.StickAndBall()
+        assert set(obj.keys()) == set(['pore.all', 'throat.all'])
+        assert len(obj.models.keys()) > 0
+        assert obj.Np == 0
+        assert obj.Nt == 0
+        assert len(obj.project) == 1
+
+    def test_initialize_GenericPhase_without_args(self):
+        obj = op.phases.GenericPhase()
+        assert obj.Np == 0
+        assert obj.Nt == 0
+        assert len(obj.project) == 1
+
+    def test_initialize_Air_without_args(self):
+        obj = op.phases.Air()
+        assert len(obj.keys()) > 4
+        assert len(obj.models.keys()) > 0
+        assert obj.Np == 0
+        assert obj.Nt == 0
+        assert len(obj.project) == 1
+
+    def test_initialize_GenericPhysics_without_args(self):
+        obj = op.physics.GenericPhysics()
+        assert set(obj.keys()) == set(['pore.all', 'throat.all'])
+        assert obj.Np == 0
+        assert obj.Nt == 0
+        assert len(obj.project) == 1
+
+
+if __name__ == '__main__':
+
+    t = ExtrasTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/physics/GenericSourceTerm.py b/tests/unit/models/physics/GenericSourceTerm.py
index 782258e6d..6a9cee6d7 100644
--- a/tests/unit/models/physics/GenericSourceTerm.py
+++ b/tests/unit/models/physics/GenericSourceTerm.py
@@ -15,260 +15,198 @@ class GenericSourceTermTest:
                                               phase=self.phase,
                                               geometry=self.geo)
         self.phys['throat.diffusive_conductance'] = 5e-8
-        self.phase['pore.mole_fraction'] = 0.
-        self.alg = op.algorithms.GenericTransport(network=self.net,
-                                                  phase=self.phase)
-        BC_pores = np.arange(20, 30)
-        self.S_pores = np.arange(55, 85)
-        self.alg.set_boundary_conditions(bctype='Dirichlet',
-                                         bcvalue=0.4,
-                                         pores=BC_pores)
+        self.phase['pore.mole_fraction'] = 0.0
+        self.BC_pores = np.arange(20, 30)
+        self.source_pores = np.arange(55, 85)
 
-#    def test_linear(self):
-#        self.phys['pore.item1'] = 0.5e-11
-#        self.phys['pore.item2'] = 1.5e-12
-#        self.phys.add_model(propname='pore.source1',
-#                            model=pm.generic_source_term.linear,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            x='mole_fraction',
-#                            return_rate=False,
-#                            regen_mode='on_demand')
-#        self.phys.add_model(propname='pore.source2',
-#                            model=pm.generic_source_term.linear,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            x='mole_fraction',
-#                            return_rate=True,
-#                            regen_mode='on_demand')
-#        self.alg.set_source_term(source_name='pore.source1',
-#                                 pores=self.S_pores,
-#                                 mode='overwrite')
-#        self.alg.run(conductance='throat.diffusive_conductance',
-#                     quantity='pore.mole_fraction',
-#                     super_pore_conductance=None)
-#        self.alg.return_results()
-#        self.phys.regenerate(props='pore.source1')
-#        self.phys.regenerate(props='pore.source2')
-#        X = self.phase['pore.mole_fraction']
-#        r1 = np.round(np.sum(0.5e-11 * X[self.S_pores] + 1.5e-12), 20)
-#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+    def test_linear(self):
+        self.phys['pore.item1'] = 0.5e-11
+        self.phys['pore.item2'] = 1.5e-12
+        self.phys.add_model(propname='pore.source1',
+                            model=pm.generic_source_term.linear,
+                            A1='pore.item1',
+                            A2='pore.item2',
+                            quantity='pore.mole_fraction',
+                            regen_mode='normal')
+        self.alg = op.algorithms.ReactiveTransport(network=self.net,
+                                                   phase=self.phase)
+        self.alg.settings.update({'conductance': 'throat.diffusive_conductance',
+                                  'quantity': 'pore.mole_fraction'})
+        self.alg.set_dirichlet_BC(values=0.4, pores=self.BC_pores)
+        self.alg.set_source(propname='pore.source1',
+                            pores=self.source_pores)
+        self.alg.run()
+        self.phase.update(self.alg.results())
+        self.phys.regenerate_models(propnames='pore.source1')
+        X = self.phase['pore.mole_fraction']
+        r1 = np.round(np.sum(0.5e-11 * X[self.source_pores] + 1.5e-12), 20)
+        r2 = np.round(np.sum(self.phys['pore.source1.rate'][self.source_pores]), 20)
 #        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
-#        assert r1 == r2
+        assert r1 == r2
 #        assert r2 == -r3
-#
-#    def test_power_law(self):
-#        self.phys['pore.item1'] = 0.5e-12
-#        self.phys['pore.item2'] = 2.5
-#        self.phys['pore.item3'] = -1.4e-11
-#        self.phys.add_model(propname='pore.source1',
-#                            model=pm.generic_source_term.power_law,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            A3='pore.item3',
-#                            x='mole_fraction',
-#                            return_rate=False,
-#                            regen_mode='on_demand')
-#        self.phys.add_model(propname='pore.source2',
-#                            model=pm.generic_source_term.power_law,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            A3='pore.item3',
-#                            x='mole_fraction',
-#                            return_rate=True,
-#                            regen_mode='on_demand')
-#        self.alg.set_source_term(source_name='pore.source1',
-#                                 pores=self.S_pores,
-#                                 mode='overwrite')
-#        self.alg.run(conductance='throat.diffusive_conductance',
-#                     quantity='pore.mole_fraction',
-#                     super_pore_conductance=None)
-#        self.alg.return_results()
-#        self.phys.regenerate(props='pore.source1')
-#        self.phys.regenerate(props='pore.source2')
-#        X = self.phase['pore.mole_fraction']
-#        r1 = np.round(np.sum(0.5e-12 * X[self.S_pores] ** 2.5 - 1.4e-11), 20)
-#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+
+    def test_power_law(self):
+        self.phys['pore.item1'] = 0.5e-12
+        self.phys['pore.item2'] = 2.5
+        self.phys['pore.item3'] = -1.4e-11
+        self.phys.add_model(propname='pore.source1',
+                            model=pm.generic_source_term.power_law,
+                            A1='pore.item1',
+                            A2='pore.item2',
+                            A3='pore.item3',
+                            quantity='pore.mole_fraction',
+                            regen_mode='normal')
+        self.alg = op.algorithms.ReactiveTransport(network=self.net,
+                                                   phase=self.phase)
+        self.alg.set_dirichlet_BC(values=0.4, pores=self.BC_pores)
+        self.alg.set_source(propname='pore.source1',
+                            pores=self.source_pores)
+        self.alg.settings.update({'conductance': 'throat.diffusive_conductance',
+                                  'quantity': 'pore.mole_fraction'})
+        self.alg.run()
+        self.phase.update(self.alg.results())
+        self.phys.regenerate_models(propnames='pore.source1')
+        X = self.phase['pore.mole_fraction']
+        r1 = np.sum(0.5e-12 * X[self.source_pores]**2.5 - 1.4e-11)
+        r2 = np.sum(self.phys['pore.source1.rate'][self.source_pores])
 #        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
-#        assert r1 == r2
+        assert r1 == r2
 #        assert r2 == -r3
-#
-#    def test_exponential(self):
-#        self.phys['pore.item1'] = 0.8e-11
-#        self.phys['pore.item2'] = 3
-#        self.phys['pore.item3'] = 0.5
-#        self.phys['pore.item4'] = 2
-#        self.phys['pore.item5'] = -0.34
-#        self.phys['pore.item6'] = 2e-14
-#        self.phys.add_model(propname='pore.source1',
-#                            model=pm.generic_source_term.exponential,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            A3='pore.item3',
-#                            A4='pore.item4',
-#                            A5='pore.item5',
-#                            A6='pore.item6',
-#                            x='mole_fraction',
-#                            return_rate=False,
-#                            regen_mode='on_demand')
-#        self.phys.add_model(propname='pore.source2',
-#                            model=pm.generic_source_term.exponential,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            A3='pore.item3',
-#                            A4='pore.item4',
-#                            A5='pore.item5',
-#                            A6='pore.item6',
-#                            x='mole_fraction',
-#                            return_rate=True,
-#                            regen_mode='on_demand')
-#        self.alg.set_source_term(source_name='pore.source1',
-#                                 pores=self.S_pores,
-#                                 mode='overwrite')
-#        self.alg.run(conductance='throat.diffusive_conductance',
-#                     quantity='pore.mole_fraction',
-#                     super_pore_conductance=None)
-#        self.alg.return_results()
-#        self.phys.regenerate(props='pore.source1')
-#        self.phys.regenerate(props='pore.source2')
-#        X = self.phase['pore.mole_fraction']
-#        r1 = np.round(np.sum(0.8e-11 * 3 ** (0.5 * X[self.S_pores] ** 2 -
-#                      0.34) + 2e-14), 20)
-#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+
+    def test_exponential(self):
+        self.phys['pore.item1'] = 0.8e-11
+        self.phys['pore.item2'] = 3
+        self.phys['pore.item3'] = 0.5
+        self.phys['pore.item4'] = 2
+        self.phys['pore.item5'] = -0.34
+        self.phys['pore.item6'] = 2e-14
+        self.phys.add_model(propname='pore.source1',
+                            model=pm.generic_source_term.exponential,
+                            A1='pore.item1',
+                            A2='pore.item2',
+                            A3='pore.item3',
+                            A4='pore.item4',
+                            A5='pore.item5',
+                            A6='pore.item6',
+                            quantity='pore.mole_fraction',
+                            regen_mode='normal')
+        self.alg = op.algorithms.ReactiveTransport(network=self.net,
+                                                   phase=self.phase)
+        self.alg.set_dirichlet_BC(values=0.4, pores=self.BC_pores)
+        self.alg.set_source(propname='pore.source1',
+                            pores=self.source_pores)
+        self.alg.settings.update({'conductance': 'throat.diffusive_conductance',
+                                  'quantity': 'pore.mole_fraction'})
+        self.alg.run()
+        self.phase.update(self.alg.results())
+        self.phys.regenerate_models(propnames='pore.source1')
+        X = self.phase['pore.mole_fraction']
+        r1 = np.round(np.sum(0.8e-11 * 3 ** (0.5 * X[self.source_pores]**2 -
+                      0.34) + 2e-14), 20)
+        r2 = np.round(np.sum(self.phys['pore.source1.rate'][self.source_pores]), 20)
 #        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
-#        assert r1 == r2
+        assert r1 == r2
 #        assert r2 == -r3
-#
-#    def test_natural_exponential(self):
-#        self.phys['pore.item1'] = 0.8e-11
-#        self.phys['pore.item2'] = 0.5
-#        self.phys['pore.item3'] = 2
-#        self.phys['pore.item4'] = -0.34
-#        self.phys['pore.item5'] = 2e-14
-#        self.phys.add_model(propname='pore.source1',
-#                            model=pm.generic_source_term.natural_exponential,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            A3='pore.item3',
-#                            A4='pore.item4',
-#                            A5='pore.item5',
-#                            x='mole_fraction',
-#                            return_rate=False,
-#                            regen_mode='on_demand')
-#        self.phys.add_model(propname='pore.source2',
-#                            model=pm.generic_source_term.natural_exponential,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            A3='pore.item3',
-#                            A4='pore.item4',
-#                            A5='pore.item5',
-#                            x='mole_fraction',
-#                            return_rate=True,
-#                            regen_mode='on_demand')
-#        self.alg.set_source_term(source_name='pore.source1',
-#                                 pores=self.S_pores,
-#                                 mode='overwrite')
-#        self.alg.run(conductance='throat.diffusive_conductance',
-#                     quantity='pore.mole_fraction',
-#                     super_pore_conductance=None)
-#        self.alg.return_results()
-#        self.phys.regenerate(props='pore.source1')
-#        self.phys.regenerate(props='pore.source2')
-#        X = self.phase['pore.mole_fraction']
-#        r1 = np.round(np.sum(0.8e-11 * np.exp(0.5 * X[self.S_pores] ** 2 -
-#                      0.34) + 2e-14), 20)
-#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+
+    def test_natural_exponential(self):
+        self.phys['pore.item1'] = 0.8e-11
+        self.phys['pore.item2'] = 0.5
+        self.phys['pore.item3'] = 2
+        self.phys['pore.item4'] = -0.34
+        self.phys['pore.item5'] = 2e-14
+        self.phys.add_model(propname='pore.source1',
+                            model=pm.generic_source_term.natural_exponential,
+                            A1='pore.item1',
+                            A2='pore.item2',
+                            A3='pore.item3',
+                            A4='pore.item4',
+                            A5='pore.item5',
+                            quantity='pore.mole_fraction',
+                            regen_mode='normal')
+        self.alg = op.algorithms.ReactiveTransport(network=self.net,
+                                                   phase=self.phase)
+        self.alg.set_dirichlet_BC(values=0.4, pores=self.BC_pores)
+        self.alg.set_source(propname='pore.source1',
+                            pores=self.source_pores)
+        self.alg.settings.update({'conductance': 'throat.diffusive_conductance',
+                                  'quantity': 'pore.mole_fraction'})
+        self.alg.run()
+        self.phase.update(self.alg.results())
+        self.phys.regenerate_models(propnames='pore.source1')
+        X = self.phase['pore.mole_fraction']
+        r1 = np.round(np.sum(0.8e-11 * np.exp(0.5 * X[self.source_pores]**2 -
+                      0.34) + 2e-14), 20)
+        r2 = np.round(np.sum(self.phys['pore.source1.rate'][self.source_pores]), 20)
 #        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
-#        assert r1 == r2
+        assert r1 == r2
 #        assert r2 == -r3
-#
-#    def test_logarithm(self):
-#        self.phys['pore.item1'] = 0.16e-13
-#        self.phys['pore.item2'] = 10
-#        self.phys['pore.item3'] = 4
-#        self.phys['pore.item4'] = 1.4
-#        self.phys['pore.item5'] = 0.133
-#        self.phys['pore.item6'] = -5.1e-13
-#        self.phys.add_model(propname='pore.source1',
-#                            model=pm.generic_source_term.logarithm,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            A3='pore.item3',
-#                            A4='pore.item4',
-#                            A5='pore.item5',
-#                            A6='pore.item6',
-#                            x='mole_fraction',
-#                            return_rate=False,
-#                            regen_mode='on_demand')
-#        self.phys.add_model(propname='pore.source2',
-#                            model=pm.generic_source_term.logarithm,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            A3='pore.item3',
-#                            A4='pore.item4',
-#                            A5='pore.item5',
-#                            A6='pore.item6',
-#                            x='mole_fraction',
-#                            return_rate=True,
-#                            regen_mode='on_demand')
-#        self.alg.set_source_term(source_name='pore.source1',
-#                                 pores=self.S_pores,
-#                                 mode='overwrite')
-#        self.alg.run(conductance='throat.diffusive_conductance',
-#                     quantity='pore.mole_fraction',
-#                     super_pore_conductance=None)
-#        self.alg.return_results()
-#        self.phys.regenerate(props='pore.source1')
-#        self.phys.regenerate(props='pore.source2')
-#        X = self.phase['pore.mole_fraction']
-#        r1 = np.round(np.sum(0.16e-13 * np.log(4 * X[self.S_pores] ** (1.4) +
-#                             0.133) / np.log(10) - 5.1e-13), 20)
-#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+
+    def test_logarithm(self):
+        self.phys['pore.item1'] = 0.16e-13
+        self.phys['pore.item2'] = 10
+        self.phys['pore.item3'] = 4
+        self.phys['pore.item4'] = 1.4
+        self.phys['pore.item5'] = 0.133
+        self.phys['pore.item6'] = -5.1e-13
+        self.phys.add_model(propname='pore.source1',
+                            model=pm.generic_source_term.logarithm,
+                            A1='pore.item1',
+                            A2='pore.item2',
+                            A3='pore.item3',
+                            A4='pore.item4',
+                            A5='pore.item5',
+                            A6='pore.item6',
+                            quantity='pore.mole_fraction',
+                            regen_mode='normal')
+        self.alg = op.algorithms.ReactiveTransport(network=self.net,
+                                                   phase=self.phase)
+        self.alg.set_dirichlet_BC(values=0.4, pores=self.BC_pores)
+        self.alg.set_source(propname='pore.source1',
+                            pores=self.source_pores)
+        self.alg.settings.update({'conductance': 'throat.diffusive_conductance',
+                                  'quantity': 'pore.mole_fraction'})
+        self.alg.run()
+        self.phase.update(self.alg.results())
+        self.phys.regenerate_models(propnames='pore.source1')
+        X = self.phase['pore.mole_fraction']
+        r1 = np.round(np.sum(0.16e-13*np.log(4*X[self.source_pores]**(1.4) +
+                             0.133) / np.log(10) - 5.1e-13), 20)
+        r2 = np.round(np.sum(self.phys['pore.source1.rate'][self.source_pores]), 20)
 #        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
-#        assert r1 == r2
+        assert r1 == r2
 #        assert r2 == -r3
-#
-#    def test_natural_logarithm(self):
-#        self.phys['pore.item1'] = 0.16e-14
-#        self.phys['pore.item2'] = 4
-#        self.phys['pore.item3'] = 1.4
-#        self.phys['pore.item4'] = 0.133
-#        self.phys['pore.item5'] = -5.1e-14
-#        self.phys.add_model(propname='pore.source1',
-#                            model=pm.generic_source_term.natural_logarithm,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            A3='pore.item3',
-#                            A4='pore.item4',
-#                            A5='pore.item5',
-#                            x='mole_fraction',
-#                            return_rate=False,
-#                            regen_mode='on_demand')
-#        self.phys.add_model(propname='pore.source2',
-#                            model=pm.generic_source_term.natural_logarithm,
-#                            A1='pore.item1',
-#                            A2='pore.item2',
-#                            A3='pore.item3',
-#                            A4='pore.item4',
-#                            A5='pore.item5',
-#                            x='mole_fraction',
-#                            return_rate=True,
-#                            regen_mode='on_demand')
-#        self.alg.set_source_term(source_name='pore.source1',
-#                                 pores=self.S_pores,
-#                                 mode='overwrite')
-#        self.alg.run(conductance='throat.diffusive_conductance',
-#                     quantity='pore.mole_fraction',
-#                     super_pore_conductance=None)
-#        self.alg.return_results()
-#        self.phys.regenerate(props='pore.source1')
-#        self.phys.regenerate(props='pore.source2')
-#        X = self.phase['pore.mole_fraction']
-#        r1 = np.round(np.sum(0.16e-14 * np.log(4 * X[self.S_pores] ** (1.4) +
-#                             0.133) - 5.1e-14), 20)
-#        r2 = np.round(np.sum(self.phys['pore.source2'][self.S_pores]), 20)
+
+    def test_natural_logarithm(self):
+        self.phys['pore.item1'] = 0.16e-14
+        self.phys['pore.item2'] = 4
+        self.phys['pore.item3'] = 1.4
+        self.phys['pore.item4'] = 0.133
+        self.phys['pore.item5'] = -5.1e-14
+        self.phys.add_model(propname='pore.source1',
+                            model=pm.generic_source_term.natural_logarithm,
+                            A1='pore.item1',
+                            A2='pore.item2',
+                            A3='pore.item3',
+                            A4='pore.item4',
+                            A5='pore.item5',
+                            quantity='pore.mole_fraction',
+                            regen_mode='on_demand')
+        self.alg = op.algorithms.ReactiveTransport(network=self.net,
+                                                   phase=self.phase)
+        self.alg.set_dirichlet_BC(values=0.4, pores=self.BC_pores)
+        self.alg.set_source(propname='pore.source1',
+                            pores=self.source_pores)
+        self.alg.settings.update({'conductance': 'throat.diffusive_conductance',
+                                  'quantity': 'pore.mole_fraction'})
+        self.alg.run()
+        self.phase.update(self.alg.results())
+        self.phys.regenerate_models(propnames='pore.source1')
+        X = self.phase['pore.mole_fraction']
+        r1 = np.round(np.sum(0.16e-14*np.log(4*X[self.source_pores]**(1.4) +
+                             0.133) - 5.1e-14), 20)
+        r2 = np.round(np.sum(self.phys['pore.source1.rate'][self.source_pores]), 20)
 #        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
-#        assert r1 == r2
+        assert r1 == r2
 #        assert r2 == -r3