Issue 665 jac class

.gitignore

@@ -13,6 +13,7 @@
 *.pyc
 *.py~
 .~lock*
+*.swp
 
 # Python cache
 *__pycache__/
@@ -68,3 +69,11 @@ pyvenv.cfg
 sundials
 sundials4
 SuiteSparse
+SuiteSparse/.gitignore
+
+# cmakefiles
+CMakeFiles
+Makefile
+*.cmake
+*.so
+

CHANGELOG.md

@@ -2,20 +2,21 @@
 
 ## Features
 
-- Add method to evaluate parameters more easily ([#669](https://github.com/pybamm-team/PyBaMM/pull/669))
-- Add `Interpolant` class to interpolate experimental data (e.g. OCP curves) ([#661](https://github.com/pybamm-team/PyBaMM/pull/661))
-- Allow parameters to be set by material or by specifying a particular paper ([#647](https://github.com/pybamm-team/PyBaMM/pull/647))
-- Set relative and absolute tolerances independently in solvers ([#645](https://github.com/pybamm-team/PyBaMM/pull/645))
-- Add some non-uniform meshes in 1D and 2D ([#617](https://github.com/pybamm-team/PyBaMM/pull/617))
+-   Add method to evaluate parameters more easily ([#669](https://github.com/pybamm-team/PyBaMM/pull/669))
+-   Add `Jacobian` class to reuse known Jacobians of expressions ([#665](https://github.com/pybamm-team/PyBaMM/pull/670))
+-   Add `Interpolant` class to interpolate experimental data (e.g. OCP curves) ([#661](https://github.com/pybamm-team/PyBaMM/pull/661))
+-   Allow parameters to be set by material or by specifying a particular paper ([#647](https://github.com/pybamm-team/PyBaMM/pull/647))
+-   Set relative and absolute tolerances independently in solvers ([#645](https://github.com/pybamm-team/PyBaMM/pull/645))
+-   Add some non-uniform meshes in 1D and 2D ([#617](https://github.com/pybamm-team/PyBaMM/pull/617))
 
 ## Optimizations
 
-- Avoid re-checking size when making a copy of an `Index` object ([#656](https://github.com/pybamm-team/PyBaMM/pull/656))
-- Avoid recalculating `_evaluation_array` when making a copy of a `StateVector` object ([#653](https://github.com/pybamm-team/PyBaMM/pull/653))
+-   Avoid re-checking size when making a copy of an `Index` object ([#656](https://github.com/pybamm-team/PyBaMM/pull/656))
+-   Avoid recalculating `_evaluation_array` when making a copy of a `StateVector` object ([#653](https://github.com/pybamm-team/PyBaMM/pull/653))
 
 ## Bug fixes
 
-- Improve the way `ProcessedVariable` objects are created in higher dimensions ([#581](https://github.com/pybamm-team/PyBaMM/pull/581))
+-   Improve the way `ProcessedVariable` objects are created in higher dimensions ([#581](https://github.com/pybamm-team/PyBaMM/pull/581))
 
 # [v0.1.0](https://github.com/pybamm-team/PyBaMM/tree/v0.1.0) - 2019-10-08
 

docs/source/expression_tree/index.rst

@@ -15,7 +15,8 @@ Expression Tree
   unary_operator
   concatenations
   broadcasts
-  simplify
   functions
   interpolant
   evaluate
+  simplify
+  jacobian

docs/source/expression_tree/jacobian.rst

@@ -0,0 +1,5 @@
+Jacobian
+========
+
+.. autoclass:: pybamm.Jacobian
+  :members:

examples/notebooks/change-input-current.ipynb

@@ -372,7 +372,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.7"
+   "version": "3.6.8"
   }
  },
  "nbformat": 4,

examples/notebooks/change-settings.ipynb

@@ -508,7 +508,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.7"
+   "version": "3.6.8"
   }
  },
  "nbformat": 4,

examples/notebooks/expression_tree/expression-tree.ipynb

@@ -220,7 +220,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.7"
+   "version": "3.6.8"
   }
  },
  "nbformat": 4,

examples/notebooks/models/compare-lithium-ion.ipynb

@@ -463,7 +463,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.5.7"
+   "version": "3.6.8"
   }
  },
  "nbformat": 4,

examples/notebooks/parameter-values.ipynb

@@ -482,7 +482,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.7"
+   "version": "3.6.8"
   }
  },
  "nbformat": 4,

pybamm/__init__.py

@@ -156,6 +156,7 @@ def version(formatted=False):
     simplify_addition_subtraction,
     simplify_multiplication_division,
 )
+from .expression_tree.jacobian import Jacobian
 from .expression_tree.evaluate import (
     find_symbols,
     id_to_python_variable,

pybamm/discretisations/discretisation.py

@@ -93,6 +93,7 @@ def process_model(self, model, inplace=True):
             If an empty model is passed (`model.rhs = {}` and `model.algebraic={}`)
 
         """
+
         pybamm.logger.info("Start discretising {}".format(model.name))
 
         # Make sure model isn't empty
@@ -424,6 +425,9 @@ def process_rhs_and_algebraic(self, model):
         # Discretise and concatenate algebraic equations
         processed_algebraic = self.process_dict(model.algebraic)
 
+        # Concatenate algebraic into a single state vector
+        # Need to concatenate in order as the ordering of equations could be different
+        # in processed_algebraic and model.algebraic
         processed_concatenated_algebraic = self._concatenate_in_order(
             processed_algebraic
         )
@@ -497,15 +501,73 @@ def create_mass_matrix(self, model):
 
         return pybamm.Matrix(mass_matrix)
 
+    def create_jacobian(self, model):
+        """Creates Jacobian of the discretised model.
+        Note that the model is assumed to be of the form M*y_dot = f(t,y), where
+        M is the (possibly singular) mass matrix. The Jacobian is df/dy.
+
+        Note: At present, calculation of the Jacobian is deferred until after
+        simplification, since it is much faster to compute the Jacobian of the
+        simplified model. However, in some use cases (e.g. running the same
+        model multiple times but with different parameters) it may be more
+        efficient to compute the Jacobian once, before simplification, so that
+        parameters in the Jacobian can be updated (see PR #670).
+
+        Parameters
+        ----------
+        model : :class:`pybamm.BaseModel`
+            Discretised model. Must have attributes rhs, initial_conditions and
+            boundary_conditions (all dicts of {variable: equation})
+
+        Returns
+        -------
+        :class:`pybamm.Concatenation`
+            The expression trees corresponding to the Jacobian of the model
+        """
+        # create state vector to differentiate with respect to
+        y = pybamm.StateVector(slice(0, np.size(model.concatenated_initial_conditions)))
+        # set up Jacobian object, for re-use of dict
+        jacobian = pybamm.Jacobian()
+
+        # calculate Jacobian of rhs by equation
+        jac_rhs_eqn_dict = {}
+        for eqn_key, eqn in model.rhs.items():
+            pybamm.logger.debug(
+                "Calculating block of Jacobian for {!r}".format(eqn_key.name)
+            )
+            jac_rhs_eqn_dict[eqn_key] = jacobian.jac(eqn, y)
+        jac_rhs = self._concatenate_in_order(jac_rhs_eqn_dict, sparse=True)
+
+        # calculate Jacobian of algebraic by equation
+        jac_algebraic_eqn_dict = {}
+        for eqn_key, eqn in model.algebraic.items():
+            pybamm.logger.debug(
+                "Calculating block of Jacobian for {!r}".format(eqn_key.name)
+            )
+            jac_algebraic_eqn_dict[eqn_key] = jacobian.jac(eqn, y)
+        jac_algebraic = self._concatenate_in_order(jac_algebraic_eqn_dict, sparse=True)
+
+        # full Jacobian
+        if model.rhs.keys() and model.algebraic.keys():
+            jac = pybamm.SparseStack(jac_rhs, jac_algebraic)
+        elif not model.algebraic.keys():
+            jac = jac_rhs
+        else:
+            jac = jac_algebraic
+
+        return jac, jac_rhs, jac_algebraic
+
     def process_dict(self, var_eqn_dict):
         """Discretise a dictionary of {variable: equation}, broadcasting if necessary
-        (can be model.rhs, model.initial_conditions or model.variables).
+        (can be model.rhs, model.algebraic, model.initial_conditions or
+        model.variables).
 
         Parameters
         ----------
         var_eqn_dict : dict
             Equations ({variable: equation} dict) to dicretise
-            (can be model.rhs, model.initial_conditions or model.variables)
+            (can be model.rhs, model.algebraic, model.initial_conditions or
+            model.variables)
 
         Returns
         -------
@@ -679,10 +741,13 @@ def _process_symbol(self, symbol):
                     "Cannot discretise symbol of type '{}'".format(type(symbol))
                 )
 
-    def concatenate(self, *symbols):
-        return pybamm.NumpyConcatenation(*symbols)
+    def concatenate(self, *symbols, sparse=False):
+        if sparse:
+            return pybamm.SparseStack(*symbols)
+        else:
+            return pybamm.NumpyConcatenation(*symbols)
 
-    def _concatenate_in_order(self, var_eqn_dict, check_complete=False):
+    def _concatenate_in_order(self, var_eqn_dict, check_complete=False, sparse=False):
         """
         Concatenate a dictionary of {variable: equation} using self.y_slices
 
@@ -695,8 +760,15 @@ def _concatenate_in_order(self, var_eqn_dict, check_complete=False):
         ----------
         var_eqn_dict : dict
             Equations ({variable: equation} dict) to dicretise
+        check_complete : bool, optional
+            Whether to check keys in var_eqn_dict against self.y_slices. Default
+            is False
+        sparse : bool, optional
+            If True the concatenation will be a :class:`pybamm.SparseStack`. If
+            False the concatenation will be a :class:`pybamm.NumpyConcatenation`.
+            Default is False
 
-                Returns
+        Returns
         -------
         var_eqn_dict : dict
             Discretised right-hand side equations
@@ -735,7 +807,7 @@ def _concatenate_in_order(self, var_eqn_dict, check_complete=False):
         # sort equations according to slices
         sorted_equations = [eq for _, eq in sorted(zip(slices, equations))]
 
-        return self.concatenate(*sorted_equations)
+        return self.concatenate(*sorted_equations, sparse=sparse)
 
     def check_model(self, model):
         """ Perform some basic checks to make sure the discretised model makes sense."""

pybamm/expression_tree/array.py

@@ -3,7 +3,7 @@
 #
 import numpy as np
 import pybamm
-from scipy.sparse import issparse
+from scipy.sparse import issparse, csr_matrix
 
 
 class Array(pybamm.Symbol):
@@ -82,6 +82,12 @@ def set_id(self):
             (self.__class__, self.name, self.entries_string) + tuple(self.domain)
         )
 
+    def _jac(self, variable):
+        """ See :meth:`pybamm.Symbol._jac()`. """
+        # Return zeros of correct size
+        jac = csr_matrix((self.size, variable.evaluation_array.count(True)))
+        return pybamm.Matrix(jac)
+
     def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
         return self.__class__(

pybamm/expression_tree/concatenations.py

@@ -78,6 +78,10 @@ def _concatenation_new_copy(self, children):
         new_symbol = self.__class__(*children)
         return new_symbol
 
+    def _concatenation_jac(self, children_jacs):
+        """ Calculate the jacobian of a concatenation """
+        return NotImplementedError
+
     def _concatenation_simplify(self, children):
         """ See :meth:`pybamm.Symbol.simplify()`. """
         new_symbol = self.__class__(*children)
@@ -127,13 +131,13 @@ def __init__(self, *children):
             concat_fun=np.concatenate
         )
 
-    def _jac(self, variable):
-        """ See :meth:`pybamm.Symbol._jac()`. """
+    def _concatenation_jac(self, children_jacs):
+        """ See :meth:`pybamm.Concatenation.concatenation_jac()`. """
         children = self.cached_children
         if len(children) == 0:
             return pybamm.Scalar(0)
         else:
-            return SparseStack(*[child.jac(variable) for child in children])
+            return SparseStack(*children_jacs)
 
     def _concatenation_simplify(self, children):
         """ See :meth:`pybamm.Symbol.simplify()`. """
@@ -253,14 +257,13 @@ def _concatenation_evaluate(self, children_eval):
 
         return vector
 
-    def _jac(self, variable):
-        """ See :meth:`pybamm.Symbol._jac()`. """
+    def _concatenation_jac(self, children_jacs):
+        """ See :meth:`pybamm.Concatenation.concatenation_jac()`. """
         # note that this assumes that the children are in the right order and only have
         # one domain each
         jacs = []
-        child_jacs = [child.jac(variable) for child in self.cached_children]
         for i in range(self.secondary_dimensions_npts):
-            for child_jac, slices in zip(child_jacs, self._children_slices):
+            for child_jac, slices in zip(children_jacs, self._children_slices):
                 if len(slices) > 1:
                     raise NotImplementedError(
                         """jacobian only implemented for when each child has

pybamm/expression_tree/functions.py

@@ -106,22 +106,20 @@ def _diff(self, children):
                 self.function.derivative(), *children, derivative="derivative"
             )
 
-    def _jac(self, variable):
-        """ See :meth:`pybamm.Symbol._jac()`. """
+    def _function_jac(self, children_jacs):
+        """ Calculate the jacobian of a function. """
 
         if all(child.evaluates_to_number() for child in self.children):
             jacobian = pybamm.Scalar(0)
         else:
-
             # if at least one child contains variable dependence, then
             # calculate the required partial jacobians and add them
             jacobian = None
             children = self.orphans
-            for child in children:
+            for i, child in enumerate(children):
                 if not child.evaluates_to_number():
-                    jac_fun = self._diff(children) * child.jac(variable)
-
-                    jac_fun.domain = self.domain
+                    jac_fun = self._diff(children) * children_jacs[i]
+                    jac_fun.domain = []
                     if jacobian is None:
                         jacobian = jac_fun
                     else:

pybamm/expression_tree/independent_variable.py

@@ -31,6 +31,10 @@ def evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
         return pybamm.evaluate_for_shape_using_domain(self.domain)
 
+    def _jac(self, variable):
+        """ See :meth:`pybamm.Symbol._jac()`. """
+        return pybamm.Scalar(0)
+
 
 class Time(IndependentVariable):
     """A node in the expression tree representing time

pybamm/expression_tree/jacobian.py

@@ -0,0 +1,88 @@
+#
+# Calculate the Jacobian of a symbol
+#
+import pybamm
+
+
+class Jacobian(object):
+    def __init__(self, known_jacs=None):
+        self._known_jacs = known_jacs or {}
+
+    def jac(self, symbol, variable):
+        """
+        This function recurses down the tree, computing the Jacobian using
+        the Jacobians defined in classes derived from pybamm.Symbol. E.g. the
+        Jacobian of a 'pybamm.Multiplication' is computed via the product rule.
+        If the Jacobian of a symbol has already been calculated, the stored value
+        is returned.
+        Note: The Jacobian is the derivative of a symbol with respect to a (slice of)
+        a State Vector.
+
+        Parameters
+        ----------
+        symbol : :class:`pybamm.Symbol`
+            The symbol to calculate the Jacobian of
+        variable : :class:`pybamm.Symbol`
+            The variable with respect to which to differentiate
+
+        Returns
+        -------
+        :class:`pybamm.Symbol`
+            Symbol representing the Jacobian
+        """
+
+        try:
+            return self._known_jacs[symbol.id]
+        except KeyError:
+            jac = self._jac(symbol, variable)
+            self._known_jacs[symbol.id] = jac
+            return jac
+
+    def _jac(self, symbol, variable):
+        """ See :meth:`Jacobian.jac()`. """
+
+        if isinstance(symbol, pybamm.BinaryOperator):
+            left, right = symbol.children
+            # process children
+            left_jac = self.jac(left, variable)
+            right_jac = self.jac(right, variable)
+            # Need to treat outer differently. If the left child of an Outer
+            # evaluates to number then we need to return a matrix of zeros
+            # of the correct size, which requires variable.evaluation_array
+            if isinstance(symbol, pybamm.Outer):
+                # _outer_jac defined in pybamm.Outer
+                jac = symbol._outer_jac(left_jac, right_jac, variable)
+            else:
+                # _binary_jac defined in derived classes for specific rules
+                jac = symbol._binary_jac(left_jac, right_jac)
+
+        elif isinstance(symbol, pybamm.UnaryOperator):
+            child_jac = self.jac(symbol.child, variable)
+            # _unary_jac defined in derived classes for specific rules
+            jac = symbol._unary_jac(child_jac)
+
+        elif isinstance(symbol, pybamm.Function):
+            children_jacs = [None] * len(symbol.children)
+            for i, child in enumerate(symbol.children):
+                children_jacs[i] = self.jac(child, variable)
+            # _function_jac defined in function class
+            jac = symbol._function_jac(children_jacs)
+
+        elif isinstance(symbol, pybamm.Concatenation):
+            children_jacs = [child.jac(variable) for child in symbol.cached_children]
+            jac = symbol._concatenation_jac(children_jacs)
+
+        else:
+            try:
+                jac = symbol._jac(variable)
+            except NotImplementedError:
+                raise NotImplementedError(
+                    "Cannot calculate Jacobian of symbol of type '{}'".format(
+                        type(symbol)
+                    )
+                )
+
+        # jacobian removes the domain(s)
+        jac.domain = []
+        jac.auxiliary_domains = {}
+        return jac

pybamm/expression_tree/scalar.py

@@ -54,8 +54,8 @@ def _base_evaluate(self, t=None, y=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         return self._value
 
-    def jac(self, variable):
-        """ See :meth:`pybamm.Symbol.jac()`. """
+    def _jac(self, variable):
+        """ See :meth:`pybamm.Symbol._jac()`. """
         return pybamm.Scalar(0)
 
     def new_copy(self):

pybamm/expression_tree/simplify.py

@@ -593,12 +593,14 @@ def _simplify(self, symbol):
 
         elif isinstance(symbol, pybamm.UnaryOperator):
             new_child = self.simplify(symbol.child)
+            # _unary_simplify defined in derived classes for specific rules
             new_symbol = symbol._unary_simplify(new_child)
 
         elif isinstance(symbol, pybamm.Function):
             simplified_children = [None] * len(symbol.children)
             for i, child in enumerate(symbol.children):
                 simplified_children[i] = self.simplify(child)
+            # _function_simplify defined in function class
             new_symbol = symbol._function_simplify(simplified_children)
 
         elif isinstance(symbol, pybamm.Concatenation):

pybamm/expression_tree/state_vector.py

@@ -117,7 +117,7 @@ def _base_evaluate(self, t=None, y=None):
                 out = out[:, np.newaxis]
             return out
 
-    def jac(self, variable):
+    def _jac(self, variable):
         """
         Differentiate a slice of a StateVector of size m with respect to another
         slice of a StateVector of size n. This returns a (sparse) matrix of size

pybamm/expression_tree/symbol.py

@@ -441,33 +441,19 @@ def _diff(self, variable):
         "Default behaviour for differentiation, overriden by Binary and Unary Operators"
         raise NotImplementedError
 
-    def jac(self, variable):
+    def jac(self, variable, known_jacs=None):
         """
         Differentiate a symbol with respect to a (slice of) a State Vector.
-        Default behaviour is to return `1` if differentiating with respect to
-        yourself and zero otherwise. Binary and Unary Operators override this.
-
-        Parameters
-        ----------
-        variable : :class:`pybamm.Symbol`
-            The variable with respect to which to differentiate
-
+        See :class:`pybamm.Jacobian`.
         """
-        if variable.id == self.id:
-            return pybamm.Scalar(1)
-        else:
-            jac = self._jac(variable)
-            # jacobian removes the domain(s)
-            jac.domain = []
-            jac.auxiliary_domains = {}
-            return jac
+        return pybamm.Jacobian(known_jacs).jac(self, variable)
 
     def _jac(self, variable):
-        "Default behaviour for jacobian, overriden by Binary and Unary Operators"
-        if variable.id == self.id:
-            return pybamm.Scalar(1)
-        else:
-            return pybamm.Scalar(0)
+        """
+        Default behaviour for jacobian, will raise a ``NotImplementedError``
+        if this member function has not been defined for the node.
+        """
+        raise NotImplementedError
 
     def _base_evaluate(self, t=None, y=None):
         """evaluate expression tree

pybamm/expression_tree/unary_operators.py

@@ -47,6 +47,10 @@ def _unary_new_copy(self, child):
 
         return self.__class__(child)
 
+    def _unary_jac(self, child_jac):
+        """ Calculate the jacobian of a unary operator. """
+        raise NotImplementedError
+
     def _unary_simplify(self, simplified_child):
         """
         Simplify a unary operator. Default behaviour is to make a new copy, with
@@ -100,9 +104,9 @@ def _diff(self, variable):
         """ See :meth:`pybamm.Symbol._diff()`. """
         return -self.child.diff(variable)
 
-    def _jac(self, variable):
-        """ See :meth:`pybamm.Symbol._jac()`. """
-        return -self.child.jac(variable)
+    def _unary_jac(self, child_jac):
+        """ See :meth:`pybamm.UnaryOperator._unary_jac()`. """
+        return -child_jac
 
     def _unary_evaluate(self, child):
         """ See :meth:`UnaryOperator._unary_evaluate()`. """
@@ -126,8 +130,8 @@ def diff(self, variable):
             "Derivative of absolute function is not defined"
         )
 
-    def jac(self, variable):
-        """ See :meth:`pybamm.Symbol.jac()`. """
+    def _unary_jac(self, child_jac):
+        """ See :meth:`pybamm.UnaryOperator._unary_jac()`. """
         # Derivative is not well-defined
         raise pybamm.UndefinedOperationError(
             "Derivative of absolute function is not defined"
@@ -189,18 +193,17 @@ def __init__(self, child, index, name=None, check_size=True):
         if isinstance(index, int):
             self.domain = []
 
-    def _jac(self, variable):
-        """ See :meth:`pybamm.Symbol._jac()`. """
+    def _unary_jac(self, child_jac):
+        """ See :meth:`pybamm.UnaryOperator._unary_jac()`. """
 
         # if child.jac returns a matrix of zeros, this subsequently gives a bug
         # when trying to simplify the node Index(child_jac). Instead, search the
         # tree for StateVectors and return a matrix of zeros of the correct size
         # if none are found.
         if all([not (isinstance(n, pybamm.StateVector)) for n in self.pre_order()]):
-            jac = csr_matrix((1, variable.evaluation_array.count(True)))
+            jac = csr_matrix((1, child_jac.shape[1]))
             return pybamm.Matrix(jac)
         else:
-            child_jac = self.child.jac(variable)
             return Index(child_jac, self.index)
 
     def set_id(self):
@@ -262,10 +265,6 @@ def diff(self, variable):
         # We shouldn't need this
         raise NotImplementedError
 
-    def jac(self, variable):
-        """ See :meth:`pybamm.Symbol.jac()`. """
-        raise NotImplementedError
-
     def _unary_simplify(self, child):
         """ See :meth:`pybamm.UnaryOperator.simplify()`. """
 

pybamm/expression_tree/vector.py

@@ -4,7 +4,6 @@
 import pybamm
 
 import numpy as np
-from scipy.sparse import csr_matrix
 
 
 class Vector(pybamm.Array):
@@ -37,9 +36,3 @@ def __init__(
             name = "Column vector of length {!s}".format(entries.shape[0])
 
         super().__init__(entries, name, domain, auxiliary_domains, entries_string)
-
-    def _jac(self, variable):
-        """ See :meth:`pybamm.Symbol._jac()`. """
-        # Return zeros of correct size
-        jac = csr_matrix((self.size, variable.evaluation_array.count(True)))
-        return pybamm.Matrix(jac)

pybamm/models/base_model.py

@@ -52,7 +52,16 @@ class BaseModel(object):
         automatically
     jacobian : :class:`pybamm.Concatenation`
         Contains the Jacobian for the model. If model.use_jacobian is True, the
-        Jacobian is computed automatically during the set up in solve
+        Jacobian is computed automatically during solver set up
+    jacobian_rhs : :class:`pybamm.Concatenation`
+        Contains the Jacobian for the part of the model which contains time derivatives.
+        If model.use_jacobian is True, the Jacobian is computed automatically during
+        solver set up
+    jacobian_algebraic : :class:`pybamm.Concatenation`
+        Contains the Jacobian for the algebraic part of the model. This may be used
+        by the solver when calculating consistent initial conditions. If
+        model.use_jacobian is True, the Jacobian is computed automatically during
+        solver set up
     use_jacobian : bool
         Whether to use the Jacobian when solving the model (default is True)
     use_simplify : bool
@@ -82,6 +91,7 @@ def __init__(self, name="Unnamed model"):
         self._concatenated_initial_conditions = None
         self._mass_matrix = None
         self._jacobian = None
+        self._jacobian_algebraic = None
 
         # Default behaviour is to use the jacobian and simplify
         self.use_jacobian = True
@@ -226,6 +236,22 @@ def jacobian(self):
     def jacobian(self, jacobian):
         self._jacobian = jacobian
 
+    @property
+    def jacobian_rhs(self):
+        return self._jacobian_rhs
+
+    @jacobian_rhs.setter
+    def jacobian_rhs(self, jacobian_rhs):
+        self._jacobian_rhs = jacobian_rhs
+
+    @property
+    def jacobian_algebraic(self):
+        return self._jacobian_algebraic
+
+    @jacobian_algebraic.setter
+    def jacobian_algebraic(self, jacobian_algebraic):
+        self._jacobian_algebraic = jacobian_algebraic
+
     @property
     def set_of_parameters(self):
         return self._set_of_parameters

pybamm/solvers/algebraic_solver.py

@@ -189,15 +189,20 @@ def set_up(self, model):
             concatenated_algebraic = simp.simplify(concatenated_algebraic)
 
         if model.use_jacobian:
+            # Create Jacobian from concatenated algebraic
             y = pybamm.StateVector(
                 slice(0, np.size(model.concatenated_initial_conditions))
             )
+            # set up Jacobian object, for re-use of dict
+            jacobian = pybamm.Jacobian()
             pybamm.logger.info("Calculating jacobian")
-            jac = concatenated_algebraic.jac(y)
+            jac = jacobian.jac(concatenated_algebraic, y)
+            model.jacobian = jac
+            model.jacobian_algebraic = jac
 
             if model.use_simplify:
                 pybamm.logger.info("Simplifying jacobian")
-                jac = jac.simplify()
+                jac = simp.simplify(jac)
 
             if model.use_to_python:
                 pybamm.logger.info("Converting jacobian to python")

pybamm/solvers/dae_solver.py

@@ -19,7 +19,7 @@ class DaeSolver(pybamm.BaseSolver):
     root_method : str, optional
         The method to use to find initial conditions (default is "lm")
     root_tol : float, optional
-        The tolerance for the initial-condition solver (default is 1e-8).
+        The tolerance for the initial-condition solver (default is 1e-6).
     max_steps: int, optional
         The maximum number of steps the solver will take before terminating
         (default is 1000).
@@ -125,15 +125,19 @@ def set_up(self, model):
             events = {name: simp.simplify(event) for name, event in events.items()}
 
         if model.use_jacobian:
-            # Create Jacobian from simplified rhs
+            # Create Jacobian from concatenated rhs and algebraic
             y = pybamm.StateVector(
                 slice(0, np.size(model.concatenated_initial_conditions))
             )
+            # set up Jacobian object, for re-use of dict
+            jacobian = pybamm.Jacobian()
             pybamm.logger.info("Calculating jacobian")
-            jac_rhs = concatenated_rhs.jac(y)
-            jac_algebraic = concatenated_algebraic.jac(y)
+            jac_rhs = jacobian.jac(concatenated_rhs, y)
+            jac_algebraic = jacobian.jac(concatenated_algebraic, y)
             jac = pybamm.SparseStack(jac_rhs, jac_algebraic)
             model.jacobian = jac
+            model.jacobian_rhs = jac_rhs
+            model.jacobian_algebraic = jac_algebraic
 
             if model.use_simplify:
                 pybamm.logger.info("Simplifying jacobian")
@@ -212,6 +216,8 @@ def jacobian(t, y):
             jacobian = None
 
         # Add the solver attributes
+        # Note: these are the (possibly) converted to python version rhs, algebraic
+        # etc. The expression tree versions of these are attributes of the model
         self.y0 = y0
         self.rhs = rhs
         self.algebraic = algebraic

pybamm/solvers/ode_solver.py

@@ -89,11 +89,14 @@ def set_up(self, model):
         y0 = model.concatenated_initial_conditions[:, 0]
 
         if model.use_jacobian:
-            # Create Jacobian from simplified rhs
+            # Create Jacobian from concatenated rhs
             y = pybamm.StateVector(slice(0, np.size(y0)))
+            # set up Jacobian object, for re-use of dict
+            jacobian = pybamm.Jacobian()
             pybamm.logger.info("Calculating jacobian")
-            jac_rhs = concatenated_rhs.jac(y)
+            jac_rhs = jacobian.jac(concatenated_rhs, y)
             model.jacobian = jac_rhs
+            model.jacobian_rhs = jac_rhs
 
             if model.use_simplify:
                 pybamm.logger.info("Simplifying jacobian")
@@ -139,6 +142,8 @@ def jacobian(t, y):
             jacobian = None
 
         # Add the solver attributes
+        # Note: these are the (possibly) converted to python version rhs, algebraic
+        # etc. The expression tree versions of these are attributes of the model
         self.y0 = y0
         self.dydt = dydt
         self.events = events

pybamm/solvers/scikits_dae_solver.py

@@ -29,7 +29,7 @@ class ScikitsDaeSolver(pybamm.DaeSolver):
     root_method : str, optional
         The method to use to find initial conditions (default is "lm")
     root_tol : float, optional
-        The tolerance for the initial-condition solver (default is 1e-8).
+        The tolerance for the initial-condition solver (default is 1e-6).
     max_steps: int, optional
         The maximum number of steps the solver will take before terminating
         (default is 1000).

tests/unit/test_discretisations/test_discretisation.py

@@ -496,6 +496,11 @@ def test_process_model_ode(self):
         jacobian = model.concatenated_rhs.jac(y).evaluate(0, y0)
         np.testing.assert_array_equal(np.eye(np.size(y0)), jacobian.toarray())
 
+        # test jacobian by eqn gives same as jacobian of concatenated rhs
+        model.jacobian, _, _ = disc.create_jacobian(model)
+        model_jacobian = model.jacobian.evaluate(0, y0)
+        np.testing.assert_array_equal(model_jacobian.toarray(), jacobian.toarray())
+
         # test that not enough initial conditions raises an error
         model = pybamm.BaseModel()
         model.rhs = {c: pybamm.div(N), T: pybamm.div(q), S: pybamm.div(p)}
@@ -590,6 +595,11 @@ def test_process_model_dae(self):
         )
         np.testing.assert_array_equal(jacobian_actual, jacobian.toarray())
 
+        # test jacobian by eqn gives same as jacobian of concatenated rhs & algebraic
+        model.jacobian, _, _ = disc.create_jacobian(model)
+        model_jacobian = model.jacobian.evaluate(0, y0)
+        np.testing.assert_array_equal(model_jacobian.toarray(), jacobian.toarray())
+
         # test known_evals
         expr = pybamm.SparseStack(jac_rhs, jac_algebraic)
         jacobian, known_evals = expr.evaluate(0, y0, known_evals={})

tests/unit/test_expression_tree/test_binary_operators.py

@@ -102,8 +102,9 @@ def test_kron(self):
         with self.assertRaises(NotImplementedError):
             kron.diff(None)
 
+        y = pybamm.StateVector(slice(0, 2))
         with self.assertRaises(NotImplementedError):
-            kron.jac(None)
+            kron.jac(y)
 
     def test_known_eval(self):
         # Scalars

tests/unit/test_expression_tree/test_independent_variable.py

@@ -31,6 +31,8 @@ def test_time(self):
         with self.assertRaises(ValueError):
             t.evaluate(None)
 
+        self.assertEqual(t.evaluate_for_shape(), 0)
+
     def test_spatial_variable(self):
         x = pybamm.SpatialVariable("x", "negative electrode")
         self.assertEqual(x.name, "x")

tests/unit/test_expression_tree/test_jac.py

@@ -80,9 +80,8 @@ def test_linear(self):
 
         # test jac of outer if left evaluates to number
         func = pybamm.Outer(pybamm.Scalar(1), pybamm.Scalar(4))
-        jacobian = np.zeros((1, 4))
         dfunc_dy = func.jac(y).evaluate(y=y0)
-        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+        np.testing.assert_array_equal(0, dfunc_dy.toarray())
 
     def test_nonlinear(self):
         y = pybamm.StateVector(slice(0, 4))
@@ -208,28 +207,6 @@ def test_index(self):
         jac = ind.jac(vec).evaluate(y=np.linspace(0, 2, 5)).toarray()
         np.testing.assert_array_equal(jac, np.array([[0, 0, 0, 0, 0]]))
 
-    def test_jac_of_self(self):
-        "Jacobian of variable with respect to itself should be one."
-        a = pybamm.Variable("a")
-        b = pybamm.Variable("b")
-
-        self.assertEqual(a.jac(a).evaluate(), 1)
-
-        add = a + b
-        self.assertEqual(add.jac(add).evaluate(), 1)
-
-        subtract = a - b
-        self.assertEqual(subtract.jac(subtract).evaluate(), 1)
-
-        multiply = a * b
-        self.assertEqual(multiply.jac(multiply).evaluate(), 1)
-
-        divide = a / b
-        self.assertEqual(divide.jac(divide).evaluate(), 1)
-
-        power = a ** b
-        self.assertEqual(power.jac(power).evaluate(), 1)
-
     def test_jac_of_number(self):
         "Jacobian of a number should be zero"
         a = pybamm.Scalar(1)
@@ -257,15 +234,26 @@ def test_jac_of_number(self):
     def test_jac_of_symbol(self):
         a = pybamm.Symbol("a")
         y = pybamm.StateVector(slice(0, 1))
-
-        self.assertEqual(a.jac(y).evaluate(), 0)
+        with self.assertRaises(NotImplementedError):
+            a.jac(y)
 
     def test_spatial_operator(self):
         a = pybamm.Variable("a")
         b = pybamm.SpatialOperator("Operator", a)
         with self.assertRaises(NotImplementedError):
             b.jac(None)
 
+    def test_jac_of_unary_operator(self):
+        a = pybamm.Scalar(1)
+        b = pybamm.UnaryOperator("Operator", a)
+        with self.assertRaises(NotImplementedError):
+            b.jac(None)
+
+    def test_jac_of_independent_variable(self):
+        a = pybamm.IndependentVariable("Variable")
+        y = pybamm.StateVector(slice(0, 1))
+        self.assertEqual(a.jac(y).evaluate(), 0)
+
     def test_jac_of_inner(self):
         a = pybamm.Scalar(1)
         b = pybamm.Scalar(2)

tests/unit/test_expression_tree/test_simplify.py

@@ -309,9 +309,13 @@ def test_matrix_simplifications(self):
         # matrix multiplication
         A = pybamm.Matrix(np.array([[1, 0], [0, 1]]))
         self.assertIsInstance((a @ A).simplify(), pybamm.Matrix)
-        np.testing.assert_array_equal((a @ A).simplify().evaluate(), 0)
+        np.testing.assert_array_equal(
+            (a @ A).simplify().evaluate().toarray(), np.zeros((2, 2))
+        )
         self.assertIsInstance((A @ a).simplify(), pybamm.Matrix)
-        np.testing.assert_array_equal((A @ a).simplify().evaluate(), 0)
+        np.testing.assert_array_equal(
+            (A @ a).simplify().evaluate().toarray(), np.zeros((2, 2))
+        )
 
         # matrix * matrix
         m1 = pybamm.Matrix(np.array([[2, 0], [0, 2]]))
@@ -495,19 +499,26 @@ def test_matrix_simplifications(self):
             )
 
     def test_matrix_divide_simplify(self):
-        m = pybamm.Matrix(np.random.rand(300, 200))
+        m = pybamm.Matrix(np.random.rand(30, 20))
         zero = pybamm.Scalar(0)
 
         expr1 = (zero / m).simplify()
         self.assertIsInstance(expr1, pybamm.Matrix)
         self.assertEqual(expr1.shape, m.shape)
-        np.testing.assert_array_equal(expr1.evaluate(), 0)
+        np.testing.assert_array_equal(expr1.evaluate().toarray(), np.zeros((30, 20)))
 
         expr2 = (m / zero).simplify()
         self.assertIsInstance(expr2, pybamm.Matrix)
         self.assertEqual(expr2.shape, m.shape)
         np.testing.assert_array_equal(expr2.evaluate(), np.inf)
 
+        m = pybamm.Matrix(np.zeros((10, 10)))
+        a = pybamm.Scalar(7)
+        expr3 = (m / a).simplify()
+        self.assertIsInstance(expr3, pybamm.Matrix)
+        self.assertEqual(expr3.shape, m.shape)
+        np.testing.assert_array_equal(expr3.evaluate(), np.zeros((10, 10)))
+
     def test_domain_concatenation_simplify(self):
         # create discretisation
         disc = get_discretisation_for_testing()
@@ -595,11 +606,11 @@ def test_simplify_inner(self):
         a3 = pybamm.Scalar(3)
 
         np.testing.assert_array_equal(
-            pybamm.inner(a1, M2).simplify().evaluate(), M1.entries
+            pybamm.inner(a1, M2).simplify().evaluate().toarray(), M1.entries
         )
         self.assertEqual(pybamm.inner(a1, a2).simplify().evaluate(), 0)
         np.testing.assert_array_equal(
-            pybamm.inner(M2, a1).simplify().evaluate(), M1.entries
+            pybamm.inner(M2, a1).simplify().evaluate().toarray(), M1.entries
         )
         self.assertEqual(pybamm.inner(a2, a1).simplify().evaluate(), 0)
         np.testing.assert_array_equal(

tests/unit/test_solvers/test_scikits_solvers.py

@@ -471,8 +471,6 @@ def test_model_solver_ode_jacobian(self):
         def jacobian(t, y):
             return J
 
-        model.jacobian = jacobian
-
         # Solve
         solver = pybamm.ScikitsOdeSolver(rtol=1e-9, atol=1e-9)
         t_eval = np.linspace(0, 1, 100)
@@ -588,8 +586,6 @@ def jacobian(t, y):
                 ]
             )
 
-        model.jacobian = jacobian
-
         # Solve
         solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
         t_eval = np.linspace(0, 1, 100)

tests/unit/test_solvers/test_scipy_solver.py

@@ -216,8 +216,6 @@ def test_model_solver_ode_with_jacobian(self):
         def jacobian(t, y):
             return J
 
-        model.jacobian = jacobian
-
         # Solve
         solver = pybamm.ScipySolver(rtol=1e-9, atol=1e-9)
         t_eval = np.linspace(0, 1, 100)