From 4a661f6d5a690cd496039ef1dd4af979acd6dbb4 Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Fri, 22 Nov 2019 22:09:12 -0500
Subject: [PATCH 01/18] add input parameter node

---
 docs/source/expression_tree/index.rst         |  1 +
 .../expression_tree/input_parameter.rst       |  5 +++
 pybamm/__init__.py                            |  1 +
 pybamm/expression_tree/binary_operators.py    | 10 ++---
 pybamm/expression_tree/concatenations.py      |  6 +--
 pybamm/expression_tree/functions.py           |  8 ++--
 pybamm/expression_tree/input_parameter.py     | 40 +++++++++++++++++++
 pybamm/expression_tree/operations/evaluate.py |  2 +-
 pybamm/expression_tree/symbol.py              |  2 +-
 pybamm/expression_tree/unary_operators.py     |  6 +--
 .../test_input_parameter.py                   | 28 +++++++++++++
 .../test_operations/test_copy.py              |  1 +
 12 files changed, 93 insertions(+), 17 deletions(-)
 create mode 100644 docs/source/expression_tree/input_parameter.rst
 create mode 100644 pybamm/expression_tree/input_parameter.py
 create mode 100644 tests/unit/test_expression_tree/test_input_parameter.py

diff --git a/docs/source/expression_tree/index.rst b/docs/source/expression_tree/index.rst
index 3eedd7b622..3efd97fb3f 100644
--- a/docs/source/expression_tree/index.rst
+++ b/docs/source/expression_tree/index.rst
@@ -16,5 +16,6 @@ Expression Tree
   concatenations
   broadcasts
   functions
+  input_parameter
   interpolant
   operations/index
diff --git a/docs/source/expression_tree/input_parameter.rst b/docs/source/expression_tree/input_parameter.rst
new file mode 100644
index 0000000000..ebbc13e7e5
--- /dev/null
+++ b/docs/source/expression_tree/input_parameter.rst
@@ -0,0 +1,5 @@
+Input Parameter
+===============
+
+.. autoclass:: pybamm.InputParameter
+  :members:
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 23f55aa76f..5c5a18353b 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -98,6 +98,7 @@ def version(formatted=False):
 from .expression_tree.unary_operators import *
 from .expression_tree.functions import *
 from .expression_tree.interpolant import Interpolant
+from .expression_tree.input_parameter import InputParameter
 from .expression_tree.parameter import Parameter, FunctionParameter
 from .expression_tree.broadcasts import Broadcast, PrimaryBroadcast, FullBroadcast
 from .expression_tree.scalar import Scalar
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index cec1dc1f05..0e85a2f3bc 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -148,21 +148,21 @@ def _binary_new_copy(self, left, right):
         "Default behaviour for new_copy"
         return self.__class__(left, right)
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             id = self.id
             try:
                 return known_evals[id], known_evals
             except KeyError:
-                left, known_evals = self.left.evaluate(t, y, known_evals)
-                right, known_evals = self.right.evaluate(t, y, known_evals)
+                left, known_evals = self.left.evaluate(t, y, u, known_evals)
+                right, known_evals = self.right.evaluate(t, y, u, known_evals)
                 value = self._binary_evaluate(left, right)
                 known_evals[id] = value
                 return value, known_evals
         else:
-            left = self.left.evaluate(t, y)
-            right = self.right.evaluate(t, y)
+            left = self.left.evaluate(t, y, u)
+            right = self.right.evaluate(t, y, u)
             return self._binary_evaluate(left, right)
 
     def evaluate_for_shape(self):
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index c90d2874e0..4b2a96e269 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -52,20 +52,20 @@ def _concatenation_evaluate(self, children_eval):
         else:
             return self.concatenation_function(children_eval)
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         children = self.cached_children
         if known_evals is not None:
             if self.id not in known_evals:
                 children_eval = [None] * len(children)
                 for idx, child in enumerate(children):
-                    children_eval[idx], known_evals = child.evaluate(t, y, known_evals)
+                    children_eval[idx], known_evals = child.evaluate(t, y, u, known_evals)
                 known_evals[self.id] = self._concatenation_evaluate(children_eval)
             return known_evals[self.id], known_evals
         else:
             children_eval = [None] * len(children)
             for idx, child in enumerate(children):
-                children_eval[idx] = child.evaluate(t, y)
+                children_eval[idx] = child.evaluate(t, y, u)
             return self._concatenation_evaluate(children_eval)
 
     def new_copy(self):
diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index 8ce840815f..c2dbc8a95a 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -14,8 +14,8 @@ class Function(pybamm.Symbol):
     ----------
     function : method
         A function can have 0 or many inputs. If no inputs are given, self.evaluate()
-        simply returns func(). Otherwise, self.evaluate(t, y) returns
-        func(child0.evaluate(t, y), child1.evaluate(t, y), etc).
+        simply returns func(). Otherwise, self.evaluate(t, y, u) returns
+        func(child0.evaluate(t, y, u), child1.evaluate(t, y, u), etc).
     children : :class:`pybamm.Symbol`
         The children nodes to apply the function to
     derivative : str, optional
@@ -152,7 +152,7 @@ def _function_jac(self, children_jacs):
 
         return jacobian
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             if self.id not in known_evals:
@@ -164,7 +164,7 @@ def evaluate(self, t=None, y=None, known_evals=None):
                 known_evals[self.id] = self._function_evaluate(evaluated_children)
             return known_evals[self.id], known_evals
         else:
-            evaluated_children = [child.evaluate(t, y) for child in self.children]
+            evaluated_children = [child.evaluate(t, y, u) for child in self.children]
             return self._function_evaluate(evaluated_children)
 
     def evaluate_for_shape(self):
diff --git a/pybamm/expression_tree/input_parameter.py b/pybamm/expression_tree/input_parameter.py
new file mode 100644
index 0000000000..4fa23012dc
--- /dev/null
+++ b/pybamm/expression_tree/input_parameter.py
@@ -0,0 +1,40 @@
+#
+# Parameter classes
+#
+import numpy as np
+import pybamm
+
+
+class InputParameter(pybamm.Symbol):
+    """A node in the expression tree representing an input parameter
+
+    This node's value can be set at the point of solving, allowing parameter estimation
+    and control
+
+    Parameters
+    ----------
+    name : str
+        name of the node
+
+    """
+
+    def __init__(self, name):
+        super().__init__(name)
+
+    def new_copy(self):
+        """ See :meth:`pybamm.Symbol.new_copy()`. """
+        return InputParameter(self.name)
+
+    def evaluate_for_shape(self):
+        """
+        Returns the scalar 'NaN' to represent the shape of a parameter.
+        See :meth:`pybamm.Symbol.evaluate_for_shape()`
+        """
+        return np.nan
+
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+        # u should be a dictionary
+        if not isinstance(u, dict):
+            raise TypeError("inputs u should be a dictionary")
+        return u[self.name]
+
diff --git a/pybamm/expression_tree/operations/evaluate.py b/pybamm/expression_tree/operations/evaluate.py
index 43065a91b8..caffad205b 100644
--- a/pybamm/expression_tree/operations/evaluate.py
+++ b/pybamm/expression_tree/operations/evaluate.py
@@ -270,7 +270,7 @@ def __init__(self, symbol):
             self._result_var, "return" + self._result_var, "eval"
         )
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """
         Acts as a drop-in replacement for :func:`pybamm.Symbol.evaluate`
         """
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index a5af685ab2..5f8c8125fc 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -459,7 +459,7 @@ def _base_evaluate(self, t=None, y=None):
             )
         )
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """Evaluate expression tree (wrapper to allow using dict of known values).
         If the dict 'known_evals' is provided, the dict is searched for self.id; if
         self.id is in the keys, return that value; otherwise, evaluate using
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 5ce164498a..0fcd87c066 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -63,15 +63,15 @@ def _unary_evaluate(self, child):
         """Perform unary operation on a child. """
         raise NotImplementedError
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             if self.id not in known_evals:
-                child, known_evals = self.child.evaluate(t, y, known_evals)
+                child, known_evals = self.child.evaluate(t, y, u, known_evals)
                 known_evals[self.id] = self._unary_evaluate(child)
             return known_evals[self.id], known_evals
         else:
-            child = self.child.evaluate(t, y)
+            child = self.child.evaluate(t, y, u)
             return self._unary_evaluate(child)
 
     def evaluate_for_shape(self):
diff --git a/tests/unit/test_expression_tree/test_input_parameter.py b/tests/unit/test_expression_tree/test_input_parameter.py
new file mode 100644
index 0000000000..ce4c1868f6
--- /dev/null
+++ b/tests/unit/test_expression_tree/test_input_parameter.py
@@ -0,0 +1,28 @@
+#
+# Tests for the InputParameter class
+#
+import numbers
+import pybamm
+import unittest
+
+
+class TestInputParameter(unittest.TestCase):
+    def test_input_parameter_init(self):
+        a = pybamm.InputParameter("a")
+        self.assertEqual(a.name, "a")
+        self.assertEqual(a.evaluate(u={"a": 1}), 1)
+        self.assertEqual(a.evaluate(u={"a": 5}), 5)
+
+    def test_evaluate_for_shape(self):
+        a = pybamm.InputParameter("a")
+        self.assertIsInstance(a.evaluate_for_shape(), numbers.Number)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()
diff --git a/tests/unit/test_expression_tree/test_operations/test_copy.py b/tests/unit/test_expression_tree/test_operations/test_copy.py
index 3161967d20..c969e89928 100644
--- a/tests/unit/test_expression_tree/test_operations/test_copy.py
+++ b/tests/unit/test_expression_tree/test_operations/test_copy.py
@@ -37,6 +37,7 @@ def test_symbol_new_copy(self):
             pybamm.NumpyConcatenation(a, b, v_s),
             pybamm.DomainConcatenation([v_n, v_s], mesh),
             pybamm.Parameter("param"),
+            pybamm.InputParameter("param"),
             pybamm.StateVector(slice(0, 56)),
             pybamm.Matrix(np.ones((50, 40))),
             pybamm.SpatialVariable("x", ["negative electrode"]),

From 7b814244f9c991cdbeb58188f51fc14de2e45347 Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Sun, 24 Nov 2019 18:10:11 -0500
Subject: [PATCH 02/18] process input parameter values

---
 pybamm/parameters/parameter_values.py           | 11 ++++++++---
 .../test_parameters/test_parameter_values.py    | 17 +++++++++++++++++
 2 files changed, 25 insertions(+), 3 deletions(-)

diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 8566ece726..b3f98aeb1e 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -175,6 +175,8 @@ def update(self, values, check_conflict=False, path=""):
                     # Special case (hacky) for zero current
                     elif value == "[zero]":
                         super().__setitem__(name, 0)
+                    elif value == "[input]":
+                        super().__setitem__(name, pybamm.InputParameter(name))
                     # Anything else should be a converted to a float
                     else:
                         super().__setitem__(name, float(value))
@@ -419,9 +421,12 @@ def _process_symbol(self, symbol):
 
         if isinstance(symbol, pybamm.Parameter):
             value = self[symbol.name]
-            # Scalar inherits name (for updating parameters) and domain (for Broadcast)
-            return pybamm.Scalar(value, name=symbol.name, domain=symbol.domain)
-
+            if isinstance(value, numbers.Number):
+                # Scalar inherits name (for updating parameters) and domain (for Broadcast)
+                return pybamm.Scalar(value, name=symbol.name, domain=symbol.domain)
+            elif isinstance(value, pybamm.InputParameter):
+                return value
+                
         elif isinstance(symbol, pybamm.FunctionParameter):
             new_children = [self.process_symbol(child) for child in symbol.children]
             function_name = self[symbol.name]
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index 05cd23b7da..c5182c3c8c 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -228,6 +228,23 @@ def test_process_symbol(self):
         with self.assertRaises(NotImplementedError):
             parameter_values.process_symbol(sym)
 
+    def test_process_input_parameter(self):
+        parameter_values = pybamm.ParameterValues({"a": "[input]", "b": 3})
+        # process input parameter
+        a = pybamm.Parameter("a")
+        processed_a = parameter_values.process_symbol(a)
+        self.assertIsInstance(processed_a, pybamm.InputParameter)
+        self.assertEqual(processed_a.evaluate(u={"a": 5}), 5)
+
+        # process binary operation
+        b = pybamm.Parameter("b")
+        add = a + b
+        processed_add = parameter_values.process_symbol(add)
+        self.assertIsInstance(processed_add, pybamm.Addition)
+        self.assertIsInstance(processed_add.children[0], pybamm.InputParameter)
+        self.assertIsInstance(processed_add.children[1], pybamm.Scalar)
+        self.assertEqual(processed_add.evaluate(u={"a": 4}), 7)
+
     def test_process_function_parameter(self):
         parameter_values = pybamm.ParameterValues(
             {

From 4a760de143c1e35a71c80175e6e02d386742d471 Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Mon, 25 Nov 2019 20:58:48 -0500
Subject: [PATCH 03/18] test scipy solver

---
 pybamm/expression_tree/input_parameter.py     |   9 +-
 .../operations/convert_to_casadi.py           |  41 ++++---
 pybamm/expression_tree/symbol.py              |  21 +++-
 pybamm/simulation.py                          |   6 +-
 pybamm/solvers/base_solver.py                 |  19 ++-
 pybamm/solvers/dae_solver.py                  |  55 ++++-----
 pybamm/solvers/ode_solver.py                  | 112 ++++++++++++------
 pybamm/solvers/scipy_solver.py                |   1 +
 .../test_input_parameter.py                   |   7 ++
 .../test_operations/test_convert_to_casadi.py |  34 ++++++
 tests/unit/test_solvers/test_scipy_solver.py  |  25 ++++
 11 files changed, 232 insertions(+), 98 deletions(-)

diff --git a/pybamm/expression_tree/input_parameter.py b/pybamm/expression_tree/input_parameter.py
index 4fa23012dc..a8371ebb2e 100644
--- a/pybamm/expression_tree/input_parameter.py
+++ b/pybamm/expression_tree/input_parameter.py
@@ -35,6 +35,13 @@ def evaluate_for_shape(self):
     def evaluate(self, t=None, y=None, u=None, known_evals=None):
         # u should be a dictionary
         if not isinstance(u, dict):
+            # if the special input "shape test" is passed, just return 1
+            if u == "shape test":
+                return 1
             raise TypeError("inputs u should be a dictionary")
-        return u[self.name]
+        try:
+            return u[self.name]
+        # raise more informative error if can't find name in dict
+        except KeyError:
+            raise KeyError("Input parameter '{}' not found".format(self.name))
 
diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py
index e6e0758410..12a91a31a1 100644
--- a/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -11,36 +11,43 @@ class CasadiConverter(object):
     def __init__(self, casadi_symbols=None):
         self._casadi_symbols = casadi_symbols or {}
 
-    def convert(self, symbol, t=None, y=None):
+    def convert(self, symbol, t=None, y=None, u=None):
         """
-        This function recurses down the tree, applying any simplifications defined in
-        classes derived from pybamm.Symbol. E.g. any expression multiplied by a
-        pybamm.Scalar(0) will be simplified to a pybamm.Scalar(0).
-        If a symbol has already been simplified, the stored value is returned.
+        This function recurses down the tree, converting the PyBaMM expression tree to
+        a CasADi expression tree
 
         Parameters
         ----------
         symbol : :class:`pybamm.Symbol`
             The symbol to convert
+        t : :class:`casadi.MX`
+            A casadi symbol representing time
+        y : :class:`casadi.MX`
+            A casadi symbol representing state vectors
+        u : dict
+            A dictionary of casadi symbols representing inputs
 
         Returns
         -------
-        CasADi symbol
-            The convert symbol
+        :class:`casadi.MX`
+            The converted symbol
         """
-
         try:
             return self._casadi_symbols[symbol.id]
         except KeyError:
-            casadi_symbol = self._convert(symbol, t, y)
+            # Change u to empty dictionary if it's None
+            u = u or {}
+            casadi_symbol = self._convert(symbol, t, y, u)
             self._casadi_symbols[symbol.id] = casadi_symbol
 
             return casadi_symbol
 
-    def _convert(self, symbol, t, y):
+    def _convert(self, symbol, t, y, u):
         """ See :meth:`CasadiConverter.convert()`. """
-        if isinstance(symbol, (pybamm.Scalar, pybamm.Array, pybamm.Time)):
-            return casadi.MX(symbol.evaluate(t, y))
+        if isinstance(
+            symbol, (pybamm.Scalar, pybamm.Array, pybamm.Time, pybamm.InputParameter)
+        ):
+            return casadi.MX(symbol.evaluate(t, y, u))
 
         elif isinstance(symbol, pybamm.StateVector):
             if y is None:
@@ -50,8 +57,8 @@ def _convert(self, symbol, t, y):
         elif isinstance(symbol, pybamm.BinaryOperator):
             left, right = symbol.children
             # process children
-            converted_left = self.convert(left, t, y)
-            converted_right = self.convert(right, t, y)
+            converted_left = self.convert(left, t, y, u)
+            converted_right = self.convert(right, t, y, u)
             if isinstance(symbol, pybamm.Outer):
                 return casadi.kron(converted_left, converted_right)
             else:
@@ -59,14 +66,14 @@ def _convert(self, symbol, t, y):
                 return symbol._binary_evaluate(converted_left, converted_right)
 
         elif isinstance(symbol, pybamm.UnaryOperator):
-            converted_child = self.convert(symbol.child, t, y)
+            converted_child = self.convert(symbol.child, t, y, u)
             if isinstance(symbol, pybamm.AbsoluteValue):
                 return casadi.fabs(converted_child)
             return symbol._unary_evaluate(converted_child)
 
         elif isinstance(symbol, pybamm.Function):
             converted_children = [
-                self.convert(child, t, y) for child in symbol.children
+                self.convert(child, t, y, u) for child in symbol.children
             ]
             # Special functions
             if symbol.function == np.min:
@@ -97,7 +104,7 @@ def _convert(self, symbol, t, y):
                 return symbol._function_evaluate(converted_children)
         elif isinstance(symbol, pybamm.Concatenation):
             converted_children = [
-                self.convert(child, t, y) for child in symbol.children
+                self.convert(child, t, y, u) for child in symbol.children
             ]
             if isinstance(symbol, (pybamm.NumpyConcatenation, pybamm.SparseStack)):
                 return casadi.vertcat(*converted_children)
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 5f8c8125fc..dbda1aead0 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -498,6 +498,7 @@ def evaluate_for_shape(self):
 
     def is_constant(self):
         """returns true if evaluating the expression is not dependent on `t` or `y`
+        or `u`
 
         See Also
         --------
@@ -505,8 +506,13 @@ def is_constant(self):
 
         """
         # if any of the nodes are instances of any of these types, then the whole
-        # expression depends on either t or y
-        search_types = (pybamm.Variable, pybamm.StateVector, pybamm.IndependentVariable)
+        # expression depends on either t or y or u
+        search_types = (
+            pybamm.Variable,
+            pybamm.StateVector,
+            pybamm.Time,
+            pybamm.InputParameter,
+        )
 
         # do the search, return true if no relevent nodes are found
         return not any((isinstance(n, search_types)) for n in self.pre_order())
@@ -531,7 +537,10 @@ def evaluate_ignoring_errors(self):
         except TypeError as error:
             # return false if specific TypeError is raised
             # (there is a e.g. StateVector in the tree)
-            if error.args[0] == "StateVector cannot evaluate input 'y=None'":
+            if error.args[0] in [
+                "StateVector cannot evaluate input 'y=None'",
+                "inputs u should be a dictionary",
+            ]:
                 return None
             else:
                 raise error
@@ -579,12 +588,12 @@ def simplify(self, simplified_symbols=None):
         """ Simplify the expression tree. See :class:`pybamm.Simplification`. """
         return pybamm.Simplification(simplified_symbols).simplify(self)
 
-    def to_casadi(self, t=None, y=None, casadi_symbols=None):
+    def to_casadi(self, t=None, y=None, u=None, casadi_symbols=None):
         """
         Convert the expression tree to a CasADi expression tree.
         See :class:`pybamm.CasadiConverter`.
         """
-        return pybamm.CasadiConverter(casadi_symbols).convert(self, t, y)
+        return pybamm.CasadiConverter(casadi_symbols).convert(self, t, y, u)
 
     def new_copy(self):
         """
@@ -614,7 +623,7 @@ def shape(self):
         # Try with some large y, to avoid having to use pre_order (slow)
         try:
             y = np.linspace(0.1, 0.9, int(1e4))
-            evaluated_self = self.evaluate(0, y)
+            evaluated_self = self.evaluate(0, y, u="shape test")
         # If that fails, fall back to calculating how big y should really be
         except ValueError:
             state_vectors_in_node = [
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 41c383166e..4ae6af837c 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -144,7 +144,7 @@ def build(self):
         self._disc = pybamm.Discretisation(self._mesh, self._spatial_methods)
         self._built_model = self._disc.process_model(self._model, inplace=False)
 
-    def solve(self, t_eval=None, solver=None):
+    def solve(self, t_eval=None, solver=None, inputs=None):
         """
         A method to solve the model. This method will automatically build
         and set the model parameters if not already done so.
@@ -158,6 +158,8 @@ def solve(self, t_eval=None, solver=None):
             non-dimensional time of 1.
         solver : :class:`pybamm.BaseSolver`
             The solver to use to solve the model.
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
         """
         self.build()
 
@@ -174,7 +176,7 @@ def solve(self, t_eval=None, solver=None):
         if solver is None:
             solver = self.solver
 
-        self._solution = solver.solve(self.built_model, t_eval)
+        self._solution = solver.solve(self.built_model, t_eval, inputs=inputs)
 
     def step(self, dt, solver=None, external_variables=None, save=True):
         """
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 5ef78b690e..de9d72a307 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -49,7 +49,7 @@ def atol(self):
     def atol(self, value):
         self._atol = value
 
-    def solve(self, model, t_eval):
+    def solve(self, model, t_eval, inputs=None):
         """
         Execute the solver setup and calculate the solution of the model at
         specified times.
@@ -61,6 +61,8 @@ def solve(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         Raises
         ------
@@ -77,14 +79,17 @@ def solve(self, model, t_eval):
         # Set up
         timer = pybamm.Timer()
         start_time = timer.time()
+        inputs = inputs or {}
         if model.convert_to_format == "casadi" or isinstance(self, pybamm.CasadiSolver):
-            self.set_up_casadi(model)
+            self.set_up_casadi(model, inputs)
         else:
             self.set_up(model)
         set_up_time = timer.time() - start_time
 
         # Solve
-        solution, solve_time, termination = self.compute_solution(model, t_eval)
+        solution, solve_time, termination = self.compute_solution(
+            model, t_eval, inputs=inputs
+        )
 
         # Assign times
         solution.solve_time = solve_time
@@ -100,7 +105,7 @@ def solve(self, model, t_eval):
         )
         return solution
 
-    def step(self, model, dt, npts=2, log=True, external_variables=None):
+    def step(self, model, dt, npts=2, log=True, external_variables=None, inputs=None):
         """
         Step the solution of the model forward by a given time increment. The
         first time this method is called it executes the necessary setup by
@@ -119,6 +124,9 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
         external_variables : dict
             A dictionary of external variables and their corresponding
             values at the current time
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
+
 
         Raises
         ------
@@ -132,6 +140,7 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
 
         # Set timer
         timer = pybamm.Timer()
+        inputs = inputs or {}
 
         if not hasattr(self, "y0"):
             # create a y_pad vector of the correct size:
@@ -162,7 +171,7 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
 
         # Step
         t_eval = np.linspace(self.t, self.t + dt, npts)
-        solution, solve_time, termination = self.compute_solution(model, t_eval)
+        solution, solve_time, termination = self.compute_solution(model, t_eval, inputs)
 
         # Set self.t and self.y0 to their values at the final step
         self.t = solution.t[-1]
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 33d312de49..d7b431c31c 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -67,7 +67,7 @@ def max_steps(self):
     def max_steps(self, max_steps):
         self._max_steps = max_steps
 
-    def compute_solution(self, model, t_eval):
+    def compute_solution(self, model, t_eval, inputs=None):
         """Calculate the solution of the model at specified times.
 
         Parameters
@@ -77,7 +77,8 @@ def compute_solution(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
-
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
         """
         timer = pybamm.Timer()
 
@@ -100,6 +101,7 @@ def compute_solution(self, model, t_eval):
             mass_matrix=model.mass_matrix.entries,
             jacobian=self.jacobian,
             model=model,
+            inputs=inputs,
         )
 
         solve_time = timer.time() - solve_start_time
@@ -439,18 +441,31 @@ def integrate(
         raise NotImplementedError
 
 
-class Rhs:
-    "Returns information about rhs at time t and state y"
-
-    def __init__(self, concatenated_rhs_fn):
-        self.concatenated_rhs_fn = concatenated_rhs_fn
-        self.y_pad = None
-        self.y_ext = None
+class SolverCallable:
+    "A class that will be called by the solver when integrating"
+    y_pad = None
+    y_ext = None
+    _inputs = {}
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
         self.y_ext = y_ext
 
+    @property
+    def inputs(self):
+        return self._inputs
+
+    @inputs.setter
+    def inputs(self, value):
+        self._inputs = value
+
+
+class Rhs(SolverCallable):
+    "Returns information about rhs at time t and state y"
+
+    def __init__(self, concatenated_rhs_fn):
+        self.concatenated_rhs_fn = concatenated_rhs_fn
+
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
@@ -466,17 +481,11 @@ def __call__(self, t, y):
         return self.concatenated_rhs_fn(t, y).full()[:, 0]
 
 
-class Algebraic:
+class Algebraic(SolverCallable):
     "Returns information about algebraic equations at time t and state y"
 
     def __init__(self, concatenated_algebraic_fn):
         self.concatenated_algebraic_fn = concatenated_algebraic_fn
-        self.y_pad = None
-        self.y_ext = None
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def __call__(self, t, y):
         y = y[:, np.newaxis]
@@ -493,7 +502,7 @@ def __call__(self, t, y):
         return self.concatenated_algebraic_fn(t, y).full()[:, 0]
 
 
-class Residuals:
+class Residuals(SolverCallable):
     "Returns information about residuals at time t and state y"
 
     def __init__(self, model, concatenated_rhs_fn, concatenated_algebraic_fn):
@@ -501,12 +510,6 @@ def __init__(self, model, concatenated_rhs_fn, concatenated_algebraic_fn):
         self.concatenated_rhs_fn = concatenated_rhs_fn
         self.concatenated_algebraic_fn = concatenated_algebraic_fn
         self.mass_matrix = model.mass_matrix.entries
-        self.y_pad = None
-        self.y_ext = None
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def __call__(self, t, y, ydot):
         pybamm.logger.debug(
@@ -530,8 +533,6 @@ def __init__(self, model, all_states_fn):
         self.model = model
         self.all_states_fn = all_states_fn
         self.mass_matrix = model.mass_matrix.entries
-        self.y_pad = None
-        self.y_ext = None
 
     def __call__(self, t, y, ydot):
         pybamm.logger.debug(
@@ -548,8 +549,6 @@ class EvalEvent:
 
     def __init__(self, event_fn):
         self.event_fn = event_fn
-        self.y_pad = None
-        self.y_ext = None
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -566,8 +565,6 @@ class Jacobian:
 
     def __init__(self, jac_fn):
         self.jac_fn = jac_fn
-        self.y_pad = None
-        self.y_ext = None
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 934e521ef5..1b86a3d66f 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -23,7 +23,7 @@ def __init__(self, method=None, rtol=1e-6, atol=1e-6):
         super().__init__(method, rtol, atol)
         self.name = "Base ODE solver"
 
-    def compute_solution(self, model, t_eval):
+    def compute_solution(self, model, t_eval, inputs=None):
         """Calculate the solution of the model at specified times.
 
         Parameters
@@ -33,16 +33,16 @@ def compute_solution(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         """
         timer = pybamm.Timer()
 
-        self.dydt.set_pad_ext(self.y_pad, self.y_ext)
-        for evnt in self.event_funs:
-            evnt.set_pad_ext(self.y_pad, self.y_ext)
-        if self.jacobian:
-            self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+        # Set inputs and external
+        self.set_inputs_and_external(inputs)
 
+        # Solve
         solve_start_time = timer.time()
         pybamm.logger.info("Calling ODE solver")
         solution = self.integrate(
@@ -146,7 +146,7 @@ def get_event_class(event):
         self.event_funs = [get_event_class(event) for event in events.values()]
         self.jacobian = jacobian
 
-    def set_up_casadi(self, model):
+    def set_up_casadi(self, model, inputs):
         """Convert model to casadi format and use their inbuilt functionalities.
 
         Parameters
@@ -154,6 +154,8 @@ def set_up_casadi(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         Raises
         ------
@@ -172,6 +174,7 @@ def set_up_casadi(self, model):
 
         t_casadi = casadi.MX.sym("t")
         y_casadi = casadi.MX.sym("y", len(y0))
+        u_casadi = {name: casadi.MX.sym(name) for name in inputs.keys()}
 
         if self.y_pad is not None:
             y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
@@ -180,31 +183,34 @@ def set_up_casadi(self, model):
             y_casadi_w_ext = y_casadi
 
         pybamm.logger.info("Converting RHS to CasADi")
-        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi_w_ext)
+        concatenated_rhs = model.concatenated_rhs.to_casadi(
+            t_casadi, y_casadi_w_ext, u_casadi
+        )
         pybamm.logger.info("Converting events to CasADi")
         casadi_events = {
-            name: event.to_casadi(t_casadi, y_casadi_w_ext)
+            name: event.to_casadi(t_casadi, y_casadi_w_ext, u_casadi)
             for name, event in model.events.items()
         }
 
         # Create function to evaluate rhs
+        u_casadi_stacked = casadi.vertcat(*[u for u in u_casadi.values()])
         concatenated_rhs_fn = casadi.Function(
-            "rhs", [t_casadi, y_casadi_w_ext], [concatenated_rhs]
+            "rhs", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [concatenated_rhs]
         )
 
         # Create event-dependent function to evaluate events
         def get_event_class(event):
             casadi_event_fn = casadi.Function(
-                "event", [t_casadi, y_casadi_w_ext], [event]
+                "event", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [event]
             )
-            return EvalEvent(casadi_event_fn)
+            return EvalEventCasadi(casadi_event_fn)
 
         # Create function to evaluate jacobian
         if model.use_jacobian:
             pybamm.logger.info("Calculating jacobian")
             casadi_jac = casadi.jacobian(concatenated_rhs, y_casadi)
             casadi_jac_fn = casadi.Function(
-                "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac]
+                "jacobian", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [casadi_jac]
             )
 
             jacobian = JacobianCasadi(casadi_jac_fn)
@@ -219,6 +225,26 @@ def get_event_class(event):
         self.event_funs = [get_event_class(event) for event in casadi_events.values()]
         self.jacobian = jacobian
 
+    def set_inputs_and_external(self, inputs):
+        """
+        Set values that are controlled externally, such as external variables and input
+        parameters
+
+        Parameters
+        ----------
+        inputs : dict
+            Any input parameters to pass to the model when solving
+            
+        """
+        self.dydt.set_pad_ext(self.y_pad, self.y_ext)
+        self.dydt.set_inputs(inputs)
+        for evnt in self.event_funs:
+            evnt.set_pad_ext(self.y_pad, self.y_ext)
+            evnt.set_inputs(inputs)
+        if self.jacobian:
+            self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+            self.jacobian.set_inputs(inputs)
+
     def integrate(
         self, derivs, y0, t_eval, events=None, mass_matrix=None, jacobian=None
     ):
@@ -244,25 +270,35 @@ def integrate(
         raise NotImplementedError
 
 
+class SolverCallable:
+    "A class that will be called by the solver when integrating"
+    y_pad = None
+    y_ext = None
+    inputs = {}
+    inputs_casadi = casadi.DM()
+
+    def set_pad_ext(self, y_pad, y_ext):
+        self.y_pad = y_pad
+        self.y_ext = y_ext
+
+    def set_inputs(self, inputs):
+        self.inputs = inputs
+        self.inputs_casadi = casadi.vertcat(*[x for x in inputs.values()])
+
+
 # Set up caller classes outside of the solver object to allow pickling
-class Dydt:
+class Dydt(SolverCallable):
     "Returns information about time derivatives at time t and state y"
 
     def __init__(self, model, concatenated_rhs_fn):
         self.model = model
         self.concatenated_rhs_fn = concatenated_rhs_fn
-        self.y_pad = None
-        self.y_ext = None
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def __call__(self, t, y):
         pybamm.logger.debug("Evaluating RHS for {} at t={}".format(self.model.name, t))
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        dy = self.concatenated_rhs_fn(t, y, known_evals={})[0]
+        dy = self.concatenated_rhs_fn(t, y, self.inputs, known_evals={})[0]
         return dy[:, 0]
 
 
@@ -273,44 +309,44 @@ def __call__(self, t, y):
         pybamm.logger.debug("Evaluating RHS for {} at t={}".format(self.model.name, t))
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        dy = self.concatenated_rhs_fn(t, y).full()
+        dy = self.concatenated_rhs_fn(t, y, self.inputs_casadi).full()
         return dy[:, 0]
 
 
-class EvalEvent:
+class EvalEvent(SolverCallable):
     "Returns information about events at time t and state y"
 
     def __init__(self, event_fn):
         self.event_fn = event_fn
-        self.y_pad = None
-        self.y_ext = None
 
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
+    def __call__(self, t, y):
+        y = y[:, np.newaxis]
+        y = add_external(y, self.y_pad, self.y_ext)
+        return self.event_fn(t, y, self.inputs)
+
+
+class EvalEventCasadi(EvalEvent):
+    "Returns information about events at time t and state y"
+
+    def __init__(self, event_fn):
+        self.event_fn = event_fn
 
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.event_fn(t, y)
+        return self.event_fn(t, y, self.inputs_casadi)
 
 
-class Jacobian:
+class Jacobian(SolverCallable):
     "Returns information about the jacobian at time t and state y"
 
     def __init__(self, jac_fn):
         self.jac_fn = jac_fn
-        self.y_pad = None
-        self.y_ext = None
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.jac_fn(t, y, known_evals={})[0]
+        return self.jac_fn(t, y, self.inputs, known_evals={})[0]
 
 
 class JacobianCasadi(Jacobian):
@@ -319,4 +355,4 @@ class JacobianCasadi(Jacobian):
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.jac_fn(t, y)
+        return self.jac_fn(t, y, self.inputs)
diff --git a/pybamm/solvers/scipy_solver.py b/pybamm/solvers/scipy_solver.py
index fd2af8d763..1a99eabd33 100644
--- a/pybamm/solvers/scipy_solver.py
+++ b/pybamm/solvers/scipy_solver.py
@@ -47,6 +47,7 @@ def integrate(
         jacobian : method, optional
             A function that takes in t and y and returns the Jacobian. If
             None, the solver will approximate the Jacobian.
+            
         Returns
         -------
         object
diff --git a/tests/unit/test_expression_tree/test_input_parameter.py b/tests/unit/test_expression_tree/test_input_parameter.py
index ce4c1868f6..b64d8d8d30 100644
--- a/tests/unit/test_expression_tree/test_input_parameter.py
+++ b/tests/unit/test_expression_tree/test_input_parameter.py
@@ -17,6 +17,13 @@ def test_evaluate_for_shape(self):
         a = pybamm.InputParameter("a")
         self.assertIsInstance(a.evaluate_for_shape(), numbers.Number)
 
+    def test_errors(self):
+        a = pybamm.InputParameter("a")
+        with self.assertRaises(TypeError):
+            a.evaluate(u="not a dictionary")
+        with self.assertRaises(KeyError):
+            a.evaluate(u={"bad param": 5})
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
index ab7d50d9b9..d0ee4fcf1e 100644
--- a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
+++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
@@ -157,6 +157,40 @@ def myfunction(x, y):
         f = pybamm.Function(myfunction, a, b).diff(b)
         self.assert_casadi_equal(f.to_casadi(), casadi.MX(3), evalf=True)
 
+    def test_convert_input_parameter(self):
+        # Arrays
+        a = np.array([1, 2, 3, 4, 5])
+        pybamm_a = pybamm.Array(a)
+        self.assert_casadi_equal(pybamm_a.to_casadi(), casadi.MX(a))
+
+        casadi_t = casadi.MX.sym("t")
+        casadi_y = casadi.MX.sym("y", 10)
+        casadi_us = {
+            "Input 1": casadi.MX.sym("Input 1"),
+            "Input 2": casadi.MX.sym("Input 2"),
+        }
+
+        pybamm_y = pybamm.StateVector(slice(0, 10))
+        pybamm_u1 = pybamm.InputParameter("Input 1")
+        pybamm_u2 = pybamm.InputParameter("Input 2")
+
+        # Input only
+        self.assert_casadi_equal(
+            pybamm_u1.to_casadi(casadi_t, casadi_y, casadi_us), casadi_us["Input 1"]
+        )
+
+        # More complex
+        expr = pybamm_u1 + pybamm_y
+        self.assert_casadi_equal(
+            expr.to_casadi(casadi_t, casadi_y, casadi_us),
+            casadi_us["Input 1"] + casadi_y,
+        )
+        expr = pybamm_u2 * pybamm_y
+        self.assert_casadi_equal(
+            expr.to_casadi(casadi_t, casadi_y, casadi_us),
+            casadi_us["Input 2"] * casadi_y,
+        )
+
     def test_errors(self):
         y = pybamm.StateVector(slice(0, 10))
         with self.assertRaisesRegex(
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 68ecd1de05..2ae309617a 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -299,6 +299,31 @@ def test_model_solver_with_event_with_casadi(self):
             np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
             np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
 
+    def test_model_solver_with_inputs_with_casadi(self):
+        # Create model
+        model = pybamm.BaseModel()
+        model.convert_to_format = "casadi"
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=domain)
+        model.rhs = {var: -pybamm.InputParameter("rate") * var}
+        model.initial_conditions = {var: 1}
+        model.events = {"var=0.5": pybamm.min(var - 0.5)}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+        # Solve
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+        t_eval = np.linspace(0, 10, 100)
+        solution = solver.solve(model, t_eval, inputs={"rate": 0.1})
+        self.assertLess(len(solution.t), len(t_eval))
+        np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+        np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

From cfcc83ee69592ea5bc59fce16392ec017059fa9d Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Mon, 25 Nov 2019 21:03:18 -0500
Subject: [PATCH 04/18] #743 flake8

---
 pybamm/expression_tree/concatenations.py       |  4 +++-
 pybamm/expression_tree/input_parameter.py      |  1 -
 pybamm/logger.py                               |  1 -
 pybamm/meshes/meshes.py                        |  4 +---
 pybamm/models/base_model.py                    |  1 -
 .../full_battery_models/base_battery_model.py  |  2 +-
 .../models/submodels/porosity/base_porosity.py |  1 -
 .../models/submodels/thermal/base_thermal.py   |  4 ++--
 pybamm/parameters/parameter_values.py          |  6 ++++--
 pybamm/solvers/casadi_solver.py                |  1 -
 pybamm/solvers/dae_solver.py                   |  1 -
 pybamm/solvers/ode_solver.py                   |  2 +-
 pybamm/solvers/scipy_solver.py                 |  2 +-
 pybamm/spatial_methods/finite_volume.py        | 18 +++++++-----------
 14 files changed, 20 insertions(+), 28 deletions(-)

diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index 4b2a96e269..2508b4cc8b 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -59,7 +59,9 @@ def evaluate(self, t=None, y=None, u=None, known_evals=None):
             if self.id not in known_evals:
                 children_eval = [None] * len(children)
                 for idx, child in enumerate(children):
-                    children_eval[idx], known_evals = child.evaluate(t, y, u, known_evals)
+                    children_eval[idx], known_evals = child.evaluate(
+                        t, y, u, known_evals
+                    )
                 known_evals[self.id] = self._concatenation_evaluate(children_eval)
             return known_evals[self.id], known_evals
         else:
diff --git a/pybamm/expression_tree/input_parameter.py b/pybamm/expression_tree/input_parameter.py
index a8371ebb2e..dab5725356 100644
--- a/pybamm/expression_tree/input_parameter.py
+++ b/pybamm/expression_tree/input_parameter.py
@@ -44,4 +44,3 @@ def evaluate(self, t=None, y=None, u=None, known_evals=None):
         # raise more informative error if can't find name in dict
         except KeyError:
             raise KeyError("Input parameter '{}' not found".format(self.name))
-
diff --git a/pybamm/logger.py b/pybamm/logger.py
index 27853fb298..2be256297f 100644
--- a/pybamm/logger.py
+++ b/pybamm/logger.py
@@ -18,4 +18,3 @@ def set_logging_level(level):
 # Create a custom logger
 logger = logging.getLogger(__name__)
 set_logging_level("WARNING")
-
diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py
index 084b159689..30c5b3be77 100644
--- a/pybamm/meshes/meshes.py
+++ b/pybamm/meshes/meshes.py
@@ -235,6 +235,4 @@ def __call__(self, lims, npts, tabs=None):
         return self.submesh_type(lims, npts, tabs, **self.submesh_params)
 
     def __repr__(self):
-        return "Generator for {}".format(
-            self.submesh_type.__name__
-        )
+        return "Generator for {}".format(self.submesh_type.__name__)
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 78625704d4..38a66107a5 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -538,4 +538,3 @@ def default_solver(self):
             return pybamm.IDAKLUSolver()
         else:
             return pybamm.CasadiSolver(mode="safe")
-
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index b16e88a637..50b2bb35a6 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -156,7 +156,7 @@ def options(self, extra_options):
             "particle": "Fickian diffusion",
             "thermal": "isothermal",
             "thermal current collector": False,
-            "external submodels": []
+            "external submodels": [],
         }
         options = default_options
         # any extra options overwrite the default options
diff --git a/pybamm/models/submodels/porosity/base_porosity.py b/pybamm/models/submodels/porosity/base_porosity.py
index 72012656ab..72a5082fce 100644
--- a/pybamm/models/submodels/porosity/base_porosity.py
+++ b/pybamm/models/submodels/porosity/base_porosity.py
@@ -100,4 +100,3 @@ def set_events(self, variables):
         self.events["Max negative electrode porosity cut-off"] = pybamm.max(eps_n) - 1
         self.events["Zero positive electrode porosity cut-off"] = pybamm.min(eps_p)
         self.events["Max positive electrode porosity cut-off"] = pybamm.max(eps_p) - 1
-
diff --git a/pybamm/models/submodels/thermal/base_thermal.py b/pybamm/models/submodels/thermal/base_thermal.py
index 2f7186b187..030a631297 100644
--- a/pybamm/models/submodels/thermal/base_thermal.py
+++ b/pybamm/models/submodels/thermal/base_thermal.py
@@ -173,8 +173,8 @@ def _get_standard_coupled_variables(self, variables):
                 "X-averaged irreversible electrochemical heating [W.m-3]": Q_rxn_av
                 * Q_scale,
                 "Volume-averaged irreversible electrochemical heating": Q_rxn_vol_av,
-                "Volume-averaged irreversible electrochemical heating [W.m-3]":
-                Q_rxn_vol_av * Q_scale,
+                "Volume-averaged irreversible electrochemical heating "
+                + "[W.m-3]": Q_rxn_vol_av * Q_scale,
                 "Reversible heating": Q_rev,
                 "Reversible heating [W.m-3]": Q_rev * Q_scale,
                 "X-averaged reversible heating": Q_rev_av,
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index b3f98aeb1e..2e98c72707 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -422,11 +422,13 @@ def _process_symbol(self, symbol):
         if isinstance(symbol, pybamm.Parameter):
             value = self[symbol.name]
             if isinstance(value, numbers.Number):
-                # Scalar inherits name (for updating parameters) and domain (for Broadcast)
+                # Scalar inherits name (for updating parameters) and domain (for
+                # Broadcast)
                 return pybamm.Scalar(value, name=symbol.name, domain=symbol.domain)
             elif isinstance(value, pybamm.InputParameter):
+                value.domain = symbol.domain
                 return value
-                
+
         elif isinstance(symbol, pybamm.FunctionParameter):
             new_children = [self.process_symbol(child) for child in symbol.children]
             function_name = self[symbol.name]
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index b88ceed73c..b5c0b5c0c5 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -261,4 +261,3 @@ def integrate_casadi(self, rhs, algebraic, y0, t_eval, mass_matrix=None):
         except RuntimeError as e:
             # If it doesn't work raise error
             raise pybamm.SolverError(e.args[0])
-
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index d7b431c31c..ba1f24004f 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -583,4 +583,3 @@ def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
         return self.jac_fn(t, y)
-
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 1b86a3d66f..9810c7c36b 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -234,7 +234,7 @@ def set_inputs_and_external(self, inputs):
         ----------
         inputs : dict
             Any input parameters to pass to the model when solving
-            
+
         """
         self.dydt.set_pad_ext(self.y_pad, self.y_ext)
         self.dydt.set_inputs(inputs)
diff --git a/pybamm/solvers/scipy_solver.py b/pybamm/solvers/scipy_solver.py
index 1a99eabd33..1541a3b801 100644
--- a/pybamm/solvers/scipy_solver.py
+++ b/pybamm/solvers/scipy_solver.py
@@ -47,7 +47,7 @@ def integrate(
         jacobian : method, optional
             A function that takes in t and y and returns the Jacobian. If
             None, the solver will approximate the Jacobian.
-            
+
         Returns
         -------
         object
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 98f8898f38..e1947aaaf5 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -639,9 +639,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
         # https://github.com/Scottmar93/extrapolation-coefficents/tree/master
         if isinstance(symbol, pybamm.BoundaryValue):
 
-            if use_bcs and pybamm.has_bc_of_form(
-                child, symbol.side, bcs, "Dirichlet"
-            ):
+            if use_bcs and pybamm.has_bc_of_form(child, symbol.side, bcs, "Dirichlet"):
                 # just use the value from the bc: f(x*)
                 sub_matrix = csr_matrix((1, prim_pts))
                 additive = bcs[child.id][symbol.side][0]
@@ -655,7 +653,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                     if use_bcs and pybamm.has_bc_of_form(
                         child, symbol.side, bcs, "Neumann"
                     ):
-                        sub_matrix = csr_matrix(([1], ([0], [0])), shape=(1, prim_pts),)
+                        sub_matrix = csr_matrix(([1], ([0], [0])), shape=(1, prim_pts))
 
                         additive = -dx0 * bcs[child.id][symbol.side][0]
 
@@ -676,7 +674,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         alpha = -(dx0 * (dx0 + dx1)) / (2 * dx0 + dx1)
 
                         sub_matrix = csr_matrix(
-                            ([a, b], ([0, 0], [0, 1])), shape=(1, prim_pts),
+                            ([a, b], ([0, 0], [0, 1])), shape=(1, prim_pts)
                         )
                         additive = alpha * bcs[child.id][symbol.side][0]
 
@@ -686,7 +684,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         c = dx0 * (dx0 + dx1) / (dx2 * (dx1 + dx2))
 
                         sub_matrix = csr_matrix(
-                            ([a, b, c], ([0, 0, 0], [0, 1, 2])), shape=(1, prim_pts),
+                            ([a, b, c], ([0, 0, 0], [0, 1, 2])), shape=(1, prim_pts)
                         )
 
                         additive = pybamm.Scalar(0)
@@ -703,7 +701,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         # use formula:
                         # f(x*) = fN + dxN * f'(x*)
                         sub_matrix = csr_matrix(
-                            ([1], ([0], [prim_pts - 1]),), shape=(1, prim_pts),
+                            ([1], ([0], [prim_pts - 1])), shape=(1, prim_pts)
                         )
                         additive = dxN * bcs[child.id][symbol.side][0]
 
@@ -727,7 +725,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         b = -(dxN ** 2) / (2 * dxN * dxNm1 + dxNm1 ** 2)
                         alpha = dxN * (dxN + dxNm1) / (2 * dxN + dxNm1)
                         sub_matrix = csr_matrix(
-                            ([b, a], ([0, 0], [prim_pts - 2, prim_pts - 1]),),
+                            ([b, a], ([0, 0], [prim_pts - 2, prim_pts - 1])),
                             shape=(1, prim_pts),
                         )
 
@@ -755,9 +753,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
 
         elif isinstance(symbol, pybamm.BoundaryGradient):
 
-            if use_bcs and pybamm.has_bc_of_form(
-                child, symbol.side, bcs, "Neumann"
-            ):
+            if use_bcs and pybamm.has_bc_of_form(child, symbol.side, bcs, "Neumann"):
                 # just use the value from the bc: f'(x*)
                 sub_matrix = csr_matrix((1, prim_pts))
                 additive = bcs[child.id][symbol.side][0]

From 6b675f67cfd15ee55dd677222e142e0d6d3df5e7 Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Sat, 30 Nov 2019 12:18:54 -0800
Subject: [PATCH 05/18] #743 fix jacobian for input parameter

---
 pybamm/expression_tree/input_parameter.py     |  7 ++
 pybamm/expression_tree/symbol.py              | 16 ++--
 pybamm/solvers/dae_solver.py                  | 77 ++++++++++++++-----
 tests/unit/test_solvers/test_casadi_solver.py | 24 ++++++
 tests/unit/test_solvers/test_scipy_solver.py  | 25 ++++++
 5 files changed, 121 insertions(+), 28 deletions(-)

diff --git a/pybamm/expression_tree/input_parameter.py b/pybamm/expression_tree/input_parameter.py
index dab5725356..67b8a87551 100644
--- a/pybamm/expression_tree/input_parameter.py
+++ b/pybamm/expression_tree/input_parameter.py
@@ -32,8 +32,15 @@ def evaluate_for_shape(self):
         """
         return np.nan
 
+    def _jac(self, variable):
+        """ See :meth:`pybamm.Symbol._jac()`. """
+        return pybamm.Scalar(0)
+
     def evaluate(self, t=None, y=None, u=None, known_evals=None):
         # u should be a dictionary
+        # convert 'None' to empty dictionary for more informative error
+        if u is None:
+            u = {}
         if not isinstance(u, dict):
             # if the special input "shape test" is passed, just return 1
             if u == "shape test":
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index dbda1aead0..ca902f5172 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -520,8 +520,8 @@ def is_constant(self):
     def evaluate_ignoring_errors(self):
         """
         Evaluates the expression. If a node exists in the tree that cannot be evaluated
-        as a scalar or vectr (e.g. Parameter, Variable, StateVector), then None is
-        returned. Otherwise the result of the evaluation is given
+        as a scalar or vector (e.g. Parameter, Variable, StateVector, InputParameter),
+        then None is returned. Otherwise the result of the evaluation is given
 
         See Also
         --------
@@ -529,22 +529,18 @@ def evaluate_ignoring_errors(self):
 
         """
         try:
-            result = self.evaluate(t=0)
+            result = self.evaluate(t=0, u="shape test")
         except NotImplementedError:
-            # return false if NotImplementedError is raised
+            # return None if NotImplementedError is raised
             # (there is a e.g. Parameter, Variable, ... in the tree)
             return None
         except TypeError as error:
-            # return false if specific TypeError is raised
+            # return None if specific TypeError is raised
             # (there is a e.g. StateVector in the tree)
-            if error.args[0] in [
-                "StateVector cannot evaluate input 'y=None'",
-                "inputs u should be a dictionary",
-            ]:
+            if error.args[0] == "StateVector cannot evaluate input 'y=None'":
                 return None
             else:
                 raise error
-
         return result
 
     def evaluates_to_number(self):
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index ba1f24004f..182eb2b25e 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -211,7 +211,7 @@ def get_event_class(event):
         self.event_funs = [get_event_class(event) for event in events.values()]
         self.jacobian = jacobian
 
-    def set_up_casadi(self, model):
+    def set_up_casadi(self, model, inputs):
         """Convert model to casadi format and use their inbuilt functionalities.
 
         Parameters
@@ -219,6 +219,9 @@ def set_up_casadi(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
+
         """
         # Convert model attributes to casadi
         t_casadi = casadi.MX.sym("t")
@@ -235,39 +238,49 @@ def set_up_casadi(self, model):
             y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
         else:
             y_casadi_w_ext = y_casadi
+        u_casadi = {name: casadi.MX.sym(name) for name in inputs.keys()}
 
         pybamm.logger.info("Converting RHS to CasADi")
-        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi_w_ext)
+        concatenated_rhs = model.concatenated_rhs.to_casadi(
+            t_casadi, y_casadi_w_ext, u_casadi
+        )
         pybamm.logger.info("Converting algebraic to CasADi")
         concatenated_algebraic = model.concatenated_algebraic.to_casadi(
-            t_casadi, y_casadi_w_ext
+            t_casadi, y_casadi_w_ext, u_casadi
         )
         all_states = casadi.vertcat(concatenated_rhs, concatenated_algebraic)
         pybamm.logger.info("Converting events to CasADi")
         casadi_events = {
-            name: event.to_casadi(t_casadi, y_casadi_w_ext)
+            name: event.to_casadi(t_casadi, y_casadi_w_ext, u_casadi)
             for name, event in model.events.items()
         }
 
         # Create functions to evaluate rhs and algebraic
+        u_casadi_stacked = casadi.vertcat(*[u for u in u_casadi.values()])
         concatenated_rhs_fn = casadi.Function(
-            "rhs", [t_casadi, y_casadi_w_ext], [concatenated_rhs]
+            "rhs", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [concatenated_rhs]
         )
         concatenated_algebraic_fn = casadi.Function(
-            "algebraic", [t_casadi, y_casadi_w_ext], [concatenated_algebraic]
+            "algebraic",
+            [t_casadi, y_casadi_w_ext, u_casadi_stacked],
+            [concatenated_algebraic],
+        )
+        all_states_fn = casadi.Function(
+            "all", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [all_states]
         )
-        all_states_fn = casadi.Function("all", [t_casadi, y_casadi_w_ext], [all_states])
 
         if model.use_jacobian:
 
             pybamm.logger.info("Calculating jacobian")
             casadi_jac = casadi.jacobian(all_states, y_casadi)
             casadi_jac_fn = casadi.Function(
-                "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac]
+                "jacobian", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [casadi_jac]
             )
             casadi_jac_alg = casadi.jacobian(concatenated_algebraic, y_casadi)
             casadi_jac_alg_fn = casadi.Function(
-                "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac_alg]
+                "jacobian",
+                [t_casadi, y_casadi_w_ext, u_casadi_stacked],
+                [casadi_jac_alg],
             )
 
             jacobian = JacobianCasadi(casadi_jac_fn)
@@ -302,7 +315,7 @@ def jacobian_alg(t, y):
         # Create event-dependent function to evaluate events
         def get_event_class(event):
             casadi_event_fn = casadi.Function(
-                "event", [t_casadi, y_casadi_w_ext], [event]
+                "event", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [event]
             )
             return EvalEvent(casadi_event_fn)
 
@@ -321,6 +334,26 @@ def get_event_class(event):
         self.casadi_rhs = concatenated_rhs_fn
         self.casadi_algebraic = concatenated_algebraic_fn
 
+    def set_inputs_and_external(self, inputs):
+        """
+        Set values that are controlled externally, such as external variables and input
+        parameters
+
+        Parameters
+        ----------
+        inputs : dict
+            Any input parameters to pass to the model when solving
+
+        """
+        self.dydt.set_pad_ext(self.y_pad, self.y_ext)
+        self.dydt.set_inputs(inputs)
+        for evnt in self.event_funs:
+            evnt.set_pad_ext(self.y_pad, self.y_ext)
+            evnt.set_inputs(inputs)
+        if self.jacobian:
+            self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+            self.jacobian.set_inputs(inputs)
+
     def calculate_consistent_initial_conditions(
         self, rhs, algebraic, y0_guess, jac=None
     ):
@@ -445,19 +478,16 @@ class SolverCallable:
     "A class that will be called by the solver when integrating"
     y_pad = None
     y_ext = None
-    _inputs = {}
+    inputs = {}
+    inputs_casadi = casadi.DM()
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
         self.y_ext = y_ext
 
-    @property
-    def inputs(self):
-        return self._inputs
-
-    @inputs.setter
-    def inputs(self, value):
-        self._inputs = value
+    def set_inputs(self, inputs):
+        self.inputs = inputs
+        self.inputs_casadi = casadi.vertcat(*[x for x in inputs.values()])
 
 
 class Rhs(SolverCallable):
@@ -560,6 +590,17 @@ def __call__(self, t, y):
         return self.event_fn(t, y)
 
 
+class EvalEventCasadi(EvalEvent):
+    "Returns information about events at time t and state y"
+
+    def __init__(self, event_fn):
+        self.event_fn = event_fn
+
+    def __call__(self, t, y):
+        y = y[:, np.newaxis]
+        y = add_external(y, self.y_pad, self.y_ext)
+        return self.event_fn(t, y, self.inputs_casadi)
+
 class Jacobian:
     "Returns information about the jacobian at time t and state y"
 
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 13d66ff0d2..3bf42901df 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -184,6 +184,30 @@ def test_model_step(self):
         solution = solver.solve(model, t_eval)
         np.testing.assert_allclose(solution.y[0], step_sol.y[0])
 
+    def test_model_solver_with_inputs(self):
+        # Create model
+        model = pybamm.BaseModel()
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=domain)
+        model.rhs = {var: -pybamm.InputParameter("rate") * var}
+        model.initial_conditions = {var: 1}
+        model.events = {"var=0.5": pybamm.min(var - 0.5)}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+        # Solve
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+        t_eval = np.linspace(0, 10, 100)
+        solution = solver.solve(model, t_eval, inputs={"rate": 0.1})
+        self.assertLess(len(solution.t), len(t_eval))
+        np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+        np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 2ae309617a..ff10eadc8f 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -272,6 +272,31 @@ def test_model_step_python(self):
         solution = solver.solve(model, t_eval)
         np.testing.assert_allclose(solution.y[0], step_sol.y[0])
 
+    def test_model_solver_with_inputs(self):
+        # Create model
+        model = pybamm.BaseModel()
+        model.convert_to_format = "python"
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=domain)
+        model.rhs = {var: -pybamm.InputParameter("rate") * var}
+        model.initial_conditions = {var: 1}
+        model.events = {"var=0.5": pybamm.min(var - 0.5)}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+        # Solve
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+        t_eval = np.linspace(0, 10, 100)
+        solution = solver.solve(model, t_eval, inputs={"rate": 0.1})
+        self.assertLess(len(solution.t), len(t_eval))
+        np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+        np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
+
     def test_model_solver_with_event_with_casadi(self):
         # Create model
         model = pybamm.BaseModel()

From c80f38a4f9d247f829052fcbe7b35304d2e10d43 Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Sat, 30 Nov 2019 12:28:28 -0800
Subject: [PATCH 06/18] #743 get scipy solver working with python

---
 pybamm/expression_tree/operations/evaluate.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/pybamm/expression_tree/operations/evaluate.py b/pybamm/expression_tree/operations/evaluate.py
index caffad205b..6558357c33 100644
--- a/pybamm/expression_tree/operations/evaluate.py
+++ b/pybamm/expression_tree/operations/evaluate.py
@@ -175,6 +175,9 @@ def find_symbols(symbol, constant_symbols, variable_symbols):
     elif isinstance(symbol, pybamm.Time):
         symbol_str = "t"
 
+    elif isinstance(symbol, pybamm.InputParameter):
+        symbol_str = "u['{}']".format(symbol.name)
+
     else:
         raise NotImplementedError(
             "Not implemented for a symbol of type '{}'".format(type(symbol))

From 544c0dda73a4756cc4271bab73f6d5783a64f6eb Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Sat, 30 Nov 2019 12:56:11 -0800
Subject: [PATCH 07/18] #743 fix scikits solvers without inputs

---
 pybamm/solvers/base_solver.py   |  2 +-
 pybamm/solvers/casadi_solver.py | 20 ++++---
 pybamm/solvers/dae_solver.py    | 96 ++++++++++++++++++---------------
 pybamm/solvers/ode_solver.py    |  2 +-
 4 files changed, 70 insertions(+), 50 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index de9d72a307..9db1573756 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -157,7 +157,7 @@ def step(self, model, dt, npts=2, log=True, external_variables=None, inputs=None
             if model.convert_to_format == "casadi" or isinstance(
                 self, pybamm.CasadiSolver
             ):
-                self.set_up_casadi(model)
+                self.set_up_casadi(model, inputs)
             else:
                 pybamm.logger.debug(
                     "Start stepping {} with {}".format(model.name, self.name)
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index b5c0b5c0c5..413eca8069 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -68,7 +68,7 @@ def __init__(
         self.extra_options = extra_options
         self.name = "CasADi solver ({}) with '{}' mode".format(method, mode)
 
-    def solve(self, model, t_eval):
+    def solve(self, model, t_eval, inputs=None):
         """
         Execute the solver setup and calculate the solution of the model at
         specified times.
@@ -80,7 +80,9 @@ def solve(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
-
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
+        
         Raises
         ------
         :class:`pybamm.ValueError`
@@ -99,7 +101,8 @@ def solve(self, model, t_eval):
         elif self.mode == "safe":
             # Step-and-check
             timer = pybamm.Timer()
-            self.set_up_casadi(model)
+            inputs = inputs or {}
+            self.set_up_casadi(model, inputs)
             set_up_time = timer.time()
             init_event_signs = np.sign(
                 np.concatenate([event(0, self.y0) for event in self.event_funs])
@@ -175,7 +178,7 @@ def solve(self, model, t_eval):
             )
             return solution
 
-    def compute_solution(self, model, t_eval):
+    def compute_solution(self, model, t_eval, inputs):
         """Calculate the solution of the model at specified times. In this class, we
         overwrite the behaviour of :class:`pybamm.DaeSolver`, since CasADi requires
         slightly different syntax.
@@ -187,7 +190,9 @@ def compute_solution(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
-
+        inputs : dict
+            Any input parameters to pass to the model when solving
+        
         """
         timer = pybamm.Timer()
 
@@ -198,6 +203,7 @@ def compute_solution(self, model, t_eval):
             self.casadi_algebraic,
             self.y0,
             t_eval,
+            inputs,
             mass_matrix=model.mass_matrix.entries,
         )
         solve_time = timer.time() - solve_start_time
@@ -207,7 +213,7 @@ def compute_solution(self, model, t_eval):
 
         return solution, solve_time, termination
 
-    def integrate_casadi(self, rhs, algebraic, y0, t_eval, mass_matrix=None):
+    def integrate_casadi(self, rhs, algebraic, y0, t_eval, inputs, mass_matrix=None):
         """
         Solve a DAE model defined by residuals with initial conditions y0.
 
@@ -220,6 +226,8 @@ def integrate_casadi(self, rhs, algebraic, y0, t_eval, mass_matrix=None):
             The initial conditions
         t_eval : numeric type
             The times at which to compute the solution
+        inputs : dict
+            Any input parameters to pass to the model when solving
         mass_matrix : array_like, optional
             The (sparse) mass matrix for the chosen spatial method. This is only passed
             to check that the mass matrix is diagonal with 1s for the odes and 0s for
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 182eb2b25e..b3d5d0f730 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -82,14 +82,8 @@ def compute_solution(self, model, t_eval, inputs=None):
         """
         timer = pybamm.Timer()
 
-        # update y_pad and y_ext
-        self.rhs.set_pad_ext(self.y_pad, self.y_ext)
-        self.algebraic.set_pad_ext(self.y_pad, self.y_ext)
-        self.residuals.set_pad_ext(self.y_pad, self.y_ext)
-        for evnt in self.event_funs:
-            evnt.set_pad_ext(self.y_pad, self.y_ext)
-        if self.jacobian:
-            self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+        # Set inputs and external
+        self.set_inputs_and_external(inputs)
 
         solve_start_time = timer.time()
         pybamm.logger.info("Calling DAE solver")
@@ -101,7 +95,6 @@ def compute_solution(self, model, t_eval, inputs=None):
             mass_matrix=model.mass_matrix.entries,
             jacobian=self.jacobian,
             model=model,
-            inputs=inputs,
         )
 
         solve_time = timer.time() - solve_start_time
@@ -161,11 +154,7 @@ def set_up(self, model):
                 jac = pybamm.EvaluatorPython(jac)
 
             jacobian = Jacobian(jac.evaluate)
-
-            def jacobian_alg(t, y):
-                y = y[:, np.newaxis]
-                y = add_external(y, self.y_pad, self.y_ext)
-                return jac_algebraic.evaluate(t, y, known_evals={})[0]
+            jacobian_alg = Jacobian(jac_algebraic.evaluate)
 
         else:
             jacobian = None
@@ -284,11 +273,7 @@ def set_up_casadi(self, model, inputs):
             )
 
             jacobian = JacobianCasadi(casadi_jac_fn)
-
-            def jacobian_alg(t, y):
-                y = y[:, np.newaxis]
-                y = add_external(y, self.y_pad, self.y_ext)
-                return casadi_jac_alg_fn(t, y)
+            jacobian_alg = JacobianAlgebraicCasadi(casadi_jac_alg_fn)
 
         else:
             jacobian = None
@@ -317,7 +302,7 @@ def get_event_class(event):
             casadi_event_fn = casadi.Function(
                 "event", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [event]
             )
-            return EvalEvent(casadi_event_fn)
+            return EvalEventCasadi(casadi_event_fn)
 
         # Add the solver attributes
         # Note: these are the (possibly) converted to python version rhs, algebraic
@@ -345,8 +330,12 @@ def set_inputs_and_external(self, inputs):
             Any input parameters to pass to the model when solving
 
         """
-        self.dydt.set_pad_ext(self.y_pad, self.y_ext)
-        self.dydt.set_inputs(inputs)
+        self.rhs.set_pad_ext(self.y_pad, self.y_ext)
+        self.rhs.set_inputs(inputs)
+        self.algebraic.set_pad_ext(self.y_pad, self.y_ext)
+        self.algebraic.set_inputs(inputs)
+        self.residuals.set_pad_ext(self.y_pad, self.y_ext)
+        self.residuals.set_inputs(inputs)
         for evnt in self.event_funs:
             evnt.set_pad_ext(self.y_pad, self.y_ext)
             evnt.set_inputs(inputs)
@@ -499,7 +488,7 @@ def __init__(self, concatenated_rhs_fn):
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.concatenated_rhs_fn(t, y, known_evals={})[0][:, 0]
+        return self.concatenated_rhs_fn(t, y, self.inputs, known_evals={})[0][:, 0]
 
 
 class RhsCasadi(Rhs):
@@ -508,7 +497,7 @@ class RhsCasadi(Rhs):
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.concatenated_rhs_fn(t, y).full()[:, 0]
+        return self.concatenated_rhs_fn(t, y, self.inputs_casadi).full()[:, 0]
 
 
 class Algebraic(SolverCallable):
@@ -520,7 +509,9 @@ def __init__(self, concatenated_algebraic_fn):
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.concatenated_algebraic_fn(t, y, known_evals={})[0][:, 0]
+        return self.concatenated_algebraic_fn(t, y, self.inputs, known_evals={})[0][
+            :, 0
+        ]
 
 
 class AlgebraicCasadi(Algebraic):
@@ -529,7 +520,7 @@ class AlgebraicCasadi(Algebraic):
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.concatenated_algebraic_fn(t, y).full()[:, 0]
+        return self.concatenated_algebraic_fn(t, y, self.inputs_casadi).full()[:, 0]
 
 
 class Residuals(SolverCallable):
@@ -547,9 +538,13 @@ def __call__(self, t, y, ydot):
         )
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        rhs_eval, known_evals = self.concatenated_rhs_fn(t, y, known_evals={})
+        rhs_eval, known_evals = self.concatenated_rhs_fn(
+            t, y, self.inputs, known_evals={}
+        )
         # reuse known_evals
-        alg_eval = self.concatenated_algebraic_fn(t, y, known_evals=known_evals)[0]
+        alg_eval = self.concatenated_algebraic_fn(
+            t, y, self.inputs, known_evals=known_evals
+        )[0]
         # turn into 1D arrays
         rhs_eval = rhs_eval[:, 0]
         alg_eval = alg_eval[:, 0]
@@ -570,24 +565,20 @@ def __call__(self, t, y, ydot):
         )
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        states_eval = self.all_states_fn(t, y).full()[:, 0]
+        states_eval = self.all_states_fn(t, y, self.inputs_casadi).full()[:, 0]
         return states_eval - self.mass_matrix @ ydot
 
 
-class EvalEvent:
+class EvalEvent(SolverCallable):
     "Returns information about events at time t and state y"
 
     def __init__(self, event_fn):
         self.event_fn = event_fn
 
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.event_fn(t, y)
+        return self.event_fn(t, y, self.inputs)
 
 
 class EvalEventCasadi(EvalEvent):
@@ -601,20 +592,17 @@ def __call__(self, t, y):
         y = add_external(y, self.y_pad, self.y_ext)
         return self.event_fn(t, y, self.inputs_casadi)
 
-class Jacobian:
+
+class Jacobian(SolverCallable):
     "Returns information about the jacobian at time t and state y"
 
     def __init__(self, jac_fn):
         self.jac_fn = jac_fn
 
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.jac_fn(t, y, known_evals={})[0]
+        return self.jac_fn(t, y, self.inputs, known_evals={})[0]
 
 
 class JacobianCasadi(Jacobian):
@@ -623,4 +611,28 @@ class JacobianCasadi(Jacobian):
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.jac_fn(t, y)
+        return self.jac_fn(t, y, self.inputs_casadi)
+
+
+class JacobianAlgebraic(SolverCallable):
+    "Returns information about the algebraic part of the jacobian at time t and state y"
+
+    def __init__(self, jac_alg_fn):
+        self.jac_fn = jac_alg_fn
+
+    def __call__(self, t, y):
+        y = y[:, np.newaxis]
+        y = add_external(y, self.y_pad, self.y_ext)
+        return self.jac_fn(t, y, self.inputs, known_evals={})[0]
+
+
+class JacobianAlgebraicCasadi(JacobianAlgebraic):
+    """
+    Returns information about the algebraic part of the jacobian at time t and
+    state y, with CasADi
+    """
+
+    def __call__(self, t, y):
+        y = y[:, np.newaxis]
+        y = add_external(y, self.y_pad, self.y_ext)
+        return self.jac_fn(t, y, self.inputs_casadi)
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 9810c7c36b..137e6c82db 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -355,4 +355,4 @@ class JacobianCasadi(Jacobian):
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.jac_fn(t, y, self.inputs)
+        return self.jac_fn(t, y, self.inputs_casadi)

From 64ffdd8b8aa1d95db3f6098a04370c5a4cb6a240 Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Sat, 30 Nov 2019 16:27:26 -0800
Subject: [PATCH 08/18] #743 scikits solvers working with inputs

---
 pybamm/solvers/base_solver.py                 |  4 +--
 pybamm/solvers/dae_solver.py                  | 22 ++++++++++++---
 pybamm/solvers/ode_solver.py                  |  5 +++-
 .../unit/test_solvers/test_scikits_solvers.py | 28 ++++++++++++++++++-
 4 files changed, 51 insertions(+), 8 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 9db1573756..c13c8ddc4d 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -83,7 +83,7 @@ def solve(self, model, t_eval, inputs=None):
         if model.convert_to_format == "casadi" or isinstance(self, pybamm.CasadiSolver):
             self.set_up_casadi(model, inputs)
         else:
-            self.set_up(model)
+            self.set_up(model, inputs)
         set_up_time = timer.time() - start_time
 
         # Solve
@@ -162,7 +162,7 @@ def step(self, model, dt, npts=2, log=True, external_variables=None, inputs=None
                 pybamm.logger.debug(
                     "Start stepping {} with {}".format(model.name, self.name)
                 )
-                self.set_up(model)
+                self.set_up(model, inputs)
             self.t = 0.0
             set_up_time = timer.time()
 
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index b3d5d0f730..bc24ff77f5 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -104,7 +104,7 @@ def compute_solution(self, model, t_eval, inputs=None):
 
         return solution, solve_time, termination
 
-    def set_up(self, model):
+    def set_up(self, model, inputs):
         """Unpack model, perform checks, simplify and calculate jacobian.
 
         Parameters
@@ -112,6 +112,9 @@ def set_up(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict
+            Any input parameters to pass to the model when solving
+
         """
         # create simplified rhs, algebraic and event expressions
         concatenated_rhs = model.concatenated_rhs
@@ -155,6 +158,8 @@ def set_up(self, model):
 
             jacobian = Jacobian(jac.evaluate)
             jacobian_alg = Jacobian(jac_algebraic.evaluate)
+            jacobian_alg.set_pad_ext(self.y_pad, self.y_ext)
+            jacobian_alg.set_inputs(inputs)
 
         else:
             jacobian = None
@@ -174,6 +179,11 @@ def set_up(self, model):
         rhs = Rhs(concatenated_rhs.evaluate)
         algebraic = Algebraic(concatenated_algebraic.evaluate)
 
+        rhs.set_pad_ext(self.y_pad, self.y_ext)
+        rhs.set_inputs(inputs)
+        algebraic.set_pad_ext(self.y_pad, self.y_ext)
+        algebraic.set_inputs(inputs)
+        
         if len(model.algebraic) > 0:
             y0 = self.calculate_consistent_initial_conditions(
                 rhs,
@@ -208,7 +218,7 @@ def set_up_casadi(self, model, inputs):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
-        inputs : dict, optional
+        inputs : dict
             Any input parameters to pass to the model when solving
 
         """
@@ -217,9 +227,11 @@ def set_up_casadi(self, model, inputs):
         y0 = model.concatenated_initial_conditions
         y0 = add_external(y0, self.y_pad, self.y_ext)
 
-        y_diff = casadi.MX.sym("y_diff", len(model.concatenated_rhs.evaluate(0, y0)))
+        y_diff = casadi.MX.sym(
+            "y_diff", len(model.concatenated_rhs.evaluate(0, y0, inputs))
+        )
         y_alg = casadi.MX.sym(
-            "y_alg", len(model.concatenated_algebraic.evaluate(0, y0))
+            "y_alg", len(model.concatenated_algebraic.evaluate(0, y0, inputs))
         )
         y_casadi = casadi.vertcat(y_diff, y_alg)
         if self.y_pad is not None:
@@ -283,7 +295,9 @@ def set_up_casadi(self, model, inputs):
         algebraic = AlgebraicCasadi(concatenated_algebraic_fn)
 
         rhs.set_pad_ext(self.y_pad, self.y_ext)
+        rhs.set_inputs(inputs)
         algebraic.set_pad_ext(self.y_pad, self.y_ext)
+        algebraic.set_inputs(inputs)
 
         if len(model.algebraic) > 0:
 
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 137e6c82db..55c06bb9c3 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -61,7 +61,7 @@ def compute_solution(self, model, t_eval, inputs=None):
 
         return solution, solve_time, termination
 
-    def set_up(self, model):
+    def set_up(self, model, inputs):
         """Unpack model, perform checks, simplify and calculate jacobian.
 
         Parameters
@@ -69,6 +69,9 @@ def set_up(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
+
 
         Raises
         ------
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 7ac16b5af3..5a4c9b377f 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -547,7 +547,6 @@ def test_model_solver_dae_events_python(self):
         np.testing.assert_array_less(solution.y[-1], 2.5)
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
         np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
-
     def test_model_solver_dae_with_jacobian_python(self):
         model = pybamm.BaseModel()
         model.convert_to_format = "python"
@@ -717,6 +716,33 @@ def test_model_solver_dae_events_casadi(self):
             np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
             np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
+    def test_model_solver_dae_inputs_events(self):
+        # Create model
+        for format in ["python", "casadi"]:
+            model = pybamm.BaseModel()
+            model.convert_to_format = format
+            whole_cell = ["negative electrode", "separator", "positive electrode"]
+            var1 = pybamm.Variable("var1", domain=whole_cell)
+            var2 = pybamm.Variable("var2", domain=whole_cell)
+            model.rhs = {var1: pybamm.InputParameter("rate 1") * var1}
+            model.algebraic = {var2: pybamm.InputParameter("rate 2") * var1 - var2}
+            model.initial_conditions = {var1: 1, var2: 2}
+            model.events = {
+                "var1 = 1.5": pybamm.min(var1 - 1.5),
+                "var2 = 2.5": pybamm.min(var2 - 2.5),
+            }
+            disc = get_discretisation_for_testing()
+            disc.process_model(model)
+
+            # Solve
+            solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+            t_eval = np.linspace(0, 5, 100)
+            solution = solver.solve(model, t_eval, inputs={"rate 1": 0.1, "rate 2": 2})
+            np.testing.assert_array_less(solution.y[0], 1.5)
+            np.testing.assert_array_less(solution.y[-1], 2.5)
+            np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
+            np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
+
     def test_solve_ode_model_with_dae_solver_casadi(self):
         model = pybamm.BaseModel()
         model.convert_to_format = "casadi"

From 166681ebe69c685d4c250addfc616f9aa9b30bb3 Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Sat, 30 Nov 2019 16:44:27 -0800
Subject: [PATCH 09/18] #743 casadi solver works without events

---
 pybamm/solvers/casadi_solver.py               | 23 +++++++++----------
 tests/unit/test_solvers/test_casadi_solver.py | 20 +++++++++++++---
 2 files changed, 28 insertions(+), 15 deletions(-)

diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index 413eca8069..237a77a259 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -213,7 +213,9 @@ def compute_solution(self, model, t_eval, inputs):
 
         return solution, solve_time, termination
 
-    def integrate_casadi(self, rhs, algebraic, y0, t_eval, inputs, mass_matrix=None):
+    def integrate_casadi(
+        self, rhs, algebraic, y0, t_eval, inputs=None, mass_matrix=None
+    ):
         """
         Solve a DAE model defined by residuals with initial conditions y0.
 
@@ -226,13 +228,14 @@ def integrate_casadi(self, rhs, algebraic, y0, t_eval, inputs, mass_matrix=None)
             The initial conditions
         t_eval : numeric type
             The times at which to compute the solution
-        inputs : dict
+        inputs : dict, optional
             Any input parameters to pass to the model when solving
         mass_matrix : array_like, optional
             The (sparse) mass matrix for the chosen spatial method. This is only passed
             to check that the mass matrix is diagonal with 1s for the odes and 0s for
             the algebraic equations, as CasADi does not allow to pass mass matrices.
         """
+        inputs = inputs or {}
         options = {
             "grid": t_eval,
             "reltol": self.rtol,
@@ -246,19 +249,15 @@ def integrate_casadi(self, rhs, algebraic, y0, t_eval, inputs, mass_matrix=None)
 
         # set up and solve
         t = casadi.MX.sym("t")
-        y_diff = casadi.MX.sym("y_diff", rhs(0, y0).shape[0])
+        u = casadi.vertcat(*[x for x in inputs.values()])
+        y_diff = casadi.MX.sym("y_diff", rhs(0, y0, u).shape[0])
+        problem = {"t": t, "x": y_diff}
         if algebraic is None:
-            problem = {"t": t, "x": y_diff, "ode": rhs(t, y_diff)}
+            problem.update({"ode": rhs(t, y_diff, u)})
         else:
-            y_alg = casadi.MX.sym("y_alg", algebraic(0, y0).shape[0])
+            y_alg = casadi.MX.sym("y_alg", algebraic(0, y0, u).shape[0])
             y = casadi.vertcat(y_diff, y_alg)
-            problem = {
-                "t": t,
-                "x": y_diff,
-                "z": y_alg,
-                "ode": rhs(t, y),
-                "alg": algebraic(t, y),
-            }
+            problem.update({"z": y_alg, "ode": rhs(t, y, u), "alg": algebraic(t, y, u)})
         integrator = casadi.integrator("F", self.method, problem, options)
         try:
             # Try solving
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 3bf42901df..74cee9d218 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -16,8 +16,9 @@ def test_integrate(self):
 
         t = casadi.MX.sym("t")
         y = casadi.MX.sym("y")
+        u = casadi.MX.sym("u")
         constant_growth = casadi.MX(0.5)
-        rhs = casadi.Function("rhs", [t, y], [constant_growth])
+        rhs = casadi.Function("rhs", [t, y, u], [constant_growth])
 
         y0 = np.array([0])
         t_eval = np.linspace(0, 1, 100)
@@ -29,7 +30,7 @@ def test_integrate(self):
         solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8, method="cvodes")
 
         exponential_decay = -0.1 * y
-        rhs = casadi.Function("rhs", [t, y], [exponential_decay])
+        rhs = casadi.Function("rhs", [t, y, u], [exponential_decay])
 
         y0 = np.array([1])
         t_eval = np.linspace(0, 1, 100)
@@ -37,18 +38,31 @@ def test_integrate(self):
         np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
         self.assertEqual(solution.termination, "final time")
 
+        # Exponential decay with input
+        solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8, method="cvodes")
+
+        exponential_decay = -u * y
+        rhs = casadi.Function("rhs", [t, y, u], [exponential_decay])
+
+        y0 = np.array([1])
+        t_eval = np.linspace(0, 1, 100)
+        solution = solver.integrate_casadi(rhs, None, y0, t_eval, inputs={"u": 0.1})
+        np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
+        self.assertEqual(solution.termination, "final time")
+
     def test_integrate_failure(self):
         # Turn off warnings to ignore sqrt error
         warnings.simplefilter("ignore")
 
         t = casadi.MX.sym("t")
         y = casadi.MX.sym("y")
+        u = casadi.MX.sym("u")
         sqrt_decay = -np.sqrt(y)
 
         y0 = np.array([1])
         t_eval = np.linspace(0, 3, 100)
         solver = pybamm.CasadiSolver()
-        rhs = casadi.Function("rhs", [t, y], [sqrt_decay])
+        rhs = casadi.Function("rhs", [t, y, u], [sqrt_decay])
         # Expect solver to fail when y goes negative
         with self.assertRaises(pybamm.SolverError):
             solver.integrate_casadi(rhs, None, y0, t_eval)

From 98a76f66b55a4eeb7ca181e3ba9339793d9eaa4a Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Sat, 30 Nov 2019 16:55:56 -0800
Subject: [PATCH 10/18] #743 fixing tests

---
 pybamm/expression_tree/functions.py             |  2 +-
 pybamm/processed_variable.py                    | 10 +++++-----
 pybamm/solvers/casadi_solver.py                 |  4 ++--
 pybamm/solvers/dae_solver.py                    |  2 +-
 tests/unit/test_solvers/test_scikits_solvers.py |  1 +
 5 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index c2dbc8a95a..0014b58fbb 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -159,7 +159,7 @@ def evaluate(self, t=None, y=None, u=None, known_evals=None):
                 evaluated_children = [None] * len(self.children)
                 for i, child in enumerate(self.children):
                     evaluated_children[i], known_evals = child.evaluate(
-                        t, y, known_evals
+                        t, y, known_evals=known_evals
                     )
                 known_evals[self.id] = self._function_evaluate(evaluated_children)
             return known_evals[self.id], known_evals
diff --git a/pybamm/processed_variable.py b/pybamm/processed_variable.py
index a4db4a6462..f74d4b5011 100644
--- a/pybamm/processed_variable.py
+++ b/pybamm/processed_variable.py
@@ -88,7 +88,7 @@ def __init__(
 
         if self.known_evals:
             self.base_eval, self.known_evals[t_sol[0]] = base_variable.evaluate(
-                t_sol[0], u_sol[:, 0], self.known_evals[t_sol[0]]
+                t_sol[0], u_sol[:, 0], known_evals=self.known_evals[t_sol[0]]
             )
         else:
             self.base_eval = base_variable.evaluate(t_sol[0], u_sol[:, 0])
@@ -131,7 +131,7 @@ def initialise_1D(self):
             t = self.t_sol[idx]
             if self.known_evals:
                 entries[idx], self.known_evals[t] = self.base_variable.evaluate(
-                    t, self.u_sol[:, idx], self.known_evals[t]
+                    t, self.u_sol[:, idx], known_evals=self.known_evals[t]
                 )
             else:
                 entries[idx] = self.base_variable.evaluate(t, self.u_sol[:, idx])
@@ -158,7 +158,7 @@ def initialise_2D(self):
             u = self.u_sol[:, idx]
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, self.known_evals[t]
+                    t, u, known_evals=self.known_evals[t]
                 )
                 entries[:, idx] = eval_and_known_evals[0][:, 0]
                 self.known_evals[t] = eval_and_known_evals[1]
@@ -301,7 +301,7 @@ def initialise_3D(self):
             u = self.u_sol[:, idx]
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, self.known_evals[t]
+                    t, u, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(
                     eval_and_known_evals[0],
@@ -366,7 +366,7 @@ def initialise_3D_scikit_fem(self):
             u = self.u_sol[:, idx]
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, self.known_evals[t]
+                    t, u, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(eval_and_known_evals[0], [len_y, len_z])
                 self.known_evals[t] = eval_and_known_evals[1]
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index 237a77a259..b008ca6a71 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -82,7 +82,7 @@ def solve(self, model, t_eval, inputs=None):
             The times at which to compute the solution
         inputs : dict, optional
             Any input parameters to pass to the model when solving
-        
+
         Raises
         ------
         :class:`pybamm.ValueError`
@@ -192,7 +192,7 @@ def compute_solution(self, model, t_eval, inputs):
             The times at which to compute the solution
         inputs : dict
             Any input parameters to pass to the model when solving
-        
+
         """
         timer = pybamm.Timer()
 
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index bc24ff77f5..3441d26284 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -183,7 +183,7 @@ def set_up(self, model, inputs):
         rhs.set_inputs(inputs)
         algebraic.set_pad_ext(self.y_pad, self.y_ext)
         algebraic.set_inputs(inputs)
-        
+
         if len(model.algebraic) > 0:
             y0 = self.calculate_consistent_initial_conditions(
                 rhs,
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 5a4c9b377f..aa55b3edd9 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -547,6 +547,7 @@ def test_model_solver_dae_events_python(self):
         np.testing.assert_array_less(solution.y[-1], 2.5)
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
         np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
+
     def test_model_solver_dae_with_jacobian_python(self):
         model = pybamm.BaseModel()
         model.convert_to_format = "python"

From 40db53b10a284a3a2aafda5ba775c15d6996ffbb Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Sat, 30 Nov 2019 17:00:58 -0800
Subject: [PATCH 11/18] #743 add inputs to processed variable

---
 pybamm/processed_variable.py | 35 ++++++++++++++++++++++++-----------
 1 file changed, 24 insertions(+), 11 deletions(-)

diff --git a/pybamm/processed_variable.py b/pybamm/processed_variable.py
index f74d4b5011..5081276be6 100644
--- a/pybamm/processed_variable.py
+++ b/pybamm/processed_variable.py
@@ -7,7 +7,9 @@
 import scipy.interpolate as interp
 
 
-def post_process_variables(variables, t_sol, u_sol, mesh=None, interp_kind="linear"):
+def post_process_variables(
+    variables, t_sol, u_sol, mesh=None, inputs=None, interp_kind="linear"
+):
     """
     Post-process all variables in a model
 
@@ -23,6 +25,8 @@ def post_process_variables(variables, t_sol, u_sol, mesh=None, interp_kind="line
     mesh : :class:`pybamm.Mesh`
         The mesh used to solve, used here to calculate the reference x values for
         interpolation
+    inputs : dict, optional
+        Any input parameters to pass to the model
     interp_kind : str
         The method to use for interpolation
 
@@ -36,7 +40,7 @@ def post_process_variables(variables, t_sol, u_sol, mesh=None, interp_kind="line
     for var, eqn in variables.items():
         pybamm.logger.debug("Post-processing {}".format(var))
         processed_variables[var] = ProcessedVariable(
-            eqn, t_sol, u_sol, mesh, interp_kind, known_evals
+            eqn, t_sol, u_sol, mesh, inputs, interp_kind, known_evals
         )
 
         for t in known_evals:
@@ -64,6 +68,8 @@ class ProcessedVariable(object):
     mesh : :class:`pybamm.Mesh`
         The mesh used to solve, used here to calculate the reference x values for
         interpolation
+    inputs : dict, optional
+            Any input parameters to pass to the model
     interp_kind : str
         The method to use for interpolation
     """
@@ -74,6 +80,7 @@ def __init__(
         t_sol,
         u_sol,
         mesh=None,
+        inputs=None,
         interp_kind="linear",
         known_evals=None,
     ):
@@ -81,6 +88,7 @@ def __init__(
         self.t_sol = t_sol
         self.u_sol = u_sol
         self.mesh = mesh
+        self.inputs = inputs or {}
         self.interp_kind = interp_kind
         self.domain = base_variable.domain
         self.auxiliary_domains = base_variable.auxiliary_domains
@@ -88,7 +96,10 @@ def __init__(
 
         if self.known_evals:
             self.base_eval, self.known_evals[t_sol[0]] = base_variable.evaluate(
-                t_sol[0], u_sol[:, 0], known_evals=self.known_evals[t_sol[0]]
+                t_sol[0],
+                u_sol[:, 0],
+                self.inputs,
+                known_evals=self.known_evals[t_sol[0]],
             )
         else:
             self.base_eval = base_variable.evaluate(t_sol[0], u_sol[:, 0])
@@ -131,10 +142,12 @@ def initialise_1D(self):
             t = self.t_sol[idx]
             if self.known_evals:
                 entries[idx], self.known_evals[t] = self.base_variable.evaluate(
-                    t, self.u_sol[:, idx], known_evals=self.known_evals[t]
+                    t, self.u_sol[:, idx], self.inputs, known_evals=self.known_evals[t]
                 )
             else:
-                entries[idx] = self.base_variable.evaluate(t, self.u_sol[:, idx])
+                entries[idx] = self.base_variable.evaluate(
+                    t, self.u_sol[:, idx], self.inputs
+                )
 
         # No discretisation provided, or variable has no domain (function of t only)
         self._interpolation_function = interp.interp1d(
@@ -158,12 +171,12 @@ def initialise_2D(self):
             u = self.u_sol[:, idx]
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, known_evals=self.known_evals[t]
+                    t, u, self.inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, idx] = eval_and_known_evals[0][:, 0]
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
-                entries[:, idx] = self.base_variable.evaluate(t, u)[:, 0]
+                entries[:, idx] = self.base_variable.evaluate(t, u, self.inputs)[:, 0]
 
         # Process the discretisation to get x values
         nodes = self.mesh.combine_submeshes(*self.domain)[0].nodes
@@ -301,7 +314,7 @@ def initialise_3D(self):
             u = self.u_sol[:, idx]
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, known_evals=self.known_evals[t]
+                    t, u, self.inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(
                     eval_and_known_evals[0],
@@ -311,7 +324,7 @@ def initialise_3D(self):
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
                 entries[:, :, idx] = np.reshape(
-                    self.base_variable.evaluate(t, u),
+                    self.base_variable.evaluate(t, u, self.inputs),
                     [first_dim_size, second_dim_size],
                     order=order,
                 )
@@ -366,13 +379,13 @@ def initialise_3D_scikit_fem(self):
             u = self.u_sol[:, idx]
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, known_evals=self.known_evals[t]
+                    t, u, self.inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(eval_and_known_evals[0], [len_y, len_z])
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
                 entries[:, :, idx] = np.reshape(
-                    self.base_variable.evaluate(t, u), [len_y, len_z]
+                    self.base_variable.evaluate(t, u, self.inputs), [len_y, len_z]
                 )
 
         # assign attributes for reference

From 1c53cfc05102bfbd09161181958427d91cb41a2d Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Tue, 3 Dec 2019 19:48:55 -0500
Subject: [PATCH 12/18] #743 fix tests [ci skip]

---
 .../Tutorial 3 - Basic Plotting.ipynb         |   2 +-
 .../Tutorial 4 - Model Options.ipynb          |   6 +-
 ...orial 5 - Saving, loading and inputs.ipynb | 152 ++++++++++++++++++
 pybamm/solvers/dae_solver.py                  |  12 +-
 pybamm/solvers/ode_solver.py                  |   7 +-
 5 files changed, 170 insertions(+), 9 deletions(-)
 create mode 100644 examples/notebooks/Getting Started/Tutorial 5 - Saving, loading and inputs.ipynb

diff --git a/examples/notebooks/Getting Started/Tutorial 3 - Basic Plotting.ipynb b/examples/notebooks/Getting Started/Tutorial 3 - Basic Plotting.ipynb
index 5263289889..0f16219533 100644
--- a/examples/notebooks/Getting Started/Tutorial 3 - Basic Plotting.ipynb	
+++ b/examples/notebooks/Getting Started/Tutorial 3 - Basic Plotting.ipynb	
@@ -560,7 +560,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.7.4"
   }
  },
  "nbformat": 4,
diff --git a/examples/notebooks/Getting Started/Tutorial 4 - Model Options.ipynb b/examples/notebooks/Getting Started/Tutorial 4 - Model Options.ipynb
index 3b21989fcd..6a044c6008 100644
--- a/examples/notebooks/Getting Started/Tutorial 4 - Model Options.ipynb	
+++ b/examples/notebooks/Getting Started/Tutorial 4 - Model Options.ipynb	
@@ -98,7 +98,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In this tutorial we have seen how to adjust the model options. To see all of the options currently available in PyBaMM, please take a look at the documentation [here](https://pybamm.readthedocs.io/en/latest/source/models/base_models/base_battery_model.html)."
+    "In this tutorial we have seen how to adjust the model options. To see all of the options currently available in PyBaMM, please take a look at the documentation [here](https://pybamm.readthedocs.io/en/latest/source/models/base_models/base_battery_model.html).\n",
+    "\n",
+    "In [Tutorial 5](./Tutorial%205%20-%20Saving%2C%20loading%20and%20inputs.ipynb) we show how to save and load a simulation, and use \"inputs\" to quickly study the effect of various parameters"
    ]
   },
   {
@@ -125,7 +127,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.7.4"
   }
  },
  "nbformat": 4,
diff --git a/examples/notebooks/Getting Started/Tutorial 5 - Saving, loading and inputs.ipynb b/examples/notebooks/Getting Started/Tutorial 5 - Saving, loading and inputs.ipynb
new file mode 100644
index 0000000000..e659ac948b
--- /dev/null
+++ b/examples/notebooks/Getting Started/Tutorial 5 - Saving, loading and inputs.ipynb	
@@ -0,0 +1,152 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Tutorial 5 - Saving and loading simulations, and using inputs"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In [Tutorial 4](Tutorial%204%20-%20Model%20Options.ipynb) we saw how to change model options. We now introduce another useful aspect of PyBaMM, which is the ability to save simulations and use inputs to quickly study the effect of specific parameters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pybamm"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this tutorial, we work with the standard SPMe with default options"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model = pybamm.lithium_ion.SPMe()\n",
+    "sim = pybamm.Simulation(model)\n",
+    "\n",
+    "# Make the negative electrode diffusivity an input\n",
+    "parameter_values = sim.parameter_values\n",
+    "parameter_values.update({\"Negative electrode diffusivity [m2.s-1]\": \"[input]\"})\n",
+    "sim.specs(parameter_values=parameter_values)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Instead of solving the simulation, we can now save it first"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim.save(\"spme.pickle\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can now successively load and solve the simulation, and plot the voltage with various diffusivities"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We now plot the cell temperature and the total heating by passing these variables to the `plot` method as we saw in Tutorial 3:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "17368df63d2b44d3982c151286110d41",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=0.2014814814814815, step=0.05), Output()), _…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "sim.plot([\"Cell temperature [K]\", \"Total heating [W.m-3]\", \"Terminal voltage [V]\"])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this tutorial we have seen how to adjust the model options. To see all of the options currently available in PyBaMM, please take a look at the documentation [here](https://pybamm.readthedocs.io/en/latest/source/models/base_models/base_battery_model.html)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A key point to be aware of is that parameters that affect mesh geometry (i.e. electrode and separator widths) cannot be used as inputs, since no remeshing is performed when inputs are changed."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 3441d26284..0b639906ac 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -104,7 +104,7 @@ def compute_solution(self, model, t_eval, inputs=None):
 
         return solution, solve_time, termination
 
-    def set_up(self, model, inputs):
+    def set_up(self, model, inputs=None):
         """Unpack model, perform checks, simplify and calculate jacobian.
 
         Parameters
@@ -112,10 +112,12 @@ def set_up(self, model, inputs):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
-        inputs : dict
+        inputs : dict, optional
             Any input parameters to pass to the model when solving
 
         """
+        inputs = inputs or {}
+
         # create simplified rhs, algebraic and event expressions
         concatenated_rhs = model.concatenated_rhs
         concatenated_algebraic = model.concatenated_algebraic
@@ -210,7 +212,7 @@ def get_event_class(event):
         self.event_funs = [get_event_class(event) for event in events.values()]
         self.jacobian = jacobian
 
-    def set_up_casadi(self, model, inputs):
+    def set_up_casadi(self, model, inputs=None):
         """Convert model to casadi format and use their inbuilt functionalities.
 
         Parameters
@@ -218,10 +220,12 @@ def set_up_casadi(self, model, inputs):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
-        inputs : dict
+        inputs : dict, optional
             Any input parameters to pass to the model when solving
 
         """
+        inputs = inputs or {}
+
         # Convert model attributes to casadi
         t_casadi = casadi.MX.sym("t")
         y0 = model.concatenated_initial_conditions
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 55c06bb9c3..8faf202cdc 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -61,7 +61,7 @@ def compute_solution(self, model, t_eval, inputs=None):
 
         return solution, solve_time, termination
 
-    def set_up(self, model, inputs):
+    def set_up(self, model, inputs=None):
         """Unpack model, perform checks, simplify and calculate jacobian.
 
         Parameters
@@ -86,6 +86,8 @@ def set_up(self, model, inputs):
                 """Cannot use ODE solver to solve model with DAEs"""
             )
 
+        inputs = inputs or {}
+
         # create simplified rhs and event expressions
         concatenated_rhs = model.concatenated_rhs
         events = model.events
@@ -149,7 +151,7 @@ def get_event_class(event):
         self.event_funs = [get_event_class(event) for event in events.values()]
         self.jacobian = jacobian
 
-    def set_up_casadi(self, model, inputs):
+    def set_up_casadi(self, model, inputs=None):
         """Convert model to casadi format and use their inbuilt functionalities.
 
         Parameters
@@ -177,6 +179,7 @@ def set_up_casadi(self, model, inputs):
 
         t_casadi = casadi.MX.sym("t")
         y_casadi = casadi.MX.sym("y", len(y0))
+        inputs = inputs or {}
         u_casadi = {name: casadi.MX.sym(name) for name in inputs.keys()}
 
         if self.y_pad is not None:

From 68d960bc11d93ed5172105988a001a0786ad2410 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <vsulzer@umich.edu>
Date: Tue, 3 Dec 2019 20:17:47 -0500
Subject: [PATCH 13/18] #743 remove notebooks

---
 .../Tutorial 3 - Basic Plotting.ipynb         |   2 +-
 .../Tutorial 4 - Model Options.ipynb          |   6 +-
 ...orial 5 - Saving, loading and inputs.ipynb | 152 ------------------
 pybamm/simulation.py                          |   1 +
 4 files changed, 4 insertions(+), 157 deletions(-)
 delete mode 100644 examples/notebooks/Getting Started/Tutorial 5 - Saving, loading and inputs.ipynb

diff --git a/examples/notebooks/Getting Started/Tutorial 3 - Basic Plotting.ipynb b/examples/notebooks/Getting Started/Tutorial 3 - Basic Plotting.ipynb
index 0f16219533..5263289889 100644
--- a/examples/notebooks/Getting Started/Tutorial 3 - Basic Plotting.ipynb	
+++ b/examples/notebooks/Getting Started/Tutorial 3 - Basic Plotting.ipynb	
@@ -560,7 +560,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,
diff --git a/examples/notebooks/Getting Started/Tutorial 4 - Model Options.ipynb b/examples/notebooks/Getting Started/Tutorial 4 - Model Options.ipynb
index 6a044c6008..3b21989fcd 100644
--- a/examples/notebooks/Getting Started/Tutorial 4 - Model Options.ipynb	
+++ b/examples/notebooks/Getting Started/Tutorial 4 - Model Options.ipynb	
@@ -98,9 +98,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In this tutorial we have seen how to adjust the model options. To see all of the options currently available in PyBaMM, please take a look at the documentation [here](https://pybamm.readthedocs.io/en/latest/source/models/base_models/base_battery_model.html).\n",
-    "\n",
-    "In [Tutorial 5](./Tutorial%205%20-%20Saving%2C%20loading%20and%20inputs.ipynb) we show how to save and load a simulation, and use \"inputs\" to quickly study the effect of various parameters"
+    "In this tutorial we have seen how to adjust the model options. To see all of the options currently available in PyBaMM, please take a look at the documentation [here](https://pybamm.readthedocs.io/en/latest/source/models/base_models/base_battery_model.html)."
    ]
   },
   {
@@ -127,7 +125,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.6.8"
   }
  },
  "nbformat": 4,
diff --git a/examples/notebooks/Getting Started/Tutorial 5 - Saving, loading and inputs.ipynb b/examples/notebooks/Getting Started/Tutorial 5 - Saving, loading and inputs.ipynb
deleted file mode 100644
index e659ac948b..0000000000
--- a/examples/notebooks/Getting Started/Tutorial 5 - Saving, loading and inputs.ipynb	
+++ /dev/null
@@ -1,152 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Tutorial 5 - Saving and loading simulations, and using inputs"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "In [Tutorial 4](Tutorial%204%20-%20Model%20Options.ipynb) we saw how to change model options. We now introduce another useful aspect of PyBaMM, which is the ability to save simulations and use inputs to quickly study the effect of specific parameters."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 2,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import pybamm"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "In this tutorial, we work with the standard SPMe with default options"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 3,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "model = pybamm.lithium_ion.SPMe()\n",
-    "sim = pybamm.Simulation(model)\n",
-    "\n",
-    "# Make the negative electrode diffusivity an input\n",
-    "parameter_values = sim.parameter_values\n",
-    "parameter_values.update({\"Negative electrode diffusivity [m2.s-1]\": \"[input]\"})\n",
-    "sim.specs(parameter_values=parameter_values)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Instead of solving the simulation, we can now save it first"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sim.save(\"spme.pickle\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We can now successively load and solve the simulation, and plot the voltage with various diffusivities"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We now plot the cell temperature and the total heating by passing these variables to the `plot` method as we saw in Tutorial 3:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "17368df63d2b44d3982c151286110d41",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "interactive(children=(FloatSlider(value=0.0, description='t', max=0.2014814814814815, step=0.05), Output()), _…"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
-   "source": [
-    "sim.plot([\"Cell temperature [K]\", \"Total heating [W.m-3]\", \"Terminal voltage [V]\"])"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "In this tutorial we have seen how to adjust the model options. To see all of the options currently available in PyBaMM, please take a look at the documentation [here](https://pybamm.readthedocs.io/en/latest/source/models/base_models/base_battery_model.html)."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "A key point to be aware of is that parameters that affect mesh geometry (i.e. electrode and separator widths) cannot be used as inputs, since no remeshing is performed when inputs are changed."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.7.4"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 99b87973fd..31dc82e4b5 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -187,6 +187,7 @@ def solve(self, t_eval=None, solver=None, inputs=None, check_model=True):
         if solver is None:
             solver = self.solver
 
+        self.t_eval = t_eval
         self._solution = solver.solve(self.built_model, t_eval, inputs=inputs)
 
     def step(self, dt, solver=None, external_variables=None, save=True):

From 22c8ae826d57a281c07a5aa7d5517b79d4214f20 Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Wed, 4 Dec 2019 00:19:51 -0500
Subject: [PATCH 14/18] #743 bug fixes and coverage

---
 pybamm/solvers/dae_solver.py                    | 8 +++++---
 tests/unit/test_solvers/test_dae_solver.py      | 5 +++++
 tests/unit/test_solvers/test_scikits_solvers.py | 8 --------
 3 files changed, 10 insertions(+), 11 deletions(-)

diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 0b639906ac..d235ae0d58 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -159,9 +159,7 @@ def set_up(self, model, inputs=None):
                 jac = pybamm.EvaluatorPython(jac)
 
             jacobian = Jacobian(jac.evaluate)
-            jacobian_alg = Jacobian(jac_algebraic.evaluate)
-            jacobian_alg.set_pad_ext(self.y_pad, self.y_ext)
-            jacobian_alg.set_inputs(inputs)
+            jacobian_alg = JacobianAlgebraic(jac_algebraic.evaluate)
 
         else:
             jacobian = None
@@ -185,6 +183,8 @@ def set_up(self, model, inputs=None):
         rhs.set_inputs(inputs)
         algebraic.set_pad_ext(self.y_pad, self.y_ext)
         algebraic.set_inputs(inputs)
+        jacobian_alg.set_pad_ext(self.y_pad, self.y_ext)
+        jacobian_alg.set_inputs(inputs)
 
         if len(model.algebraic) > 0:
             y0 = self.calculate_consistent_initial_conditions(
@@ -302,6 +302,8 @@ def set_up_casadi(self, model, inputs=None):
         rhs.set_inputs(inputs)
         algebraic.set_pad_ext(self.y_pad, self.y_ext)
         algebraic.set_inputs(inputs)
+        jacobian_alg.set_pad_ext(self.y_pad, self.y_ext)
+        jacobian_alg.set_inputs(inputs)
 
         if len(model.algebraic) > 0:
 
diff --git a/tests/unit/test_solvers/test_dae_solver.py b/tests/unit/test_solvers/test_dae_solver.py
index 6e35964415..4eb334e722 100644
--- a/tests/unit/test_solvers/test_dae_solver.py
+++ b/tests/unit/test_solvers/test_dae_solver.py
@@ -77,6 +77,11 @@ def algebraic(t, y):
         ):
             solver.calculate_consistent_initial_conditions(rhs, algebraic, y0)
 
+    def test_errors(self):
+        solver = pybamm.DaeSolver()
+        with self.assertRaises(NotImplementedError):
+            solver.integrate(None, None, None)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index aa55b3edd9..f7c15fcbf7 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -568,14 +568,6 @@ def test_model_solver_dae_with_jacobian_python(self):
         )
         N = combined_submesh[0].npts
 
-        def jacobian(t, y):
-            return np.block(
-                [
-                    [0.1 * np.eye(N), np.zeros((N, N))],
-                    [2.0 * np.eye(N), -1.0 * np.eye(N)],
-                ]
-            )
-
         # Solve
         solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
         t_eval = np.linspace(0, 1, 100)

From 15e7dc3b2ba61e25bdc70be776607ee7982a5cc1 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <vsulzer@umich.edu>
Date: Wed, 4 Dec 2019 10:38:37 -0500
Subject: [PATCH 15/18] #743 more bugs and coverage

---
 pybamm/solvers/dae_solver.py                            | 8 ++++----
 tests/unit/test_expression_tree/test_input_parameter.py | 3 +++
 tests/unit/test_solvers/test_scikits_solvers.py         | 9 +++++++++
 3 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index d235ae0d58..1822074368 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -160,6 +160,8 @@ def set_up(self, model, inputs=None):
 
             jacobian = Jacobian(jac.evaluate)
             jacobian_alg = JacobianAlgebraic(jac_algebraic.evaluate)
+            jacobian_alg.set_pad_ext(self.y_pad, self.y_ext)
+            jacobian_alg.set_inputs(inputs)
 
         else:
             jacobian = None
@@ -183,8 +185,6 @@ def set_up(self, model, inputs=None):
         rhs.set_inputs(inputs)
         algebraic.set_pad_ext(self.y_pad, self.y_ext)
         algebraic.set_inputs(inputs)
-        jacobian_alg.set_pad_ext(self.y_pad, self.y_ext)
-        jacobian_alg.set_inputs(inputs)
 
         if len(model.algebraic) > 0:
             y0 = self.calculate_consistent_initial_conditions(
@@ -290,6 +290,8 @@ def set_up_casadi(self, model, inputs=None):
 
             jacobian = JacobianCasadi(casadi_jac_fn)
             jacobian_alg = JacobianAlgebraicCasadi(casadi_jac_alg_fn)
+            jacobian_alg.set_pad_ext(self.y_pad, self.y_ext)
+            jacobian_alg.set_inputs(inputs)
 
         else:
             jacobian = None
@@ -302,8 +304,6 @@ def set_up_casadi(self, model, inputs=None):
         rhs.set_inputs(inputs)
         algebraic.set_pad_ext(self.y_pad, self.y_ext)
         algebraic.set_inputs(inputs)
-        jacobian_alg.set_pad_ext(self.y_pad, self.y_ext)
-        jacobian_alg.set_inputs(inputs)
 
         if len(model.algebraic) > 0:
 
diff --git a/tests/unit/test_expression_tree/test_input_parameter.py b/tests/unit/test_expression_tree/test_input_parameter.py
index b64d8d8d30..f5c24e3909 100644
--- a/tests/unit/test_expression_tree/test_input_parameter.py
+++ b/tests/unit/test_expression_tree/test_input_parameter.py
@@ -23,6 +23,9 @@ def test_errors(self):
             a.evaluate(u="not a dictionary")
         with self.assertRaises(KeyError):
             a.evaluate(u={"bad param": 5})
+        # if u is not provided it gets turned into a dictionary and then raises KeyError
+        with self.assertRaises(KeyError):
+            a.evaluate()
 
 
 if __name__ == "__main__":
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index f7c15fcbf7..1108b8b7e6 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -568,6 +568,15 @@ def test_model_solver_dae_with_jacobian_python(self):
         )
         N = combined_submesh[0].npts
 
+        def jacobian(t, y):
+            return np.block(
+                [
+                    [0.1 * np.eye(N), np.zeros((N, N))],
+                    [2.0 * np.eye(N), -1.0 * np.eye(N)],
+                ]
+            )
+
+        model.jacobian = jacobian
         # Solve
         solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
         t_eval = np.linspace(0, 1, 100)

From 136244ca5b4c98b855a2067d4ab5a1e52e785886 Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Thu, 5 Dec 2019 09:48:20 -0500
Subject: [PATCH 16/18] #743 codacy

---
 .../expression_tree/operations/convert_to_casadi.py  |  2 +-
 pybamm/solvers/base_solver.py                        | 12 +++++++++---
 pybamm/solvers/casadi_solver.py                      |  4 ++--
 tests/unit/test_solvers/test_scikits_solvers.py      |  4 ++--
 4 files changed, 14 insertions(+), 8 deletions(-)

diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py
index 12a91a31a1..983008d706 100644
--- a/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -42,7 +42,7 @@ def convert(self, symbol, t=None, y=None, u=None):
 
             return casadi_symbol
 
-    def _convert(self, symbol, t, y, u):
+    def _convert(self, symbol, t=None, y=None, u=None):
         """ See :meth:`CasadiConverter.convert()`. """
         if isinstance(
             symbol, (pybamm.Scalar, pybamm.Array, pybamm.Time, pybamm.InputParameter)
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index c13c8ddc4d..7e93be8b58 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -230,7 +230,7 @@ def set_external_variables(self, model, external_variables):
                 )
             self.y_ext[y_slice] = var_vals
 
-    def compute_solution(self, model, t_eval):
+    def compute_solution(self, model, t_eval, inputs=None):
         """Calculate the solution of the model at specified times. Note: this
         does *not* execute the solver setup.
 
@@ -241,11 +241,13 @@ def compute_solution(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         """
         raise NotImplementedError
 
-    def set_up(self, model):
+    def set_up(self, model, inputs=None):
         """Unpack model, perform checks, simplify and calculate jacobian.
 
         Parameters
@@ -253,11 +255,13 @@ def set_up(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         """
         raise NotImplementedError
 
-    def set_up_casadi(self, model):
+    def set_up_casadi(self, model, inputs=None):
         """Convert model to casadi format and use their inbuilt functionalities.
 
         Parameters
@@ -265,6 +269,8 @@ def set_up_casadi(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         """
         raise NotImplementedError
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index b008ca6a71..7c0329d622 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -178,7 +178,7 @@ def solve(self, model, t_eval, inputs=None):
             )
             return solution
 
-    def compute_solution(self, model, t_eval, inputs):
+    def compute_solution(self, model, t_eval, inputs=None):
         """Calculate the solution of the model at specified times. In this class, we
         overwrite the behaviour of :class:`pybamm.DaeSolver`, since CasADi requires
         slightly different syntax.
@@ -190,7 +190,7 @@ def compute_solution(self, model, t_eval, inputs):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
-        inputs : dict
+        inputs : dict, optional
             Any input parameters to pass to the model when solving
 
         """
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 1108b8b7e6..c2259b6a79 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -720,9 +720,9 @@ def test_model_solver_dae_events_casadi(self):
 
     def test_model_solver_dae_inputs_events(self):
         # Create model
-        for format in ["python", "casadi"]:
+        for form in ["python", "casadi"]:
             model = pybamm.BaseModel()
-            model.convert_to_format = format
+            model.convert_to_format = form
             whole_cell = ["negative electrode", "separator", "positive electrode"]
             var1 = pybamm.Variable("var1", domain=whole_cell)
             var2 = pybamm.Variable("var2", domain=whole_cell)

From 662f40d352a0e049d1ee1c685c1e9314f5f4e438 Mon Sep 17 00:00:00 2001
From: tinosulzer <valentinsulzer@hotmail.com>
Date: Thu, 5 Dec 2019 09:49:26 -0500
Subject: [PATCH 17/18] #743 changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index ca6f825ebb..b3a3a23a06 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added `InputParameter` node for quickly changing parameter values ([#752](https://github.com/pybamm-team/PyBaMM/pull/752))
 -   Generalized importing of external variables ([#728](https://github.com/pybamm-team/PyBaMM/pull/728))
 -   Separated active and inactive material volume fractions ([#726](https://github.com/pybamm-team/PyBaMM/pull/726))
 -   Added submodels for tortuosity ([#726](https://github.com/pybamm-team/PyBaMM/pull/726))

From f9c8a261bd108910466e540b97f484fff8066808 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <vsulzer@umich.edu>
Date: Thu, 5 Dec 2019 14:04:53 -0500
Subject: [PATCH 18/18] #743 add input to simulation step

---
 pybamm/simulation.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 31dc82e4b5..4ce57ce3bd 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -190,7 +190,7 @@ def solve(self, t_eval=None, solver=None, inputs=None, check_model=True):
         self.t_eval = t_eval
         self._solution = solver.solve(self.built_model, t_eval, inputs=inputs)
 
-    def step(self, dt, solver=None, external_variables=None, save=True):
+    def step(self, dt, solver=None, external_variables=None, inputs=None, save=True):
         """
         A method to step the model forward one timestep. This method will
         automatically build and set the model parameters if not already done so.
@@ -206,6 +206,8 @@ def step(self, dt, solver=None, external_variables=None, save=True):
             values at the current time. The variables must correspond to
             the variables that would normally be found by solving the
             submodels that have been made external.
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
         save : bool
             Turn on to store the solution of all previous timesteps
         """
@@ -215,7 +217,7 @@ def step(self, dt, solver=None, external_variables=None, save=True):
             solver = self.solver
 
         solution = solver.step(
-            self.built_model, dt, external_variables=external_variables
+            self.built_model, dt, external_variables=external_variables, inputs=inputs
         )
 
         if save is False or self._made_first_step is False: