diff --git a/.requirements-docs.txt b/.requirements-docs.txt index 1aa9e60102..c7bcd6446f 100644 --- a/.requirements-docs.txt +++ b/.requirements-docs.txt @@ -5,5 +5,6 @@ pandas>=0.23 anytree>=2.4.3 autograd>=1.2 scikit-fem>=0.2.0 +casadi>=3.5.0 guzzle-sphinx-theme sphinx>=1.5 diff --git a/CHANGELOG.md b/CHANGELOG.md index 07421aef02..80b9c20626 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,10 +2,12 @@ ## Features -- Added interface (via pybind11) to sundials with the IDA KLU sparse linear solver ([#657](https://github.com/pybamm-team/PyBaMM/pull/657)) +- Add interface to CasADi solver ([#687](https://github.com/pybamm-team/PyBaMM/pull/687)) +- Add option to use CasADi's Algorithmic Differentiation framework to calculate Jacobians ([#687](https://github.com/pybamm-team/PyBaMM/pull/687)) - Add method to evaluate parameters more easily ([#669](https://github.com/pybamm-team/PyBaMM/pull/669)) - Add `Jacobian` class to reuse known Jacobians of expressions ([#665](https://github.com/pybamm-team/PyBaMM/pull/670)) - Add `Interpolant` class to interpolate experimental data (e.g. OCP curves) ([#661](https://github.com/pybamm-team/PyBaMM/pull/661)) +- Added interface (via pybind11) to sundials with the IDA KLU sparse linear solver ([#657](https://github.com/pybamm-team/PyBaMM/pull/657)) - Allow parameters to be set by material or by specifying a particular paper ([#647](https://github.com/pybamm-team/PyBaMM/pull/647)) - Set relative and absolute tolerances independently in solvers ([#645](https://github.com/pybamm-team/PyBaMM/pull/645)) - Add some non-uniform meshes in 1D and 2D ([#617](https://github.com/pybamm-team/PyBaMM/pull/617)) @@ -17,6 +19,7 @@ ## Bug fixes +- Fix differentiation of functions that have more than one argument ([#687](https://github.com/pybamm-team/PyBaMM/pull/687)) - Add warning if `ProcessedVariable` is called outisde its interpolation range ([#681](https://github.com/pybamm-team/PyBaMM/pull/681)) - Improve the way `ProcessedVariable` objects are created in higher dimensions ([#581](https://github.com/pybamm-team/PyBaMM/pull/581)) diff --git a/docs/source/expression_tree/index.rst b/docs/source/expression_tree/index.rst index c5055334dd..3eedd7b622 100644 --- a/docs/source/expression_tree/index.rst +++ b/docs/source/expression_tree/index.rst @@ -17,6 +17,4 @@ Expression Tree broadcasts functions interpolant - evaluate - simplify - jacobian + operations/index diff --git a/docs/source/expression_tree/operations/convert_to_casadi.rst b/docs/source/expression_tree/operations/convert_to_casadi.rst new file mode 100644 index 0000000000..76a6ed614e --- /dev/null +++ b/docs/source/expression_tree/operations/convert_to_casadi.rst @@ -0,0 +1,5 @@ +Convert to CasADi +================= + +.. autoclass:: pybamm.CasadiConverter + :members: diff --git a/docs/source/expression_tree/evaluate.rst b/docs/source/expression_tree/operations/evaluate.rst similarity index 100% rename from docs/source/expression_tree/evaluate.rst rename to docs/source/expression_tree/operations/evaluate.rst diff --git a/docs/source/expression_tree/operations/index.rst b/docs/source/expression_tree/operations/index.rst new file mode 100644 index 0000000000..31004a2204 --- /dev/null +++ b/docs/source/expression_tree/operations/index.rst @@ -0,0 +1,11 @@ +Operations on expression trees +============================== + +Classes and functions that operate on the expression tree + +.. toctree:: + + simplify + evaluate + jacobian + convert_to_casadi diff --git a/docs/source/expression_tree/jacobian.rst b/docs/source/expression_tree/operations/jacobian.rst similarity index 100% rename from docs/source/expression_tree/jacobian.rst rename to docs/source/expression_tree/operations/jacobian.rst diff --git a/docs/source/expression_tree/simplify.rst b/docs/source/expression_tree/operations/simplify.rst similarity index 100% rename from docs/source/expression_tree/simplify.rst rename to docs/source/expression_tree/operations/simplify.rst diff --git a/docs/source/solvers/casadi_solver.rst b/docs/source/solvers/casadi_solver.rst new file mode 100644 index 0000000000..f21e3c74cf --- /dev/null +++ b/docs/source/solvers/casadi_solver.rst @@ -0,0 +1,5 @@ +Casadi Solver +============= + +.. autoclass:: pybamm.CasadiSolver + :members: diff --git a/docs/source/solvers/index.rst b/docs/source/solvers/index.rst index 6cc6c0021d..d85865dec7 100644 --- a/docs/source/solvers/index.rst +++ b/docs/source/solvers/index.rst @@ -7,4 +7,5 @@ Solvers base_solvers scipy_solver scikits_solvers + casadi_solver solution diff --git a/examples/scripts/compare-dae-solver.py b/examples/scripts/compare-dae-solver.py index 05476790cf..e4dff68a3d 100644 --- a/examples/scripts/compare-dae-solver.py +++ b/examples/scripts/compare-dae-solver.py @@ -16,9 +16,7 @@ # set mesh var = pybamm.standard_spatial_vars - -var_pts = {var.x_n: 60, var.x_s: 100, var.x_p: 60, var.r_n: 50, var.r_p: 50} -# var_pts = model.default_var_pts +var_pts = {var.x_n: 50, var.x_s: 50, var.x_p: 50, var.r_n: 20, var.r_p: 20} mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) # discretise model @@ -26,13 +24,14 @@ disc.process_model(model) # solve model -t_eval = np.linspace(0, 0.2, 100) +t_eval = np.linspace(0, 0.17, 100) -klu_sol = pybamm.IDAKLU(atol=1e-8, rtol=1e-8).solve(model, t_eval) +casadi_sol = pybamm.CasadiSolver(atol=1e-8, rtol=1e-8).solve(model, t_eval) +klu_sol = pybamm.IDAKLUSolver(atol=1e-8, rtol=1e-8).solve(model, t_eval) scikits_sol = pybamm.ScikitsDaeSolver(atol=1e-8, rtol=1e-8).solve(model, t_eval) # plot -models = [model, model] -solutions = [scikits_sol, klu_sol] +models = [model, model, model] +solutions = [scikits_sol, klu_sol, casadi_sol] plot = pybamm.QuickPlot(models, mesh, solutions) plot.dynamic_plot() diff --git a/examples/scripts/compare_lead_acid.py b/examples/scripts/compare_lead_acid.py index dc02b534dc..028a2a1629 100644 --- a/examples/scripts/compare_lead_acid.py +++ b/examples/scripts/compare_lead_acid.py @@ -41,16 +41,18 @@ # solve model solutions = [None] * len(models) -t_eval = np.linspace(0, 3, 1000) +t_eval = np.linspace(0, 1, 1000) for i, model in enumerate(models): solution = model.default_solver.solve(model, t_eval) solutions[i] = solution # plot output_variables = [ - "Electrolyte pressure", - "Electrolyte concentration", - "Volume-averaged velocity [m.s-1]", + "Interfacial current density [A.m-2]", + "Electrolyte concentration [mol.m-3]", + "Current [A]", + "Porosity", + "Electrolyte potential [V]", "Terminal voltage [V]", ] plot = pybamm.QuickPlot(models, mesh, solutions, output_variables) diff --git a/examples/scripts/compare_lead_acid_3D.py b/examples/scripts/compare_lead_acid_3D.py index b7b8cfb66a..584d027351 100644 --- a/examples/scripts/compare_lead_acid_3D.py +++ b/examples/scripts/compare_lead_acid_3D.py @@ -21,8 +21,7 @@ # {"current collector": "potential pair", "dimensionality": 2}, name="2+1D LOQS" # ), pybamm.lead_acid.Full( - {"current collector": "potential pair", "dimensionality": 1}, - name="1+1D Full", + {"current collector": "potential pair", "dimensionality": 1}, name="1+1D Full" ), # pybamm.lead_acid.Full( # {"dimensionality": 1}, name="1+1D uniform Full" diff --git a/examples/scripts/compare_lithium_ion_3D.py b/examples/scripts/compare_lithium_ion_3D.py index e62fc7697b..f271b63cf1 100644 --- a/examples/scripts/compare_lithium_ion_3D.py +++ b/examples/scripts/compare_lithium_ion_3D.py @@ -18,12 +18,10 @@ # load models models = [ pybamm.lithium_ion.SPM( - {"current collector": "potential pair", "dimensionality": 2}, - name="2+1D SPM", + {"current collector": "potential pair", "dimensionality": 2}, name="2+1D SPM" ), pybamm.lithium_ion.SPMe( - {"current collector": "potential pair", "dimensionality": 2}, - name="2+1D SPMe", + {"current collector": "potential pair", "dimensionality": 2}, name="2+1D SPMe" ), ] diff --git a/pybamm/__init__.py b/pybamm/__init__.py index 482d6c9476..3a94c5dc99 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -150,19 +150,22 @@ def version(formatted=False): UndefinedOperationError, GeometryError, ) -from .expression_tree.simplify import ( + +# Operations +from .expression_tree.operations.simplify import ( Simplification, simplify_if_constant, simplify_addition_subtraction, simplify_multiplication_division, ) -from .expression_tree.jacobian import Jacobian -from .expression_tree.evaluate import ( +from .expression_tree.operations.evaluate import ( find_symbols, id_to_python_variable, to_python, EvaluatorPython, ) +from .expression_tree.operations.jacobian import Jacobian +from .expression_tree.operations.convert_to_casadi import CasadiConverter # # Model classes @@ -258,12 +261,12 @@ def version(formatted=False): from .solvers.base_solver import BaseSolver from .solvers.ode_solver import OdeSolver from .solvers.dae_solver import DaeSolver -from .solvers.scipy_solver import ScipySolver -from .solvers.scikits_dae_solver import ScikitsDaeSolver -from .solvers.scikits_ode_solver import ScikitsOdeSolver -from .solvers.scikits_ode_solver import have_scikits_odes from .solvers.algebraic_solver import AlgebraicSolver -from .solvers.idaklu_solver import IDAKLU, have_idaklu +from .solvers.casadi_solver import CasadiSolver +from .solvers.scikits_dae_solver import ScikitsDaeSolver +from .solvers.scikits_ode_solver import ScikitsOdeSolver, have_scikits_odes +from .solvers.scipy_solver import ScipySolver +from .solvers.idaklu_solver import IDAKLUSolver, have_idaklu # diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py index a78206b3f0..284e375cb7 100644 --- a/pybamm/discretisations/discretisation.py +++ b/pybamm/discretisations/discretisation.py @@ -128,7 +128,7 @@ def process_model(self, model, inplace=True): model_disc.options = model.options model_disc.use_jacobian = model.use_jacobian model_disc.use_simplify = model.use_simplify - model_disc.use_to_python = model.use_to_python + model_disc.convert_to_format = model.convert_to_format model_disc.bcs = self.bcs diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py index 41c2492edd..1f1d28e524 100644 --- a/pybamm/expression_tree/functions.py +++ b/pybamm/expression_tree/functions.py @@ -21,10 +21,19 @@ class Function(pybamm.Symbol): derivative : str, optional Which derivative to use when differentiating ("autograd" or "derivative"). Default is "autograd". + differentiated_function : method, optional + The function which was differentiated to obtain this one. Default is None. **Extends:** :class:`pybamm.Symbol` """ - def __init__(self, function, *children, name=None, derivative="autograd"): + def __init__( + self, + function, + *children, + name=None, + derivative="autograd", + differentiated_function=None + ): if name is not None: self.name = name @@ -39,6 +48,7 @@ def __init__(self, function, *children, name=None, derivative="autograd"): self.function = function self.derivative = derivative + self.differentiated_function = differentiated_function # hack to work out whether function takes any params # (signature doesn't work for numpy) @@ -85,7 +95,9 @@ def diff(self, variable): # if variable appears in the function,use autograd to differentiate # function, and apply chain rule if variable.id in [symbol.id for symbol in child.pre_order()]: - partial_derivatives[i] = self._diff(children) * child.diff(variable) + partial_derivatives[i] = self._function_diff( + children, i + ) * child.diff(variable) # remove None entries partial_derivatives = list(filter(None, partial_derivatives)) @@ -96,15 +108,34 @@ def diff(self, variable): return derivative - def _diff(self, children): - """ See :meth:`pybamm.Symbol._diff()`. """ + def _function_diff(self, children, idx): + """ + Derivative with respect to child number 'idx'. + See :meth:`pybamm.Symbol._diff()`. + """ + # Store differentiated function, needed in case we want to convert to CasADi if self.derivative == "autograd": - return Function(autograd.elementwise_grad(self.function), *children) - elif self.derivative == "derivative": - # keep using "derivative" as derivative - return pybamm.Function( - self.function.derivative(), *children, derivative="derivative" + return Function( + autograd.elementwise_grad(self.function, idx), + *children, + differentiated_function=self.function ) + elif self.derivative == "derivative": + if len(children) > 1: + raise ValueError( + """ + differentiation using '.derivative()' not implemented for functions + with more than one child + """ + ) + else: + # keep using "derivative" as derivative + return pybamm.Function( + self.function.derivative(), + *children, + derivative="derivative", + differentiated_function=self.function + ) def _function_jac(self, children_jacs): """ Calculate the jacobian of a function. """ @@ -118,7 +149,7 @@ def _function_jac(self, children_jacs): children = self.orphans for i, child in enumerate(children): if not child.evaluates_to_number(): - jac_fun = self._diff(children) * children_jacs[i] + jac_fun = self._function_diff(children, i) * children_jacs[i] jac_fun.domain = [] if jacobian is None: jacobian = jac_fun @@ -175,7 +206,11 @@ def _function_new_copy(self, children): A new copy of the function """ return pybamm.Function( - self.function, *children, name=self.name, derivative=self.derivative + self.function, + *children, + name=self.name, + derivative=self.derivative, + differentiated_function=self.differentiated_function ) def _function_simplify(self, simplified_children): @@ -203,7 +238,8 @@ def _function_simplify(self, simplified_children): self.function, *simplified_children, name=self.name, - derivative=self.derivative + derivative=self.derivative, + differentiated_function=self.differentiated_function ) @@ -239,8 +275,8 @@ class Cos(SpecificFunction): def __init__(self, child): super().__init__(np.cos, child) - def _diff(self, children): - """ See :meth:`pybamm.Symbol._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Symbol._function_diff()`. """ return -Sin(children[0]) @@ -255,8 +291,8 @@ class Cosh(SpecificFunction): def __init__(self, child): super().__init__(np.cosh, child) - def _diff(self, children): - """ See :meth:`pybamm.Function._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Function._function_diff()`. """ return Sinh(children[0]) @@ -271,8 +307,8 @@ class Exponential(SpecificFunction): def __init__(self, child): super().__init__(np.exp, child) - def _diff(self, children): - """ See :meth:`pybamm.Function._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Function._function_diff()`. """ return Exponential(children[0]) @@ -287,8 +323,8 @@ class Log(SpecificFunction): def __init__(self, child): super().__init__(np.log, child) - def _diff(self, children): - """ See :meth:`pybamm.Function._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Function._function_diff()`. """ return 1 / children[0] @@ -313,8 +349,8 @@ class Sin(SpecificFunction): def __init__(self, child): super().__init__(np.sin, child) - def _diff(self, children): - """ See :meth:`pybamm.Function._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Function._function_diff()`. """ return Cos(children[0]) @@ -329,8 +365,8 @@ class Sinh(SpecificFunction): def __init__(self, child): super().__init__(np.sinh, child) - def _diff(self, children): - """ See :meth:`pybamm.Function._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Function._function_diff()`. """ return Cosh(children[0]) diff --git a/pybamm/expression_tree/operations/__init__.py b/pybamm/expression_tree/operations/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py new file mode 100644 index 0000000000..d613badead --- /dev/null +++ b/pybamm/expression_tree/operations/convert_to_casadi.py @@ -0,0 +1,129 @@ +# +# Convert a PyBaMM expression tree to a CasADi expression tree +# +import pybamm +import casadi +import numpy as np + + +class CasadiConverter(object): + def __init__(self, casadi_symbols=None): + self._casadi_symbols = casadi_symbols or {} + + def convert(self, symbol, t=None, y=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. + + Parameters + ---------- + symbol : :class:`pybamm.Symbol` + The symbol to convert + + Returns + ------- + CasADi symbol + The convert symbol + """ + + try: + return self._casadi_symbols[symbol.id] + except KeyError: + casadi_symbol = self._convert(symbol, t, y) + self._casadi_symbols[symbol.id] = casadi_symbol + + return casadi_symbol + + def _convert(self, symbol, t, y): + """ See :meth:`CasadiConverter.convert()`. """ + if isinstance(symbol, (pybamm.Scalar, pybamm.Array, pybamm.Time)): + return casadi.SX(symbol.evaluate(t, y)) + + elif isinstance(symbol, pybamm.StateVector): + if y is None: + raise ValueError("Must provide a 'y' for converting state vectors") + return casadi.vertcat(*[y[y_slice] for y_slice in symbol.y_slices]) + + 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) + if isinstance(symbol, pybamm.Outer): + return casadi.kron(converted_left, converted_right) + else: + # _binary_evaluate defined in derived classes for specific rules + return symbol._binary_evaluate(converted_left, converted_right) + + elif isinstance(symbol, pybamm.UnaryOperator): + converted_child = self.convert(symbol.child, t, y) + 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 + ] + # Special functions + if symbol.function == np.min: + return casadi.mmin(*converted_children) + elif symbol.function == np.max: + return casadi.mmax(*converted_children) + elif symbol.function == np.abs: + return casadi.fabs(*converted_children) + elif not isinstance( + symbol.function, pybamm.GetCurrent + ) and symbol.function.__name__.startswith("elementwise_grad_of_"): + differentiating_child_idx = int(symbol.function.__name__[-1]) + # Create dummy symbolic variables in order to differentiate using CasADi + dummy_vars = [ + casadi.SX.sym("y_" + str(i)) for i in range(len(converted_children)) + ] + func_diff = casadi.gradient( + symbol.differentiated_function(*dummy_vars), + dummy_vars[differentiating_child_idx], + ) + # Create function and evaluate it using the children + casadi_func_diff = casadi.Function("func_diff", dummy_vars, [func_diff]) + return casadi_func_diff(*converted_children) + # Other functions + else: + return symbol._function_evaluate(converted_children) + elif isinstance(symbol, pybamm.Concatenation): + converted_children = [ + self.convert(child, t, y) for child in symbol.children + ] + if isinstance(symbol, (pybamm.NumpyConcatenation, pybamm.SparseStack)): + return casadi.vertcat(*converted_children) + # DomainConcatenation specifies a particular ordering for the concatenation, + # which we must follow + elif isinstance(symbol, pybamm.DomainConcatenation): + slice_starts = [] + all_child_vectors = [] + for i in range(symbol.secondary_dimensions_npts): + child_vectors = [] + for child_var, slices in zip( + converted_children, symbol._children_slices + ): + for child_dom, child_slice in slices.items(): + slice_starts.append(symbol._slices[child_dom][i].start) + child_vectors.append( + child_var[child_slice[i].start : child_slice[i].stop] + ) + all_child_vectors.extend( + [v for _, v in sorted(zip(slice_starts, child_vectors))] + ) + return casadi.vertcat(*all_child_vectors) + + else: + raise TypeError( + """ + Cannot convert symbol of type '{}' to CasADi. Symbols must all be + 'linear algebra' at this stage. + """.format( + type(symbol) + ) + ) diff --git a/pybamm/expression_tree/evaluate.py b/pybamm/expression_tree/operations/evaluate.py similarity index 100% rename from pybamm/expression_tree/evaluate.py rename to pybamm/expression_tree/operations/evaluate.py diff --git a/pybamm/expression_tree/jacobian.py b/pybamm/expression_tree/operations/jacobian.py similarity index 100% rename from pybamm/expression_tree/jacobian.py rename to pybamm/expression_tree/operations/jacobian.py diff --git a/pybamm/expression_tree/simplify.py b/pybamm/expression_tree/operations/simplify.py similarity index 100% rename from pybamm/expression_tree/simplify.py rename to pybamm/expression_tree/operations/simplify.py diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py index 1aa170e9d8..688f07173a 100644 --- a/pybamm/expression_tree/state_vector.py +++ b/pybamm/expression_tree/state_vector.py @@ -113,7 +113,7 @@ def _base_evaluate(self, t=None, y=None): ) else: out = (y[: len(self._evaluation_array)])[self._evaluation_array] - if out.ndim == 1: + if isinstance(out, np.ndarray) and out.ndim == 1: out = out[:, np.newaxis] return out diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py index a189e539c9..8a94f05447 100644 --- a/pybamm/expression_tree/symbol.py +++ b/pybamm/expression_tree/symbol.py @@ -599,6 +599,13 @@ 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): + """ + Convert the expression tree to a CasADi expression tree. + See :class:`pybamm.CasadiConverter`. + """ + return pybamm.CasadiConverter(casadi_symbols).convert(self, t, y) + def new_copy(self): """ Make a new copy of a symbol, to avoid Tree corruption errors while bypassing diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py index 505988e7c7..9d192af47a 100644 --- a/pybamm/models/base_model.py +++ b/pybamm/models/base_model.py @@ -68,11 +68,18 @@ class BaseModel(object): Whether to simplify the expression tress representing the rhs and algebraic equations, Jacobain (if using) and events, before solving the model (default is True) - use_to_python : bool - Whether to convert the expression tress representing the rhs and - algebraic equations, Jacobain (if using) and events into pure python code - that will calculate the result of calling `evaluate(t, y)` on the given - expression tree (default is True) + convert_to_format : str + Whether to convert the expression trees representing the rhs and + algebraic equations, Jacobain (if using) and events into a different format: + + - None: keep PyBaMM expression tree structure. + - "python": convert into pure python code that will calculate the result of \ + calling `evaluate(t, y)` on the given expression treeself. + - "casadi": convert into CasADi expression tree, which then uses CasADi's \ + algorithm to calculate the Jacobian. + + Default is "python". + """ def __init__(self, name="Unnamed model"): @@ -96,7 +103,7 @@ def __init__(self, name="Unnamed model"): # Default behaviour is to use the jacobian and simplify self.use_jacobian = True self.use_simplify = True - self.use_to_python = True + self.convert_to_format = "python" def _set_dictionary(self, dict, name): """ diff --git a/pybamm/models/submodels/electrode/ohm/composite_ohm.py b/pybamm/models/submodels/electrode/ohm/composite_ohm.py index fa6950beb3..8cd5efd504 100644 --- a/pybamm/models/submodels/electrode/ohm/composite_ohm.py +++ b/pybamm/models/submodels/electrode/ohm/composite_ohm.py @@ -96,4 +96,3 @@ def set_boundary_conditions(self, variables): rbc = (-i_boundary_cc_0 / sigma_eff_0, "Neumann") self.boundary_conditions[phi_s] = {"left": lbc, "right": rbc} - diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py index ab3e5c7f0f..c8bf288fc1 100644 --- a/pybamm/models/submodels/particle/base_particle.py +++ b/pybamm/models/submodels/particle/base_particle.py @@ -53,11 +53,11 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav): + self.domain.lower() + " particle surface concentration [mol.m-3]": c_scale * c_s_surf_av, self.domain + " electrode active volume fraction": active_volume, + self.domain + " electrode volume-averaged concentration": c_s_r_av_vol, self.domain - + " electrode volume-averaged concentration": c_s_r_av_vol, - self.domain + " electrode " + + " electrode " + "volume-averaged concentration [mol.m-3]": c_s_r_av_vol * c_scale, - self.domain + " electrode average extent of lithiation": c_s_r_av + self.domain + " electrode average extent of lithiation": c_s_r_av, } return variables diff --git a/pybamm/parameters/electrical_parameters.py b/pybamm/parameters/electrical_parameters.py index 36c4f6b215..8db66dbe29 100644 --- a/pybamm/parameters/electrical_parameters.py +++ b/pybamm/parameters/electrical_parameters.py @@ -5,13 +5,6 @@ import numpy as np -def abs_non_zero(x): - if x == 0: # pragma: no cover - return 1 - else: - return abs(x) - - # -------------------------------------------------------------------------------------- # Dimensional Parameters I_typ = pybamm.Parameter("Typical current [A]") @@ -21,7 +14,7 @@ def abs_non_zero(x): "Number of electrodes connected in parallel to make a cell" ) i_typ = pybamm.Function( - abs_non_zero, (I_typ / (n_electrodes_parallel * pybamm.geometric_parameters.A_cc)) + np.abs, I_typ / (n_electrodes_parallel * pybamm.geometric_parameters.A_cc) ) voltage_low_cut_dimensional = pybamm.Parameter("Lower voltage cut-off [V]") voltage_high_cut_dimensional = pybamm.Parameter("Upper voltage cut-off [V]") diff --git a/pybamm/parameters/parameter_sets.py b/pybamm/parameters/parameter_sets.py index 943591434d..31c09106cc 100644 --- a/pybamm/parameters/parameter_sets.py +++ b/pybamm/parameters/parameter_sets.py @@ -46,4 +46,3 @@ "electrolyte": "sulfuric_acid_Sulzer2019", "experiment": "1C_discharge_from_full", } - diff --git a/pybamm/solvers/algebraic_solver.py b/pybamm/solvers/algebraic_solver.py index fe051cc344..3f085daf06 100644 --- a/pybamm/solvers/algebraic_solver.py +++ b/pybamm/solvers/algebraic_solver.py @@ -204,14 +204,14 @@ def set_up(self, model): pybamm.logger.info("Simplifying jacobian") jac = simp.simplify(jac) - if model.use_to_python: + if model.convert_to_format == "python": pybamm.logger.info("Converting jacobian to python") jac = pybamm.EvaluatorPython(jac) else: jac = None - if model.use_to_python: + if model.convert_to_format == "python": pybamm.logger.info("Converting algebraic to python") concatenated_algebraic = pybamm.EvaluatorPython(concatenated_algebraic) diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index ad1ad2b4a8..3b95b10f73 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -73,7 +73,10 @@ def solve(self, model, t_eval): # Set up timer = pybamm.Timer() start_time = timer.time() - self.set_up(model) + if model.convert_to_format == "casadi" or isinstance(self, pybamm.CasadiSolver): + self.set_up_casadi(model) + else: + self.set_up(model) set_up_time = timer.time() - start_time # Solve @@ -127,7 +130,12 @@ def step(self, model, dt, npts=2): # Run set up on first step if not hasattr(self, "y0"): start_time = timer.time() - self.set_up(model) + if model.convert_to_format == "casadi" or isinstance( + self, pybamm.CasadiSolver + ): + self.set_up_casadi(model) + else: + self.set_up(model) self.t = 0.0 set_up_time = timer.time() - start_time else: @@ -187,11 +195,17 @@ def set_up(self, model): The model whose solution to calculate. Must have attributes rhs and initial_conditions - Raises - ------ - :class:`pybamm.SolverError` - If the model contains any algebraic equations (in which case a DAE solver - should be used instead) + """ + raise NotImplementedError + + def set_up_casadi(self, model): + """Convert model to casadi format and use their inbuilt functionalities. + + Parameters + ---------- + model : :class:`pybamm.BaseModel` + The model whose solution to calculate. Must have attributes rhs and + initial_conditions """ raise NotImplementedError diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py new file mode 100644 index 0000000000..373bbc0563 --- /dev/null +++ b/pybamm/solvers/casadi_solver.py @@ -0,0 +1,103 @@ +# +# CasADi Solver class +# +import casadi +import pybamm +import numpy as np + + +class CasadiSolver(pybamm.DaeSolver): + """Solve a discretised model, using CasADi. + + Parameters + ---------- + rtol : float, optional + The relative tolerance for the solver (default is 1e-6). + atol : float, optional + The absolute tolerance for the solver (default is 1e-6). + + **Extends**: :class:`pybamm.DaeSolver` + """ + + def __init__( + self, + method="idas", + rtol=1e-6, + atol=1e-6, + root_method="lm", + root_tol=1e-6, + max_steps=1000, + **extra_options, + ): + super().__init__(method, rtol, atol, root_method, root_tol, max_steps) + self.extra_options = extra_options + + def compute_solution(self, model, t_eval): + """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. + + Parameters + ---------- + model : :class:`pybamm.BaseModel` + The model whose solution to calculate. Must have attributes rhs and + initial_conditions + t_eval : numeric type + The times at which to compute the solution + + """ + timer = pybamm.Timer() + + solve_start_time = timer.time() + pybamm.logger.info("Calling DAE solver") + solution = self.integrate_casadi( + self.casadi_problem, self.y0, t_eval, mass_matrix=model.mass_matrix.entries + ) + solve_time = timer.time() - solve_start_time + + # Events not implemented, termination is always 'final time' + termination = "final time" + + return solution, solve_time, termination + + def integrate_casadi(self, problem, y0, t_eval, mass_matrix=None): + """ + Solve a DAE model defined by residuals with initial conditions y0. + + Parameters + ---------- + residuals : method + A function that takes in t, y and ydot and returns the residuals of the + equations + y0 : numeric type + The initial conditions + t_eval : numeric type + The times at which to compute the solution + 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. + """ + options = { + "grid": t_eval, + "reltol": self.rtol, + "abstol": self.atol, + "output_t0": True, + } + options.update(self.extra_options) + if self.method == "idas": + options["calc_ic"] = True + + # set up and solve + integrator = casadi.integrator("F", self.method, problem, options) + try: + # Try solving + len_rhs = problem["x"].size()[0] + y0_diff, y0_alg = np.split(y0, [len_rhs]) + sol = integrator(x0=y0_diff, z0=y0_alg) + y_values = np.concatenate([sol["xf"].full(), sol["zf"].full()]) + return pybamm.Solution(t_eval, y_values, None, None, "final time") + 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 ba14d26d2c..39fe18f364 100644 --- a/pybamm/solvers/dae_solver.py +++ b/pybamm/solvers/dae_solver.py @@ -1,6 +1,7 @@ # # Base solver class # +import casadi import pybamm import numpy as np from scipy import optimize @@ -102,12 +103,6 @@ def set_up(self, model): model : :class:`pybamm.BaseModel` The model whose solution to calculate. Must have attributes rhs and initial_conditions - - Raises - ------ - :class:`pybamm.SolverError` - If the model contains any algebraic equations (in which case a DAE solver - should be used instead) """ # create simplified rhs, algebraic and event expressions concatenated_rhs = model.concatenated_rhs @@ -144,19 +139,22 @@ def set_up(self, model): jac_algebraic = simp.simplify(jac_algebraic) jac = simp.simplify(jac) - if model.use_to_python: + if model.convert_to_format == "python": pybamm.logger.info("Converting jacobian to python") jac_algebraic = pybamm.EvaluatorPython(jac_algebraic) jac = pybamm.EvaluatorPython(jac) - def jac_alg_fn(t, y): - return jac_algebraic.evaluate(t, y) + def jacobian(t, y): + return jac.evaluate(t, y, known_evals={})[0] + + def jacobian_alg(t, y): + return jac_algebraic.evaluate(t, y, known_evals={})[0] else: - jac = None - jac_alg_fn = None + jacobian = None + jacobian_alg = None - if model.use_to_python: + if model.convert_to_format == "python": pybamm.logger.info("Converting RHS to python") concatenated_rhs = pybamm.EvaluatorPython(concatenated_rhs) pybamm.logger.info("Converting algebraic to python") @@ -175,7 +173,10 @@ def algebraic(t, y): if len(model.algebraic) > 0: y0 = self.calculate_consistent_initial_conditions( - rhs, algebraic, model.concatenated_initial_conditions[:, 0], jac_alg_fn + rhs, + algebraic, + model.concatenated_initial_conditions[:, 0], + jacobian_alg, ) else: # can use DAE solver to solve ODE model @@ -206,14 +207,113 @@ def eval_event(t, y): event_funs = [event_fun(event) for event in events.values()] - # Create function to evaluate jacobian - if jac is not None: + # Add the solver attributes + self.y0 = y0 + self.rhs = rhs + self.algebraic = algebraic + self.residuals = residuals + self.events = events + self.event_funs = event_funs + self.jacobian = jacobian + + def set_up_casadi(self, model): + """Convert model to casadi format and use their inbuilt functionalities. + + Parameters + ---------- + model : :class:`pybamm.BaseModel` + The model whose solution to calculate. Must have attributes rhs and + initial_conditions + """ + # Convert model attributes to casadi + t_casadi = casadi.SX.sym("t") + y0 = model.concatenated_initial_conditions + y_diff = casadi.SX.sym("y_diff", len(model.concatenated_rhs.evaluate(0, y0))) + y_alg = casadi.SX.sym( + "y_alg", len(model.concatenated_algebraic.evaluate(0, y0)) + ) + y_casadi = casadi.vertcat(y_diff, y_alg) + pybamm.logger.info("Converting RHS to CasADi") + concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi) + pybamm.logger.info("Converting algebraic to CasADi") + concatenated_algebraic = model.concatenated_algebraic.to_casadi( + t_casadi, y_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) + for name, event in model.events.items() + } + + # Create functions to evaluate rhs and algebraic + concatenated_rhs_fn = casadi.Function( + "rhs", [t_casadi, y_casadi], [concatenated_rhs] + ) + concatenated_algebraic_fn = casadi.Function( + "algebraic", [t_casadi, y_casadi], [concatenated_algebraic] + ) + all_states_fn = casadi.Function("all", [t_casadi, y_casadi], [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], [casadi_jac] + ) + casadi_jac_alg = casadi.jacobian(concatenated_algebraic, y_casadi) + casadi_jac_alg_fn = casadi.Function( + "jacobian", [t_casadi, y_casadi], [casadi_jac_alg] + ) def jacobian(t, y): - return jac.evaluate(t, y, known_evals={})[0] + return casadi_jac_fn(t, y) + + def jacobian_alg(t, y): + return casadi_jac_alg_fn(t, y) else: jacobian = None + jacobian_alg = None + + # Calculate consistent initial conditions for the algebraic equations + def rhs(t, y): + return concatenated_rhs_fn(t, y).full()[:, 0] + + def algebraic(t, y): + return concatenated_algebraic_fn(t, y).full()[:, 0] + + if len(model.algebraic) > 0: + y0 = self.calculate_consistent_initial_conditions( + rhs, + algebraic, + model.concatenated_initial_conditions[:, 0], + jacobian_alg, + ) + else: + # can use DAE solver to solve ODE model + y0 = model.concatenated_initial_conditions[:, 0] + + # Create functions to evaluate residuals + + def residuals(t, y, ydot): + pybamm.logger.debug( + "Evaluating residuals for {} at t={}".format(model.name, t) + ) + states_eval = all_states_fn(t, y).full()[:, 0] + return states_eval - model.mass_matrix.entries @ ydot + + # Create event-dependent function to evaluate events + def event_fun(event): + casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi], [event]) + + def eval_event(t, y): + return casadi_event_fn(t, y) + + return eval_event + + event_funs = [event_fun(event) for event in casadi_events.values()] # Add the solver attributes # Note: these are the (possibly) converted to python version rhs, algebraic @@ -222,10 +322,18 @@ def jacobian(t, y): self.rhs = rhs self.algebraic = algebraic self.residuals = residuals - self.events = events + self.events = model.events self.event_funs = event_funs self.jacobian = jacobian + # Create CasADi problem for the CasADi solver + self.casadi_problem = { + "x": y_diff, + "z": y_alg, + "ode": concatenated_rhs, + "alg": concatenated_algebraic, + } + def calculate_consistent_initial_conditions( self, rhs, algebraic, y0_guess, jac=None ): diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py index b055c95db5..131f0fdd8b 100644 --- a/pybamm/solvers/idaklu_solver.py +++ b/pybamm/solvers/idaklu_solver.py @@ -17,7 +17,7 @@ def have_idaklu(): return idaklu_spec is None -class IDAKLU(pybamm.DaeSolver): +class IDAKLUSolver(pybamm.DaeSolver): """Solve a discretised model, using sundials with the KLU sparse linear solver. Parameters diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py index 973c10401c..9b650241f0 100644 --- a/pybamm/solvers/ode_solver.py +++ b/pybamm/solvers/ode_solver.py @@ -1,6 +1,7 @@ # # Base solver class # +import casadi import pybamm import numpy as np @@ -102,13 +103,13 @@ def set_up(self, model): pybamm.logger.info("Simplifying jacobian") jac_rhs = simp.simplify(jac_rhs) - if model.use_to_python: + if model.convert_to_format == "python": pybamm.logger.info("Converting jacobian to python") jac_rhs = pybamm.EvaluatorPython(jac_rhs) else: jac_rhs = None - if model.use_to_python: + if model.convert_to_format == "python": pybamm.logger.info("Converting RHS to python") concatenated_rhs = pybamm.EvaluatorPython(concatenated_rhs) pybamm.logger.info("Converting events to python") @@ -150,6 +151,82 @@ def jacobian(t, y): self.event_funs = event_funs self.jacobian = jacobian + def set_up_casadi(self, model): + """Convert model to casadi format and use their inbuilt functionalities. + + Parameters + ---------- + model : :class:`pybamm.BaseModel` + The model whose solution to calculate. Must have attributes rhs and + initial_conditions + + Raises + ------ + :class:`pybamm.SolverError` + If the model contains any algebraic equations (in which case a DAE solver + should be used instead) + + """ + # Check for algebraic equations + if len(model.algebraic) > 0: + raise pybamm.SolverError( + """Cannot use ODE solver to solve model with DAEs""" + ) + + y0 = model.concatenated_initial_conditions[:, 0] + + t_casadi = casadi.SX.sym("t") + y_casadi = casadi.SX.sym("y", len(y0)) + pybamm.logger.info("Converting RHS to CasADi") + concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi) + pybamm.logger.info("Converting events to CasADi") + casadi_events = { + name: event.to_casadi(t_casadi, y_casadi) + for name, event in model.events.items() + } + + # Create function to evaluate rhs + concatenated_rhs_fn = casadi.Function( + "rhs", [t_casadi, y_casadi], [concatenated_rhs] + ) + + def dydt(t, y): + pybamm.logger.debug("Evaluating RHS for {} at t={}".format(model.name, t)) + dy = concatenated_rhs_fn(t, y).full() + return dy[:, 0] + + # Create event-dependent function to evaluate events + def event_fun(event): + casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi], [event]) + + def eval_event(t, y): + return casadi_event_fn(t, y) + + return eval_event + + event_funs = [event_fun(event) for event in casadi_events.values()] + + # 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], [casadi_jac] + ) + + def jacobian(t, y): + return casadi_jac_fn(t, y) + + else: + jacobian = None + + # Add the solver attributes + self.y0 = y0 + self.dydt = dydt + self.events = model.events + self.event_funs = event_funs + self.jacobian = jacobian + def integrate( self, derivs, y0, t_eval, events=None, mass_matrix=None, jacobian=None ): diff --git a/pybamm/solvers/scipy_solver.py b/pybamm/solvers/scipy_solver.py index fd7afbf498..23e11483c5 100644 --- a/pybamm/solvers/scipy_solver.py +++ b/pybamm/solvers/scipy_solver.py @@ -83,7 +83,7 @@ def integrate( termination = "event" t_event = [] for time in sol.t_events: - if time: + if time.size > 0: t_event = np.append(t_event, np.max(time)) t_event = np.array([np.max(t_event)]) y_event = sol.sol(t_event) diff --git a/results/2plus1D/dfn_2plus1D.py b/results/2plus1D/dfn_2plus1D.py index 9beff8d5a4..86cbad603c 100644 --- a/results/2plus1D/dfn_2plus1D.py +++ b/results/2plus1D/dfn_2plus1D.py @@ -46,7 +46,7 @@ tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) t_end = 3600 / tau.evaluate(0) t_eval = np.linspace(0, t_end, 120) -solution = pybamm.IDAKLU().solve(model, t_eval) +solution = pybamm.IDAKLUSolver().solve(model, t_eval) # TO DO: 2+1D automated plotting phi_s_cn = pybamm.ProcessedVariable( diff --git a/results/2plus1D/spm_2plus1D.py b/results/2plus1D/spm_2plus1D.py index 1b9a08411c..943a1afdd5 100644 --- a/results/2plus1D/spm_2plus1D.py +++ b/results/2plus1D/spm_2plus1D.py @@ -46,7 +46,7 @@ tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) t_end = 3600 / tau.evaluate(0) t_eval = np.linspace(0, t_end, 120) -solution = pybamm.IDAKLU().solve(model, t_eval) +solution = pybamm.IDAKLUSolver().solve(model, t_eval) # TO DO: 2+1D automated plotting phi_s_cn = pybamm.ProcessedVariable( diff --git a/results/2plus1D/spme_2plus1D.py b/results/2plus1D/spme_2plus1D.py index d356f1a62f..1cef748871 100644 --- a/results/2plus1D/spme_2plus1D.py +++ b/results/2plus1D/spme_2plus1D.py @@ -46,7 +46,7 @@ tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) t_end = 3600 / tau.evaluate(0) t_eval = np.linspace(0, t_end, 120) -solution = pybamm.IDAKLU().solve(model, t_eval) +solution = pybamm.IDAKLUSolver().solve(model, t_eval) # TO DO: 2+1D automated plotting phi_s_cn = pybamm.ProcessedVariable( diff --git a/setup.py b/setup.py index 14227d4205..8e9dac78cd 100644 --- a/setup.py +++ b/setup.py @@ -45,6 +45,7 @@ def load_version(): "anytree>=2.4.3", "autograd>=1.2", "scikit-fem>=0.2.0", + "casadi>=3.5.0", # Note: Matplotlib is loaded for debug plots, but to ensure pybamm runs # on systems without an attached display, it should never be imported # outside of plot() methods. diff --git a/tests/integration/test_models/standard_model_tests.py b/tests/integration/test_models/standard_model_tests.py index 2b5afa1244..ba827e9cd3 100644 --- a/tests/integration/test_models/standard_model_tests.py +++ b/tests/integration/test_models/standard_model_tests.py @@ -172,6 +172,7 @@ def evaluate_model(self, simplify=False, use_known_evals=False, to_python=False) def set_up_model(self, simplify=False, to_python=False): self.model.use_simplify = simplify - self.model.use_to_python = to_python + if to_python is True: + self.model.convert_to_format = "python" self.model.default_solver.set_up(self.model) return None diff --git a/tests/integration/test_solvers/test_idaklu.py b/tests/integration/test_solvers/test_idaklu.py index 7885d54433..ea78201f74 100644 --- a/tests/integration/test_solvers/test_idaklu.py +++ b/tests/integration/test_solvers/test_idaklu.py @@ -16,7 +16,7 @@ def test_on_spme(self): disc = pybamm.Discretisation(mesh, model.default_spatial_methods) disc.process_model(model) t_eval = np.linspace(0, 0.2, 100) - solution = pybamm.IDAKLU().solve(model, t_eval) + solution = pybamm.IDAKLUSolver().solve(model, t_eval) np.testing.assert_array_less(1, solution.t.size) diff --git a/tests/unit/test_expression_tree/test_functions.py b/tests/unit/test_expression_tree/test_functions.py index bb10dd2429..28d1586b36 100644 --- a/tests/unit/test_expression_tree/test_functions.py +++ b/tests/unit/test_expression_tree/test_functions.py @@ -21,6 +21,10 @@ def test_multi_var_function(arg1, arg2): return arg1 + arg2 +def test_multi_var_function_cube(arg1, arg2): + return arg1 + arg2 ** 3 + + class TestFunction(unittest.TestCase): def test_constant_functions(self): d = pybamm.Scalar(6) @@ -52,8 +56,9 @@ def test_function_of_one_variable(self): logvar.evaluate(y=y, known_evals={})[0], np.log1p(y) ) - def test_with_autograd(self): + def test_diff(self): a = pybamm.StateVector(slice(0, 1)) + b = pybamm.StateVector(slice(1, 2)) y = np.array([5]) func = pybamm.Function(test_function, a) self.assertEqual(func.diff(a).evaluate(y=y), 2) @@ -68,6 +73,21 @@ def test_with_autograd(self): # multiple variables func = pybamm.Function(test_multi_var_function, 4 * a, 3 * a) self.assertEqual(func.diff(a).evaluate(y=y), 7) + func = pybamm.Function(test_multi_var_function, 4 * a, 3 * b) + self.assertEqual(func.diff(a).evaluate(y=np.array([5, 6])), 4) + self.assertEqual(func.diff(b).evaluate(y=np.array([5, 6])), 3) + func = pybamm.Function(test_multi_var_function_cube, 4 * a, 3 * b) + self.assertEqual(func.diff(a).evaluate(y=np.array([5, 6])), 4) + self.assertEqual( + func.diff(b).evaluate(y=np.array([5, 6])), 3 * 3 * (3 * 6) ** 2 + ) + + # exceptions + func = pybamm.Function( + test_multi_var_function_cube, 4 * a, 3 * b, derivative="derivative" + ) + with self.assertRaises(ValueError): + func.diff(a) def test_function_of_multiple_variables(self): a = pybamm.Variable("a") diff --git a/tests/unit/test_expression_tree/test_operations/__init__.py b/tests/unit/test_expression_tree/test_operations/__init__.py new file mode 100644 index 0000000000..e69de29bb2 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 new file mode 100644 index 0000000000..a516cb7955 --- /dev/null +++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py @@ -0,0 +1,151 @@ +# +# Test for the Simplify class +# +import casadi +import numpy as np +import autograd.numpy as anp +import pybamm +import unittest +from tests import get_mesh_for_testing, get_1p1d_discretisation_for_testing + + +class TestCasadiConverter(unittest.TestCase): + def test_convert_scalar_symbols(self): + a = pybamm.Scalar(0) + b = pybamm.Scalar(1) + c = pybamm.Scalar(-1) + d = pybamm.Scalar(2) + + self.assertEqual(a.to_casadi(), casadi.SX(0)) + self.assertEqual(d.to_casadi(), casadi.SX(2)) + + # negate + self.assertEqual((-b).to_casadi(), casadi.SX(-1)) + # absolute value + self.assertEqual(abs(c).to_casadi(), casadi.SX(1)) + + # function + def sin(x): + return np.sin(x) + + f = pybamm.Function(sin, b) + self.assertEqual(f.to_casadi(), casadi.SX(np.sin(1))) + + def myfunction(x, y): + return x + y + + f = pybamm.Function(myfunction, b, d) + self.assertEqual(f.to_casadi(), casadi.SX(3)) + + # addition + self.assertEqual((a + b).to_casadi(), casadi.SX(1)) + # subtraction + self.assertEqual((c - d).to_casadi(), casadi.SX(-3)) + # multiplication + self.assertEqual((c * d).to_casadi(), casadi.SX(-2)) + # power + self.assertEqual((c ** d).to_casadi(), casadi.SX(1)) + # division + self.assertEqual((b / d).to_casadi(), casadi.SX(1 / 2)) + + def test_convert_array_symbols(self): + # Arrays + a = np.array([1, 2, 3, 4, 5]) + pybamm_a = pybamm.Array(a) + self.assertTrue(casadi.is_equal(pybamm_a.to_casadi(), casadi.SX(a))) + + casadi_t = casadi.SX.sym("t") + casadi_y = casadi.SX.sym("y", 10) + + pybamm_t = pybamm.Time() + pybamm_y = pybamm.StateVector(slice(0, 10)) + + # Time + self.assertEqual(pybamm_t.to_casadi(casadi_t, casadi_y), casadi_t) + + # State Vector + self.assertTrue( + casadi.is_equal(pybamm_y.to_casadi(casadi_t, casadi_y), casadi_y) + ) + + # outer product + outer = pybamm.Outer(pybamm_a, pybamm_a) + self.assertTrue(casadi.is_equal(outer.to_casadi(), casadi.SX(outer.evaluate()))) + + def test_special_functions(self): + a = pybamm.Array(np.array([1, 2, 3, 4, 5])) + self.assertEqual(pybamm.max(a).to_casadi(), casadi.SX(5)) + self.assertEqual(pybamm.min(a).to_casadi(), casadi.SX(1)) + b = pybamm.Array(np.array([-2])) + c = pybamm.Array(np.array([3])) + self.assertEqual(pybamm.Function(np.abs, b).to_casadi(), casadi.SX(2)) + self.assertEqual(pybamm.Function(np.abs, c).to_casadi(), casadi.SX(3)) + + def test_concatenations(self): + y = np.linspace(0, 1, 10)[:, np.newaxis] + a = pybamm.Vector(y) + b = pybamm.Scalar(16) + c = pybamm.Scalar(3) + conc = pybamm.NumpyConcatenation(a, b, c) + self.assertTrue(casadi.is_equal(conc.to_casadi(), casadi.SX(conc.evaluate()))) + + # Domain concatenation + mesh = get_mesh_for_testing() + a_dom = ["negative electrode"] + b_dom = ["positive electrode"] + a = 2 * pybamm.Vector(np.ones_like(mesh[a_dom[0]][0].nodes), domain=a_dom) + b = pybamm.Vector(np.ones_like(mesh[b_dom[0]][0].nodes), domain=b_dom) + conc = pybamm.DomainConcatenation([b, a], mesh) + self.assertTrue(casadi.is_equal(conc.to_casadi(), casadi.SX(conc.evaluate()))) + + # 2d + disc = get_1p1d_discretisation_for_testing() + a = pybamm.Variable("a", domain=a_dom) + b = pybamm.Variable("b", domain=b_dom) + conc = pybamm.Concatenation(a, b) + disc.set_variable_slices([conc]) + expr = disc.process_symbol(conc) + y = casadi.SX.sym("y", expr.size) + x = expr.to_casadi(None, y) + f = casadi.Function("f", [x], [x]) + y_eval = np.linspace(0, 1, expr.size) + self.assertTrue(casadi.is_equal(f(y_eval), casadi.SX(expr.evaluate(y=y_eval)))) + + def test_convert_differentiated_function(self): + a = pybamm.Scalar(0) + b = pybamm.Scalar(1) + + # function + def sin(x): + return anp.sin(x) + + f = pybamm.Function(sin, b).diff(b) + self.assertEqual(f.to_casadi(), casadi.SX(np.cos(1))) + + def myfunction(x, y): + return x + y ** 3 + + f = pybamm.Function(myfunction, a, b).diff(a) + self.assertEqual(f.to_casadi(), casadi.SX(1)) + f = pybamm.Function(myfunction, a, b).diff(b) + self.assertEqual(f.to_casadi(), casadi.SX(3)) + + def test_errors(self): + y = pybamm.StateVector(slice(0, 10)) + with self.assertRaisesRegex( + ValueError, "Must provide a 'y' for converting state vectors" + ): + y.to_casadi() + var = pybamm.Variable("var") + with self.assertRaisesRegex(TypeError, "Cannot convert symbol of type"): + var.to_casadi() + + +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_copy.py b/tests/unit/test_expression_tree/test_operations/test_copy.py similarity index 100% rename from tests/unit/test_expression_tree/test_copy.py rename to tests/unit/test_expression_tree/test_operations/test_copy.py diff --git a/tests/unit/test_expression_tree/test_evaluate.py b/tests/unit/test_expression_tree/test_operations/test_evaluate.py similarity index 100% rename from tests/unit/test_expression_tree/test_evaluate.py rename to tests/unit/test_expression_tree/test_operations/test_evaluate.py diff --git a/tests/unit/test_expression_tree/test_jac.py b/tests/unit/test_expression_tree/test_operations/test_jac.py similarity index 100% rename from tests/unit/test_expression_tree/test_jac.py rename to tests/unit/test_expression_tree/test_operations/test_jac.py diff --git a/tests/unit/test_expression_tree/test_jac_2D.py b/tests/unit/test_expression_tree/test_operations/test_jac_2D.py similarity index 100% rename from tests/unit/test_expression_tree/test_jac_2D.py rename to tests/unit/test_expression_tree/test_operations/test_jac_2D.py diff --git a/tests/unit/test_expression_tree/test_simplify.py b/tests/unit/test_expression_tree/test_operations/test_simplify.py similarity index 100% rename from tests/unit/test_expression_tree/test_simplify.py rename to tests/unit/test_expression_tree/test_operations/test_simplify.py diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py new file mode 100644 index 0000000000..6d6b45eecd --- /dev/null +++ b/tests/unit/test_solvers/test_casadi_solver.py @@ -0,0 +1,133 @@ +# +# Tests for the Casadi Solver class +# +import casadi +import pybamm +import unittest +import numpy as np +from tests import get_mesh_for_testing +import warnings + + +class TestCasadiSolver(unittest.TestCase): + def test_integrate(self): + # Constant + solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8, method="idas") + + y = casadi.SX.sym("y") + constant_growth = casadi.SX(0.5) + problem = {"x": y, "ode": constant_growth} + + y0 = np.array([0]) + t_eval = np.linspace(0, 1, 100) + solution = solver.integrate_casadi(problem, y0, t_eval) + np.testing.assert_array_equal(solution.t, t_eval) + np.testing.assert_allclose(0.5 * solution.t, solution.y[0]) + + # Exponential decay + solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8, method="cvodes") + + exponential_decay = -0.1 * y + problem = {"x": y, "ode": exponential_decay} + + y0 = np.array([1]) + t_eval = np.linspace(0, 1, 100) + solution = solver.integrate_casadi(problem, y0, t_eval) + 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") + + y = casadi.SX.sym("y") + sqrt_decay = -np.sqrt(y) + + y0 = np.array([1]) + t_eval = np.linspace(0, 3, 100) + solver = pybamm.CasadiSolver( + rtol=1e-8, + atol=1e-8, + method="idas", + disable_internal_warnings=True, + regularity_check=False, + ) + problem = {"x": y, "ode": sqrt_decay} + # Expect solver to fail when y goes negative + with self.assertRaises(pybamm.SolverError): + solver.integrate_casadi(problem, y0, t_eval) + + # Turn warnings back on + warnings.simplefilter("default") + + def test_model_solver(self): + # Create model + model = pybamm.BaseModel() + domain = ["negative electrode", "separator", "positive electrode"] + var = pybamm.Variable("var", domain=domain) + model.rhs = {var: 0.1 * var} + model.initial_conditions = {var: 1} + # 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.CasadiSolver(rtol=1e-8, atol=1e-8, method="idas") + t_eval = np.linspace(0, 1, 100) + solution = solver.solve(model, t_eval) + np.testing.assert_array_equal(solution.t, t_eval) + np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t)) + + # Test time + self.assertGreater( + solution.total_time, solution.solve_time + solution.set_up_time + ) + + def test_model_step(self): + # Create model + model = pybamm.BaseModel() + domain = ["negative electrode", "separator", "positive electrode"] + var = pybamm.Variable("var", domain=domain) + model.rhs = {var: 0.1 * var} + model.initial_conditions = {var: 1} + # 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) + + solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8, method="idas") + + # Step once + dt = 0.1 + step_sol = solver.step(model, dt) + np.testing.assert_array_equal(step_sol.t, [0, dt]) + np.testing.assert_allclose(step_sol.y[0], np.exp(0.1 * step_sol.t)) + + # Step again (return 5 points) + step_sol_2 = solver.step(model, dt, npts=5) + np.testing.assert_array_equal(step_sol_2.t, np.linspace(dt, 2 * dt, 5)) + np.testing.assert_allclose(step_sol_2.y[0], np.exp(0.1 * step_sol_2.t)) + + # append solutions + step_sol.append(step_sol_2) + + # Check steps give same solution as solve + t_eval = step_sol.t + solution = solver.solve(model, t_eval) + np.testing.assert_allclose(solution.y[0], step_sol.y[0]) + + +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_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py index a644ec7e3b..6d442e12d0 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -53,7 +53,7 @@ def rhs(t, y): def alg(t, y): return np.array([1 - y[1]]) - solver = pybamm.IDAKLU() + solver = pybamm.IDAKLUSolver() solver.residuals = res solver.rhs = rhs solver.algebraic = alg diff --git a/tests/unit/test_solvers/test_ode_solver.py b/tests/unit/test_solvers/test_ode_solver.py index 009dbcca23..5469e4f726 100644 --- a/tests/unit/test_solvers/test_ode_solver.py +++ b/tests/unit/test_solvers/test_ode_solver.py @@ -24,6 +24,10 @@ def test_wrong_solver(self): pybamm.SolverError, "Cannot use ODE solver to solve model with DAEs" ): solver.solve(model, None) + with self.assertRaisesRegex( + pybamm.SolverError, "Cannot use ODE solver to solve model with DAEs" + ): + solver.set_up_casadi(model) if __name__ == "__main__": diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py index e8b503f2fe..c0977fd709 100644 --- a/tests/unit/test_solvers/test_scikits_solvers.py +++ b/tests/unit/test_solvers/test_scikits_solvers.py @@ -676,6 +676,57 @@ def test_model_step_dae(self): np.testing.assert_allclose(solution.y[0], step_sol.y[0, :]) np.testing.assert_allclose(solution.y[-1], step_sol.y[-1, :]) + def test_model_solver_ode_events_casadi(self): + # Create model + model = pybamm.BaseModel() + model.convert_to_format = "casadi" + whole_cell = ["negative electrode", "separator", "positive electrode"] + var = pybamm.Variable("var", domain=whole_cell) + model.rhs = {var: 0.1 * var} + model.initial_conditions = {var: 1} + model.events = { + "2 * var = 2.5": pybamm.min(2 * var - 2.5), + "var = 1.5": pybamm.min(var - 1.5), + } + disc = get_discretisation_for_testing() + disc.process_model(model) + + # Solve + solver = pybamm.ScikitsOdeSolver(rtol=1e-9, atol=1e-9) + t_eval = np.linspace(0, 10, 100) + solution = solver.solve(model, t_eval) + np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t)) + np.testing.assert_array_less(solution.y[0], 1.5) + np.testing.assert_array_less(solution.y[0], 1.25) + + def test_model_solver_dae_events_casadi(self): + # Create model + model = pybamm.BaseModel() + for use_jacobian in [True, False]: + model.use_jacobian = use_jacobian + model.convert_to_format = "casadi" + whole_cell = ["negative electrode", "separator", "positive electrode"] + var1 = pybamm.Variable("var1", domain=whole_cell) + var2 = pybamm.Variable("var2", domain=whole_cell) + model.rhs = {var1: 0.1 * var1} + model.algebraic = {var2: 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) + 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)) + 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 1e5aaeebb9..f5d7e783ce 100644 --- a/tests/unit/test_solvers/test_scipy_solver.py +++ b/tests/unit/test_solvers/test_scipy_solver.py @@ -268,6 +268,33 @@ 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_event_with_casadi(self): + # Create model + model = pybamm.BaseModel() + for use_jacobian in [True, False]: + model.use_jacobian = use_jacobian + model.convert_to_format = "casadi" + domain = ["negative electrode", "separator", "positive electrode"] + var = pybamm.Variable("var", domain=domain) + model.rhs = {var: -0.1 * 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) + 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")