Repository: PMEAL/OpenPNM
Pull request #884 (2018-04-27 05:18:08)

Test utils
cleaning up the utils module and adding tests


Patch:
diff --git a/conda_requirements.txt b/conda_requirements.txt
index 3a6e441c7..e791d1b43 100644
--- a/conda_requirements.txt
+++ b/conda_requirements.txt
@@ -3,4 +3,3 @@ scipy
 matplotlib
 pandas
 scikit-image
-sympy
diff --git a/openpnm/models/physics/electrical_conductance.py b/openpnm/models/physics/electrical_conductance.py
index ca728ee31..25cf07dd9 100644
--- a/openpnm/models/physics/electrical_conductance.py
+++ b/openpnm/models/physics/electrical_conductance.py
@@ -1,3 +1,10 @@
+r"""
+===============================================================================
+Submodule -- electrical_conductance
+===============================================================================
+
+"""
+
 import scipy as _sp
 
 
@@ -5,8 +12,7 @@ def series_resistors(target, conductivity='pore.electrical_conductivity',
                      pore_area='pore.area', pore_diameter='pore.diameter',
                      throat_area='throat.area', throat_length='throat.length'):
     r"""
-    Calculates the electrical conductance of throat assuming cylindrical
-    geometry
+    Calculates the electrical conductance of throat assuming cylindrical geometry
 
     Parameters
     ----------
diff --git a/openpnm/models/physics/generic_source_term.py b/openpnm/models/physics/generic_source_term.py
index dc038be45..2ec020d3a 100644
--- a/openpnm/models/physics/generic_source_term.py
+++ b/openpnm/models/physics/generic_source_term.py
@@ -1,5 +1,4 @@
 import scipy as _sp
-import sympy as _syp
 
 
 def standard_kinetics(target, quantity, prefactor, exponent):
@@ -17,706 +16,353 @@ def standard_kinetics(target, quantity, prefactor, exponent):
     return values
 
 
-def _parse_args(target, key, default):
-    if key == '':
-        val = default
-    else:
-        val = target[key]
-    return val
-
-
-def linear(target, X, A1='', A2=''):
+def linear(target, A1='', A2='', quantity=''):
     r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of `X`:
+    For the following source term:
 
         .. math::
-            r = A_{1}   X  +  A_{2}
-
-    Parameters
-    ----------
-    A1 -> A2 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
-
-    X : string
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
-
-    Returns
-    -------
-    A dictionary containing the following three items:
-
-        **'rate'** - The value of the source term function at the given X.
-
-        **'S1'** - The slope of the source term function at the given X.
+            r = A_{1}   x  +  A_{2}
 
-        **'S2'** - The intercept of the source term function at the given X.
+    If return_rate is True, it returns the value of source term for the
+    provided x in each pore.
 
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
+    If return_rate is False, it calculates the slope and intercept for the
+    following linear form :
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
-
-    """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=0.0)
-    X = target[X]
-
-    r = A * X + B
-    S1 = A
-    S2 = B
-    values = {'S1': S1, 'S2': S2, 'rate': r}
-    return values
-
-
-def power_law(target, X, A1='', A2='', A3=''):
-    r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of *X*:
-
-        .. math::
-            r = A_{1}   x^{A_{2}}  +  A_{3}
+            r = S_{1}   x  +  S_{2}
 
     Parameters
     ----------
-    A1 -> A3 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+    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
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
+    quantity : string
+        The property name for the main quantity being produced or consumed.
 
     Returns
     -------
-    A dictionary containing the following three items:
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
 
-        **'rate'** - The value of the source term function at the given X.
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
 
-        **'S1'** - The slope of the source term function at the given X.
+    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:
 
-        **'S2'** - The intercept of the source term function at the given X.
+        obj['pore.source_term']  = d
 
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
+    is equivalent to:
 
-        .. math::
-            rate = S_{1}   X  +  S_{2}
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=1.0)
-    C = _parse_args(target=target, key=A3, default=0.0)
-    X = target[X]
-
-    r = A * X ** B + C
-    S1 = A * B * X ** (B - 1)
-    S2 = A * X ** B * (1 - B) + C
+    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 exponential(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
+def power_law(target, A1='', A2='', A3='', quantity=''):
     r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of `X`:
+    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::
-            r =  A_{1} A_{2}^{( A_{3} x^{ A_{4} } + A_{5})} + A_{6}
+            r = S_{1}   x  +  S_{2}
 
     Parameters
     ----------
-    A1 -> A6 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+    A1 -> A3 : string
+        The property name of the coefficients in the source term model
 
-    X : string
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
+    quantity : string
+        The property name for the main quantity being produced or consumed.
 
     Returns
     -------
-    A dictionary containing the following three items:
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
 
-        **'rate'** - The value of the source term function at the given X.
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
 
-        **'S1'** - The slope of the source term function at the given X.
+    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:
 
-        **'S2'** - The intercept of the source term function at the given X.
+        obj['pore.source_term']  = d
 
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
+    is equivalent to:
 
-        .. math::
-            rate = S_{1}   X  +  S_{2}
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=1.0)
-    C = _parse_args(target=target, key=A3, default=1.0)
-    D = _parse_args(target=target, key=A4, default=1.0)
-    E = _parse_args(target=target, key=A5, default=0.0)
-    F = _parse_args(target=target, key=A6, default=0.0)
-    X = target[X]
-
-    r = A * B ** (C * X ** D + E) + F
-    S1 = A * C * D * X ** (D - 1) * _sp.log(B) * B ** (C * X ** D + E)
-    S2 = A * B ** (C * X ** D + E) * (1 - C * D * _sp.log(B) * X ** D) + F
+    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 natural_exponential(target, X, A1='', A2='', A3='', A4='', A5=''):
+def exponential(target, A1='', A2='', A3='', A4='', A5='', A6='',
+                quantity='', return_rate=True):
     r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of `X`:
-
+    For the following source term:
         .. math::
-            r =   A_{1} exp( A_{2}  x^{ A_{3} } + A_{4} )+ A_{5}
-
-    Parameters
-    ----------
-    A1 -> A5 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
-
-    X : string
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
-
-    Returns
-    -------
-    A dictionary containing the following three items:
-
-        **'rate'** - The value of the source term function at the given X.
-
-        **'S1'** - The slope of the source term function at the given X.
-
-        **'S2'** - The intercept of the source term function at the given X.
-
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
-
-        .. math::
-            rate = S_{1}   X  +  S_{2}
-
-    """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=1.0)
-    C = _parse_args(target=target, key=A3, default=1.0)
-    D = _parse_args(target=target, key=A4, default=1.0)
-    E = _parse_args(target=target, key=A5, default=0.0)
-    X = target[X]
-
-    r = A * _sp.exp(B * X ** C + D) + E
-    S1 = A * B * C * X ** (C - 1) * _sp.exp(B * X ** C + D)
-    S2 = A * (1 - B * C * X ** C) * _sp.exp(B * X ** C + D) + E
-    values = {'pore.S1': S1, 'pore.S2': S2, 'pore.rate': r}
-    return values
-
+            r =  A_{1} A_{2}^{( A_{3} x^{ A_{4} } + A_{5})} + A_{6}
 
-def logarithm(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
-    r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of `X`:
+    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 =  A_{1}   Log_{ A_{2} }( A_{3} x^{ A_{4} }+ A_{5})+ A_{6}
+            r = S_{1}   x  +  S_{2}
 
     Parameters
     ----------
     A1 -> A6 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+        The property name of the coefficients in the source term model
 
-    X : string
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
+    quantity : string
+        The property name for the main quantity being produced or consumed.
 
     Returns
     -------
-    A dictionary containing the following three items:
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
 
-        **'rate'** - The value of the source term function at the given X.
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
 
-        **'S1'** - The slope of the source term function at the given X.
+    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:
 
-        **'S2'** - The intercept of the source term function at the given X.
+        obj['pore.source_term']  = d
 
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
+    is equivalent to:
 
-        .. math::
-            rate = S_{1}   X  +  S_{2}
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=1.0)
-    C = _parse_args(target=target, key=A3, default=1.0)
-    D = _parse_args(target=target, key=A4, default=1.0)
-    E = _parse_args(target=target, key=A5, default=0.0)
-    F = _parse_args(target=target, key=A6, default=0.0)
-    X = target[X]
-
-    r = (A * _sp.log(C * X ** D + E)/_sp.log(B) + F)
-    S1 = A * C * D * X ** (D - 1) / (_sp.log(B) * (C * X ** D + E))
-    S2 = A * _sp.log(C * X ** D + E) / _sp.log(B) + F - A * C * D * X ** D / \
-        (_sp.log(B) * (C * X ** D + E))
+    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:
+            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_logarithm(target, X, A1='', A2='', A3='', A4='', A5=''):
+def natural_exponential(target, A1='', A2='', A3='', A4='', A5='',
+                        quantity=''):
     r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of `X`:
+    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::
-            r =   A_{1}  Ln( A_{2} x^{ A_{3} }+ A_{4})+ A_{5}
+            r = S_{1}   x  +  S_{2}
 
     Parameters
     ----------
     A1 -> A5 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+        The property name of the coefficients in the source term model
 
-    X : string
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
+    quantity : string
+        The property name for the main quantity being produced or consumed.
 
     Returns
     -------
-    A dictionary containing the following three items:
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
 
-        **'rate'** - The value of the source term function at the given X.
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
 
-        **'S1'** - The slope of the source term function at the given X.
+    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:
 
-        **'S2'** - The intercept of the source term function at the given X.
+        obj['pore.source_term']  = d
 
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
+    is equivalent to:
 
-        .. math::
-            rate = S_{1}   X  +  S_{2}
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=1.0)
-    C = _parse_args(target=target, key=A3, default=1.0)
-    D = _parse_args(target=target, key=A4, default=1.0)
-    E = _parse_args(target=target, key=A5, default=0.0)
-    X = target[X]
-
-    r = A*_sp.log(B*X**C + D) + E
-    S1 = A*B*C*X**(C - 1) / (B * X ** C + D)
-    S2 = A*_sp.log(B*X**C + D) + E - A*B*C*X**C / (B*X**C + D)
-    values = {'pore.S1': S1, 'pore.S2': S2, 'pore.rate': r}
-    return values
-
-
-def _build_func(eq, **args):
-    r'''
-    Take a symbolic equation and return the lambdified version plus the
-    linearization of form S1 * x + S2
-    '''
-    eq_prime = eq.diff(args['x'])
-    s1 = eq_prime
-    s2 = eq - eq_prime*args['x']
-    EQ = _syp.lambdify(args.values(), expr=eq, modules='numpy')
-    S1 = _syp.lambdify(args.values(), expr=s1, modules='numpy')
-    S2 = _syp.lambdify(args.values(), expr=s2, modules='numpy')
-    return EQ, S1, S2
-
-
-def linear_sym(target, X, A1='', A2=''):
-    r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of *x*:
-
-        .. math::
-            r = A_{1}   x  +  A_{2}
-
-    Parameters
-    ----------
-    A1 -> A2 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
-
-    X : string
-        The dictionary key on the target object containing the the quantity
-        of interest
-
-    Returns
-    -------
-    A dictionary containing the following three items:
-
-        **'rate'** - The value of the source term function at the given 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:
+            a[item] = target[a[item]]
 
-        **'S1'** - The slope of the source term function at the given X.
+    r = a['1'] * _sp.exp(a['2'] * X ** a['3'] + a['4']) + a['5']
 
-        **'S2'** - The intercept of the source term function at the given X.
+    S1 = a['1'] * a['2'] * a['3'] * X ** (a['3'] - 1) * \
+        _sp.exp(a['2'] * X ** a['3'] + a['4'])
 
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
+    S2 = a['1'] * (1 - a['2'] * a['3'] * X ** a['3']) * \
+        _sp.exp(a['2'] * X ** a['3'] + a['4']) + a['5']
 
-        .. math::
-            rate = S_{1}   X  +  S_{2}
-
-    """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=0.0)
-    X = target[X]
-    # Symbols used in symbolic function
-    a, b, x = _syp.symbols('a,b,x')
-    # Equation
-    y = a*x + b
-    # Callable functions
-    r, s1, s2 = _build_func(eq=y, a=a, b=b,  x=x)
-    # Values
-    r_val = r(A, B, X)
-    s1_val = s1(A, B, X)
-    s2_val = s2(A, B, X)
-    values = {'S1': s1_val, 'S2': s2_val, 'rate': r_val}
+    values = {'S1': S1, 'S2': S2, 'rate': r}
     return values
 
 
-def power_law_sym(target, X, A1='', A2='', A3=''):
+def logarithm(target, A1='', A2='', A3='', A4='', A5='', A6='', quantity=''):
     r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of *x*:
-
-        .. math::
-            r = A_{1}   x^{A_{2}}  +  A_{3}
-
-    Parameters
-    ----------
-    A1 -> A3 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
-
-    X : string
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
-
-    Returns
-    -------
-    A dictionary containing the following three items:
-
-        **'rate'** - The value of the source term function at the given X.
-
-        **'S1'** - The slope of the source term function at the given X.
-
-        **'S2'** - The intercept of the source term function at the given X.
-
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
-
+    For the following source term:
         .. math::
-            rate = S_{1}   X  +  S_{2}
-
-    """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=1.0)
-    C = _parse_args(target=target, key=A3, default=0.0)
-    X = target[X]
-    # Symbols used in symbolic function
-    a, b, c, x = _syp.symbols('a,b,c,x')
-    # Equation
-    y = a*x**b + c
-    # Callable functions
-    r, s1, s2 = _build_func(eq=y, a=a, b=b, c=c, x=x)
-    # Values
-    r_val = r(A, B, C, X)
-    s1_val = s1(A, B, C, X)
-    s2_val = s2(A, B, C, X)
-    values = {'S1': s1_val, 'S2': s2_val, 'rate': r_val}
-    return values
-
+            r =  A_{1}   Log_{ A_{2} }( A_{3} x^{ A_{4} }+ A_{5})+ A_{6}
 
-def exponential_sym(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
-    r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of *x*:
+    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 =  A_{1} A_{2}^{( A_{3} x^{ A_{4} } + A_{5})} + A_{6}
+            r = S_{1}   x  +  S_{2}
 
     Parameters
     ----------
     A1 -> A6 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+        The property name of the coefficients in the source term model
 
-    X : string or float/int or array/list
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
+    quantity : string
+        The property name for the main quantity being produced or consumed.
 
     Returns
     -------
-    A dictionary containing the following three items:
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
 
-        **'rate'** - The value of the source term function at the given X.
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
 
-        **'S1'** - The slope of the source term function at the given X.
+    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:
 
-        **'S2'** - The intercept of the source term function at the given X.
+        obj['pore.source_term']  = d
 
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
+    is equivalent to:
 
-        .. math::
-            rate = S_{1}   X  +  S_{2}
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=1.0)
-    C = _parse_args(target=target, key=A3, default=1.0)
-    D = _parse_args(target=target, key=A4, default=1.0)
-    E = _parse_args(target=target, key=A5, default=0.0)
-    F = _parse_args(target=target, key=A6, default=0.0)
-    X = target[X]
-    # Symbols used in symbolic function
-    a, b, c, d, e, f, x = _syp.symbols('a,b,c,d,e,f,x')
-    # Equation
-    y = a*b**(c*x**d + e) + f
-    # Callable functions
-    r, s1, s2 = _build_func(eq=y, a=a, b=b, c=c, d=d, e=e, f=f, x=x)
-    # Values
-    r_val = r(A, B, C, D, E, F, X)
-    s1_val = s1(A, B, C, D, E, F, X)
-    s2_val = s2(A, B, C, D, E, F, X)
-    values = {'S1': s1_val, 'S2': s2_val, 'rate': r_val}
-    return values
-
-
-def natural_exponential_sym(target, X, A1='', A2='', A3='', A4='', A5=''):
-    r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of *x*:
-
-        .. math::
-            r =   A_{1} exp( A_{2}  x^{ A_{3} } + A_{4} )+ A_{5}
-
-    Parameters
-    ----------
-    A1 -> A6 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
-
-    X : string or float/int or array/list
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
-
-    Returns
-    -------
-    A dictionary containing the following three items:
-
-        **'rate'** - The value of the source term function at the given X.
+    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:
+            a[item] = target[a[item]]
 
-        **'S1'** - The slope of the source term function at the given X.
+    r = a['1'] * _sp.log(a['3'] * X ** a['4'] + a['5']) / _sp.log(a['2']) + \
+        a['6']
 
-        **'S2'** - The intercept of the source term function at the given X.
+    S1 = a['1'] * a['3'] * a['4'] * X ** (a['4'] - 1) / \
+        (_sp.log(a['2']) * (a['3'] * X ** a['4'] + a['5']))
 
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
+    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']))
 
-        .. math::
-            rate = S_{1}   X  +  S_{2}
-
-    """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=1.0)
-    C = _parse_args(target=target, key=A3, default=1.0)
-    D = _parse_args(target=target, key=A4, default=1.0)
-    E = _parse_args(target=target, key=A5, default=0.0)
-    X = target[X]
-    # Symbols used in symbolic function
-    a, b, c, d, e, x = _syp.symbols('a,b,c,d,e,x')
-    # Equation
-    y = a*_syp.exp(b*x**c + d) + e
-    # Callable functions
-    r, s1, s2 = _build_func(eq=y, a=a, b=b, c=c, d=d, e=e, x=x)
-    # Values
-    r_val = r(A, B, C, D, E, X)
-    s1_val = s1(A, B, C, D, E, X)
-    s2_val = s2(A, B, C, D, E, X)
-    values = {'S1': s1_val, 'S2': s2_val, 'rate': r_val}
+    values = {'S1': S1, 'S2': S2, 'rate': r}
     return values
 
 
-def logarithm_sym(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
+def natural_logarithm(target, A1='', A2='', A3='', A4='', A5='', quantity=''):
     r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of *x*:
+    For the following source term:
 
         .. math::
-            r =  A_{1}   Log_{ A_{2} }( A_{3} x^{ A_{4} }+ A_{5})+ A_{6}
-
-    Parameters
-    ----------
-    A1 -> A6 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
-
-    x : string or float/int or array/list
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
-
-    Returns
-    -------
-    A dictionary containing the following three items:
-
-        **'rate'** - The value of the source term function at the given X.
-
-        **'S1'** - The slope of the source term function at the given X.
-
-        **'S2'** - The intercept of the source term function at the given X.
+            r =   A_{1}  Ln( A_{2} x^{ A_{3} }+ A_{4})+ A_{5}
 
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
+    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::
-            rate = S_{1}   X  +  S_{2}
-
-    """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=1.0)
-    C = _parse_args(target=target, key=A3, default=1.0)
-    D = _parse_args(target=target, key=A4, default=1.0)
-    E = _parse_args(target=target, key=A5, default=0.0)
-    F = _parse_args(target=target, key=A6, default=0.0)
-    X = target[X]
-    # Symbols used in symbolic function
-    a, b, c, d, e, f, x = _syp.symbols('a,b,c,d,e,f,x')
-    # Equation
-    y = a*_syp.log((c*x**d + e), b) + f
-    # Callable functions
-    r, s1, s2 = _build_func(eq=y, a=a, b=b, c=c, d=d, e=e, f=f, x=x)
-    # Values
-    r_val = r(A, B, C, D, E, F, X)
-    s1_val = s1(A, B, C, D, E, F, X)
-    s2_val = s2(A, B, C, D, E, F, X)
-    values = {'S1': s1_val, 'S2': s2_val, 'rate': r_val}
-    return values
-
-
-def natural_logarithm_sym(target, X, A1='', A2='', A3='', A4='', A5=''):
-    r"""
-    Calculates the rate, as well as slope and intercept of the following
-    function at the given value of *x*:
+            r = S_{1}   x  +  S_{2}
 
-        .. math::
-            rate =   A_{1}  Ln( A_{2} x^{ A_{3} }+ A_{4})+ A_{5}
+    which can be used in linear solvers to produce a better guess of
+    ``quantity``, which can
 
     Parameters
     ----------
     A1 -> A5 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+        The property name of the coefficients in the source term model
 
-    X : string or float/int or array/list
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
+    quantity : string
+        The property name for the main quantity being produced or consumed.
 
     Returns
     -------
-    A dictionary containing the following three items:
+    A dictionary containing the values for *S1*, *S2*, and *rate*:
 
-        **'rate'** - The value of the source term function at the given X.
+        d = {'pore.S1': vals_1, 'pore.S2': vals_2, 'pore.rate': vals_3}
 
-        **'S1'** - The slope of the source term function at the given X.
+    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:
 
-        **'S2'** - The intercept of the source term function at the given X.
+        obj['pore.source_term']  = d
 
-    The slope and intercept provide a linearized source term equation about the
-    current value of X as follow:
+    is equivalent to:
 
-        .. math::
-            rate = S_{1}   X  +  S_{2}
+        obj['pore.source_term_S1'] = d['pore.S1']
+        obj['pore.source_term_S2'] = d['pore.S2']
+        obj['pore.source_term_rate'] = d['pore.rate']
 
     """
-    A = _parse_args(target=target, key=A1, default=1.0)
-    B = _parse_args(target=target, key=A2, default=1.0)
-    C = _parse_args(target=target, key=A3, default=1.0)
-    D = _parse_args(target=target, key=A4, default=1.0)
-    E = _parse_args(target=target, key=A5, default=0.0)
-    X = target[X]
-    # Symbols used in symbolic function
-    a, b, c, d, e, x = _syp.symbols('a,b,c,d,e,x')
-    # Equation
-    y = a*_syp.ln(b*x**c + d) + e
-    # Callable functions
-    r, s1, s2 = _build_func(eq=y, a=a, b=b, c=c, d=d, e=e, x=x)
-    # Values
-    r_val = r(A, B, C, D, E, X)
-    s1_val = s1(A, B, C, D, E, X)
-    s2_val = s2(A, B, C, D, E, X)
-    values = {'S1': s1_val, 'S2': s2_val, 'rate': r_val}
-    return values
-
-
-def general_symbolic(target, eqn=None, arg_map=None):
-    r'''
-    A general function to interpret a sympy equation and evaluate the linear
-    components of the source term.
+    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:
+            a[item] = target[a[item]]
 
-    Parameters
-    ----------
-    target : OpenPNM object
-        The OpenPNM object where the result will be applied.
+    r = a['1'] * _sp.log(a['2'] * X ** a['3'] + a['4']) + a['5']
 
-    eqn : sympy symbolic expression for the source terms
-        e.g. y = a*x**b + c
+    S1 = a['1'] * a['2'] * a['3'] * X ** (a['3'] - 1) / \
+        (a['2'] * X ** a['3'] + a['4'])
 
-    arg_map : Dict mapping the symbols in the expression to OpenPNM data
-        on the target. Must contain 'x' which is the independent variable.
-        e.g. arg_map={'a':'pore.a', 'b':'pore.b', 'c':'pore.c', 'x':'pore.x'}
+    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'])
 
-    Example
-    ----------
-    >>> import openpnm as op
-    >>> from openpnm.models.physics import generic_source_term as gst
-    >>> import scipy as sp
-    >>> import sympy as _syp
-    >>> pn = op.network.Cubic(shape=[5, 5, 5], spacing=0.0001)
-    >>> water = op.phases.Water(network=pn)
-    >>> water['pore.a'] = 1
-    >>> water['pore.b'] = 2
-    >>> water['pore.c'] = 3
-    >>> water['pore.x'] = sp.random.random(water.Np)
-    >>> a, b, c, x = _syp.symbols('a,b,c,x')
-    >>> y = a*x**b + c
-    >>> arg_map = {'a':'pore.a', 'b':'pore.b', 'c':'pore.c', 'x':'pore.x'}
-    >>> water.add_model(propname='pore.general',
-    ...                 model=gst.general_symbolic,
-    ...                 eqn=y, arg_map=arg_map,
-    ...                 regen_mode='normal')
-    >>> assert 'pore.general.rate' in water.props()
-    >>> assert 'pore.general.S1' in water.props()
-    >>> assert 'pore.general.S1' in water.props()
-    '''
-    # First make sure all the symbols have been allocated dict items
-    for arg in _syp.postorder_traversal(eqn):
-        if _syp.srepr(arg)[:6] == 'Symbol':
-            key = _syp.srepr(arg)[7:].strip('(').strip(')').strip("'")
-            if key not in arg_map.keys():
-                raise Exception('argument mapping incomplete, missing '+key)
-    if 'x' not in arg_map.keys():
-        raise Exception('argument mapping must contain "x" for the ' +
-                        'independent variable')
-    # Get the data
-    data = {}
-    args = {}
-    for key in arg_map.keys():
-        data[key] = target[arg_map[key]]
-        # Callable functions
-        args[key] = _syp.symbols(key)
-    r, s1, s2 = _build_func(eqn, **args)
-    r_val = r(*data.values())
-    s1_val = s1(*data.values())
-    s2_val = s2(*data.values())
-    values = {'S1': s1_val, 'S2': s2_val, 'rate': r_val}
+    values = {'S1': S1, 'S2': S2, 'rate': r}
     return values
diff --git a/openpnm/network/Cubic.py b/openpnm/network/Cubic.py
index 05b6e65a4..9f4f86929 100644
--- a/openpnm/network/Cubic.py
+++ b/openpnm/network/Cubic.py
@@ -145,8 +145,8 @@ class Cubic(GenericNetwork):
         Notes
         -----
         This method uses ``clone_pores`` to clone the surface pores (labeled
-        'left','right', etc), then shifts them to the periphery of the domain,
-        and gives them the label 'right_face', 'left_face', etc.
+        'left', 'right', etc), then shifts them to the periphery of the domain,
+        and gives them the label 'right_boundary', 'left_boundary', etc.
         """
         if type(labels) == str:
             labels = [labels]
@@ -188,3 +188,48 @@ class Cubic(GenericNetwork):
         return spacing
 
     spacing = property(fget=_get_spacing)
+
+    def to_array(self, values):
+        r"""
+        Converts the values to a rectangular array with the same shape as the
+        network
+
+        Parameters
+        ----------
+        values : array_like
+            An Np-long array of values to convert to
+
+        Notes
+        -----
+        This method can break on networks that have had boundaries added.  It
+        will usually work IF the given values came only from 'internal'
+        pores.
+        """
+        if sp.shape(values)[0] > self.num_pores('internal'):
+            raise Exception('The array shape does not match the network')
+        Ps = sp.array(self['pore.index'][self.pores('internal')], dtype=int)
+        arr = sp.ones(self._shape)*sp.nan
+        ind = sp.unravel_index(Ps, self._shape)
+        arr[ind[0], ind[1], ind[2]] = values
+        return arr
+
+    def from_array(self, array, propname):
+        r"""
+        Apply data to the network based on a rectangular array filled with
+        values.  Each array location corresponds to a pore in the network.
+        Parameters
+        ----------
+        array : array_like
+            The rectangular array containing the values to be added to the
+            network. This array must be the same shape as the original network.
+        propname : string
+            The name of the pore property being added.
+        """
+        array = sp.atleast_3d(array)
+        if sp.shape(array) != self._shape:
+            raise Exception('The array shape does not match the network')
+        temp = array.flatten()
+        Ps = sp.array(self['pore.index'][self.pores('internal')], dtype=int)
+        propname = 'pore.' + propname.split('.')[-1]
+        self[propname] = sp.nan
+        self[propname][self.pores('internal')] = temp[Ps]
diff --git a/openpnm/network/CubicDual.py b/openpnm/network/CubicDual.py
index 0cde5c1d0..91775329d 100644
--- a/openpnm/network/CubicDual.py
+++ b/openpnm/network/CubicDual.py
@@ -59,7 +59,7 @@ class CubicDual(GenericNetwork):
         faces = [['front', 'back'], ['left', 'right'], ['top', 'bottom']]
         faces = [faces[i] for i in sp.where(~single_dim)[0]]
         faces = sp.array(faces).flatten().tolist()
-        dual.add_boundaries(faces)
+        dual.add_boundary_pores(faces)
         # Add secondary network name as a label
         dual['pore.'+label_2] = True
         dual['throat.'+label_2] = True
@@ -92,4 +92,4 @@ class CubicDual(GenericNetwork):
         net['throat.internal'][Ts] = True
         # Transfer all dictionary items from 'net' to 'self'
         [self.update({item: net[item]}) for item in net]
-        del self.workspace[net.name]
+        self.project.remove(self)
diff --git a/openpnm/network/CubicTemplate.py b/openpnm/network/CubicTemplate.py
index a5025cc5f..34cf54dbb 100644
--- a/openpnm/network/CubicTemplate.py
+++ b/openpnm/network/CubicTemplate.py
@@ -1,27 +1,15 @@
-# -*- coding: utf-8 -*-
-"""
-===============================================================================
-CubicTemplate: Generate cubic lattice networks of arbitrary overall shape
-===============================================================================
-
-"""
 import scipy as sp
 from openpnm.network import Cubic
 from openpnm import topotools
 
-shapes = [[10, 10, 10], [20, 20, 5], [5, 10, 8]]
-spacings = [1, 0.5, [2, 1, 0.8]]
-origins = [[0, 0, 0], [0, 0, 10], [0, 0, 12.5]]
-labels = None
-
 
 class CubicTemplate(Cubic):
     r"""
     """
-    def __init__(self, template, **kwargs):
+    def __init__(self, template, spacing, **kwargs):
 
         template = sp.atleast_3d(template)
-        super().__init__(shape=template.shape, **kwargs)
+        super().__init__(shape=template.shape, spacing=spacing, **kwargs)
 
         coords = sp.unravel_index(range(template.size), template.shape)
         self['pore.template_coords'] = sp.vstack(coords).T
@@ -29,9 +17,3 @@ class CubicTemplate(Cubic):
         self['pore.drop'] = template.flatten() == 0
         topotools.trim(network=self, pores=self.pores('drop'))
         del self['pore.drop']
-
-    def toarray(self, prop):
-        pass
-
-    def fromarray(self, array):
-        pass
diff --git a/openpnm/utils/misc.py b/openpnm/utils/misc.py
index 9890fb52f..586f4f823 100644
--- a/openpnm/utils/misc.py
+++ b/openpnm/utils/misc.py
@@ -1,6 +1,6 @@
 import scipy as _sp
 import time as _time
-from collections import OrderedDict, abc
+from collections import OrderedDict
 
 
 class PrintableList(list):
@@ -9,8 +9,7 @@ class PrintableList(list):
         lines = [horizontal_rule]
         self.sort()
         for i, item in enumerate(self):
-            if '._' not in item:
-                lines.append('{0}\t: {1}'.format(i + 1, item))
+            lines.append('{0}\t: {1}'.format(i + 1, item))
         lines.append(horizontal_rule)
         return '\n'.join(lines)
 
@@ -45,6 +44,7 @@ class PrintableDict(OrderedDict):
 
 
 class SettingsDict(PrintableDict):
+
     def __missing__(self, key):
         self[key] = None
         return self[key]
@@ -159,7 +159,7 @@ def toc(quiet=False):
         else:
             return t
     else:
-        print("Toc: start time not set")
+        raise Exception('Start time not set, call tic first')
 
 
 def unique_list(input_list):
@@ -202,7 +202,7 @@ def flat_list(input_list):
 def sanitize_dict(input_dict):
     r"""
     Given a nested dictionary, ensures that all nested dicts are normal
-    Python dict.  This is necessary for pickling, or just converting
+    Python dicts.  This is necessary for pickling, or just converting
     an 'auto-vivifying' dict to something that acts normal.
     """
     plain_dict = dict()
@@ -215,73 +215,6 @@ def sanitize_dict(input_dict):
     return plain_dict
 
 
-def amalgamate_data(objs=[], delimiter='_'):
-    r"""
-    Returns a dictionary containing ALL pore data from all netowrk and/or
-    phase objects received as arguments
-
-    Parameters
-    ----------
-    obj : list of OpenPNM objects
-        The network and Phase objects whose data should be amalgamated into a
-        single dict
-
-    delimiter : string
-        The delimiter to place between the prop name and the object name.  For
-        instance \'pore.air_molar_density\' or \'pore.air|molar_density'\.  The
-        use of underscores can be problematic for reloading the data since they
-        are also used in multiple word properties.  The default is '_' for
-        backwards compatibility, but the '|' option is preferred.
-
-    Returns
-    -------
-    A standard Python dict containing all the data from the supplied OpenPNM
-    objects
-    """
-    if type(objs) is not list:
-        objs = list(objs)
-    data_amalgamated = {}
-    dlim = delimiter
-    exclusion_list = ['pore.centroid', 'pore.vertices', 'throat.centroid',
-                      'throat.offset_vertices', 'throat.vertices', 'throat.normal',
-                      'throat.perimeter', 'pore.vert_index', 'throat.vert_index']
-    for item in objs:
-        mro = item._mro()
-        # If Network object, combine Geometry and Network keys
-        if 'GenericNetwork' in mro:
-            keys = []
-            for key in list(item.keys()):
-                keys.append(key)
-            for geom in item._geometries:
-                for key in list(geom.keys()):
-                    if key not in keys:
-                        keys.append(key)
-        else:
-            if 'GenericPhase' in mro:
-                keys = []
-                for key in list(item.keys()):
-                    keys.append(key)
-                for physics in item._physics:
-                    for key in list(physics.keys()):
-                        if key not in keys:
-                            keys.append(key)
-        keys.sort()
-        for key in keys:
-            if key not in exclusion_list:
-                try:
-                    if _sp.amax(item[key]) < _sp.inf:
-                        element = key.split('.')[0]
-                        propname = key.split('.')[1]
-                        dict_name = element + '.' + item.name + dlim + propname
-                        if key in ['pore.coords', 'throat.conns',
-                                   'pore.all', 'throat.all']:
-                            dict_name = key
-                        data_amalgamated.update({dict_name: item[key]})
-                except TypeError:
-                    pass
-    return data_amalgamated
-
-
 def conduit_lengths(network, throats=None, mode='pore'):
     r"""
     Return the respective lengths of the conduit components defined by the throat


Test patch:
diff --git a/tests/unit/models/physics/GenericSourceTermTest.py b/tests/unit/models/physics/GenericSourceTermTest.py
index fdb0bf0c0..6a9cee6d7 100644
--- a/tests/unit/models/physics/GenericSourceTermTest.py
+++ b/tests/unit/models/physics/GenericSourceTermTest.py
@@ -1,8 +1,6 @@
 import openpnm as op
 import numpy as np
 import openpnm.models.physics as pm
-import sympy as syp
-import collections
 
 
 class GenericSourceTermTest:
@@ -28,13 +26,7 @@ class GenericSourceTermTest:
                             model=pm.generic_source_term.linear,
                             A1='pore.item1',
                             A2='pore.item2',
-                            X='pore.mole_fraction',
-                            regen_mode='normal')
-        self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.linear_sym,
-                            A1='pore.item1',
-                            A2='pore.item2',
-                            X='pore.mole_fraction',
+                            quantity='pore.mole_fraction',
                             regen_mode='normal')
         self.alg = op.algorithms.ReactiveTransport(network=self.net,
                                                    phase=self.phase)
@@ -46,13 +38,12 @@ class GenericSourceTermTest:
         self.alg.run()
         self.phase.update(self.alg.results())
         self.phys.regenerate_models(propnames='pore.source1')
-        self.phys.regenerate_models(propnames='pore.source2')
         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)
-        rs = np.round(np.sum(self.phys['pore.source2.rate'][self.source_pores]), 20)
+#        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
         assert r1 == r2
-        assert r1 == rs
+#        assert r2 == -r3
 
     def test_power_law(self):
         self.phys['pore.item1'] = 0.5e-12
@@ -63,14 +54,7 @@ class GenericSourceTermTest:
                             A1='pore.item1',
                             A2='pore.item2',
                             A3='pore.item3',
-                            X='pore.mole_fraction',
-                            regen_mode='normal')
-        self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.power_law_sym,
-                            A1='pore.item1',
-                            A2='pore.item2',
-                            A3='pore.item3',
-                            X='pore.mole_fraction',
+                            quantity='pore.mole_fraction',
                             regen_mode='normal')
         self.alg = op.algorithms.ReactiveTransport(network=self.net,
                                                    phase=self.phase)
@@ -82,13 +66,12 @@ class GenericSourceTermTest:
         self.alg.run()
         self.phase.update(self.alg.results())
         self.phys.regenerate_models(propnames='pore.source1')
-        self.phys.regenerate_models(propnames='pore.source2')
         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])
-        rs = np.sum(self.phys['pore.source2.rate'][self.source_pores])
+#        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
         assert r1 == r2
-        assert r1 == rs
+#        assert r2 == -r3
 
     def test_exponential(self):
         self.phys['pore.item1'] = 0.8e-11
@@ -105,17 +88,7 @@ class GenericSourceTermTest:
                             A4='pore.item4',
                             A5='pore.item5',
                             A6='pore.item6',
-                            X='pore.mole_fraction',
-                            regen_mode='normal')
-        self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.exponential_sym,
-                            A1='pore.item1',
-                            A2='pore.item2',
-                            A3='pore.item3',
-                            A4='pore.item4',
-                            A5='pore.item5',
-                            A6='pore.item6',
-                            X='pore.mole_fraction',
+                            quantity='pore.mole_fraction',
                             regen_mode='normal')
         self.alg = op.algorithms.ReactiveTransport(network=self.net,
                                                    phase=self.phase)
@@ -127,14 +100,13 @@ class GenericSourceTermTest:
         self.alg.run()
         self.phase.update(self.alg.results())
         self.phys.regenerate_models(propnames='pore.source1')
-        self.phys.regenerate_models(propnames='pore.source2')
         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)
-        rs = np.round(np.sum(self.phys['pore.source2.rate'][self.source_pores]), 20)
+#        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
         assert r1 == r2
-        assert r1== rs
+#        assert r2 == -r3
 
     def test_natural_exponential(self):
         self.phys['pore.item1'] = 0.8e-11
@@ -149,16 +121,7 @@ class GenericSourceTermTest:
                             A3='pore.item3',
                             A4='pore.item4',
                             A5='pore.item5',
-                            X='pore.mole_fraction',
-                            regen_mode='normal')
-        self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.natural_exponential_sym,
-                            A1='pore.item1',
-                            A2='pore.item2',
-                            A3='pore.item3',
-                            A4='pore.item4',
-                            A5='pore.item5',
-                            X='pore.mole_fraction',
+                            quantity='pore.mole_fraction',
                             regen_mode='normal')
         self.alg = op.algorithms.ReactiveTransport(network=self.net,
                                                    phase=self.phase)
@@ -170,14 +133,13 @@ class GenericSourceTermTest:
         self.alg.run()
         self.phase.update(self.alg.results())
         self.phys.regenerate_models(propnames='pore.source1')
-        self.phys.regenerate_models(propnames='pore.source2')
         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)
-        rs = 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 == rs
+#        assert r2 == -r3
 
     def test_logarithm(self):
         self.phys['pore.item1'] = 0.16e-13
@@ -194,17 +156,7 @@ class GenericSourceTermTest:
                             A4='pore.item4',
                             A5='pore.item5',
                             A6='pore.item6',
-                            X='pore.mole_fraction',
-                            regen_mode='normal')
-        self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.logarithm_sym,
-                            A1='pore.item1',
-                            A2='pore.item2',
-                            A3='pore.item3',
-                            A4='pore.item4',
-                            A5='pore.item5',
-                            A6='pore.item6',
-                            X='pore.mole_fraction',
+                            quantity='pore.mole_fraction',
                             regen_mode='normal')
         self.alg = op.algorithms.ReactiveTransport(network=self.net,
                                                    phase=self.phase)
@@ -216,14 +168,13 @@ class GenericSourceTermTest:
         self.alg.run()
         self.phase.update(self.alg.results())
         self.phys.regenerate_models(propnames='pore.source1')
-        self.phys.regenerate_models(propnames='pore.source2')
         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)
-        rs = np.round(np.sum(self.phys['pore.source2.rate'][self.source_pores]), 20)
+#        r3 = np.round(self.alg.rate(pores=self.S_pores)[0], 20)
         assert r1 == r2
-        assert r1 == rs
+#        assert r2 == -r3
 
     def test_natural_logarithm(self):
         self.phys['pore.item1'] = 0.16e-14
@@ -238,17 +189,8 @@ class GenericSourceTermTest:
                             A3='pore.item3',
                             A4='pore.item4',
                             A5='pore.item5',
-                            X='pore.mole_fraction',
+                            quantity='pore.mole_fraction',
                             regen_mode='on_demand')
-        self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.natural_logarithm_sym,
-                            A1='pore.item1',
-                            A2='pore.item2',
-                            A3='pore.item3',
-                            A4='pore.item4',
-                            A5='pore.item5',
-                            X='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)
@@ -259,47 +201,14 @@ class GenericSourceTermTest:
         self.alg.run()
         self.phase.update(self.alg.results())
         self.phys.regenerate_models(propnames='pore.source1')
-        self.phys.regenerate_models(propnames='pore.source2')
         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)
-        rs = 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 == rs
+#        assert r2 == -r3
 
-    def test_general_symbolic(self):
-        a, b, c, d, e, x = syp.symbols('a,b,c,d,e,x')
-        # natural log function
-        y =  a*syp.ln(b*x**c + d)+e
-        phys = self.phys
-        phys['pore.item1'] = 0.16e-14
-        phys['pore.item2'] = 4
-        phys['pore.item3'] = 1.4
-        phys['pore.item4'] = 0.133
-        phys['pore.item5'] = -5.1e-14
-        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='pore.mole_fraction',
-                       regen_mode='normal')
-        arg_map=collections.OrderedDict([('a', 'pore.item1'),
-                                         ('b', 'pore.item2'),
-                                         ('c', 'pore.item3'),
-                                         ('d', 'pore.item4'),
-                                         ('e', 'pore.item5'),
-                                         ('x', 'pore.mole_fraction')])
-        phys.add_model(propname='pore.general',
-                       model=op.models.physics.generic_source_term.general_symbolic,
-                       eqn=y, arg_map=arg_map,
-                       regen_mode='normal')
-        assert np.allclose(phys['pore.source1.rate'], phys['pore.general.rate'])
-        assert np.allclose(phys['pore.source1.S1'], phys['pore.general.S1'])
-        assert np.allclose(phys['pore.source1.S2'], phys['pore.general.S2'])
 
 if __name__ == '__main__':
 
diff --git a/tests/unit/network/CubicDualTest.py b/tests/unit/network/CubicDualTest.py
new file mode 100644
index 000000000..e27dfb6df
--- /dev/null
+++ b/tests/unit/network/CubicDualTest.py
@@ -0,0 +1,101 @@
+import openpnm as op
+
+
+class CubicDualTest:
+    def setup_class(self):
+        pass
+
+    def teardown_class(self):
+        pass
+
+    def test_generation_3D(self):
+        net = op.network.CubicDual(shape=[5, 5, 5], label_1='primary',
+                                   label_2='secondary')
+        assert net.Np == 285
+        assert net.Nt == 1436
+        assert net.num_pores('all') == 285
+        assert net.num_pores('back') == 41
+        assert net.num_pores('bottom') == 41
+        assert net.num_pores('front') == 41
+        assert net.num_pores('internal') == 91
+        assert net.num_pores('left') == 41
+        assert net.num_pores('primary') == 125
+        assert net.num_pores('right') == 41
+        assert net.num_pores('secondary') == 160
+        assert net.num_pores('surface') == 194
+        assert net.num_pores('top') == 41
+        assert net.num_throats('all') == 1436
+        assert net.num_throats('interconnect') == 896
+        assert net.num_throats('internal') == 860
+        assert net.num_throats('primary') == 300
+        assert net.num_throats('secondary') == 240
+        assert net.num_throats('surface') == 576
+        assert net.num_throats('top') == 104
+        assert net.num_throats('bottom') == 104
+        assert net.num_throats('left') == 104
+        assert net.num_throats('right') == 104
+        assert net.num_throats('front') == 104
+        assert net.num_throats('back') == 104
+
+    def test_generation_2D_XY(self):
+        net = op.network.CubicDual(shape=[5, 5, 1], label_1='primary',
+                                   label_2='secondary')
+        assert net.Np == 57
+        assert net.Nt == 176
+        assert net.num_pores('left') == 9
+        assert net.num_pores('right') == 9
+        assert net.num_pores('front') == 9
+        assert net.num_pores('back') == 9
+        assert net.num_pores('top') == 41
+        assert net.num_pores('bottom') == 41
+        assert net.num_throats('interconnect') == 96
+
+    def test_generation_2D_XZ(self):
+        net = op.network.CubicDual(shape=[5, 1, 5], label_1='primary',
+                                   label_2='secondary')
+        assert net.Np == 57
+        assert net.Nt == 176
+        assert net.num_pores('left') == 41
+        assert net.num_pores('right') == 41
+        assert net.num_pores('front') == 9
+        assert net.num_pores('back') == 9
+        assert net.num_pores('top') == 9
+        assert net.num_pores('bottom') == 9
+        assert net.num_throats('interconnect') == 96
+
+    def test_generation_2D_YZ(self):
+        net = op.network.CubicDual(shape=[1, 5, 5], label_1='primary',
+                                   label_2='secondary')
+        assert net.Np == 57
+        assert net.Nt == 176
+        assert net.num_pores('left') == 9
+        assert net.num_pores('right') == 9
+        assert net.num_pores('front') == 41
+        assert net.num_pores('back') == 41
+        assert net.num_pores('top') == 9
+        assert net.num_pores('bottom') == 9
+        assert net.num_throats('interconnect') == 96
+
+    def test_generation_2D_2_dims(self):
+        net = op.network.CubicDual(shape=[5, 5], label_1='primary',
+                                   label_2='secondary')
+        assert net.Np == 57
+        assert net.Nt == 176
+        assert net.num_pores('left') == 9
+        assert net.num_pores('right') == 9
+        assert net.num_pores('front') == 9
+        assert net.num_pores('back') == 9
+        assert net.num_pores('top') == 41
+        assert net.num_pores('bottom') == 41
+        assert net.num_throats('interconnect') == 96
+
+
+if __name__ == '__main__':
+
+    t = CubicDualTest()
+    t.setup_class()
+    self = t
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/network/CubicTemplateTest.py b/tests/unit/network/CubicTemplateTest.py
new file mode 100644
index 000000000..30ba4892c
--- /dev/null
+++ b/tests/unit/network/CubicTemplateTest.py
@@ -0,0 +1,31 @@
+import openpnm as op
+from skimage.morphology import ball, disk
+
+
+class CubicTemplateTest:
+    def setup_class(self):
+        pass
+
+    def teardown_class(self):
+        pass
+
+    def test_2D_template(self):
+        net = op.network.CubicTemplate(template=disk(10), spacing=1)
+        assert net.Np == 317
+        assert net.Nt == 592
+
+    def test_3D_template(self):
+        net = op.network.CubicTemplate(template=ball(5), spacing=1)
+        assert net.Np == 515
+        assert net.Nt == 1302
+
+
+if __name__ == '__main__':
+
+    t = CubicTemplateTest()
+    t.setup_class()
+    self = t
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()
diff --git a/tests/unit/utils/UtilsTest.py b/tests/unit/utils/UtilsTest.py
new file mode 100644
index 000000000..9b56b0f58
--- /dev/null
+++ b/tests/unit/utils/UtilsTest.py
@@ -0,0 +1,57 @@
+import openpnm as op
+import scipy as sp
+import pytest
+
+
+class UtilsTest:
+
+    def setup_class(self):
+        self.net = op.network.Cubic(shape=[3, 3, 3])
+        self.geo = op.geometry.StickAndBall(network=self.net,
+                                            pores=self.net.Ps,
+                                            throats=self.net.Ts)
+
+    def teardown_class(self):
+        ws = op.core.Workspace()
+        ws.clear()
+
+    def test_tic_toc(self):
+        with pytest.raises(Exception):
+            op.utils.toc()
+        op.utils.tic()
+        t1 = op.utils.toc()
+        assert t1 is None
+        t2 = op.utils.toc(quiet=True)
+        assert t2 >= 0
+
+    def test_nested_dict(self):
+        d = op.utils.NestedDict()
+        d['top']['middle']['bottom'] = 1
+        assert d == {'top': {'middle': {'bottom': 1}}}
+        s = d.__str__()
+        assert s == '-top\n--middle\n---bottom\n'
+        a = d.to_dict()
+        assert type(a) is dict
+
+    def test_printable_list(self):
+        L = op.utils.PrintableList(['item1', 'item2', 'item2'])
+        s = L.__str__().split('\n')
+        assert len(s) == 5
+
+    def test_printable_dict(self):
+        D = op.utils.PrintableDict({'item1': 1, 'item2': 2, 'item3': sp.array([1, 2])})
+        s = D.__str__().split('\n')
+        assert len(s) == 7
+        r = D.__repr__()
+        assert r == "{'item1': 1, 'item2': 2, 'item3': array([1, 2])}"
+
+
+if __name__ == '__main__':
+
+    t = UtilsTest()
+    self = t
+    t.setup_class()
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print('running test: '+item)
+            t.__getattribute__(item)()