Repository: PMEAL/OpenPNM
Pull request #837 (2018-03-31 03:36:50)

Adding more tests
Basically just pulling tests over from V1


Patch:
diff --git a/example_script.py b/example_script.py
index ca6933cf5..2b8341f46 100644
--- a/example_script.py
+++ b/example_script.py
@@ -39,7 +39,7 @@ 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'}
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 74df8ca51..3389e131f 100644
--- a/openpnm/algorithms/GenericTransport.py
+++ b/openpnm/algorithms/GenericTransport.py
@@ -160,9 +160,13 @@ class GenericTransport(GenericAlgorithm):
             self.A.setdiag(datadiag)
             self.A.eliminate_zeros()  # Remove 0 entries
 
-    def setup(self):
+    def setup(self, quantity=None, conductance=None):
         r"""
         """
+        if quantity:
+            self.settings['quantity'] = quantity
+        if conductance:
+            self.settings['conductance'] = conductance
         self.build_A()
         self.build_b()
         self.apply_BCs()
diff --git a/openpnm/algorithms/ReactiveTransport.py b/openpnm/algorithms/ReactiveTransport.py
index 2aa1086b2..da4e5d73d 100644
--- a/openpnm/algorithms/ReactiveTransport.py
+++ b/openpnm/algorithms/ReactiveTransport.py
@@ -32,10 +32,10 @@ class ReactiveTransport(GenericTransport):
             # 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"""
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/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')


Test patch:
diff --git a/tests/integration/AlgorithmTest.py b/tests/integration/AlgorithmTest.py
new file mode 100644
index 000000000..7864e45b3
--- /dev/null
+++ b/tests/integration/AlgorithmTest.py
@@ -0,0 +1,126 @@
+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_dirchlet_BC(values=T_in, pores=inlets)
+    alg.set_dirchlet_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():
+    pn = op.network.Cubic([5, 5, 5], spacing=1)
+    pn.add_boundary_pores()
+    Ps = pn.pores('*boundary', mode='not')
+    Ts = pn.find_neighbor_throats(pores=Ps, mode='intersection', flatten=True)
+    geom = op.geometry.StickAndBall(network=pn, pores=Ps, throats=Ts)
+    geom['pore.diameter'] = 0.999999
+    geom['throat.diameter'] = 0.999999
+    geom.regenerate(['pore.diameter', 'throat.diameter'], mode='exclude')
+    Ps = pn.pores('boundary')
+    Ts = pn.find_neighbor_throats(pores=Ps, mode='not_intersection')
+    boun = op.geometry.Boundary(network=pn, pores=Ps, throats=Ts)
+    air = op.phases.Air(network=pn)
+    Ps = pn.pores()
+    Ts = pn.throats()
+    phys_air = op.physics.Standard(network=pn, phase=air, pores=Ps, throats=Ts)
+    BC1_pores = pn.pores(labels=['top_boundary'])
+    BC2_pores = pn.pores(labels=['bottom_boundary'])
+    Diff = 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()
+    air.update(Diff.results())
+    Diff_deff = Diff.calc_eff_diffusivity()/np.mean(air['pore.diffusivity'])
+    assert np.round(Diff_deff, 3) == 1
+
+
+def test_Darcy_alg():
+    # 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_dirchlet_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_dirchlet_BC(values=P_in, pores=inlets)
+    alg2.set_dirchlet_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/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/MiscTest.py b/tests/unit/core/MiscTest.py
new file mode 100644
index 000000000..ed694bdbb
--- /dev/null
+++ b/tests/unit/core/MiscTest.py
@@ -0,0 +1,67 @@
+import openpnm as op
+import scipy as sp
+import pytest
+
+
+class MiscTest:
+
+    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 = MiscTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/models/physics/GenericSourceTerm.py b/tests/unit/models/physics/GenericSourceTerm.py
index 782258e6d..8c937159d 100644
--- a/tests/unit/models/physics/GenericSourceTerm.py
+++ b/tests/unit/models/physics/GenericSourceTerm.py
@@ -15,260 +15,204 @@ 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.setup(conductance='throat.diffusive_conductance',
+                       quantity='pore.mole_fraction')
+        self.alg.set_BC(bctype='dirichlet', bcvalues=0.4,
+                        pores=self.BC_pores)
+        self.alg.set_source_term(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_BC(bctype='dirichlet', bcvalues=0.4,
+                        pores=self.BC_pores)
+        self.alg.set_source_term(propname='pore.source1',
+                                 pores=self.source_pores)
+        self.alg.setup(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_BC(bctype='dirichlet', bcvalues=0.4,
+                        pores=self.BC_pores)
+        self.alg.set_source_term(propname='pore.source1',
+                                 pores=self.source_pores)
+        self.alg.setup(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_BC(bctype='dirichlet', bcvalues=0.4,
+                        pores=self.BC_pores)
+        self.alg.set_source_term(propname='pore.source1',
+                                 pores=self.source_pores)
+        self.alg.setup(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_BC(bctype='dirichlet', bcvalues=0.4,
+                        pores=self.BC_pores)
+        self.alg.set_source_term(propname='pore.source1',
+                                 pores=self.source_pores)
+        self.alg.setup(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_BC(bctype='dirichlet', bcvalues=0.4,
+                        pores=self.BC_pores)
+        self.alg.set_source_term(propname='pore.source1',
+                                 pores=self.source_pores)
+        self.alg.setup(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