Merge pull request #670 from pybamm-team/issue-665-jac-class

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__(