Issue 937 sensitivity

CHANGELOG.md

@@ -2,6 +2,8 @@
 
 ## Features
 
+-   Added sensitivity to `CasadiAlgebraicSolver` ([#940](https://github.com/pybamm-team/PyBaMM/pull/940))
+-   Added `ProcessedCasadiVariable` class, which can handle symbolic variables (i.e. variables for which the inputs are symbolic) ([#940](https://github.com/pybamm-team/PyBaMM/pull/940))
 -   Made `QuickPlot` compatible with Google Colab ([#935](https://github.com/pybamm-team/PyBaMM/pull/935))
 
 # [v0.2.1](https://github.com/pybamm-team/PyBaMM/tree/v0.2.1) - 2020-03-31

docs/index.rst

@@ -33,7 +33,6 @@ Contents
     source/experiments/index
     source/simulation
     source/quick_plot
-    source/processed_variable
     source/util
     source/citations
     source/parameters_cli

docs/source/solvers/index.rst

@@ -10,3 +10,5 @@ Solvers
   casadi_solver
   algebraic_solvers
   solution
+  processed_variable
+

docs/source/processed_variable.rst → docs/source/solvers/processed_variable.rst

@@ -3,3 +3,6 @@ Post-Process Variables
 
 .. autoclass:: pybamm.ProcessedVariable
   :members:
+
+.. autoclass:: pybamm.ProcessedCasadiVariable
+  :members:

docs/source/solvers/solution.rst

@@ -1,5 +1,5 @@
-Solution
-========
+Solutions
+=========
 
 .. autoclass:: pybamm._BaseSolution
   :members:

pybamm/__init__.py

@@ -198,6 +198,8 @@ def version(formatted=False):
 # Solver classes
 #
 from .solvers.solution import Solution, _BaseSolution
+from .solvers.processed_variable import ProcessedVariable
+from .solvers.processed_casadi_variable import ProcessedCasadiVariable
 from .solvers.base_solver import BaseSolver
 from .solvers.dummy_solver import DummySolver
 from .solvers.algebraic_solver import AlgebraicSolver
@@ -217,7 +219,6 @@ def version(formatted=False):
 #
 # other
 #
-from .processed_variable import ProcessedVariable
 from .quick_plot import QuickPlot, dynamic_plot, ax_min, ax_max
 
 from .simulation import Simulation, load_sim, is_notebook

pybamm/solvers/base_solver.py

@@ -161,14 +161,26 @@ def set_up(self, model, inputs=None):
             # set up Jacobian object, for re-use of dict
             jacobian = pybamm.Jacobian()
         else:
+            # Create placeholder inputs for evaluating rhs and algebraic sizes
+            placeholder_inputs = {}
+            for k, v in inputs.items():
+                if isinstance(v, casadi.MX):
+                    placeholder_inputs[k] = 0
+                else:
+                    placeholder_inputs[k] = v
             # Convert model attributes to casadi
             t_casadi = casadi.MX.sym("t")
             y_diff = casadi.MX.sym(
-                "y_diff", len(model.concatenated_rhs.evaluate(0, y0, inputs=inputs))
+                "y_diff",
+                len(model.concatenated_rhs.evaluate(0, y0, inputs=placeholder_inputs)),
             )
             y_alg = casadi.MX.sym(
                 "y_alg",
-                len(model.concatenated_algebraic.evaluate(0, y0, inputs=inputs)),
+                len(
+                    model.concatenated_algebraic.evaluate(
+                        0, y0, inputs=placeholder_inputs
+                    )
+                ),
             )
             y_casadi = casadi.vertcat(y_diff, y_alg)
             p_casadi = {}
@@ -522,15 +534,30 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         if (np.diff(t_eval) < 0).any():
             raise pybamm.SolverError("t_eval must increase monotonically")
 
-        # Non-dimensionalise t_eval
-
-        # Set up
-        timer = pybamm.Timer()
-
         # Set up external variables and inputs
         external_variables = external_variables or {}
         inputs = inputs or {}
         ext_and_inputs = {**external_variables, **inputs}
+        # Only allow symbolic inputs for CasadiAlgebraicSolver
+        if not isinstance(self, pybamm.CasadiAlgebraicSolver) and any(
+            isinstance(v, str) and v.startswith("[sym]") for v in inputs.values()
+        ):
+            raise pybamm.SolverError(
+                "Only CasadiAlgebraicSolver can have symbolic inputs"
+            )
+        # Make symbolic inputs
+        for k, v in ext_and_inputs.items():
+            if isinstance(v, str) and v.startswith("[sym]"):
+                # If only [sym] is specified then symbolic input has size 1
+                if v == "[sym]":
+                    size = 1
+                # Otherwise read symbolic input size n from '[sym]n'
+                else:
+                    size = int(v[5:])
+                ext_and_inputs[k] = casadi.MX.sym(k, size)
+
+        # Set up
+        timer = pybamm.Timer()
 
         # Set up (if not done already)
         if model not in self.models_set_up:

pybamm/solvers/casadi_algebraic_solver.py

@@ -46,13 +46,20 @@ def _integrate(self, model, t_eval, inputs=None):
         t_eval : :class:`numpy.array`, size (k,)
             The times at which to compute the solution
         inputs : dict, optional
-            Any input parameters to pass to the model when solving
+            Any input parameters to pass to the model when solving. If this contains
+            symbolic variables (e.g. via the keyword [sym]), then the solution will
+            consist of `ProcessedCasadiVariable` objects, whose sensitivity can be
+            evaluated.
         """
         y0 = model.y0
 
-        y = np.empty((len(y0), len(t_eval)))
+        y = None
 
         # Set up
+        # Record whether there are any symbolic inputs
+        has_symbolic_inputs = any(isinstance(v, casadi.MX) for v in inputs.values())
+
+        # Create casadi objects for the root-finder
         inputs = casadi.vertcat(*[x for x in inputs.values()])
         t_sym = casadi.MX.sym("t")
         y_sym = casadi.MX.sym("y_alg", y0.shape[0])
@@ -71,17 +78,23 @@ def _integrate(self, model, t_eval, inputs=None):
         for idx, t in enumerate(t_eval):
             # Evaluate algebraic with new t and previous y0, if it's already close
             # enough then keep it
-            if np.all(abs(model.algebraic_eval(t, y0, inputs)) < self.tol):
+            # We can't do this if there are symbolic inputs
+            if has_symbolic_inputs is False and np.all(
+                abs(model.casadi_algebraic(t, y0, inputs).full()) < self.tol
+            ):
                 pybamm.logger.debug(
                     "Keeping same solution at t={}".format(t * model.timescale_eval)
                 )
-                y[:, idx] = y0
-            # Otherwise calculate new y0
+                if y is None:
+                    y = y0
+                else:
+                    y = casadi.horzcat(y, y0)
+            # Otherwise calculate new y_sol
             else:
                 t_inputs = casadi.vertcat(t, inputs)
                 # Solve
                 try:
-                    y_sol = roots(y0, t_inputs).full().flatten()
+                    y_sol = roots(y0, t_inputs)
                     success = True
                     message = None
                     # Check final output
@@ -91,11 +104,18 @@ def _integrate(self, model, t_eval, inputs=None):
                     message = err.args[0]
                     fun = None
 
-                if success and np.all(casadi.fabs(fun) < self.tol):
+                # If there are no symbolic inputs, check the function is below the tol
+                # Skip this check if there are symbolic inputs
+                if success and (
+                    has_symbolic_inputs is True or np.all(casadi.fabs(fun) < self.tol)
+                ):
                     # update initial guess for the next iteration
                     y0 = y_sol
                     # update solution array
-                    y[:, idx] = y_sol
+                    if y is None:
+                        y = y_sol
+                    else:
+                        y = casadi.horzcat(y, y_sol)
                 elif not success:
                     raise pybamm.SolverError(
                         "Could not find acceptable solution: {}".format(message)