Issue 664 casadi

.requirements-docs.txt

@@ -5,5 +5,6 @@ pandas>=0.23
 anytree>=2.4.3
 autograd>=1.2
 scikit-fem>=0.2.0
+casadi>=3.5.0
 guzzle-sphinx-theme
 sphinx>=1.5

CHANGELOG.md

@@ -2,10 +2,12 @@
 
 ## Features
 
--   Added interface (via pybind11) to sundials with the IDA KLU sparse linear solver ([#657](https://github.com/pybamm-team/PyBaMM/pull/657))
+-   Add interface to CasADi solver ([#687](https://github.com/pybamm-team/PyBaMM/pull/687))
+-   Add option to use CasADi's Algorithmic Differentiation framework to calculate Jacobians ([#687](https://github.com/pybamm-team/PyBaMM/pull/687))
 -   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))
+-   Added interface (via pybind11) to sundials with the IDA KLU sparse linear solver ([#657](https://github.com/pybamm-team/PyBaMM/pull/657))
 -   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))
@@ -17,6 +19,7 @@
 
 ## Bug fixes
 
+-   Fix differentiation of functions that have more than one argument ([#687](https://github.com/pybamm-team/PyBaMM/pull/687))
 -   Add warning if `ProcessedVariable` is called outisde its interpolation range ([#681](https://github.com/pybamm-team/PyBaMM/pull/681))
 -   Improve the way `ProcessedVariable` objects are created in higher dimensions ([#581](https://github.com/pybamm-team/PyBaMM/pull/581))
 

docs/source/expression_tree/index.rst

@@ -17,6 +17,4 @@ Expression Tree
   broadcasts
   functions
   interpolant
-  evaluate
-  simplify
-  jacobian
+  operations/index

docs/source/expression_tree/operations/convert_to_casadi.rst

@@ -0,0 +1,5 @@
+Convert to CasADi
+=================
+
+.. autoclass:: pybamm.CasadiConverter
+  :members:

docs/source/expression_tree/operations/index.rst

@@ -0,0 +1,11 @@
+Operations on expression trees
+==============================
+
+Classes and functions that operate on the expression tree
+
+.. toctree::
+
+  simplify
+  evaluate
+  jacobian
+  convert_to_casadi

docs/source/solvers/casadi_solver.rst

@@ -0,0 +1,5 @@
+Casadi Solver
+=============
+
+.. autoclass:: pybamm.CasadiSolver
+  :members:

docs/source/solvers/index.rst

@@ -7,4 +7,5 @@ Solvers
   base_solvers
   scipy_solver
   scikits_solvers
+  casadi_solver
   solution

examples/scripts/compare-dae-solver.py

@@ -16,23 +16,22 @@
 
 # set mesh
 var = pybamm.standard_spatial_vars
-
-var_pts = {var.x_n: 60, var.x_s: 100, var.x_p: 60, var.r_n: 50, var.r_p: 50}
-# var_pts = model.default_var_pts
+var_pts = {var.x_n: 50, var.x_s: 50, var.x_p: 50, var.r_n: 20, var.r_p: 20}
 mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
 
 # discretise model
 disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
 disc.process_model(model)
 
 # solve model
-t_eval = np.linspace(0, 0.2, 100)
+t_eval = np.linspace(0, 0.17, 100)
 
-klu_sol = pybamm.IDAKLU(atol=1e-8, rtol=1e-8).solve(model, t_eval)
+casadi_sol = pybamm.CasadiSolver(atol=1e-8, rtol=1e-8).solve(model, t_eval)
+klu_sol = pybamm.IDAKLUSolver(atol=1e-8, rtol=1e-8).solve(model, t_eval)
 scikits_sol = pybamm.ScikitsDaeSolver(atol=1e-8, rtol=1e-8).solve(model, t_eval)
 
 # plot
-models = [model, model]
-solutions = [scikits_sol, klu_sol]
+models = [model, model, model]
+solutions = [scikits_sol, klu_sol, casadi_sol]
 plot = pybamm.QuickPlot(models, mesh, solutions)
 plot.dynamic_plot()

examples/scripts/compare_lead_acid.py

@@ -41,16 +41,18 @@
 
 # solve model
 solutions = [None] * len(models)
-t_eval = np.linspace(0, 3, 1000)
+t_eval = np.linspace(0, 1, 1000)
 for i, model in enumerate(models):
     solution = model.default_solver.solve(model, t_eval)
     solutions[i] = solution
 
 # plot
 output_variables = [
-    "Electrolyte pressure",
-    "Electrolyte concentration",
-    "Volume-averaged velocity [m.s-1]",
+    "Interfacial current density [A.m-2]",
+    "Electrolyte concentration [mol.m-3]",
+    "Current [A]",
+    "Porosity",
+    "Electrolyte potential [V]",
     "Terminal voltage [V]",
 ]
 plot = pybamm.QuickPlot(models, mesh, solutions, output_variables)

examples/scripts/compare_lead_acid_3D.py

@@ -21,8 +21,7 @@
     #     {"current collector": "potential pair", "dimensionality": 2}, name="2+1D LOQS"
     # ),
     pybamm.lead_acid.Full(
-        {"current collector": "potential pair", "dimensionality": 1},
-        name="1+1D Full",
+        {"current collector": "potential pair", "dimensionality": 1}, name="1+1D Full"
     ),
     # pybamm.lead_acid.Full(
     #     {"dimensionality": 1}, name="1+1D uniform Full"

examples/scripts/compare_lithium_ion_3D.py

@@ -18,12 +18,10 @@
 # load models
 models = [
     pybamm.lithium_ion.SPM(
-        {"current collector": "potential pair", "dimensionality": 2},
-        name="2+1D SPM",
+        {"current collector": "potential pair", "dimensionality": 2}, name="2+1D SPM"
     ),
     pybamm.lithium_ion.SPMe(
-        {"current collector": "potential pair", "dimensionality": 2},
-        name="2+1D SPMe",
+        {"current collector": "potential pair", "dimensionality": 2}, name="2+1D SPMe"
     ),
 ]
 

pybamm/__init__.py

@@ -150,19 +150,22 @@ def version(formatted=False):
     UndefinedOperationError,
     GeometryError,
 )
-from .expression_tree.simplify import (
+
+# Operations
+from .expression_tree.operations.simplify import (
     Simplification,
     simplify_if_constant,
     simplify_addition_subtraction,
     simplify_multiplication_division,
 )
-from .expression_tree.jacobian import Jacobian
-from .expression_tree.evaluate import (
+from .expression_tree.operations.evaluate import (
     find_symbols,
     id_to_python_variable,
     to_python,
     EvaluatorPython,
 )
+from .expression_tree.operations.jacobian import Jacobian
+from .expression_tree.operations.convert_to_casadi import CasadiConverter
 
 #
 # Model classes
@@ -258,12 +261,12 @@ def version(formatted=False):
 from .solvers.base_solver import BaseSolver
 from .solvers.ode_solver import OdeSolver
 from .solvers.dae_solver import DaeSolver
-from .solvers.scipy_solver import ScipySolver
-from .solvers.scikits_dae_solver import ScikitsDaeSolver
-from .solvers.scikits_ode_solver import ScikitsOdeSolver
-from .solvers.scikits_ode_solver import have_scikits_odes
 from .solvers.algebraic_solver import AlgebraicSolver
-from .solvers.idaklu_solver import IDAKLU, have_idaklu
+from .solvers.casadi_solver import CasadiSolver
+from .solvers.scikits_dae_solver import ScikitsDaeSolver
+from .solvers.scikits_ode_solver import ScikitsOdeSolver, have_scikits_odes
+from .solvers.scipy_solver import ScipySolver
+from .solvers.idaklu_solver import IDAKLUSolver, have_idaklu
 
 
 #

pybamm/discretisations/discretisation.py

@@ -128,7 +128,7 @@ def process_model(self, model, inplace=True):
             model_disc.options = model.options
             model_disc.use_jacobian = model.use_jacobian
             model_disc.use_simplify = model.use_simplify
-            model_disc.use_to_python = model.use_to_python
+            model_disc.convert_to_format = model.convert_to_format
 
         model_disc.bcs = self.bcs
 

pybamm/expression_tree/functions.py

@@ -21,10 +21,19 @@ class Function(pybamm.Symbol):
     derivative : str, optional
         Which derivative to use when differentiating ("autograd" or "derivative").
         Default is "autograd".
+    differentiated_function : method, optional
+        The function which was differentiated to obtain this one. Default is None.
     **Extends:** :class:`pybamm.Symbol`
     """
 
-    def __init__(self, function, *children, name=None, derivative="autograd"):
+    def __init__(
+        self,
+        function,
+        *children,
+        name=None,
+        derivative="autograd",
+        differentiated_function=None
+    ):
 
         if name is not None:
             self.name = name
@@ -39,6 +48,7 @@ def __init__(self, function, *children, name=None, derivative="autograd"):
 
         self.function = function
         self.derivative = derivative
+        self.differentiated_function = differentiated_function
 
         # hack to work out whether function takes any params
         # (signature doesn't work for numpy)
@@ -85,7 +95,9 @@ def diff(self, variable):
                 # if variable appears in the function,use autograd to differentiate
                 # function, and apply chain rule
                 if variable.id in [symbol.id for symbol in child.pre_order()]:
-                    partial_derivatives[i] = self._diff(children) * child.diff(variable)
+                    partial_derivatives[i] = self._function_diff(
+                        children, i
+                    ) * child.diff(variable)
 
             # remove None entries
             partial_derivatives = list(filter(None, partial_derivatives))
@@ -96,15 +108,34 @@ def diff(self, variable):
 
             return derivative
 
-    def _diff(self, children):
-        """ See :meth:`pybamm.Symbol._diff()`. """
+    def _function_diff(self, children, idx):
+        """
+        Derivative with respect to child number 'idx'.
+        See :meth:`pybamm.Symbol._diff()`.
+        """
+        # Store differentiated function, needed in case we want to convert to CasADi
         if self.derivative == "autograd":
-            return Function(autograd.elementwise_grad(self.function), *children)
-        elif self.derivative == "derivative":
-            # keep using "derivative" as derivative
-            return pybamm.Function(
-                self.function.derivative(), *children, derivative="derivative"
+            return Function(
+                autograd.elementwise_grad(self.function, idx),
+                *children,
+                differentiated_function=self.function
             )
+        elif self.derivative == "derivative":
+            if len(children) > 1:
+                raise ValueError(
+                    """
+                    differentiation using '.derivative()' not implemented for functions
+                    with more than one child
+                    """
+                )
+            else:
+                # keep using "derivative" as derivative
+                return pybamm.Function(
+                    self.function.derivative(),
+                    *children,
+                    derivative="derivative",
+                    differentiated_function=self.function
+                )
 
     def _function_jac(self, children_jacs):
         """ Calculate the jacobian of a function. """
@@ -118,7 +149,7 @@ def _function_jac(self, children_jacs):
             children = self.orphans
             for i, child in enumerate(children):
                 if not child.evaluates_to_number():
-                    jac_fun = self._diff(children) * children_jacs[i]
+                    jac_fun = self._function_diff(children, i) * children_jacs[i]
                     jac_fun.domain = []
                     if jacobian is None:
                         jacobian = jac_fun
@@ -175,7 +206,11 @@ def _function_new_copy(self, children):
             A new copy of the function
         """
         return pybamm.Function(
-            self.function, *children, name=self.name, derivative=self.derivative
+            self.function,
+            *children,
+            name=self.name,
+            derivative=self.derivative,
+            differentiated_function=self.differentiated_function
         )
 
     def _function_simplify(self, simplified_children):
@@ -203,7 +238,8 @@ def _function_simplify(self, simplified_children):
                 self.function,
                 *simplified_children,
                 name=self.name,
-                derivative=self.derivative
+                derivative=self.derivative,
+                differentiated_function=self.differentiated_function
             )
 
 
@@ -239,8 +275,8 @@ class Cos(SpecificFunction):
     def __init__(self, child):
         super().__init__(np.cos, child)
 
-    def _diff(self, children):
-        """ See :meth:`pybamm.Symbol._diff()`. """
+    def _function_diff(self, children, idx):
+        """ See :meth:`pybamm.Symbol._function_diff()`. """
         return -Sin(children[0])
 
 
@@ -255,8 +291,8 @@ class Cosh(SpecificFunction):
     def __init__(self, child):
         super().__init__(np.cosh, child)
 
-    def _diff(self, children):
-        """ See :meth:`pybamm.Function._diff()`. """
+    def _function_diff(self, children, idx):
+        """ See :meth:`pybamm.Function._function_diff()`. """
         return Sinh(children[0])
 
 
@@ -271,8 +307,8 @@ class Exponential(SpecificFunction):
     def __init__(self, child):
         super().__init__(np.exp, child)
 
-    def _diff(self, children):
-        """ See :meth:`pybamm.Function._diff()`. """
+    def _function_diff(self, children, idx):
+        """ See :meth:`pybamm.Function._function_diff()`. """
         return Exponential(children[0])
 
 
@@ -287,8 +323,8 @@ class Log(SpecificFunction):
     def __init__(self, child):
         super().__init__(np.log, child)
 
-    def _diff(self, children):
-        """ See :meth:`pybamm.Function._diff()`. """
+    def _function_diff(self, children, idx):
+        """ See :meth:`pybamm.Function._function_diff()`. """
         return 1 / children[0]
 
 
@@ -313,8 +349,8 @@ class Sin(SpecificFunction):
     def __init__(self, child):
         super().__init__(np.sin, child)
 
-    def _diff(self, children):
-        """ See :meth:`pybamm.Function._diff()`. """
+    def _function_diff(self, children, idx):
+        """ See :meth:`pybamm.Function._function_diff()`. """
         return Cos(children[0])
 
 
@@ -329,8 +365,8 @@ class Sinh(SpecificFunction):
     def __init__(self, child):
         super().__init__(np.sinh, child)
 
-    def _diff(self, children):
-        """ See :meth:`pybamm.Function._diff()`. """
+    def _function_diff(self, children, idx):
+        """ See :meth:`pybamm.Function._function_diff()`. """
         return Cosh(children[0])