From 802cfd0ca50a50b09436b8356c9c9250c09e5e95 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Tue, 15 Oct 2019 09:41:14 -0400 Subject: [PATCH 01/23] #664 convert scalars and operations to casadi --- .requirements-docs.txt | 1 + docs/source/expression_tree/index.rst | 3 +- .../operations/convert_to_casadi.rst | 5 ++ .../{ => operations}/evaluate.rst | 0 .../expression_tree/operations/index.rst | 10 +++ .../{ => operations}/simplify.rst | 0 docs/source/solvers/casadi_solver.rst | 5 ++ docs/source/solvers/index.rst | 1 + pybamm/__init__.py | 7 +- pybamm/expression_tree/operations/__init__.py | 0 .../operations/convert_to_casadi.py | 75 +++++++++++++++++++ .../{ => operations}/evaluate.py | 0 .../{ => operations}/simplify.py | 0 pybamm/expression_tree/symbol.py | 7 ++ pybamm/solvers/casadi_solver.py | 9 +++ setup.py | 1 + .../test_operations/__init__.py | 0 .../test_operations/test_convert_to_casadi.py | 62 +++++++++++++++ .../{ => test_operations}/test_copy.py | 0 .../{ => test_operations}/test_evaluate.py | 0 .../{ => test_operations}/test_jac.py | 0 .../{ => test_operations}/test_jac_2D.py | 0 .../{ => test_operations}/test_simplify.py | 0 23 files changed, 182 insertions(+), 4 deletions(-) create mode 100644 docs/source/expression_tree/operations/convert_to_casadi.rst rename docs/source/expression_tree/{ => operations}/evaluate.rst (100%) create mode 100644 docs/source/expression_tree/operations/index.rst rename docs/source/expression_tree/{ => operations}/simplify.rst (100%) create mode 100644 docs/source/solvers/casadi_solver.rst create mode 100644 pybamm/expression_tree/operations/__init__.py create mode 100644 pybamm/expression_tree/operations/convert_to_casadi.py rename pybamm/expression_tree/{ => operations}/evaluate.py (100%) rename pybamm/expression_tree/{ => operations}/simplify.py (100%) create mode 100644 pybamm/solvers/casadi_solver.py create mode 100644 tests/unit/test_expression_tree/test_operations/__init__.py create mode 100644 tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py rename tests/unit/test_expression_tree/{ => test_operations}/test_copy.py (100%) rename tests/unit/test_expression_tree/{ => test_operations}/test_evaluate.py (100%) rename tests/unit/test_expression_tree/{ => test_operations}/test_jac.py (100%) rename tests/unit/test_expression_tree/{ => test_operations}/test_jac_2D.py (100%) rename tests/unit/test_expression_tree/{ => test_operations}/test_simplify.py (100%) 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/docs/source/expression_tree/index.rst b/docs/source/expression_tree/index.rst index 2e95d07b16..3eedd7b622 100644 --- a/docs/source/expression_tree/index.rst +++ b/docs/source/expression_tree/index.rst @@ -15,7 +15,6 @@ Expression Tree unary_operator concatenations broadcasts - simplify functions interpolant - evaluate + 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..108c6eeedf --- /dev/null +++ b/docs/source/expression_tree/operations/index.rst @@ -0,0 +1,10 @@ +Operations on expression trees +============================== + +Classes and functions that operate on the expression tree + +.. toctree:: + + simplify + evaluate + convert_to_casadi 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/pybamm/__init__.py b/pybamm/__init__.py index 7723eb8742..dda95eea37 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -150,18 +150,21 @@ 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.evaluate import ( +from .expression_tree.operations.evaluate import ( find_symbols, id_to_python_variable, to_python, EvaluatorPython, ) +from .expression_tree.operations.convert_to_casadi import CasadiConverter # # Model classes 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..41a7251c33 --- /dev/null +++ b/pybamm/expression_tree/operations/convert_to_casadi.py @@ -0,0 +1,75 @@ +# +# Convert a PyBaMM expression tree to a CasADi expression tree +# +import pybamm +import casadi + + +class CasadiConverter(object): + def __init__(self, casadi_symbols=None): + self._casadi_symbols = casadi_symbols or {} + + def convert(self, symbol): + """ + 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) + self._casadi_symbols[symbol.id] = casadi_symbol + + return casadi_symbol + + def _convert(self, symbol): + """ See :meth:`Simplification.convert()`. """ + if isinstance(symbol, pybamm.Scalar): + return casadi.SX(symbol.evaluate()) + + if isinstance(symbol, pybamm.BinaryOperator): + left, right = symbol.children + # process children + converted_left = self.convert(left) + converted_right = self.convert(right) + # _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) + if isinstance(symbol, pybamm.AbsoluteValue): + return casadi.fabs(converted_child) + return symbol._unary_evaluate(converted_child) + + elif isinstance(symbol, pybamm.Function): + converted_children = [None] * len(symbol.children) + for i, child in enumerate(symbol.children): + converted_children[i] = self.convert(child) + return symbol._function_evaluate(converted_children) + + elif isinstance(symbol, pybamm.Concatenation): + converted_children = [self.convert(child) for child in symbol.children] + return symbol._concatenation_evaluate(converted_children) + + 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/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/symbol.py b/pybamm/expression_tree/symbol.py index 8bbbb8db99..633875c492 100644 --- a/pybamm/expression_tree/symbol.py +++ b/pybamm/expression_tree/symbol.py @@ -584,6 +584,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, casadi_symbols=None): + """ + Convert the expression tree to a CasADi expression tree. + See :class:`pybamm.CasadiConverter`. + """ + return pybamm.CasadiConverter(casadi_symbols).convert(self) + def new_copy(self): """ Make a new copy of a symbol, to avoid Tree corruption errors while bypassing diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py new file mode 100644 index 0000000000..dc0e0acc9f --- /dev/null +++ b/pybamm/solvers/casadi_solver.py @@ -0,0 +1,9 @@ +# +# Wrap CasADi +# +import pybamm +import casadi + + +class CasadiSolver(pybamm.DaeSolver): + pass 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/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..7de38d55f1 --- /dev/null +++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py @@ -0,0 +1,62 @@ +# +# Test for the Simplify class +# +import casadi +import math +import numpy as np +import pybamm +import unittest +from tests import get_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): + pass + + +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 From e9b4b109d784ec61c7725fdcdc7e0b048583445e Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Tue, 15 Oct 2019 23:57:37 -0400 Subject: [PATCH 02/23] #664 convert more symbols to casadi --- examples/scripts/compare_lithium_ion.py | 6 +- pybamm/__init__.py | 10 +- .../operations/convert_to_casadi.py | 38 ++++-- pybamm/expression_tree/state_vector.py | 2 +- pybamm/expression_tree/symbol.py | 4 +- pybamm/geometry/geometry.py | 13 +- pybamm/models/standard_variables.py | 6 +- .../submodels/electrode/ohm/composite_ohm.py | 1 - .../base_electrolyte_conductivity.py | 3 +- .../submodels/particle/base_particle.py | 6 +- pybamm/parameters/electrical_parameters.py | 6 +- pybamm/parameters/parameter_sets.py | 1 - pybamm/processed_variable.py | 6 +- pybamm/solvers/casadi_solver.py | 127 +++++++++++++++++- .../test_operations/test_convert_to_casadi.py | 19 ++- 15 files changed, 198 insertions(+), 50 deletions(-) diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py index 1aca4c7816..7f2d3a283d 100644 --- a/examples/scripts/compare_lithium_ion.py +++ b/examples/scripts/compare_lithium_ion.py @@ -19,8 +19,8 @@ options = {"thermal": "isothermal"} models = [ pybamm.lithium_ion.SPM(options), - pybamm.lithium_ion.SPMe(options), - pybamm.lithium_ion.DFN(options), + # pybamm.lithium_ion.SPMe(options), + # pybamm.lithium_ion.DFN(options), ] @@ -47,7 +47,7 @@ solutions = [None] * len(models) t_eval = np.linspace(0, 0.17, 100) for i, model in enumerate(models): - solutions[i] = model.default_solver.solve(model, t_eval) + solutions[i] = pybamm.CasadiSolver().solve(model, t_eval) # plot plot = pybamm.QuickPlot(models, mesh, solutions) diff --git a/pybamm/__init__.py b/pybamm/__init__.py index dda95eea37..e525fdd640 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -25,7 +25,7 @@ def _load_version_int(): __version_int__ = _load_version_int() __version__ = ".".join([str(x) for x in __version_int__]) if sys.version_info[0] < 3: - del (x) # Before Python3, list comprehension iterators leaked + del x # Before Python3, list comprehension iterators leaked # # Expose PyBaMM version @@ -248,11 +248,11 @@ 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.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 # # Current profiles diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py index 41a7251c33..98b514d754 100644 --- a/pybamm/expression_tree/operations/convert_to_casadi.py +++ b/pybamm/expression_tree/operations/convert_to_casadi.py @@ -9,7 +9,7 @@ class CasadiConverter(object): def __init__(self, casadi_symbols=None): self._casadi_symbols = casadi_symbols or {} - def convert(self, symbol): + 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 @@ -30,26 +30,34 @@ def convert(self, symbol): try: return self._casadi_symbols[symbol.id] except KeyError: - casadi_symbol = self._convert(symbol) + casadi_symbol = self._convert(symbol, t, y) self._casadi_symbols[symbol.id] = casadi_symbol return casadi_symbol - def _convert(self, symbol): - """ See :meth:`Simplification.convert()`. """ - if isinstance(symbol, pybamm.Scalar): - return casadi.SX(symbol.evaluate()) + 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)) - if isinstance(symbol, pybamm.BinaryOperator): + 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) - converted_right = self.convert(right) - # _binary_evaluate defined in derived classes for specific rules - return symbol._binary_evaluate(converted_left, converted_right) + converted_left = self.convert(left, t, y) + converted_right = self.convert(right, t, y) + if isinstance(symbol, pybamm.Outer): + return casadi.outer_prod(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) + converted_child = self.convert(symbol.child, t, y) if isinstance(symbol, pybamm.AbsoluteValue): return casadi.fabs(converted_child) return symbol._unary_evaluate(converted_child) @@ -57,11 +65,13 @@ def _convert(self, symbol): elif isinstance(symbol, pybamm.Function): converted_children = [None] * len(symbol.children) for i, child in enumerate(symbol.children): - converted_children[i] = self.convert(child) + converted_children[i] = self.convert(child, t, y) return symbol._function_evaluate(converted_children) elif isinstance(symbol, pybamm.Concatenation): - converted_children = [self.convert(child) for child in symbol.children] + converted_children = [ + self.convert(child, t, y) for child in symbol.children + ] return symbol._concatenation_evaluate(converted_children) else: diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py index dbdd7c4251..c8f91b686c 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 633875c492..1351be141f 100644 --- a/pybamm/expression_tree/symbol.py +++ b/pybamm/expression_tree/symbol.py @@ -584,12 +584,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, casadi_symbols=None): + 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) + return pybamm.CasadiConverter(casadi_symbols).convert(self, t, y) def new_copy(self): """ diff --git a/pybamm/geometry/geometry.py b/pybamm/geometry/geometry.py index bb7389b8c3..4d38bd3bb9 100644 --- a/pybamm/geometry/geometry.py +++ b/pybamm/geometry/geometry.py @@ -113,8 +113,8 @@ def add_domain(self, name, geometry): for k, v in geometry.items(): if k not in ["primary", "secondary", "tabs"]: raise ValueError( - "keys of geometry must be either \"primary\", \"secondary\" or " - "\"tabs\"" + 'keys of geometry must be either "primary", "secondary" or ' + '"tabs"' ) if k != "tabs": for variable, rnge in v.items(): @@ -135,15 +135,12 @@ def add_domain(self, name, geometry): else: for region, params in v.items(): if region not in ["negative", "positive"]: - raise ValueError( - "tabs region must be \"negative\" or \"positive\"" - - ) + raise ValueError('tabs region must be "negative" or "positive"') for pname in params.keys(): if pname not in ["y_centre", "z_centre", "width"]: raise ValueError( - "tabs region params must be \"y_centre\", " - "\"z_centre\" or \"width\"" + 'tabs region params must be "y_centre", ' + '"z_centre" or "width"' ) self.update({name: geometry}) diff --git a/pybamm/models/standard_variables.py b/pybamm/models/standard_variables.py index 277aeebe73..64636416c5 100644 --- a/pybamm/models/standard_variables.py +++ b/pybamm/models/standard_variables.py @@ -135,12 +135,10 @@ auxiliary_domains={"secondary": "current collector"}, ) c_s_n_surf_xav = pybamm.Variable( - "X-averaged negative particle surface concentration", - "current collector", + "X-averaged negative particle surface concentration", "current collector" ) c_s_p_surf_xav = pybamm.Variable( - "X-averaged positive particle surface concentration", - "current collector", + "X-averaged positive particle surface concentration", "current collector" ) 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/electrolyte/base_electrolyte_conductivity.py b/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py index 47c6f69890..7d685dc7fe 100644 --- a/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py +++ b/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py @@ -61,8 +61,7 @@ def _get_standard_potential_variables(self, phi_e, phi_e_av): "Positive electrolyte potential": phi_e_p, "Positive electrolyte potential [V]": -param.U_n_ref + pot_scale * phi_e_p, "Electrolyte potential": phi_e, - "Electrolyte potential [V]": -param.U_n_ref - + pot_scale * phi_e, + "Electrolyte potential [V]": -param.U_n_ref + pot_scale * phi_e, "X-averaged electrolyte potential": phi_e_av, "X-averaged electrolyte potential [V]": -param.U_n_ref + pot_scale * phi_e_av, 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..f9260935c0 100644 --- a/pybamm/parameters/electrical_parameters.py +++ b/pybamm/parameters/electrical_parameters.py @@ -1,6 +1,7 @@ # # Standard electrical parameters # +import casadi import pybamm import numpy as np @@ -9,7 +10,10 @@ def abs_non_zero(x): if x == 0: # pragma: no cover return 1 else: - return abs(x) + if isinstance(x, casadi.SX): + return casadi.fabs(x) + else: + return abs(x) # -------------------------------------------------------------------------------------- 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/processed_variable.py b/pybamm/processed_variable.py index c7842ea10d..d2c78339c7 100644 --- a/pybamm/processed_variable.py +++ b/pybamm/processed_variable.py @@ -94,8 +94,10 @@ def __init__( self.base_eval = base_variable.evaluate(t_sol[0], u_sol[:, 0]) # handle 2D (in space) finite element variables differently - if mesh and "current collector" in self.domain and isinstance( - self.mesh[self.domain[0]][0], pybamm.Scikit2DSubMesh + if ( + mesh + and "current collector" in self.domain + and isinstance(self.mesh[self.domain[0]][0], pybamm.Scikit2DSubMesh) ): if len(self.t_sol) == 1: # space only (steady solution) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index dc0e0acc9f..228e6195d0 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -5,5 +5,128 @@ import casadi -class CasadiSolver(pybamm.DaeSolver): - pass +class CasadiSolver(pybamm.BaseSolver): + """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). + """ + + def __init__(self, method=None, rtol=1e-6, atol=1e-6): + super().__init__(method, rtol, atol) + + def compute_solution(self, model, t_eval): + """Calculate the solution of the model at specified times. + + 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 ODE solver") + solution = self.integrate( + self.dydt, + self.y0, + t_eval, + events=self.event_funs, + mass_matrix=model.mass_matrix.entries, + jacobian=self.jacobian, + ) + solve_time = timer.time() - solve_start_time + + # Identify the event that caused termination + termination = self.get_termination_reason(solution, self.events) + + return solution, solve_time, termination + + def set_up(self, model): + """Unpack model, perform checks, simplify and calculate jacobian. + + 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) + + """ + y0 = model.concatenated_initial_conditions[:, 0] + + t = casadi.SX.sym("t") + y_diff = casadi.SX.sym("y_diff", len(model.concatenated_rhs.evaluate(0, y0))) + # y_alg = casadi.SX.sym("y_alg") + # create simplified rhs and event expressions + concatenated_rhs = model.concatenated_rhs.to_casadi(t, y_diff) + events = model.events.to_casadi(t, y_diff) + + # Create function to evaluate rhs + def dydt(t, y): + pybamm.logger.debug("Evaluating RHS for {} at t={}".format(model.name, t)) + y = y[:, np.newaxis] + dy = concatenated_rhs.evaluate(t, y, known_evals={})[0] + return dy[:, 0] + + # Create event-dependent function to evaluate events + def event_fun(event): + def eval_event(t, y): + return event.evaluate(t, y) + + return eval_event + + event_funs = [event_fun(event) for event in events.values()] + + # Create function to evaluate jacobian + if jac_rhs is not None: + + def jacobian(t, y): + return jac_rhs.evaluate(t, y, known_evals={})[0] + + else: + jacobian = None + + # Add the solver attributes + self.y0 = y0 + self.dydt = dydt + self.events = events + self.event_funs = event_funs + self.jacobian = jacobian + + def integrate( + self, derivs, y0, t_eval, events=None, mass_matrix=None, jacobian=None + ): + """ + Solve a model defined by dydt with initial conditions y0. + + Parameters + ---------- + derivs : method + A function that takes in t and y and returns the time-derivative dydt + y0 : numeric type + The initial conditions + t_eval : numeric type + The times at which to compute the solution + events : method, optional + A function that takes in t and y and returns conditions for the solver to + stop + mass_matrix : array_like, optional + The (sparse) mass matrix for the chosen spatial method. + jacobian : method, optional + A function that takes in t and y and returns the Jacobian + """ + raise NotImplementedError 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 7de38d55f1..11a0b80952 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 @@ -49,7 +49,24 @@ def myfunction(x, y): self.assertEqual((b / d).to_casadi(), casadi.SX(1 / 2)) def test_convert_array_symbols(self): - pass + # 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) + ) if __name__ == "__main__": From 69aac748a9ed2fd7209aaccb64e89c0f2d74c220 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Thu, 17 Oct 2019 22:16:18 -0400 Subject: [PATCH 03/23] #664 add some more functions --- .flake8 | 3 +- examples/scripts/compare_lithium_ion.py | 4 +- .../operations/convert_to_casadi.py | 38 ++++++++++++++--- pybamm/solvers/casadi_solver.py | 10 +++-- .../test_operations/test_convert_to_casadi.py | 41 ++++++++++++++++++- 5 files changed, 83 insertions(+), 13 deletions(-) diff --git a/.flake8 b/.flake8 index 9acc35b757..f7b9fdd18b 100644 --- a/.flake8 +++ b/.flake8 @@ -10,7 +10,8 @@ exclude= lib, lib64, share, - pyvenv.cfg + pyvenv.cfg, + third-party ignore= # False positive for white space before ':' on list slice # black should format these correctly diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py index 7f2d3a283d..41ad311860 100644 --- a/examples/scripts/compare_lithium_ion.py +++ b/examples/scripts/compare_lithium_ion.py @@ -18,9 +18,9 @@ # load models options = {"thermal": "isothermal"} models = [ - pybamm.lithium_ion.SPM(options), + # pybamm.lithium_ion.SPM(options), # pybamm.lithium_ion.SPMe(options), - # pybamm.lithium_ion.DFN(options), + pybamm.lithium_ion.DFN(options) ] diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py index 98b514d754..0f1a0d02bf 100644 --- a/pybamm/expression_tree/operations/convert_to_casadi.py +++ b/pybamm/expression_tree/operations/convert_to_casadi.py @@ -3,6 +3,7 @@ # import pybamm import casadi +import numpy as np class CasadiConverter(object): @@ -51,7 +52,7 @@ def _convert(self, symbol, t, y): converted_left = self.convert(left, t, y) converted_right = self.convert(right, t, y) if isinstance(symbol, pybamm.Outer): - return casadi.outer_prod(converted_left, converted_right) + 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) @@ -63,16 +64,41 @@ def _convert(self, symbol, t, y): return symbol._unary_evaluate(converted_child) elif isinstance(symbol, pybamm.Function): - converted_children = [None] * len(symbol.children) - for i, child in enumerate(symbol.children): - converted_children[i] = self.convert(child, t, y) - return symbol._function_evaluate(converted_children) + converted_children = [ + self.convert(child, t, y) for child in symbol.children + ] + if symbol.function == np.min: + return casadi.mmin(*converted_children) + elif symbol.function == np.max: + return casadi.mmax(*converted_children) + else: + return symbol._function_evaluate(converted_children) elif isinstance(symbol, pybamm.Concatenation): converted_children = [ self.convert(child, t, y) for child in symbol.children ] - return symbol._concatenation_evaluate(converted_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( diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 228e6195d0..d8a6b02c92 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -70,10 +70,14 @@ def set_up(self, model): t = casadi.SX.sym("t") y_diff = casadi.SX.sym("y_diff", len(model.concatenated_rhs.evaluate(0, y0))) - # y_alg = casadi.SX.sym("y_alg") + y_alg = casadi.SX.sym( + "y_alg", len(model.concatenated_algebraic.evaluate(0, y0)) + ) + y = casadi.vertcat(y_diff, y_alg) # create simplified rhs and event expressions - concatenated_rhs = model.concatenated_rhs.to_casadi(t, y_diff) - events = model.events.to_casadi(t, y_diff) + concatenated_rhs = model.concatenated_rhs.to_casadi(t, y) + concatenated_algebraic = model.concatenated_algebraic.to_casadi(t, y) + events = {name: event.to_casadi(t, y) for name, event in model.events.items()} # Create function to evaluate rhs def dydt(t, y): 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 11a0b80952..0eb7e34762 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 @@ -6,7 +6,7 @@ import numpy as np import pybamm import unittest -from tests import get_discretisation_for_testing +from tests import get_mesh_for_testing, get_1p1d_discretisation_for_testing class TestCasadiConverter(unittest.TestCase): @@ -68,6 +68,45 @@ def test_convert_array_symbols(self): 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 = np.array([1, 2, 3, 4, 5]) + pybamm_a = pybamm.Array(a) + self.assertEqual(pybamm.min(pybamm_a).to_casadi(), casadi.SX(1)) + + 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)))) + if __name__ == "__main__": print("Add -v for more debug output") From 9572e8e9ff646bad4de53391d7745c8745dd694e Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Fri, 18 Oct 2019 21:08:59 -0400 Subject: [PATCH 04/23] #664 add casadi set up to solver --- examples/scripts/compare_lithium_ion.py | 6 +- pybamm/discretisations/discretisation.py | 2 +- pybamm/models/base_model.py | 18 +++-- pybamm/solvers/algebraic_solver.py | 4 +- pybamm/solvers/base_solver.py | 26 ++++-- pybamm/solvers/dae_solver.py | 4 +- pybamm/solvers/ode_solver.py | 81 ++++++++++++++++++- .../test_models/standard_model_tests.py | 3 +- 8 files changed, 120 insertions(+), 24 deletions(-) diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py index 41ad311860..7fcfae8788 100644 --- a/examples/scripts/compare_lithium_ion.py +++ b/examples/scripts/compare_lithium_ion.py @@ -18,9 +18,9 @@ # load models options = {"thermal": "isothermal"} models = [ - # pybamm.lithium_ion.SPM(options), - # pybamm.lithium_ion.SPMe(options), - pybamm.lithium_ion.DFN(options) + pybamm.lithium_ion.SPM(options), + pybamm.lithium_ion.SPMe(options), + pybamm.lithium_ion.DFN(options), ] diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py index a00e8f6aec..b96f04f7ba 100644 --- a/pybamm/discretisations/discretisation.py +++ b/pybamm/discretisations/discretisation.py @@ -127,7 +127,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/models/base_model.py b/pybamm/models/base_model.py index 0470d86773..c2e4ab60bc 100644 --- a/pybamm/models/base_model.py +++ b/pybamm/models/base_model.py @@ -59,11 +59,17 @@ 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"): @@ -86,7 +92,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/solvers/algebraic_solver.py b/pybamm/solvers/algebraic_solver.py index baea9ebb5d..f3be30d360 100644 --- a/pybamm/solvers/algebraic_solver.py +++ b/pybamm/solvers/algebraic_solver.py @@ -199,14 +199,14 @@ def set_up(self, model): pybamm.logger.info("Simplifying jacobian") jac = jac.simplify() - 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..1e1a0d0f4a 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": + self.set_up_casadi(model) + else: + self.set_up(model) set_up_time = timer.time() - start_time # Solve @@ -127,7 +130,10 @@ 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": + self.set_up_casadi(model) + else: + self.set_up(model) self.t = 0.0 set_up_time = timer.time() - start_time else: @@ -187,11 +193,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/dae_solver.py b/pybamm/solvers/dae_solver.py index 003403254f..656413c9d0 100644 --- a/pybamm/solvers/dae_solver.py +++ b/pybamm/solvers/dae_solver.py @@ -140,7 +140,7 @@ 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) @@ -152,7 +152,7 @@ def jac_alg_fn(t, y): jac = None jac_alg_fn = 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") diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py index db279b2577..cbe2514c99 100644 --- a/pybamm/solvers/ode_solver.py +++ b/pybamm/solvers/ode_solver.py @@ -99,13 +99,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") @@ -145,6 +145,83 @@ 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""" + ) + + # create simplified rhs and event expressions + concatenated_rhs = model.concatenated_rhs + events = model.events + + 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 = concatenated_rhs.to_casadi(t_casadi, y_casadi) + pybamm.logger.info("Converting events to CasADi") + events = { + name: event.to_casadi(t_casadi, y_casadi) for name, event in 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): + event_fn = casadi.Function("event", [t_casadi, y_casadi], [event]) + + def eval_event(t, y): + return event_fun(t, y) + + return eval_event + + event_funs = [event_fun(event) for event in events.values()] + + # Create function to evaluate jacobian + if model.use_jacobian: + + casadi_jac = casadi.jacobian() + casadi_jac_fn = casadi.Function("jac", [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 = 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/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 From 25ac27678c459ad3bb3ecd7e141ff798d4c89939 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 21 Oct 2019 13:50:22 -0400 Subject: [PATCH 05/23] #664 integrating casadi into solvers --- examples/scripts/compare_lead_acid.py | 7 +- examples/scripts/compare_lithium_ion.py | 7 +- pybamm/solvers/base_solver.py | 1 + pybamm/solvers/dae_solver.py | 127 ++++++++++++++++++- pybamm/solvers/ode_solver.py | 15 ++- pybamm/solvers/scipy_solver.py | 2 +- tests/unit/test_solvers/test_scipy_solver.py | 24 ++++ 7 files changed, 164 insertions(+), 19 deletions(-) diff --git a/examples/scripts/compare_lead_acid.py b/examples/scripts/compare_lead_acid.py index dc02b534dc..468d816bd0 100644 --- a/examples/scripts/compare_lead_acid.py +++ b/examples/scripts/compare_lead_acid.py @@ -18,9 +18,9 @@ # load models models = [ pybamm.lead_acid.LOQS(), - pybamm.lead_acid.FOQS(), - pybamm.lead_acid.Composite(), - pybamm.lead_acid.Full(), + # pybamm.lead_acid.FOQS(), + # pybamm.lead_acid.Composite(), + # pybamm.lead_acid.Full(), ] # load parameter values and process models and geometry @@ -43,6 +43,7 @@ solutions = [None] * len(models) t_eval = np.linspace(0, 3, 1000) for i, model in enumerate(models): + model.convert_to_format = "casadi" solution = model.default_solver.solve(model, t_eval) solutions[i] = solution diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py index 7fcfae8788..42056e5c5b 100644 --- a/examples/scripts/compare_lithium_ion.py +++ b/examples/scripts/compare_lithium_ion.py @@ -20,7 +20,7 @@ models = [ pybamm.lithium_ion.SPM(options), pybamm.lithium_ion.SPMe(options), - pybamm.lithium_ion.DFN(options), + # pybamm.lithium_ion.DFN(options), ] @@ -32,7 +32,7 @@ # set mesh var = pybamm.standard_spatial_vars -var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.r_n: 5, var.r_p: 5} +var_pts = {var.x_n: 100, var.x_s: 100, var.x_p: 100, var.r_n: 50, var.r_p: 50} # discretise models for model in models: @@ -47,7 +47,8 @@ solutions = [None] * len(models) t_eval = np.linspace(0, 0.17, 100) for i, model in enumerate(models): - solutions[i] = pybamm.CasadiSolver().solve(model, t_eval) + # model.convert_to_format = "casadi" + solutions[i] = model.default_solver.solve(model, t_eval) # plot plot = pybamm.QuickPlot(models, mesh, solutions) diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index 1e1a0d0f4a..6674544bd8 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -1,6 +1,7 @@ # # Base solver class # +import casadi import pybamm import numpy as np diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py index 656413c9d0..5d19357b6e 100644 --- a/pybamm/solvers/dae_solver.py +++ b/pybamm/solvers/dae_solver.py @@ -102,12 +102,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 @@ -220,6 +214,127 @@ 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 + """ + # create simplified rhs, algebraic and event expressions + concatenated_rhs = model.concatenated_rhs + concatenated_algebraic = model.concatenated_algebraic + events = model.events + + t_casadi = casadi.SX.sym("t") + y_casadi = casadi.SX.sym("y", len(y0)) + pybamm.logger.info("Converting RHS to CasADi") + concatenated_rhs = concatenated_rhs.to_casadi(t_casadi, y_casadi) + concatenated_algebraic = concatenated_algebraic.to_casadi(t_casadi, y_casadi) + if model.use_simplify: + simp = pybamm.Simplification() + pybamm.logger.info("Simplifying events") + events = {name: simp.simplify(event) for name, event in events.items()} + # pybamm.logger.info("Converting events to python") + # events = {name: pybamm.EvaluatorPython(event) for name, event in 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] + ) + + 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] + + if model.use_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 + jacobian_alg = None + + 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") + concatenated_algebraic = pybamm.EvaluatorPython(concatenated_algebraic) + pybamm.logger.info("Converting events to python") + events = { + name: pybamm.EvaluatorPython(event) for name, event in events.items() + } + + # Calculate consistent initial conditions for the algebraic equations + def rhs(t, y): + return concatenated_rhs.evaluate(t, y, known_evals={})[0][:, 0] + + def algebraic(t, y): + return concatenated_algebraic.evaluate(t, y, known_evals={})[0][:, 0] + + if len(model.algebraic) > 0: + y0 = self.calculate_consistent_initial_conditions( + rhs, algebraic, model.concatenated_initial_conditions[:, 0], jac_alg_fn + ) + 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) + ) + y = y[:, np.newaxis] + rhs_eval, known_evals = concatenated_rhs.evaluate(t, y, known_evals={}) + # reuse known_evals + alg_eval = concatenated_algebraic.evaluate(t, y, known_evals=known_evals)[0] + # turn into 1D arrays + rhs_eval = rhs_eval[:, 0] + alg_eval = alg_eval[:, 0] + return ( + np.concatenate((rhs_eval, alg_eval)) - model.mass_matrix.entries @ ydot + ) + + # Create event-dependent function to evaluate events + def event_fun(event): + def eval_event(t, y): + return event.evaluate(t, y) + + return eval_event + + event_funs = [event_fun(event) for event in events.values()] + + # Create function to evaluate jacobian + if jac is not None: + + def jacobian(t, y): + return jac.evaluate(t, y, known_evals={})[0] + + else: + jacobian = 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 calculate_consistent_initial_conditions( self, rhs, algebraic, y0_guess, jac=None ): diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py index cbe2514c99..455f06ed5b 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 @@ -178,7 +179,7 @@ def set_up_casadi(self, model): pybamm.logger.info("Converting RHS to CasADi") concatenated_rhs = concatenated_rhs.to_casadi(t_casadi, y_casadi) pybamm.logger.info("Converting events to CasADi") - events = { + casadi_events = { name: event.to_casadi(t_casadi, y_casadi) for name, event in events.items() } @@ -194,20 +195,22 @@ def dydt(t, y): # Create event-dependent function to evaluate events def event_fun(event): - event_fn = casadi.Function("event", [t_casadi, y_casadi], [event]) + casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi], [event]) def eval_event(t, y): - return event_fun(t, y) + return casadi_event_fn(t, y) return eval_event - event_funs = [event_fun(event) for event in events.values()] + event_funs = [event_fun(event) for event in casadi_events.values()] # Create function to evaluate jacobian if model.use_jacobian: - casadi_jac = casadi.jacobian() - casadi_jac_fn = casadi.Function("jac", [t_casadi, y_casadi], [casadi_jac]) + 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) 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/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py index cef29c65e0..58ad303efb 100644 --- a/tests/unit/test_solvers/test_scipy_solver.py +++ b/tests/unit/test_solvers/test_scipy_solver.py @@ -270,6 +270,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_event_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: -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") From 5a1922e55e5f25896018cb583876039cb61b2882 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 21 Oct 2019 21:57:42 -0400 Subject: [PATCH 06/23] #664 set up dae solver --- examples/scripts/compare_lithium_ion.py | 6 +- pybamm/solvers/dae_solver.py | 117 +++++++----------- pybamm/solvers/ode_solver.py | 11 +- .../unit/test_solvers/test_scikits_solvers.py | 49 ++++++++ 4 files changed, 104 insertions(+), 79 deletions(-) diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py index 42056e5c5b..45c08bd699 100644 --- a/examples/scripts/compare_lithium_ion.py +++ b/examples/scripts/compare_lithium_ion.py @@ -20,7 +20,7 @@ models = [ pybamm.lithium_ion.SPM(options), pybamm.lithium_ion.SPMe(options), - # pybamm.lithium_ion.DFN(options), + pybamm.lithium_ion.DFN(options), ] @@ -32,7 +32,7 @@ # set mesh var = pybamm.standard_spatial_vars -var_pts = {var.x_n: 100, var.x_s: 100, var.x_p: 100, var.r_n: 50, var.r_p: 50} +var_pts = {var.x_n: 20, var.x_s: 20, var.x_p: 20, var.r_n: 10, var.r_p: 10} # discretise models for model in models: @@ -47,7 +47,7 @@ solutions = [None] * len(models) t_eval = np.linspace(0, 0.17, 100) for i, model in enumerate(models): - # model.convert_to_format = "casadi" + model.convert_to_format = "casadi" solutions[i] = model.default_solver.solve(model, t_eval) # plot diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py index 5d19357b6e..846abf64ae 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 @@ -139,12 +140,15 @@ def set_up(self, model): 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.convert_to_format == "python": pybamm.logger.info("Converting RHS to python") @@ -165,7 +169,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 @@ -196,15 +203,6 @@ 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: - - def jacobian(t, y): - return jac.evaluate(t, y, known_evals={})[0] - - else: - jacobian = None - # Add the solver attributes self.y0 = y0 self.rhs = rhs @@ -223,22 +221,21 @@ def set_up_casadi(self, model): The model whose solution to calculate. Must have attributes rhs and initial_conditions """ - # create simplified rhs, algebraic and event expressions - concatenated_rhs = model.concatenated_rhs - concatenated_algebraic = model.concatenated_algebraic - events = model.events - + # Convert model attributes to casadi t_casadi = casadi.SX.sym("t") - y_casadi = casadi.SX.sym("y", len(y0)) + y_casadi = casadi.SX.sym("y", len(model.concatenated_initial_conditions)) pybamm.logger.info("Converting RHS to CasADi") - concatenated_rhs = concatenated_rhs.to_casadi(t_casadi, y_casadi) - concatenated_algebraic = concatenated_algebraic.to_casadi(t_casadi, y_casadi) - if model.use_simplify: - simp = pybamm.Simplification() - pybamm.logger.info("Simplifying events") - events = {name: simp.simplify(event) for name, event in events.items()} - # pybamm.logger.info("Converting events to python") - # events = {name: pybamm.EvaluatorPython(event) for name, event in events.items()} + 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( @@ -247,91 +244,73 @@ def set_up_casadi(self, model): concatenated_algebraic_fn = casadi.Function( "algebraic", [t_casadi, y_casadi], [concatenated_algebraic] ) - - 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] + all_states_fn = casadi.Function("all", [t_casadi, y_casadi], [all_states]) if model.use_jacobian: - casadi_jac = casadi.jacobian(concatenated_rhs, y_casadi) + 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 casadi_jac_fn(t, y) + def jacobian_alg(t, y): + return casadi_jac_alg_fn(t, y) + else: jacobian = None jacobian_alg = None - 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") - concatenated_algebraic = pybamm.EvaluatorPython(concatenated_algebraic) - pybamm.logger.info("Converting events to python") - events = { - name: pybamm.EvaluatorPython(event) for name, event in events.items() - } - # Calculate consistent initial conditions for the algebraic equations def rhs(t, y): - return concatenated_rhs.evaluate(t, y, known_evals={})[0][:, 0] + return concatenated_rhs_fn(t, y).full()[:, 0] def algebraic(t, y): - return concatenated_algebraic.evaluate(t, y, known_evals={})[0][:, 0] + 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], jac_alg_fn + 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) ) - y = y[:, np.newaxis] - rhs_eval, known_evals = concatenated_rhs.evaluate(t, y, known_evals={}) - # reuse known_evals - alg_eval = concatenated_algebraic.evaluate(t, y, known_evals=known_evals)[0] - # turn into 1D arrays - rhs_eval = rhs_eval[:, 0] - alg_eval = alg_eval[:, 0] - return ( - np.concatenate((rhs_eval, alg_eval)) - model.mass_matrix.entries @ ydot - ) + 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 event.evaluate(t, y) + return casadi_event_fn(t, y) return eval_event - event_funs = [event_fun(event) for event in events.values()] - - # Create function to evaluate jacobian - if jac is not None: - - def jacobian(t, y): - return jac.evaluate(t, y, known_evals={})[0] - - else: - jacobian = None + event_funs = [event_fun(event) for event in casadi_events.values()] # Add the solver attributes self.y0 = y0 self.rhs = rhs self.algebraic = algebraic self.residuals = residuals - self.events = events + self.events = model.events self.event_funs = event_funs self.jacobian = jacobian diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py index 455f06ed5b..386c6b4c7a 100644 --- a/pybamm/solvers/ode_solver.py +++ b/pybamm/solvers/ode_solver.py @@ -168,19 +168,16 @@ def set_up_casadi(self, model): """Cannot use ODE solver to solve model with DAEs""" ) - # create simplified rhs and event expressions - concatenated_rhs = model.concatenated_rhs - events = model.events - 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 = concatenated_rhs.to_casadi(t_casadi, y_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 events.items() + name: event.to_casadi(t_casadi, y_casadi) + for name, event in model.events.items() } # Create function to evaluate rhs @@ -221,7 +218,7 @@ def jacobian(t, y): # Add the solver attributes self.y0 = y0 self.dydt = dydt - self.events = events + self.events = model.events self.event_funs = event_funs self.jacobian = jacobian diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py index 36b3cb21f8..ce4d384acd 100644 --- a/tests/unit/test_solvers/test_scikits_solvers.py +++ b/tests/unit/test_solvers/test_scikits_solvers.py @@ -680,6 +680,55 @@ 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() + 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") From c4bda5645db0c04443a498de7490fced6629b1c1 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 21 Oct 2019 23:19:01 -0400 Subject: [PATCH 07/23] #664 set up casadi solver for odes and add tests --- .../full_battery_models/lithium_ion/spm.py | 3 +- .../full_battery_models/lithium_ion/spme.py | 3 +- pybamm/solvers/casadi_solver.py | 151 ++++++++---------- pybamm/solvers/dae_solver.py | 3 + tests/unit/test_solvers/test_casadi_solver.py | 129 +++++++++++++++ 5 files changed, 199 insertions(+), 90 deletions(-) create mode 100644 tests/unit/test_solvers/test_casadi_solver.py diff --git a/pybamm/models/full_battery_models/lithium_ion/spm.py b/pybamm/models/full_battery_models/lithium_ion/spm.py index c8a9549127..8d4ce0302d 100644 --- a/pybamm/models/full_battery_models/lithium_ion/spm.py +++ b/pybamm/models/full_battery_models/lithium_ion/spm.py @@ -120,6 +120,7 @@ def default_solver(self): # Different solver depending on whether we solve ODEs or DAEs dimensionality = self.options["dimensionality"] if dimensionality == 0: - return pybamm.ScipySolver() + return pybamm.CasadiSolver() + # return pybamm.ScipySolver() else: return pybamm.ScikitsDaeSolver() diff --git a/pybamm/models/full_battery_models/lithium_ion/spme.py b/pybamm/models/full_battery_models/lithium_ion/spme.py index e5e6167592..7812f87856 100644 --- a/pybamm/models/full_battery_models/lithium_ion/spme.py +++ b/pybamm/models/full_battery_models/lithium_ion/spme.py @@ -124,6 +124,7 @@ def default_solver(self): # Different solver depending on whether we solve ODEs or DAEs dimensionality = self.options["dimensionality"] if dimensionality == 0: - return pybamm.ScipySolver() + return pybamm.CasadiSolver() + # return pybamm.ScipySolver() else: return pybamm.ScikitsDaeSolver() diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index d8a6b02c92..d4b828ead5 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -1,11 +1,12 @@ # -# Wrap CasADi +# CasADi Solver class # -import pybamm import casadi +import pybamm +import numpy as np -class CasadiSolver(pybamm.BaseSolver): +class CasadiSolver(pybamm.DaeSolver): """Solve a discretised model, using CasADi. Parameters @@ -14,13 +15,25 @@ class CasadiSolver(pybamm.BaseSolver): 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=None, rtol=1e-6, atol=1e-6): - super().__init__(method, rtol, atol) + def __init__( + self, + method="idas", + rtol=1e-6, + atol=1e-6, + root_method="lm", + root_tol=1e-6, + max_steps=1000, + ): + super().__init__(method, rtol, atol, root_method, root_tol, max_steps) def compute_solution(self, model, t_eval): - """Calculate the solution of the model at specified times. + """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 ---------- @@ -31,17 +44,19 @@ def compute_solution(self, model, t_eval): The times at which to compute the solution """ + # Convert to CasADi if not already done + if not hasattr(self, "casadi_problem"): + pybamm.logger.info( + "Converting model to CasADi format, required for CasADi solver" + ) + self.set_up_casadi(model) + timer = pybamm.Timer() solve_start_time = timer.time() - pybamm.logger.info("Calling ODE solver") - solution = self.integrate( - self.dydt, - self.y0, - t_eval, - events=self.event_funs, - mass_matrix=model.mass_matrix.entries, - jacobian=self.jacobian, + 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 @@ -50,87 +65,47 @@ def compute_solution(self, model, t_eval): return solution, solve_time, termination - def set_up(self, model): - """Unpack model, perform checks, simplify and calculate jacobian. - - 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) - - """ - y0 = model.concatenated_initial_conditions[:, 0] - - t = casadi.SX.sym("t") - 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.vertcat(y_diff, y_alg) - # create simplified rhs and event expressions - concatenated_rhs = model.concatenated_rhs.to_casadi(t, y) - concatenated_algebraic = model.concatenated_algebraic.to_casadi(t, y) - events = {name: event.to_casadi(t, y) for name, event in model.events.items()} - - # Create function to evaluate rhs - def dydt(t, y): - pybamm.logger.debug("Evaluating RHS for {} at t={}".format(model.name, t)) - y = y[:, np.newaxis] - dy = concatenated_rhs.evaluate(t, y, known_evals={})[0] - return dy[:, 0] - - # Create event-dependent function to evaluate events - def event_fun(event): - def eval_event(t, y): - return event.evaluate(t, y) - - return eval_event - - event_funs = [event_fun(event) for event in events.values()] - - # Create function to evaluate jacobian - if jac_rhs is not None: - - def jacobian(t, y): - return jac_rhs.evaluate(t, y, known_evals={})[0] - - else: - jacobian = None - - # Add the solver attributes - self.y0 = y0 - self.dydt = dydt - self.events = events - self.event_funs = event_funs - self.jacobian = jacobian - - def integrate( - self, derivs, y0, t_eval, events=None, mass_matrix=None, jacobian=None - ): + def integrate_casadi(self, problem, y0, t_eval, mass_matrix=None): """ - Solve a model defined by dydt with initial conditions y0. + Solve a DAE model defined by residuals with initial conditions y0. Parameters ---------- - derivs : method - A function that takes in t and y and returns the time-derivative dydt + 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 - events : method, optional - A function that takes in t and y and returns conditions for the solver to - stop mass_matrix : array_like, optional - The (sparse) mass matrix for the chosen spatial method. - jacobian : method, optional - A function that takes in t and y and returns the Jacobian + 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. """ - raise NotImplementedError + options = { + "grid": t_eval, + "reltol": self.rtol, + "abstol": self.atol, + "output_t0": True, + } + if self.method == "idas": + options["calc_ic"] = True + + # set up and solve + integrator = casadi.integrator("F", self.method, problem, options) + try: + # return solution + sol = integrator(x0=y0) + 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: + raise pybamm.SolverError(e.args[0]) + + def set_up(self, model): + "Skip classic set up with this solver" + pass + + # def calculate_consistent_initial_conditions(self): + # "No need to calculate initial conditions separately with this solver" + # pass diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py index 846abf64ae..efd5431df2 100644 --- a/pybamm/solvers/dae_solver.py +++ b/pybamm/solvers/dae_solver.py @@ -314,6 +314,9 @@ def eval_event(t, y): self.event_funs = event_funs self.jacobian = jacobian + # Create CasADi problem for the CasADi solver + self.casadi_problem = {"x": y_casadi, "ode": concatenated_rhs} + def calculate_consistent_initial_conditions( self, rhs, algebraic, y0_guess, jac=None ): 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..cfed566600 --- /dev/null +++ b/tests/unit/test_solvers/test_casadi_solver.py @@ -0,0 +1,129 @@ +# +# 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") + 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() From cbf97b7e6baec7d782f00cee0cee3065b5138a01 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Tue, 22 Oct 2019 23:06:26 -0400 Subject: [PATCH 08/23] #664 fix bugs with function diff and converting autograd to casadi --- CHANGELOG.md | 17 +++-- examples/scripts/compare_lead_acid.py | 10 +-- examples/scripts/compare_lithium_ion.py | 2 +- examples/scripts/compare_lithium_ion_3D.py | 7 +- pybamm/expression_tree/functions.py | 65 ++++++++++++----- .../operations/convert_to_casadi.py | 18 ++++- .../full_battery_models/lithium_ion/spm.py | 3 +- .../full_battery_models/lithium_ion/spme.py | 3 +- pybamm/solvers/casadi_solver.py | 20 ++++-- pybamm/solvers/dae_solver.py | 15 +++- pybamm/solvers/ode_solver.py | 1 + .../test_expression_tree/test_functions.py | 22 +++++- .../test_operations/test_convert_to_casadi.py | 24 ++++++- tests/unit/test_solvers/test_casadi_solver.py | 72 +++++++++---------- tests/unit/test_solvers/test_scipy_solver.py | 3 + 15 files changed, 196 insertions(+), 86 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index eb312ad6ee..352a1a4a38 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,18 +2,23 @@ ## Features -- Add `Interpolant` class to interpolate experimental data (e.g. OCP curves) (#661) -- Allow parameters to be set by material or by specifying a particular paper (#647) -- Set relative and absolute tolerances independently in solvers (#645) -- Add some non-uniform meshes in 1D and 2D (#617) +- Add interface to CasADi solver +- Add option to use CasADi's Algorithmic Differentiation framework to calculate Jacobians +- Add `Interpolant` class to interpolate experimental data (e.g. OCP curves) ([#661](https://github.com/pybamm-team/PyBaMM/pull/661)) +- Allow parameters to be set by material or by specifying a particular paper ([#647](https://github.com/pybamm-team/PyBaMM/pull/647)) +- Set relative and absolute tolerances independently in solvers ([#645](https://github.com/pybamm-team/PyBaMM/pull/645)) +- Add some non-uniform meshes in 1D and 2D ([#617](https://github.com/pybamm-team/PyBaMM/pull/617)) ## Optimizations -- Avoid re-checking size when making a copy of an `Index` object (#656) -- Avoid recalculating `_evaluation_array` when making a copy of a `StateVector` object (#653) +- Avoid re-checking size when making a copy of an `Index` object ([#656](https://github.com/pybamm-team/PyBaMM/pull/656)) +- Avoid recalculating `_evaluation_array` when making a copy of a `StateVector` object ([#653](https://github.com/pybamm-team/PyBaMM/pull/653)) ## Bug fixes +- Fix differentiation of functions that have more than one argument + + # [v0.1.0](https://github.com/pybamm-team/PyBaMM/tree/v0.1.0) - 2019-10-08 This is the first official version of PyBaMM. diff --git a/examples/scripts/compare_lead_acid.py b/examples/scripts/compare_lead_acid.py index 468d816bd0..a10669af3f 100644 --- a/examples/scripts/compare_lead_acid.py +++ b/examples/scripts/compare_lead_acid.py @@ -18,9 +18,9 @@ # load models models = [ pybamm.lead_acid.LOQS(), - # pybamm.lead_acid.FOQS(), - # pybamm.lead_acid.Composite(), - # pybamm.lead_acid.Full(), + pybamm.lead_acid.FOQS(), + pybamm.lead_acid.Composite(), + pybamm.lead_acid.Full(), ] # load parameter values and process models and geometry @@ -41,10 +41,10 @@ # solve model solutions = [None] * len(models) -t_eval = np.linspace(0, 3, 1000) +t_eval = np.linspace(0, 0.66, 1000) for i, model in enumerate(models): model.convert_to_format = "casadi" - solution = model.default_solver.solve(model, t_eval) + solution = pybamm.CasadiSolver().solve(model, t_eval) solutions[i] = solution # plot diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py index 45c08bd699..fc06298566 100644 --- a/examples/scripts/compare_lithium_ion.py +++ b/examples/scripts/compare_lithium_ion.py @@ -32,7 +32,7 @@ # set mesh var = pybamm.standard_spatial_vars -var_pts = {var.x_n: 20, var.x_s: 20, var.x_p: 20, var.r_n: 10, var.r_p: 10} +var_pts = {var.x_n: 50, var.x_s: 50, var.x_p: 50, var.r_n: 50, var.r_p: 50} # discretise models for model in models: diff --git a/examples/scripts/compare_lithium_ion_3D.py b/examples/scripts/compare_lithium_ion_3D.py index e62fc7697b..9449f2e521 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" ), ] @@ -54,6 +52,7 @@ solutions = [None] * len(models) t_eval = np.linspace(0, 1, 1000) for i, model in enumerate(models): + model.convert_to_format = "casadi" solution = model.default_solver.solve(model, t_eval) solutions[i] = solution diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py index 20f6ae46e4..3bba9f8f6d 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 @@ -38,6 +47,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) @@ -79,7 +89,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._diff(children, i) * child.diff( + variable + ) # remove None entries partial_derivatives = list(filter(None, partial_derivatives)) @@ -90,15 +102,31 @@ def diff(self, variable): return derivative - def _diff(self, children): + def _diff(self, children, 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 _jac(self, variable): """ See :meth:`pybamm.Symbol._jac()`. """ @@ -171,7 +199,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): @@ -199,7 +231,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 ) @@ -235,7 +268,7 @@ class Cos(SpecificFunction): def __init__(self, child): super().__init__(np.cos, child) - def _diff(self, children): + def _diff(self, children, idx): """ See :meth:`pybamm.Symbol._diff()`. """ return -Sin(children[0]) @@ -251,7 +284,7 @@ class Cosh(SpecificFunction): def __init__(self, child): super().__init__(np.cosh, child) - def _diff(self, children): + def _diff(self, children, idx): """ See :meth:`pybamm.Function._diff()`. """ return Sinh(children[0]) @@ -267,7 +300,7 @@ class Exponential(SpecificFunction): def __init__(self, child): super().__init__(np.exp, child) - def _diff(self, children): + def _diff(self, children, idx): """ See :meth:`pybamm.Function._diff()`. """ return Exponential(children[0]) @@ -283,7 +316,7 @@ class Log(SpecificFunction): def __init__(self, child): super().__init__(np.log, child) - def _diff(self, children): + def _diff(self, children, idx): """ See :meth:`pybamm.Function._diff()`. """ return 1 / children[0] @@ -309,7 +342,7 @@ class Sin(SpecificFunction): def __init__(self, child): super().__init__(np.sin, child) - def _diff(self, children): + def _diff(self, children, idx): """ See :meth:`pybamm.Function._diff()`. """ return Cos(children[0]) @@ -325,7 +358,7 @@ class Sinh(SpecificFunction): def __init__(self, child): super().__init__(np.sinh, child) - def _diff(self, children): + def _diff(self, children, idx): """ See :meth:`pybamm.Function._diff()`. """ return Cosh(children[0]) diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py index 0f1a0d02bf..395bb6ec78 100644 --- a/pybamm/expression_tree/operations/convert_to_casadi.py +++ b/pybamm/expression_tree/operations/convert_to_casadi.py @@ -67,13 +67,29 @@ def _convert(self, symbol, t, y): 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 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 diff --git a/pybamm/models/full_battery_models/lithium_ion/spm.py b/pybamm/models/full_battery_models/lithium_ion/spm.py index 8d4ce0302d..c8a9549127 100644 --- a/pybamm/models/full_battery_models/lithium_ion/spm.py +++ b/pybamm/models/full_battery_models/lithium_ion/spm.py @@ -120,7 +120,6 @@ def default_solver(self): # Different solver depending on whether we solve ODEs or DAEs dimensionality = self.options["dimensionality"] if dimensionality == 0: - return pybamm.CasadiSolver() - # return pybamm.ScipySolver() + return pybamm.ScipySolver() else: return pybamm.ScikitsDaeSolver() diff --git a/pybamm/models/full_battery_models/lithium_ion/spme.py b/pybamm/models/full_battery_models/lithium_ion/spme.py index 7812f87856..e5e6167592 100644 --- a/pybamm/models/full_battery_models/lithium_ion/spme.py +++ b/pybamm/models/full_battery_models/lithium_ion/spme.py @@ -124,7 +124,6 @@ def default_solver(self): # Different solver depending on whether we solve ODEs or DAEs dimensionality = self.options["dimensionality"] if dimensionality == 0: - return pybamm.CasadiSolver() - # return pybamm.ScipySolver() + return pybamm.ScipySolver() else: return pybamm.ScikitsDaeSolver() diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index d4b828ead5..2840474775 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -88,6 +88,7 @@ def integrate_casadi(self, problem, y0, t_eval, mass_matrix=None): "reltol": self.rtol, "abstol": self.atol, "output_t0": True, + "disable_internal_warnings": True, } if self.method == "idas": options["calc_ic"] = True @@ -95,17 +96,22 @@ def integrate_casadi(self, problem, y0, t_eval, mass_matrix=None): # set up and solve integrator = casadi.integrator("F", self.method, problem, options) try: - # return solution - sol = integrator(x0=y0) + # 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]) def set_up(self, model): - "Skip classic set up with this solver" - pass + "Skip classic set up with this solver, just reset initial conditions" + self.y0 = model.concatenated_initial_conditions - # def calculate_consistent_initial_conditions(self): - # "No need to calculate initial conditions separately with this solver" - # pass + def calculate_consistent_initial_conditions( + self, rhs, algebraic, y0_guess, jac=None + ): + "No need to calculate initial conditions separately with this solver" + return y0_guess diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py index efd5431df2..5b60d048f3 100644 --- a/pybamm/solvers/dae_solver.py +++ b/pybamm/solvers/dae_solver.py @@ -223,7 +223,12 @@ def set_up_casadi(self, model): """ # Convert model attributes to casadi t_casadi = casadi.SX.sym("t") - y_casadi = casadi.SX.sym("y", len(model.concatenated_initial_conditions)) + 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") @@ -248,6 +253,7 @@ def set_up_casadi(self, model): 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] @@ -315,7 +321,12 @@ def eval_event(t, y): self.jacobian = jacobian # Create CasADi problem for the CasADi solver - self.casadi_problem = {"x": y_casadi, "ode": concatenated_rhs} + 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/ode_solver.py b/pybamm/solvers/ode_solver.py index 386c6b4c7a..ecd9965f87 100644 --- a/pybamm/solvers/ode_solver.py +++ b/pybamm/solvers/ode_solver.py @@ -204,6 +204,7 @@ def eval_event(t, y): # 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] 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/test_convert_to_casadi.py b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py index 0eb7e34762..d21966eb94 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 @@ -4,6 +4,7 @@ import casadi import math 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 @@ -29,13 +30,13 @@ def sin(x): return np.sin(x) f = pybamm.Function(sin, b) - self.assertEqual((f).to_casadi(), casadi.SX(np.sin(1))) + 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)) + self.assertEqual(f.to_casadi(), casadi.SX(3)) # addition self.assertEqual((a + b).to_casadi(), casadi.SX(1)) @@ -107,6 +108,25 @@ def test_concatenations(self): 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)) + if __name__ == "__main__": print("Add -v for more debug output") diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py index cfed566600..b2c75802ed 100644 --- a/tests/unit/test_solvers/test_casadi_solver.py +++ b/tests/unit/test_solvers/test_casadi_solver.py @@ -80,43 +80,41 @@ def test_model_solver(self): 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]) - # + 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__": diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py index 58ad303efb..e0bed5251f 100644 --- a/tests/unit/test_solvers/test_scipy_solver.py +++ b/tests/unit/test_solvers/test_scipy_solver.py @@ -269,6 +269,9 @@ def test_model_step(self): t_eval = step_sol.t solution = solver.solve(model, t_eval) np.testing.assert_allclose(solution.y[0], step_sol.y[0]) + import ipdb + + ipdb.set_trace() def test_model_solver_with_event_with_casadi(self): # Create model From df58218e5a6615bdfc9034e632f8d2ab67027715 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Tue, 22 Oct 2019 23:15:00 -0400 Subject: [PATCH 09/23] #664 testing lead acid example --- examples/scripts/compare_lead_acid.py | 13 +++++++------ pybamm/expression_tree/functions.py | 4 ++-- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/examples/scripts/compare_lead_acid.py b/examples/scripts/compare_lead_acid.py index a10669af3f..028a2a1629 100644 --- a/examples/scripts/compare_lead_acid.py +++ b/examples/scripts/compare_lead_acid.py @@ -41,17 +41,18 @@ # solve model solutions = [None] * len(models) -t_eval = np.linspace(0, 0.66, 1000) +t_eval = np.linspace(0, 1, 1000) for i, model in enumerate(models): - model.convert_to_format = "casadi" - solution = pybamm.CasadiSolver().solve(model, t_eval) + 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/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py index 3bba9f8f6d..650d4a6767 100644 --- a/pybamm/expression_tree/functions.py +++ b/pybamm/expression_tree/functions.py @@ -139,9 +139,9 @@ def _jac(self, variable): # calculate the required partial jacobians and add them jacobian = None children = self.orphans - for child in children: + for i, child in enumerate(children): if not child.evaluates_to_number(): - jac_fun = self._diff(children) * child.jac(variable) + jac_fun = self._diff(children, i) * child.jac(variable) jac_fun.domain = self.domain if jacobian is None: From 51b577828c0775b1457ba6ab1464526350afd755 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Tue, 22 Oct 2019 23:28:28 -0400 Subject: [PATCH 10/23] #664 testing 3d lead-acid --- examples/scripts/compare_lead_acid_3D.py | 3 +-- pybamm/solvers/base_solver.py | 6 ++++-- pybamm/solvers/casadi_solver.py | 11 ----------- 3 files changed, 5 insertions(+), 15 deletions(-) diff --git a/examples/scripts/compare_lead_acid_3D.py b/examples/scripts/compare_lead_acid_3D.py index 5fa825328d..8d441092ac 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/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index 6674544bd8..e0a4758ea9 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -74,7 +74,7 @@ def solve(self, model, t_eval): # Set up timer = pybamm.Timer() start_time = timer.time() - if model.convert_to_format == "casadi": + if model.convert_to_format == "casadi" or isinstance(self, pybamm.CasadiSolver): self.set_up_casadi(model) else: self.set_up(model) @@ -131,7 +131,9 @@ def step(self, model, dt, npts=2): # Run set up on first step if not hasattr(self, "y0"): start_time = timer.time() - if model.convert_to_format == "casadi": + if model.convert_to_format == "casadi" or isinstance( + self, pybamm.CasadiSolver + ): self.set_up_casadi(model) else: self.set_up(model) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 2840474775..42080444fd 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -44,13 +44,6 @@ def compute_solution(self, model, t_eval): The times at which to compute the solution """ - # Convert to CasADi if not already done - if not hasattr(self, "casadi_problem"): - pybamm.logger.info( - "Converting model to CasADi format, required for CasADi solver" - ) - self.set_up_casadi(model) - timer = pybamm.Timer() solve_start_time = timer.time() @@ -106,10 +99,6 @@ def integrate_casadi(self, problem, y0, t_eval, mass_matrix=None): # If it doesn't work raise error raise pybamm.SolverError(e.args[0]) - def set_up(self, model): - "Skip classic set up with this solver, just reset initial conditions" - self.y0 = model.concatenated_initial_conditions - def calculate_consistent_initial_conditions( self, rhs, algebraic, y0_guess, jac=None ): From 9df5a90e36fcbdceeec93e98cfddc453f4c91da1 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Wed, 23 Oct 2019 08:39:28 -0400 Subject: [PATCH 11/23] #664 fix tests and style --- pybamm/solvers/base_solver.py | 1 - .../test_operations/test_convert_to_casadi.py | 1 - tests/unit/test_solvers/test_scipy_solver.py | 3 --- 3 files changed, 5 deletions(-) diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index e0a4758ea9..3b95b10f73 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -1,7 +1,6 @@ # # Base solver class # -import casadi import pybamm import numpy as np 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 d21966eb94..01de9b159d 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 @@ -2,7 +2,6 @@ # Test for the Simplify class # import casadi -import math import numpy as np import autograd.numpy as anp import pybamm diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py index e0bed5251f..58ad303efb 100644 --- a/tests/unit/test_solvers/test_scipy_solver.py +++ b/tests/unit/test_solvers/test_scipy_solver.py @@ -269,9 +269,6 @@ def test_model_step(self): t_eval = step_sol.t solution = solver.solve(model, t_eval) np.testing.assert_allclose(solution.y[0], step_sol.y[0]) - import ipdb - - ipdb.set_trace() def test_model_solver_with_event_with_casadi(self): # Create model From 28eed5537877a2e79b93e3d165264f2f3389340a Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Thu, 24 Oct 2019 17:22:21 -0400 Subject: [PATCH 12/23] #664 pass extra options to casadi --- pybamm/solvers/casadi_solver.py | 4 +++- tests/unit/test_solvers/test_casadi_solver.py | 8 +++++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 42080444fd..5f373246f5 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -27,8 +27,10 @@ def __init__( 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 @@ -81,8 +83,8 @@ def integrate_casadi(self, problem, y0, t_eval, mass_matrix=None): "reltol": self.rtol, "abstol": self.atol, "output_t0": True, - "disable_internal_warnings": True, } + options.update(self.extra_options) if self.method == "idas": options["calc_ic"] = True diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py index b2c75802ed..6d6b45eecd 100644 --- a/tests/unit/test_solvers/test_casadi_solver.py +++ b/tests/unit/test_solvers/test_casadi_solver.py @@ -45,7 +45,13 @@ def test_integrate_failure(self): y0 = np.array([1]) t_eval = np.linspace(0, 3, 100) - solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8, method="idas") + 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): From 0a7d5f781c65dae23999d34977f18a07d7476715 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Thu, 24 Oct 2019 22:52:19 -0400 Subject: [PATCH 13/23] #664 fix docs --- docs/source/expression_tree/index.rst | 6 ------ pybamm/models/base_model.py | 5 +++-- 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/docs/source/expression_tree/index.rst b/docs/source/expression_tree/index.rst index 7049d4a275..3eedd7b622 100644 --- a/docs/source/expression_tree/index.rst +++ b/docs/source/expression_tree/index.rst @@ -17,10 +17,4 @@ Expression Tree broadcasts functions interpolant -<<<<<<< HEAD operations/index -======= - evaluate - simplify - jacobian ->>>>>>> master diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py index aaf8474eb6..9d192af47a 100644 --- a/pybamm/models/base_model.py +++ b/pybamm/models/base_model.py @@ -73,12 +73,13 @@ class BaseModel(object): 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 + - "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 + - "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"): From d1be2b9afaa491e1d4fc829130eb684b31bdfbc7 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Fri, 25 Oct 2019 17:23:09 -0400 Subject: [PATCH 14/23] #664 now getting ida conv fail :( --- examples/scripts/compare-dae-solver.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/scripts/compare-dae-solver.py b/examples/scripts/compare-dae-solver.py index 05476790cf..dd9c15e7da 100644 --- a/examples/scripts/compare-dae-solver.py +++ b/examples/scripts/compare-dae-solver.py @@ -17,7 +17,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 = {var.x_n: 50, var.x_s: 10, var.x_p: 50, var.r_n: 20, var.r_p: 20} # var_pts = model.default_var_pts mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) @@ -26,13 +26,14 @@ disc.process_model(model) # solve model -t_eval = np.linspace(0, 0.2, 100) +t_eval = np.linspace(0, 0.17, 100) +casadi_sol = pybamm.CasadiSolver(atol=1e-8, rtol=1e-8).solve(model, t_eval) klu_sol = pybamm.IDAKLU(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() From 372fa0cfb1e93658ec7651622ba84cd55d03a69b Mon Sep 17 00:00:00 2001 From: Travis CI Date: Fri, 25 Oct 2019 21:46:03 -0400 Subject: [PATCH 15/23] #664 fix casadi solver --- examples/scripts/compare-dae-solver.py | 4 +--- examples/scripts/compare_lithium_ion.py | 2 +- pybamm/solvers/casadi_solver.py | 5 ----- 3 files changed, 2 insertions(+), 9 deletions(-) diff --git a/examples/scripts/compare-dae-solver.py b/examples/scripts/compare-dae-solver.py index dd9c15e7da..99ee230333 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: 50, var.x_s: 10, var.x_p: 50, var.r_n: 20, var.r_p: 20} -# 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 diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py index fc06298566..2d37f5ced3 100644 --- a/examples/scripts/compare_lithium_ion.py +++ b/examples/scripts/compare_lithium_ion.py @@ -48,7 +48,7 @@ t_eval = np.linspace(0, 0.17, 100) for i, model in enumerate(models): model.convert_to_format = "casadi" - solutions[i] = model.default_solver.solve(model, t_eval) + solutions[i] = pybamm.CasadiSolver().solve(model, t_eval) # plot plot = pybamm.QuickPlot(models, mesh, solutions) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 5f373246f5..b12d36f17d 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -101,8 +101,3 @@ def integrate_casadi(self, problem, y0, t_eval, mass_matrix=None): # If it doesn't work raise error raise pybamm.SolverError(e.args[0]) - def calculate_consistent_initial_conditions( - self, rhs, algebraic, y0_guess, jac=None - ): - "No need to calculate initial conditions separately with this solver" - return y0_guess From 933354fc529e3371f47d993941124c48ff6fd37e Mon Sep 17 00:00:00 2001 From: Travis CI Date: Fri, 25 Oct 2019 21:47:02 -0400 Subject: [PATCH 16/23] #664 rename IDAKLU to IDAKLUSolver --- examples/scripts/compare-dae-solver.py | 2 +- pybamm/__init__.py | 2 +- pybamm/solvers/idaklu_solver.py | 2 +- results/2plus1D/dfn_2plus1D.py | 2 +- results/2plus1D/spm_2plus1D.py | 2 +- results/2plus1D/spme_2plus1D.py | 2 +- tests/integration/test_solvers/test_idaklu.py | 2 +- tests/unit/test_solvers/test_idaklu_solver.py | 2 +- 8 files changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/scripts/compare-dae-solver.py b/examples/scripts/compare-dae-solver.py index 99ee230333..e4dff68a3d 100644 --- a/examples/scripts/compare-dae-solver.py +++ b/examples/scripts/compare-dae-solver.py @@ -27,7 +27,7 @@ t_eval = np.linspace(0, 0.17, 100) casadi_sol = pybamm.CasadiSolver(atol=1e-8, rtol=1e-8).solve(model, t_eval) -klu_sol = pybamm.IDAKLU(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 diff --git a/pybamm/__init__.py b/pybamm/__init__.py index bb13078333..3a94c5dc99 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -266,7 +266,7 @@ def version(formatted=False): 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 IDAKLU, have_idaklu +from .solvers.idaklu_solver import IDAKLUSolver, have_idaklu # 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/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/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_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 From b1626991e345279bdec46c8ebb3e7e38dd61f372 Mon Sep 17 00:00:00 2001 From: tinosulzer Date: Fri, 25 Oct 2019 22:06:07 -0400 Subject: [PATCH 17/23] #664 codacy --- pybamm/expression_tree/functions.py | 39 ++++++++++++++++------------- 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py index 0464bb02ab..1f1d28e524 100644 --- a/pybamm/expression_tree/functions.py +++ b/pybamm/expression_tree/functions.py @@ -95,9 +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, i) * child.diff( - variable - ) + partial_derivatives[i] = self._function_diff( + children, i + ) * child.diff(variable) # remove None entries partial_derivatives = list(filter(None, partial_derivatives)) @@ -108,8 +108,11 @@ def diff(self, variable): return derivative - def _diff(self, children, idx): - """ 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( @@ -146,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, i) * children_jacs[i] + jac_fun = self._function_diff(children, i) * children_jacs[i] jac_fun.domain = [] if jacobian is None: jacobian = jac_fun @@ -272,8 +275,8 @@ class Cos(SpecificFunction): def __init__(self, child): super().__init__(np.cos, child) - def _diff(self, children, idx): - """ See :meth:`pybamm.Symbol._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Symbol._function_diff()`. """ return -Sin(children[0]) @@ -288,8 +291,8 @@ class Cosh(SpecificFunction): def __init__(self, child): super().__init__(np.cosh, child) - def _diff(self, children, idx): - """ See :meth:`pybamm.Function._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Function._function_diff()`. """ return Sinh(children[0]) @@ -304,8 +307,8 @@ class Exponential(SpecificFunction): def __init__(self, child): super().__init__(np.exp, child) - def _diff(self, children, idx): - """ See :meth:`pybamm.Function._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Function._function_diff()`. """ return Exponential(children[0]) @@ -320,8 +323,8 @@ class Log(SpecificFunction): def __init__(self, child): super().__init__(np.log, child) - def _diff(self, children, idx): - """ See :meth:`pybamm.Function._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Function._function_diff()`. """ return 1 / children[0] @@ -346,8 +349,8 @@ class Sin(SpecificFunction): def __init__(self, child): super().__init__(np.sin, child) - def _diff(self, children, idx): - """ See :meth:`pybamm.Function._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Function._function_diff()`. """ return Cos(children[0]) @@ -362,8 +365,8 @@ class Sinh(SpecificFunction): def __init__(self, child): super().__init__(np.sinh, child) - def _diff(self, children, idx): - """ See :meth:`pybamm.Function._diff()`. """ + def _function_diff(self, children, idx): + """ See :meth:`pybamm.Function._function_diff()`. """ return Cosh(children[0]) From b46512962a9991e5b5031e0174d5b1c33808dc0c Mon Sep 17 00:00:00 2001 From: tinosulzer Date: Sat, 26 Oct 2019 21:25:12 -0400 Subject: [PATCH 18/23] #664 coverage --- pybamm/parameters/electrical_parameters.py | 15 +----- pybamm/solvers/ode_solver.py | 1 - .../test_operations/test_convert_to_casadi.py | 11 +++++ .../test_electrical_parameters.py | 1 + tests/unit/test_solvers/test_ode_solver.py | 4 ++ .../unit/test_solvers/test_scikits_solvers.py | 46 ++++++++++--------- tests/unit/test_solvers/test_scipy_solver.py | 43 +++++++++-------- 7 files changed, 64 insertions(+), 57 deletions(-) diff --git a/pybamm/parameters/electrical_parameters.py b/pybamm/parameters/electrical_parameters.py index f9260935c0..853823c55a 100644 --- a/pybamm/parameters/electrical_parameters.py +++ b/pybamm/parameters/electrical_parameters.py @@ -1,21 +1,10 @@ # # Standard electrical parameters # -import casadi import pybamm import numpy as np -def abs_non_zero(x): - if x == 0: # pragma: no cover - return 1 - else: - if isinstance(x, casadi.SX): - return casadi.fabs(x) - else: - return abs(x) - - # -------------------------------------------------------------------------------------- # Dimensional Parameters I_typ = pybamm.Parameter("Typical current [A]") @@ -24,9 +13,7 @@ def abs_non_zero(x): n_electrodes_parallel = pybamm.Parameter( "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)) -) +i_typ = 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/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py index af416073ec..9b650241f0 100644 --- a/pybamm/solvers/ode_solver.py +++ b/pybamm/solvers/ode_solver.py @@ -208,7 +208,6 @@ def eval_event(t, y): # 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( 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 01de9b159d..f861bb09fd 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 @@ -76,6 +76,7 @@ def test_special_functions(self): a = np.array([1, 2, 3, 4, 5]) pybamm_a = pybamm.Array(a) self.assertEqual(pybamm.min(pybamm_a).to_casadi(), casadi.SX(1)) + self.assertEqual(pybamm.max(pybamm_a).to_casadi(), casadi.SX(5)) def test_concatenations(self): y = np.linspace(0, 1, 10)[:, np.newaxis] @@ -126,6 +127,16 @@ def myfunction(x, y): 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") diff --git a/tests/unit/test_parameters/test_electrical_parameters.py b/tests/unit/test_parameters/test_electrical_parameters.py index 8f8f7cad79..fb5d97ad8e 100644 --- a/tests/unit/test_parameters/test_electrical_parameters.py +++ b/tests/unit/test_parameters/test_electrical_parameters.py @@ -1,6 +1,7 @@ # # Tests for the electrical parameters # +import casadi import pybamm import unittest 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 a46e977925..c0977fd709 100644 --- a/tests/unit/test_solvers/test_scikits_solvers.py +++ b/tests/unit/test_solvers/test_scikits_solvers.py @@ -702,28 +702,30 @@ def test_model_solver_ode_events_casadi(self): def test_model_solver_dae_events_casadi(self): # Create model model = pybamm.BaseModel() - 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)) + 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__": diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py index 3315db7199..127a751873 100644 --- a/tests/unit/test_solvers/test_scipy_solver.py +++ b/tests/unit/test_solvers/test_scipy_solver.py @@ -271,26 +271,29 @@ def test_model_step(self): def test_model_solver_with_event_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: -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)) + 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__": From 9a0699f7a0dd927495f83f1f1874e1e8f02aa4da Mon Sep 17 00:00:00 2001 From: tinosulzer Date: Sat, 26 Oct 2019 22:52:21 -0400 Subject: [PATCH 19/23] #664 implement abs --- .../expression_tree/operations/convert_to_casadi.py | 2 ++ pybamm/parameters/electrical_parameters.py | 4 +++- .../test_operations/test_convert_to_casadi.py | 11 +++++++---- 3 files changed, 12 insertions(+), 5 deletions(-) diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py index 395bb6ec78..d613badead 100644 --- a/pybamm/expression_tree/operations/convert_to_casadi.py +++ b/pybamm/expression_tree/operations/convert_to_casadi.py @@ -72,6 +72,8 @@ def _convert(self, symbol, t, y): 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_"): diff --git a/pybamm/parameters/electrical_parameters.py b/pybamm/parameters/electrical_parameters.py index 853823c55a..8db66dbe29 100644 --- a/pybamm/parameters/electrical_parameters.py +++ b/pybamm/parameters/electrical_parameters.py @@ -13,7 +13,9 @@ n_electrodes_parallel = pybamm.Parameter( "Number of electrodes connected in parallel to make a cell" ) -i_typ = I_typ / (n_electrodes_parallel * pybamm.geometric_parameters.A_cc) +i_typ = pybamm.Function( + 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/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 f861bb09fd..a516cb7955 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 @@ -73,10 +73,13 @@ def test_convert_array_symbols(self): self.assertTrue(casadi.is_equal(outer.to_casadi(), casadi.SX(outer.evaluate()))) def test_special_functions(self): - a = np.array([1, 2, 3, 4, 5]) - pybamm_a = pybamm.Array(a) - self.assertEqual(pybamm.min(pybamm_a).to_casadi(), casadi.SX(1)) - self.assertEqual(pybamm.max(pybamm_a).to_casadi(), casadi.SX(5)) + 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] From b9d160cb25e5eefb40ef477f18162c4e8c270294 Mon Sep 17 00:00:00 2001 From: tinosulzer Date: Sun, 27 Oct 2019 09:14:21 -0400 Subject: [PATCH 20/23] #664 flake8 --- tests/unit/test_parameters/test_electrical_parameters.py | 1 - tests/unit/test_solvers/test_scipy_solver.py | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/unit/test_parameters/test_electrical_parameters.py b/tests/unit/test_parameters/test_electrical_parameters.py index fb5d97ad8e..8f8f7cad79 100644 --- a/tests/unit/test_parameters/test_electrical_parameters.py +++ b/tests/unit/test_parameters/test_electrical_parameters.py @@ -1,7 +1,6 @@ # # Tests for the electrical parameters # -import casadi import pybamm import unittest diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py index 127a751873..f5d7e783ce 100644 --- a/tests/unit/test_solvers/test_scipy_solver.py +++ b/tests/unit/test_solvers/test_scipy_solver.py @@ -279,7 +279,7 @@ def test_model_solver_with_event_with_casadi(self): 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 + # No need to set parameters; can use base discretisation (no spatial # operators) # create discretisation From 2402649d82e70c80f51ee227c766cee3936da8a5 Mon Sep 17 00:00:00 2001 From: tinosulzer Date: Sun, 27 Oct 2019 09:16:12 -0400 Subject: [PATCH 21/23] #664 changelog --- CHANGELOG.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index cbc2094ea7..80b9c20626 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,8 +2,8 @@ ## Features -- Add interface to CasADi solver -- Add option to use CasADi's Algorithmic Differentiation framework to calculate Jacobians +- 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)) @@ -19,7 +19,7 @@ ## Bug fixes -- Fix differentiation of functions that have more than one argument +- 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)) From a0ad975bdc778c885b4a1698d7cde67874bb5221 Mon Sep 17 00:00:00 2001 From: tinosulzer Date: Sun, 27 Oct 2019 13:19:48 -0400 Subject: [PATCH 22/23] #664 restore examples --- examples/scripts/compare_lithium_ion.py | 5 ++--- examples/scripts/compare_lithium_ion_3D.py | 1 - 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py index 2d37f5ced3..1aca4c7816 100644 --- a/examples/scripts/compare_lithium_ion.py +++ b/examples/scripts/compare_lithium_ion.py @@ -32,7 +32,7 @@ # set mesh var = pybamm.standard_spatial_vars -var_pts = {var.x_n: 50, var.x_s: 50, var.x_p: 50, var.r_n: 50, var.r_p: 50} +var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.r_n: 5, var.r_p: 5} # discretise models for model in models: @@ -47,8 +47,7 @@ solutions = [None] * len(models) t_eval = np.linspace(0, 0.17, 100) for i, model in enumerate(models): - model.convert_to_format = "casadi" - solutions[i] = pybamm.CasadiSolver().solve(model, t_eval) + solutions[i] = model.default_solver.solve(model, t_eval) # plot plot = pybamm.QuickPlot(models, mesh, solutions) diff --git a/examples/scripts/compare_lithium_ion_3D.py b/examples/scripts/compare_lithium_ion_3D.py index 9449f2e521..f271b63cf1 100644 --- a/examples/scripts/compare_lithium_ion_3D.py +++ b/examples/scripts/compare_lithium_ion_3D.py @@ -52,7 +52,6 @@ solutions = [None] * len(models) t_eval = np.linspace(0, 1, 1000) for i, model in enumerate(models): - model.convert_to_format = "casadi" solution = model.default_solver.solve(model, t_eval) solutions[i] = solution From 466a4c4dee9ebf945b8d96198e3b1c04becd8d81 Mon Sep 17 00:00:00 2001 From: tinosulzer Date: Mon, 28 Oct 2019 08:10:45 -0400 Subject: [PATCH 23/23] #664 remove termination from casadi solver --- pybamm/solvers/casadi_solver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index b12d36f17d..373bbc0563 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -55,8 +55,8 @@ def compute_solution(self, model, t_eval): ) solve_time = timer.time() - solve_start_time - # Identify the event that caused termination - termination = self.get_termination_reason(solution, self.events) + # Events not implemented, termination is always 'final time' + termination = "final time" return solution, solve_time, termination