From 1dd5ecd21604a0412027a0e3f506a0f4704766d8 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Thu, 20 Feb 2020 22:16:03 -0500 Subject: [PATCH 01/10] #846 tackle some of the bugs --- .../expression_tree/binary_operator.rst | 6 +++ pybamm/__init__.py | 2 + pybamm/expression_tree/binary_operators.py | 51 +++++++++++-------- pybamm/expression_tree/symbol.py | 16 +++--- pybamm/parameters/parameter_values.py | 5 +- pybamm/quick_plot.py | 2 +- .../test_binary_operators.py | 2 - .../test_finite_volume/test_finite_volume.py | 16 ++++++ 8 files changed, 67 insertions(+), 33 deletions(-) diff --git a/docs/source/expression_tree/binary_operator.rst b/docs/source/expression_tree/binary_operator.rst index d8082d8c4c..b46908ee89 100644 --- a/docs/source/expression_tree/binary_operator.rst +++ b/docs/source/expression_tree/binary_operator.rst @@ -28,4 +28,10 @@ Binary Operators .. autoclass:: pybamm.Heaviside :members: +.. autoclass:: pybamm.EqualHeaviside + :members: + +.. autoclass:: pybamm.NotEqualHeaviside + :members: + .. autofunction:: pybamm.source diff --git a/pybamm/__init__.py b/pybamm/__init__.py index 022ac64370..7c8fc9a96c 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -82,6 +82,8 @@ def version(formatted=False): Inner, inner, Heaviside, + EqualHeaviside, + NotEqualHeaviside, source, ) from .expression_tree.concatenations import ( diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py index 65c459926b..ab27513123 100644 --- a/pybamm/expression_tree/binary_operators.py +++ b/pybamm/expression_tree/binary_operators.py @@ -657,23 +657,10 @@ class Heaviside(BinaryOperator): **Extends:** :class:`BinaryOperator` """ - def __init__(self, left, right, equal): + def __init__(self, name, left, right): """ See :meth:`pybamm.BinaryOperator.__init__()`. """ - # 'equal' determines whether to return 1 or 0 when left = right - self.equal = equal - if equal is True: - name = "<=" - else: - name = "<" super().__init__(name, left, right) - def __str__(self): - """ See :meth:`pybamm.Symbol.__str__()`. """ - if self.equal is True: - return "{!s} <= {!s}".format(self.left, self.right) - else: - return "{!s} < {!s}".format(self.left, self.right) - def diff(self, variable): """ See :meth:`pybamm.Symbol.diff()`. """ # Heaviside should always be multiplied by something else so hopefully don't @@ -686,18 +673,40 @@ def _binary_jac(self, left_jac, right_jac): # need to worry about shape return pybamm.Scalar(0) + +class EqualHeaviside(Heaviside): + "A heaviside function with equality (return 1 when left = right)" + + def __init__(self, left, right): + """ See :meth:`pybamm.BinaryOperator.__init__()`. """ + super().__init__("<=", left, right) + + def __str__(self): + """ See :meth:`pybamm.Symbol.__str__()`. """ + return "{!s} <= {!s}".format(self.left, self.right) + def _binary_evaluate(self, left, right): """ See :meth:`pybamm.BinaryOperator._binary_evaluate()`. """ # don't raise RuntimeWarning for NaNs with np.errstate(invalid="ignore"): - if self.equal is True: - return left <= right - else: - return left < right + return left <= right - def _binary_new_copy(self, left, right): - """ See :meth:`pybamm.BinaryOperator._binary_new_copy()`. """ - return Heaviside(left, right, self.equal) + +class NotEqualHeaviside(Heaviside): + "A heaviside function without equality (return 0 when left = right)" + + def __init__(self, left, right): + super().__init__("<", left, right) + + def __str__(self): + """ See :meth:`pybamm.Symbol.__str__()`. """ + return "{!s} < {!s}".format(self.left, self.right) + + def _binary_evaluate(self, left, right): + """ See :meth:`pybamm.BinaryOperator._binary_evaluate()`. """ + # don't raise RuntimeWarning for NaNs + with np.errstate(invalid="ignore"): + return left < right def source(left, right, boundary=False): diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py index 3bc390cc26..b61efc36e1 100644 --- a/pybamm/expression_tree/symbol.py +++ b/pybamm/expression_tree/symbol.py @@ -431,27 +431,27 @@ def __rpow__(self, other): return pybamm.simplify_if_constant(pybamm.Power(other, self), keep_domains=True) def __lt__(self, other): - """return a :class:`Heaviside` object""" + """return a :class:`NotEqualHeaviside` object""" return pybamm.simplify_if_constant( - pybamm.Heaviside(self, other, equal=False), keep_domains=True + pybamm.NotEqualHeaviside(self, other), keep_domains=True ) def __le__(self, other): - """return a :class:`Heaviside` object""" + """return a :class:`EqualHeaviside` object""" return pybamm.simplify_if_constant( - pybamm.Heaviside(self, other, equal=True), keep_domains=True + pybamm.EqualHeaviside(self, other), keep_domains=True ) def __gt__(self, other): - """return a :class:`Heaviside` object""" + """return a :class:`NotEqualHeaviside` object""" return pybamm.simplify_if_constant( - pybamm.Heaviside(other, self, equal=False), keep_domains=True + pybamm.NotEqualHeaviside(other, self), keep_domains=True ) def __ge__(self, other): - """return a :class:`Heaviside` object""" + """return a :class:`EqualHeaviside` object""" return pybamm.simplify_if_constant( - pybamm.Heaviside(other, self, equal=True), keep_domains=True + pybamm.EqualHeaviside(other, self), keep_domains=True ) def __neg__(self): diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py index 784dc74db7..0df4c5206b 100644 --- a/pybamm/parameters/parameter_values.py +++ b/pybamm/parameters/parameter_values.py @@ -72,9 +72,12 @@ def __init__(self, values=None, chemistry=None): if values is not None: # If base_parameters is a filename, load from that filename if isinstance(values, str): + path = os.path.split(values)[0] values = self.read_parameters_csv(values) + else: + path = None # Don't check parameter already exists when first creating it - self.update(values, check_already_exists=False) + self.update(values, check_already_exists=False, path=path) # Initialise empty _processed_symbols dict (for caching) self._processed_symbols = {} diff --git a/pybamm/quick_plot.py b/pybamm/quick_plot.py index 4d44fb8d77..b604d8dcb7 100644 --- a/pybamm/quick_plot.py +++ b/pybamm/quick_plot.py @@ -370,7 +370,7 @@ def dynamic_plot(self, testing=False): axcolor = "lightgoldenrodyellow" axfreq = plt.axes([0.315, 0.02, 0.37, 0.03], facecolor=axcolor) - self.sfreq = Slider(axfreq, "Time", 0, self.max_t, valinit=0) + self.sfreq = Slider(axfreq, "Time [h]", 0, self.max_t, valinit=0) self.sfreq.on_changed(self.update) # ignore the warning about tight layout diff --git a/tests/unit/test_expression_tree/test_binary_operators.py b/tests/unit/test_expression_tree/test_binary_operators.py index a35a4df290..78d8ef94a0 100644 --- a/tests/unit/test_expression_tree/test_binary_operators.py +++ b/tests/unit/test_expression_tree/test_binary_operators.py @@ -292,14 +292,12 @@ def test_heaviside(self): a = pybamm.Scalar(1) b = pybamm.StateVector(slice(0, 1)) heav = a < b - self.assertFalse(heav.equal) self.assertEqual(heav.evaluate(y=np.array([2])), 1) self.assertEqual(heav.evaluate(y=np.array([1])), 0) self.assertEqual(heav.evaluate(y=np.array([0])), 0) self.assertEqual(str(heav), "1.0 < y[0:1]") heav = a >= b - self.assertTrue(heav.equal) self.assertEqual(heav.evaluate(y=np.array([2])), 0) self.assertEqual(heav.evaluate(y=np.array([1])), 1) self.assertEqual(heav.evaluate(y=np.array([0])), 1) diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py index af36bdcfc5..cc74227b3b 100644 --- a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py +++ b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py @@ -1117,6 +1117,22 @@ def test_delta_function(self): np.sum(delta_fn_int_disc.evaluate(y=y)), ) + def test_heaviside(self): + mesh = get_mesh_for_testing() + spatial_methods = {"macroscale": pybamm.FiniteVolume()} + disc = pybamm.Discretisation(mesh, spatial_methods) + + var = pybamm.Variable("var", domain="negative electrode") + heav = var > 1 + + disc.set_variable_slices([var]) + # process_binary_operators should work with heaviside + disc_heav = disc.process_symbol(heav * var) + nodes = mesh["negative electrode"][0].nodes + self.assertEqual(disc_heav.size, nodes.size) + np.testing.assert_array_equal(disc_heav.evaluate(y=2 * np.ones_like(nodes)), 2) + np.testing.assert_array_equal(disc_heav.evaluate(y=-2 * np.ones_like(nodes)), 0) + def test_grad_div_with_bcs_on_tab(self): # 2d macroscale mesh = get_1p1d_mesh_for_testing() From 716f8ea6226b754eefff664e0d11d98d487f1ee0 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Fri, 21 Feb 2020 17:02:03 -0500 Subject: [PATCH 02/10] #735 add minimum and maximum --- .../expression_tree/binary_operator.rst | 10 +++ examples/scripts/DFN.py | 4 +- pybamm/__init__.py | 4 ++ pybamm/expression_tree/binary_operators.py | 72 +++++++++++++++++++ pybamm/expression_tree/functions.py | 10 ++- .../operations/convert_to_casadi.py | 4 ++ pybamm/expression_tree/operations/evaluate.py | 4 ++ .../full_stefan_maxwell_diffusion.py | 2 +- pybamm/solvers/base_solver.py | 13 +++- .../test_binary_operators.py | 15 ++++ .../test_operations/test_convert_to_casadi.py | 15 ++-- .../test_operations/test_evaluate.py | 16 ++++- .../test_operations/test_jac.py | 12 ++++ 13 files changed, 168 insertions(+), 13 deletions(-) diff --git a/docs/source/expression_tree/binary_operator.rst b/docs/source/expression_tree/binary_operator.rst index b46908ee89..826956a3ef 100644 --- a/docs/source/expression_tree/binary_operator.rst +++ b/docs/source/expression_tree/binary_operator.rst @@ -34,4 +34,14 @@ Binary Operators .. autoclass:: pybamm.NotEqualHeaviside :members: +.. autoclass:: pybamm.Minimum + :members: + +.. autoclass:: pybamm.Maximum + :members: + +.. autofunction:: pybamm.minimum + +.. autofunction:: pybamm.maximum + .. autofunction:: pybamm.source diff --git a/examples/scripts/DFN.py b/examples/scripts/DFN.py index 12eda03ae6..b78ec6ff61 100644 --- a/examples/scripts/DFN.py +++ b/examples/scripts/DFN.py @@ -10,7 +10,7 @@ # load model model = pybamm.lithium_ion.DFN() - +model.convert_to_format = "python" # create geometry geometry = model.default_geometry @@ -30,7 +30,7 @@ # solve model t_eval = np.linspace(0, 3600, 100) -solver = model.default_solver +solver = pybamm.IDAKLUSolver(root_method="lm") solver.rtol = 1e-3 solver.atol = 1e-6 solution = solver.solve(model, t_eval) diff --git a/pybamm/__init__.py b/pybamm/__init__.py index 7c8fc9a96c..f00e5b978e 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -84,6 +84,10 @@ def version(formatted=False): Heaviside, EqualHeaviside, NotEqualHeaviside, + Minimum, + minimum, + Maximum, + maximum, source, ) from .expression_tree.concatenations import ( diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py index ab27513123..392c1641ab 100644 --- a/pybamm/expression_tree/binary_operators.py +++ b/pybamm/expression_tree/binary_operators.py @@ -709,6 +709,78 @@ def _binary_evaluate(self, left, right): return left < right +class Minimum(BinaryOperator): + " Returns the smaller of two objects " + + def __init__(self, left, right): + super().__init__("minimum", left, right) + + def __str__(self): + """ See :meth:`pybamm.Symbol.__str__()`. """ + return "minimum({!s}, {!s})".format(self.left, self.right) + + def _diff(self, variable): + """ See :meth:`pybamm.Symbol._diff()`. """ + left, right = self.orphans + return (left <= right) * left.diff(variable) + (left > right) * right.diff( + variable + ) + + def _binary_jac(self, left_jac, right_jac): + """ See :meth:`pybamm.BinaryOperator._binary_jac()`. """ + left, right = self.orphans + return (left <= right) * left_jac + (left > right) * right_jac + + def _binary_evaluate(self, left, right): + """ See :meth:`pybamm.BinaryOperator._binary_evaluate()`. """ + # don't raise RuntimeWarning for NaNs + return np.minimum(left, right) + + +class Maximum(BinaryOperator): + " Returns the smaller of two objects " + + def __init__(self, left, right): + super().__init__("maximum", left, right) + + def __str__(self): + """ See :meth:`pybamm.Symbol.__str__()`. """ + return "maximum({!s}, {!s})".format(self.left, self.right) + + def _diff(self, variable): + """ See :meth:`pybamm.Symbol._diff()`. """ + left, right = self.orphans + return (left >= right) * left.diff(variable) + (left < right) * right.diff( + variable + ) + + def _binary_jac(self, left_jac, right_jac): + """ See :meth:`pybamm.BinaryOperator._binary_jac()`. """ + left, right = self.orphans + return (left >= right) * left_jac + (left < right) * right_jac + + def _binary_evaluate(self, left, right): + """ See :meth:`pybamm.BinaryOperator._binary_evaluate()`. """ + # don't raise RuntimeWarning for NaNs + return np.maximum(left, right) + + +def minimum(left, right): + """ + Returns the smaller of two objects. Not to be confused with :meth:`pybamm.min`, + which returns min function of child. + """ + return pybamm.simplify_if_constant(Minimum(left, right), keep_domains=True) + + +def maximum(left, right): + """ + Returns the larger of two objects. Not to be confused with :meth:`pybamm.max`, + which returns max function of child. + """ + return pybamm.simplify_if_constant(Maximum(left, right), keep_domains=True) + + def source(left, right, boundary=False): """A convinience function for creating (part of) an expression tree representing a source term. This is necessary for spatial methods where the mass matrix diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py index 2601034b98..cd6b50bf01 100644 --- a/pybamm/expression_tree/functions.py +++ b/pybamm/expression_tree/functions.py @@ -341,12 +341,18 @@ def log10(child): def max(child): - " Returns max function of child. " + """ + Returns max function of child. Not to be confused with :meth:`pybamm.maximum`, which + returns the larger of two objects. + """ return pybamm.simplify_if_constant(Function(np.max, child), keep_domains=True) def min(child): - " Returns min function of child. " + """ + Returns min function of child. Not to be confused with :meth:`pybamm.minimum`, which + returns the smaller of two objects. + """ return pybamm.simplify_if_constant(Function(np.min, child), keep_domains=True) diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py index a04aec5b0a..9525780e73 100644 --- a/pybamm/expression_tree/operations/convert_to_casadi.py +++ b/pybamm/expression_tree/operations/convert_to_casadi.py @@ -66,6 +66,10 @@ def _convert(self, symbol, t=None, y=None, u=None): # process children converted_left = self.convert(left, t, y, u) converted_right = self.convert(right, t, y, u) + if isinstance(symbol, pybamm.Minimum): + return casadi.fmin(converted_left, converted_right) + if isinstance(symbol, pybamm.Maximum): + return casadi.fmax(converted_left, converted_right) # _binary_evaluate defined in derived classes for specific rules return symbol._binary_evaluate(converted_left, converted_right) diff --git a/pybamm/expression_tree/operations/evaluate.py b/pybamm/expression_tree/operations/evaluate.py index a11769b065..1e1436fe2f 100644 --- a/pybamm/expression_tree/operations/evaluate.py +++ b/pybamm/expression_tree/operations/evaluate.py @@ -92,6 +92,10 @@ def find_symbols(symbol, constant_symbols, variable_symbols): "if scipy.sparse.issparse({1}) else " "{0} * {1}".format(children_vars[0], children_vars[1]) ) + elif isinstance(symbol, pybamm.Minimum): + symbol_str = "np.minimum({},{})".format(children_vars[0], children_vars[1]) + elif isinstance(symbol, pybamm.Maximum): + symbol_str = "np.maximum({},{})".format(children_vars[0], children_vars[1]) else: symbol_str = children_vars[0] + " " + symbol.name + " " + children_vars[1] diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py index e007a5e6f2..07c1d5ea82 100644 --- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py +++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py @@ -45,7 +45,7 @@ def get_coupled_variables(self, variables): # N_e = N_e_diffusion + N_e_migration + N_e_convection - N_e = N_e_diffusion + c_e * v_box + N_e = N_e_diffusion + v_box variables.update(self._get_standard_flux_variables(N_e)) diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index 806015e900..59043fba4a 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -131,10 +131,19 @@ def set_up(self, model, inputs=None): if self.ode_solver is True: self.root_method = None if ( - isinstance(self, pybamm.CasadiSolver) or self.root_method == "casadi" + isinstance(self, pybamm.CasadiSolver) ) and model.convert_to_format != "casadi": pybamm.logger.warning( - f"Converting {model.name} to CasADi for solving with CasADi solver" + "Converting {} to CasADi for solving with CasADi solver".format( + model.name + ) + ) + model.convert_to_format = "casadi" + if self.root_method == "casadi" and model.convert_to_format != "casadi": + pybamm.logger.warning( + "Converting {} to CasADi for calculating ICs with CasADi".format( + model.name + ) ) model.convert_to_format = "casadi" diff --git a/tests/unit/test_expression_tree/test_binary_operators.py b/tests/unit/test_expression_tree/test_binary_operators.py index 78d8ef94a0..e6c8dc26b0 100644 --- a/tests/unit/test_expression_tree/test_binary_operators.py +++ b/tests/unit/test_expression_tree/test_binary_operators.py @@ -303,6 +303,21 @@ def test_heaviside(self): self.assertEqual(heav.evaluate(y=np.array([0])), 1) self.assertEqual(str(heav), "y[0:1] <= 1.0") + def test_minimum_maximum(self): + a = pybamm.Scalar(1) + b = pybamm.StateVector(slice(0, 1)) + minimum = pybamm.minimum(a, b) + self.assertEqual(minimum.evaluate(y=np.array([2])), 1) + self.assertEqual(minimum.evaluate(y=np.array([1])), 1) + self.assertEqual(minimum.evaluate(y=np.array([0])), 0) + self.assertEqual(str(minimum), "minimum(1.0, y[0:1])") + + maximum = pybamm.maximum(a, b) + self.assertEqual(maximum.evaluate(y=np.array([2])), 2) + self.assertEqual(maximum.evaluate(y=np.array([1])), 1) + self.assertEqual(maximum.evaluate(y=np.array([0])), 1) + self.assertEqual(str(maximum), "maximum(1.0, y[0:1])") + class TestIsZero(unittest.TestCase): def test_is_scalar_zero(self): 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 017cae7f7a..915a449c31 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 @@ -42,16 +42,21 @@ def myfunction(x, y): f = pybamm.Function(myfunction, b, d) self.assertEqual(f.to_casadi(), casadi.MX(3)) + # use classes to avoid simplification # addition - self.assertEqual((a + b).to_casadi(), casadi.MX(1)) + self.assertEqual((pybamm.Addition(a, b)).to_casadi(), casadi.MX(1)) # subtraction - self.assertEqual((c - d).to_casadi(), casadi.MX(-3)) + self.assertEqual(pybamm.Subtraction(c, d).to_casadi(), casadi.MX(-3)) # multiplication - self.assertEqual((c * d).to_casadi(), casadi.MX(-2)) + self.assertEqual(pybamm.Multiplication(c, d).to_casadi(), casadi.MX(-2)) # power - self.assertEqual((c ** d).to_casadi(), casadi.MX(1)) + self.assertEqual(pybamm.Power(c, d).to_casadi(), casadi.MX(1)) # division - self.assertEqual((b / d).to_casadi(), casadi.MX(1 / 2)) + self.assertEqual(pybamm.Division(b, d).to_casadi(), casadi.MX(1 / 2)) + + # minimum and maximum + self.assertEqual(pybamm.Minimum(a, b).to_casadi(), casadi.MX(0)) + self.assertEqual(pybamm.Maximum(a, b).to_casadi(), casadi.MX(1)) def test_convert_array_symbols(self): # Arrays diff --git a/tests/unit/test_expression_tree/test_operations/test_evaluate.py b/tests/unit/test_expression_tree/test_operations/test_evaluate.py index 88bc911e1f..ba4f642ed9 100644 --- a/tests/unit/test_expression_tree/test_operations/test_evaluate.py +++ b/tests/unit/test_expression_tree/test_operations/test_evaluate.py @@ -19,7 +19,7 @@ def test_find_symbols(self): a = pybamm.StateVector(slice(0, 1)) b = pybamm.StateVector(slice(1, 2)) - # test a * b + # test a + b constant_symbols = OrderedDict() variable_symbols = OrderedDict() expr = a + b @@ -356,6 +356,20 @@ def test_evaluator_python(self): result = evaluator.evaluate(t=t, y=y) np.testing.assert_allclose(result, expr.evaluate(t=t, y=y)) + # test something with a minimum or maximum + a = pybamm.Vector(np.array([1, 2])) + expr = pybamm.minimum(a, pybamm.StateVector(slice(0, 2))) + evaluator = pybamm.EvaluatorPython(expr) + for t, y in zip(t_tests, y_tests): + result = evaluator.evaluate(t=t, y=y) + np.testing.assert_allclose(result, expr.evaluate(t=t, y=y)) + + expr = pybamm.maximum(a, pybamm.StateVector(slice(0, 2))) + evaluator = pybamm.EvaluatorPython(expr) + for t, y in zip(t_tests, y_tests): + result = evaluator.evaluate(t=t, y=y) + np.testing.assert_allclose(result, expr.evaluate(t=t, y=y)) + # test something with an index expr = pybamm.Index(A @ pybamm.StateVector(slice(0, 2)), 0) evaluator = pybamm.EvaluatorPython(expr) diff --git a/tests/unit/test_expression_tree/test_operations/test_jac.py b/tests/unit/test_expression_tree/test_operations/test_jac.py index 270b030c44..a33791306f 100644 --- a/tests/unit/test_expression_tree/test_operations/test_jac.py +++ b/tests/unit/test_expression_tree/test_operations/test_jac.py @@ -248,6 +248,18 @@ def test_jac_of_heaviside(self): ((a < y) * y ** 2).jac(y).evaluate(y=-5 * np.ones(5)), 0 ) + def test_jac_of_minimum_maximum(self): + y = pybamm.StateVector(slice(0, 10)) + y_test = np.linspace(0, 2, 10) + np.testing.assert_array_equal( + np.diag(pybamm.minimum(1, y ** 2).jac(y).evaluate(y=y_test)), + 2 * y_test * (y_test < 1), + ) + np.testing.assert_array_equal( + np.diag(pybamm.maximum(1, y ** 2).jac(y).evaluate(y=y_test)), + 2 * y_test * (y_test > 1), + ) + def test_jac_of_domain_concatenation(self): # create mesh mesh = get_mesh_for_testing() From 4537550d7449556c8b893925bfd5104c0aef7e1b Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Fri, 21 Feb 2020 18:13:28 -0500 Subject: [PATCH 03/10] #742 add derivative of abs node --- .../source/expression_tree/unary_operator.rst | 11 ++++ pybamm/__init__.py | 63 ++----------------- pybamm/expression_tree/exceptions.py | 8 --- pybamm/expression_tree/unary_operators.py | 45 ++++++++++--- .../test_operations/test_jac.py | 20 ++++-- .../test_symbolic_diff.py | 14 +++++ .../test_unary_operators.py | 17 +++-- 7 files changed, 96 insertions(+), 82 deletions(-) diff --git a/docs/source/expression_tree/unary_operator.rst b/docs/source/expression_tree/unary_operator.rst index ed0bc19af3..d481f51bd3 100644 --- a/docs/source/expression_tree/unary_operator.rst +++ b/docs/source/expression_tree/unary_operator.rst @@ -10,6 +10,9 @@ Unary Operators .. autoclass:: pybamm.AbsoluteValue :members: +.. autoclass:: pybamm.Sign + :members: + .. autoclass:: pybamm.Index :members: @@ -67,4 +70,12 @@ Unary Operators .. autofunction:: pybamm.x_average +.. autofunction:: pybamm.r_average + +.. autofunction:: pybamm.z_average + +.. autofunction:: pybamm.yz_average + .. autofunction:: pybamm.boundary_value + +.. autofunction:: pybamm.sign diff --git a/pybamm/__init__.py b/pybamm/__init__.py index f00e5b978e..72624f66dd 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -63,39 +63,9 @@ def version(formatted=False): # # Classes for the Expression Tree # -from .expression_tree.symbol import ( - Symbol, - domain_size, - create_object_of_size, - evaluate_for_shape_using_domain, -) -from .expression_tree.binary_operators import ( - is_scalar_zero, - is_matrix_zero, - BinaryOperator, - Addition, - Power, - Subtraction, - Multiplication, - MatrixMultiplication, - Division, - Inner, - inner, - Heaviside, - EqualHeaviside, - NotEqualHeaviside, - Minimum, - minimum, - Maximum, - maximum, - source, -) -from .expression_tree.concatenations import ( - Concatenation, - NumpyConcatenation, - DomainConcatenation, - SparseStack, -) +from .expression_tree.symbol import * +from .expression_tree.binary_operators import * +from .expression_tree.concatenations import * from .expression_tree.array import Array from .expression_tree.matrix import Matrix from .expression_tree.unary_operators import * @@ -103,36 +73,15 @@ def version(formatted=False): from .expression_tree.interpolant import Interpolant from .expression_tree.input_parameter import InputParameter from .expression_tree.parameter import Parameter, FunctionParameter -from .expression_tree.broadcasts import ( - Broadcast, - PrimaryBroadcast, - SecondaryBroadcast, - FullBroadcast, - ones_like, -) +from .expression_tree.broadcasts import * from .expression_tree.scalar import Scalar from .expression_tree.variable import Variable, ExternalVariable -from .expression_tree.independent_variable import ( - IndependentVariable, - Time, - SpatialVariable, -) +from .expression_tree.independent_variable import * from .expression_tree.independent_variable import t from .expression_tree.vector import Vector from .expression_tree.state_vector import StateVector -from .expression_tree.exceptions import ( - DomainError, - OptionError, - ModelError, - SolverError, - SolverWarning, - ShapeError, - ModelWarning, - UndefinedOperationError, - GeometryError, - InputError, -) +from .expression_tree.exceptions import * # Operations from .expression_tree.operations.simplify import ( diff --git a/pybamm/expression_tree/exceptions.py b/pybamm/expression_tree/exceptions.py index a71172cc48..d0f4341a20 100644 --- a/pybamm/expression_tree/exceptions.py +++ b/pybamm/expression_tree/exceptions.py @@ -61,14 +61,6 @@ class ModelWarning(UserWarning): pass -class UndefinedOperationError(Exception): - """ - Undefined operation: Raised when a mathematical operation is not well-defined - """ - - pass - - class InputError(Exception): """ An external variable has been input incorrectly into PyBaMM diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py index 0b4f68e22d..cf7f94faf4 100644 --- a/pybamm/expression_tree/unary_operators.py +++ b/pybamm/expression_tree/unary_operators.py @@ -3,7 +3,7 @@ # import numpy as np import pybamm -from scipy.sparse import csr_matrix +from scipy.sparse import issparse, csr_matrix class UnaryOperator(pybamm.Symbol): @@ -125,23 +125,45 @@ def __init__(self, child): def diff(self, variable): """ See :meth:`pybamm.Symbol.diff()`. """ - # Derivative is not well-defined - raise pybamm.UndefinedOperationError( - "Derivative of absolute function is not defined" - ) + child = self.child.new_copy() + return Sign(child) * child.diff(variable) def _unary_jac(self, child_jac): """ See :meth:`pybamm.UnaryOperator._unary_jac()`. """ - # Derivative is not well-defined - raise pybamm.UndefinedOperationError( - "Derivative of absolute function is not defined" - ) + child = self.child.new_copy() + return Sign(child) * child_jac def _unary_evaluate(self, child): """ See :meth:`UnaryOperator._unary_evaluate()`. """ return np.abs(child) +class Sign(UnaryOperator): + """A node in the expression tree representing a `sign` operator + + **Extends:** :class:`UnaryOperator` + """ + + def __init__(self, child): + """ See :meth:`pybamm.UnaryOperator.__init__()`. """ + super().__init__("sign", child) + + def diff(self, variable): + """ See :meth:`pybamm.Symbol.diff()`. """ + return pybamm.Scalar(0) + + def _unary_jac(self, child_jac): + """ See :meth:`pybamm.UnaryOperator._unary_jac()`. """ + return pybamm.Scalar(0) + + def _unary_evaluate(self, child): + """ See :meth:`UnaryOperator._unary_evaluate()`. """ + if issparse(child): + return csr_matrix.sign(child) + else: + return np.sign(child) + + class Index(UnaryOperator): """A node in the expression tree, which stores the index that should be extracted from its child after the child has been evaluated. @@ -1067,3 +1089,8 @@ def r_average(symbol): pybamm.Scalar(1), symbol.domain, symbol.auxiliary_domains ) return Integral(symbol, r) / Integral(v, r) + + +def sign(symbol): + " Returns a :class:`Sign` object. " + return Sign(symbol) diff --git a/tests/unit/test_expression_tree/test_operations/test_jac.py b/tests/unit/test_expression_tree/test_operations/test_jac.py index a33791306f..2711e06d29 100644 --- a/tests/unit/test_expression_tree/test_operations/test_jac.py +++ b/tests/unit/test_expression_tree/test_operations/test_jac.py @@ -105,10 +105,6 @@ def test_nonlinear(self): dfunc_dy = func.jac(y).evaluate(y=y0) np.testing.assert_array_equal(jacobian, dfunc_dy.toarray()) - func = pybamm.AbsoluteValue(v) - with self.assertRaises(pybamm.UndefinedOperationError): - func.jac(y) - def test_functions(self): y = pybamm.StateVector(slice(0, 4)) u = pybamm.StateVector(slice(0, 2)) @@ -260,6 +256,22 @@ def test_jac_of_minimum_maximum(self): 2 * y_test * (y_test > 1), ) + def test_jac_of_abs(self): + y = pybamm.StateVector(slice(0, 10)) + absy = abs(y) + jac = absy.jac(y) + y_test = np.linspace(-2, 2, 10) + np.testing.assert_array_equal( + np.diag(jac.evaluate(y=y_test).toarray()), np.sign(y_test) + ) + + def test_jac_of_sign(self): + y = pybamm.StateVector(slice(0, 10)) + func = pybamm.sign(y) * y + jac = func.jac(y) + y_test = np.linspace(-2, 2, 10) + np.testing.assert_array_equal(np.diag(jac.evaluate(y=y_test)), np.sign(y_test)) + def test_jac_of_domain_concatenation(self): # create mesh mesh = get_mesh_for_testing() diff --git a/tests/unit/test_expression_tree/test_symbolic_diff.py b/tests/unit/test_expression_tree/test_symbolic_diff.py index 1c0e6d1e55..68370ddfd3 100644 --- a/tests/unit/test_expression_tree/test_symbolic_diff.py +++ b/tests/unit/test_expression_tree/test_symbolic_diff.py @@ -68,6 +68,20 @@ def test_diff_heaviside(self): self.assertEqual(func.diff(b).evaluate(y=np.array([2])), 2) self.assertEqual(func.diff(b).evaluate(y=np.array([-2])), 0) + def test_diff_maximum_minimum(self): + a = pybamm.Scalar(1) + b = pybamm.StateVector(slice(0, 1)) + + func = pybamm.minimum(a, b ** 3) + self.assertEqual(func.diff(b).evaluate(y=np.array([10])), 0) + self.assertEqual(func.diff(b).evaluate(y=np.array([2])), 0) + self.assertEqual(func.diff(b).evaluate(y=np.array([-2])), 3 * (-2) ** 2) + + func = pybamm.maximum(a, b ** 3) + self.assertEqual(func.diff(b).evaluate(y=np.array([10])), 3 * 10 ** 2) + self.assertEqual(func.diff(b).evaluate(y=np.array([2])), 3 * 2 ** 2) + self.assertEqual(func.diff(b).evaluate(y=np.array([-2])), 0) + def test_exceptions(self): a = pybamm.Symbol("a") b = pybamm.Symbol("b") diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py index 4ae1eb7740..4944afcbf8 100644 --- a/tests/unit/test_expression_tree/test_unary_operators.py +++ b/tests/unit/test_expression_tree/test_unary_operators.py @@ -38,6 +38,11 @@ def test_absolute(self): absb = pybamm.AbsoluteValue(b) self.assertEqual(absb.evaluate(), 4) + def test_sign(self): + b = pybamm.Scalar(-4) + signb = pybamm.sign(b) + self.assertEqual(signb.evaluate(), -1) + def test_gradient(self): a = pybamm.Symbol("a") grad = pybamm.Gradient(a) @@ -158,10 +163,14 @@ def test_diff(self): self.assertEqual((-a).diff(a).evaluate(y=y), -1) self.assertEqual((-a).diff(-a).evaluate(), 1) - # absolute value (not implemented) - absa = abs(a) - with self.assertRaises(pybamm.UndefinedOperationError): - absa.diff(a) + # absolute value + self.assertEqual((a ** 3).diff(a).evaluate(y=y), 3 * 5 ** 2) + self.assertEqual((abs(a ** 3)).diff(a).evaluate(y=y), 3 * 5 ** 2) + self.assertEqual((a ** 3).diff(a).evaluate(y=-y), 3 * 5 ** 2) + self.assertEqual((abs(a ** 3)).diff(a).evaluate(y=-y), -3 * 5 ** 2) + + # sign + self.assertEqual((pybamm.sign(a)).diff(a).evaluate(y=y), 0) # spatial operator (not implemented) spatial_a = pybamm.SpatialOperator("name", a) From 1f53aac69198b0469a69459c448d98f93e710d0e Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Fri, 21 Feb 2020 18:45:26 -0500 Subject: [PATCH 04/10] #846 add more informative error for discretisation failure --- pybamm/discretisations/discretisation.py | 20 +++++++++++++++++-- .../test_discretisation.py | 17 ++++++++++++++++ 2 files changed, 35 insertions(+), 2 deletions(-) diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py index b74ac08a72..8f4e4ea3a5 100644 --- a/pybamm/discretisations/discretisation.py +++ b/pybamm/discretisations/discretisation.py @@ -234,7 +234,8 @@ def set_variable_slices(self, variables): y_slices[variable.id].append(slice(start, end)) start = end - self.y_slices = y_slices + # Convert y_slices back to normal dictionary + self.y_slices = dict(y_slices) # reset discretised_symbols self._discretised_symbols = {} @@ -885,8 +886,23 @@ def _process_symbol(self, symbol): return out else: + # add a try except block for a more informative error if a variable + # can't be found. This should usually be caught earlier by + # model.check_well_posedness, but won't be if debug_mode is False + try: + y_slices = self.y_slices[symbol.id] + except KeyError: + raise pybamm.ModelError( + """ + No key set for variable '{}'. Make sure it is included in either + model.rhs, model.algebraic, or model.external_variables in an + unmodified form (e.g. not Broadcasted) + """.format( + symbol.name + ) + ) return pybamm.StateVector( - *self.y_slices[symbol.id], + *y_slices, domain=symbol.domain, auxiliary_domains=symbol.auxiliary_domains ) diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py index d30b6f11f5..da38f07782 100644 --- a/tests/unit/test_discretisations/test_discretisation.py +++ b/tests/unit/test_discretisations/test_discretisation.py @@ -693,6 +693,23 @@ def test_process_model_ode(self): with self.assertRaises(pybamm.ModelError): disc.process_model(model) + def test_process_model_fail(self): + # one equation + c = pybamm.Variable("c") + d = pybamm.Variable("d") + model = pybamm.BaseModel() + model.rhs = {c: -c} + model.initial_conditions = {c: pybamm.Scalar(3)} + model.variables = {"c": c, "d": d} + + disc = pybamm.Discretisation() + # turn debug mode off to not check well posedness + debug_mode = pybamm.settings.debug_mode + pybamm.settings.debug_mode = False + with self.assertRaisesRegex(pybamm.ModelError, "No key set for variable"): + disc.process_model(model) + pybamm.settings.debug_mode = debug_mode + def test_process_model_dae(self): # one rhs equation and one algebraic whole_cell = ["negative electrode", "separator", "positive electrode"] From 98bab54fba04661ff5615dc5671ab2e8779caf6f Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 24 Feb 2020 10:23:44 -0500 Subject: [PATCH 05/10] #846 fix accidentally broken convection model --- .../stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py index 07c1d5ea82..e007a5e6f2 100644 --- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py +++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py @@ -45,7 +45,7 @@ def get_coupled_variables(self, variables): # N_e = N_e_diffusion + N_e_migration + N_e_convection - N_e = N_e_diffusion + v_box + N_e = N_e_diffusion + c_e * v_box variables.update(self._get_standard_flux_variables(N_e)) From c03af6217cd1e4607a4dd7506252d7432fbd2562 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 24 Feb 2020 12:40:51 -0500 Subject: [PATCH 06/10] #846 fix dfn example --- examples/scripts/DFN.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/scripts/DFN.py b/examples/scripts/DFN.py index b78ec6ff61..12eda03ae6 100644 --- a/examples/scripts/DFN.py +++ b/examples/scripts/DFN.py @@ -10,7 +10,7 @@ # load model model = pybamm.lithium_ion.DFN() -model.convert_to_format = "python" + # create geometry geometry = model.default_geometry @@ -30,7 +30,7 @@ # solve model t_eval = np.linspace(0, 3600, 100) -solver = pybamm.IDAKLUSolver(root_method="lm") +solver = model.default_solver solver.rtol = 1e-3 solver.atol = 1e-6 solution = solver.solve(model, t_eval) From da91cd97780c080e9ec0e9fefe1117c8eae326ad Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 24 Feb 2020 19:16:54 -0500 Subject: [PATCH 07/10] #846 windows test and skipping non-time discontinuities --- pybamm/solvers/base_solver.py | 19 ++++++++++++++----- .../test_parameters/test_parameter_values.py | 10 ++++++++-- .../unit/test_solvers/test_scikits_solvers.py | 12 +++++------- 3 files changed, 27 insertions(+), 14 deletions(-) diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index 59043fba4a..e4435c23c8 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -238,22 +238,31 @@ def report(string): model.concatenated_algebraic.pre_order(), ): if isinstance(symbol, pybamm.Heaviside): + found_t = False # Dimensionless if symbol.right.id == pybamm.t.id: expr = symbol.left + found_t = True elif symbol.left.id == pybamm.t.id: expr = symbol.right + found_t = True # Dimensional elif symbol.right.id == (pybamm.t * model.timescale).id: expr = symbol.left.new_copy() / symbol.right.right.new_copy() + found_t = True elif symbol.left.id == (pybamm.t * model.timescale).id: expr = symbol.right.new_copy() / symbol.left.right.new_copy() - - model.events.append( - pybamm.Event( - str(symbol), expr.new_copy(), pybamm.EventType.DISCONTINUITY + found_t = True + + # Update the events if the heaviside function depended on t + if found_t: + model.events.append( + pybamm.Event( + str(symbol), + expr.new_copy(), + pybamm.EventType.DISCONTINUITY, + ) ) - ) # Process rhs, algebraic and event expressions rhs, rhs_eval, jac_rhs = process(model.concatenated_rhs, "RHS") diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py index da3adfccf6..7b6993769c 100644 --- a/tests/unit/test_parameters/test_parameter_values.py +++ b/tests/unit/test_parameters/test_parameter_values.py @@ -26,8 +26,14 @@ def test_init(self): # from file param = pybamm.ParameterValues( - values="input/parameters/lithium-ion/cathodes/lico2_Marquis2019/" - + "parameters.csv" + values=os.path.join( + "input", + "parameters", + "lithium-ion", + "cathodes", + "lico2_Marquis2019", + "parameters.csv", + ) ) self.assertEqual(param["Reference temperature [K]"], 298.15) diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py index 905507c692..956c65f18c 100644 --- a/tests/unit/test_solvers/test_scikits_solvers.py +++ b/tests/unit/test_solvers/test_scikits_solvers.py @@ -262,7 +262,9 @@ def nonsmooth_mult(t): rate = pybamm.Function(nonsmooth_rate, pybamm.t) mult = pybamm.Function(nonsmooth_mult, pybamm.t) - model.rhs = {var1: rate * var1} + # put in an extra heaviside with no time dependence, this should be ignored by + # the solver i.e. no extra discontinuities added + model.rhs = {var1: rate * var1 + (var1 < 0)} model.algebraic = {var2: mult * var1 - var2} model.initial_conditions = {var1: 1, var2: 2} model.events = [ @@ -630,9 +632,7 @@ def nonsmooth_rate(t): model2.rhs = {var1: (0.1 * (pybamm.t < discontinuity) + 0.1) * var1} model2.algebraic = {var2: var2} model2.initial_conditions = {var1: 1, var2: 0} - model2.events = [ - pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)), - ] + model2.events = [pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5))] # third model implicitly adds a discontinuity event via another heaviside # function @@ -640,9 +640,7 @@ def nonsmooth_rate(t): model3.rhs = {var1: (-0.1 * (discontinuity < pybamm.t) + 0.2) * var1} model3.algebraic = {var2: var2} model3.initial_conditions = {var1: 1, var2: 0} - model3.events = [ - pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)), - ] + model3.events = [pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5))] for model in [model1, model2, model3]: From 10ddc82d92db9d5b10293804760b690ddf8a95b0 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 24 Feb 2020 20:14:02 -0500 Subject: [PATCH 08/10] #846 generalized x_average --- pybamm/__init__.py | 2 +- pybamm/expression_tree/independent_variable.py | 5 ----- pybamm/expression_tree/unary_operators.py | 7 +++---- .../test_independent_variable.py | 2 -- .../test_unary_operators.py | 17 +++++++++-------- 5 files changed, 13 insertions(+), 20 deletions(-) diff --git a/pybamm/__init__.py b/pybamm/__init__.py index 72624f66dd..d32b096f60 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -146,7 +146,7 @@ def version(formatted=False): Geometry2DCurrentCollector, ) -from .expression_tree.independent_variable import KNOWN_SPATIAL_VARS, KNOWN_COORD_SYS +from .expression_tree.independent_variable import KNOWN_COORD_SYS from .geometry import standard_spatial_vars # diff --git a/pybamm/expression_tree/independent_variable.py b/pybamm/expression_tree/independent_variable.py index 9dea22f5a7..8c11288d98 100644 --- a/pybamm/expression_tree/independent_variable.py +++ b/pybamm/expression_tree/independent_variable.py @@ -4,9 +4,6 @@ import pybamm KNOWN_COORD_SYS = ["cartesian", "spherical polar"] -KNOWN_SPATIAL_VARS = ["x", "y", "z", "r", "x_n", "x_s", "x_p", "r_n", "r_p"] -KNOWN_SPATIAL_VARS_EXTENDED = [v + "_edge" for v in KNOWN_SPATIAL_VARS] -KNOWN_SPATIAL_VARS.extend(KNOWN_SPATIAL_VARS_EXTENDED) class IndependentVariable(pybamm.Symbol): @@ -84,8 +81,6 @@ def __init__(self, name, domain=None, auxiliary_domains=None, coord_sys=None): super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains) domain = self.domain - if name not in KNOWN_SPATIAL_VARS: - raise ValueError(f"name must be in {KNOWN_SPATIAL_VARS} but is '{name}'") if domain == []: raise ValueError("domain must be provided") diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py index cf7f94faf4..f8fbc9c98e 100644 --- a/pybamm/expression_tree/unary_operators.py +++ b/pybamm/expression_tree/unary_operators.py @@ -453,8 +453,6 @@ def integration_variable(self): def set_id(self): """ See :meth:`pybamm.Symbol.set_id()` """ - if not isinstance(self.integration_variable, list): - self.integration_variable = [self.integration_variable] self._id = hash( (self.__class__, self.name) + tuple( @@ -946,8 +944,9 @@ def x_average(symbol): x = pybamm.standard_spatial_vars.x_p l = pybamm.geometric_parameters.l_p else: - raise pybamm.DomainError("domain '{}' not recognised".format(symbol.domain)) - + x = pybamm.SpatialVariable("x", domain=symbol.domain) + v = pybamm.ones_like(symbol) + l = pybamm.Integral(v, x) return Integral(symbol, x) / l diff --git a/tests/unit/test_expression_tree/test_independent_variable.py b/tests/unit/test_expression_tree/test_independent_variable.py index 734d839ec9..72f8f40620 100644 --- a/tests/unit/test_expression_tree/test_independent_variable.py +++ b/tests/unit/test_expression_tree/test_independent_variable.py @@ -45,8 +45,6 @@ def test_spatial_variable(self): with self.assertRaises(NotImplementedError): x.evaluate() - with self.assertRaisesRegex(ValueError, "name must be"): - pybamm.SpatialVariable("not a variable", ["negative electrode"]) with self.assertRaisesRegex(ValueError, "domain must be"): pybamm.SpatialVariable("x", []) with self.assertRaises(pybamm.DomainError): diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py index 4944afcbf8..42edd7fedf 100644 --- a/tests/unit/test_expression_tree/test_unary_operators.py +++ b/tests/unit/test_expression_tree/test_unary_operators.py @@ -291,12 +291,17 @@ def test_average(self): self.assertIsInstance(av_a, pybamm.Division) self.assertIsInstance(av_a.children[0], pybamm.Integral) self.assertEqual(av_a.children[0].integration_variable[0].domain, x.domain) - # electrode domains go to current collector when averaged self.assertEqual(av_a.domain, []) - a = pybamm.Symbol("a", domain="bad domain") - with self.assertRaises(pybamm.DomainError): - pybamm.x_average(a) + a = pybamm.Symbol("a", domain="new domain") + av_a = pybamm.x_average(a) + self.assertEqual(av_a.domain, []) + self.assertIsInstance(av_a, pybamm.Division) + self.assertIsInstance(av_a.children[0], pybamm.Integral) + self.assertEqual(av_a.children[0].integration_variable[0].domain, a.domain) + self.assertIsInstance(av_a.children[1], pybamm.Integral) + self.assertEqual(av_a.children[1].integration_variable[0].domain, a.domain) + self.assertEqual(av_a.children[1].children[0].id, pybamm.ones_like(a).id) def test_r_average(self): a = pybamm.Scalar(1) @@ -318,10 +323,6 @@ def test_r_average(self): # electrode domains go to current collector when averaged self.assertEqual(av_a.domain, []) - a = pybamm.Symbol("a", domain="bad domain") - with self.assertRaises(pybamm.DomainError): - pybamm.x_average(a) - def test_yz_average(self): a = pybamm.Scalar(1) z_average_a = pybamm.z_average(a) From 21d6ad38645dbd9ae2dd1524108e9991eca200a1 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 2 Mar 2020 13:33:55 -0500 Subject: [PATCH 09/10] #846 changelog --- CHANGELOG.md | 8 ++++++++ pybamm/solvers/casadi_algebraic_solver.py | 6 ++++-- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fd837ede8a..3cc9485bef 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,13 @@ # [Unreleased](https://github.com/pybamm-team/PyBaMM) +## Features + +- Added `Minimum`, `Maximum` and `Sign` operators + +## Bug fixes + +- Some bug fixes to generalize specifying models that aren't battery models, see [#846](https://github.com/pybamm-team/PyBaMM/issues/846) + ## Breaking changes - Removed "set external temperature" and "set external potential" options. Use "external submodels" option instead ([#862](https://github.com/pybamm-team/PyBaMM/pull/862)) diff --git a/pybamm/solvers/casadi_algebraic_solver.py b/pybamm/solvers/casadi_algebraic_solver.py index e766a4e33c..013499de51 100644 --- a/pybamm/solvers/casadi_algebraic_solver.py +++ b/pybamm/solvers/casadi_algebraic_solver.py @@ -62,13 +62,15 @@ def _integrate(self, model, t_eval, inputs=None): # Set up rootfinder roots = casadi.rootfinder( - "roots", "newton", dict(x=y_sym, p=t_u_sym, g=alg), {"abstol": self.tol}, + "roots", "newton", dict(x=y_sym, p=t_u_sym, g=alg), {"abstol": self.tol} ) for idx, t in enumerate(t_eval): # Evaluate algebraic with new t and previous y0, if it's already close # enough then keep it if np.all(abs(model.algebraic_eval(t, y0)) < self.tol): - pybamm.logger.debug("Keeping same solution at t={}".format(t)) + pybamm.logger.debug( + "Keeping same solution at t={}".format(t * model.timescale_eval) + ) y[:, idx] = y0 # Otherwise calculate new y0 else: From 3d7bc4abb443e6fded5d4d7537916fc9d9584e0d Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Tue, 10 Mar 2020 11:23:41 -0400 Subject: [PATCH 10/10] #846 coverage --- pybamm/solvers/base_solver.py | 1 - .../test_unary_operators.py | 8 ++++++++ tests/unit/test_solvers/test_base_solver.py | 17 +++++++++++++++++ 3 files changed, 25 insertions(+), 1 deletion(-) diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index fd93e1cbe9..f849d8b3c5 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -139,7 +139,6 @@ def set_up(self, model, inputs=None): self.root_method = None if ( isinstance(self, (pybamm.CasadiSolver, pybamm.CasadiAlgebraicSolver)) - or self.root_method == "casadi" ) and model.convert_to_format != "casadi": pybamm.logger.warning( "Converting {} to CasADi for solving with CasADi solver".format( diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py index 42edd7fedf..83ea1e81a0 100644 --- a/tests/unit/test_expression_tree/test_unary_operators.py +++ b/tests/unit/test_expression_tree/test_unary_operators.py @@ -5,6 +5,7 @@ import unittest import numpy as np +from scipy.sparse import diags class TestUnaryOperators(unittest.TestCase): @@ -43,6 +44,13 @@ def test_sign(self): signb = pybamm.sign(b) self.assertEqual(signb.evaluate(), -1) + A = diags(np.linspace(-1, 1, 5)) + b = pybamm.Matrix(A) + signb = pybamm.sign(b) + np.testing.assert_array_equal( + np.diag(signb.evaluate().toarray()), [-1, -1, 0, 1, 1] + ) + def test_gradient(self): a = pybamm.Symbol("a") grad = pybamm.Gradient(a) diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py index 7100d18b55..cd2004908f 100644 --- a/tests/unit/test_solvers/test_base_solver.py +++ b/tests/unit/test_solvers/test_base_solver.py @@ -176,6 +176,23 @@ def algebraic_eval(self, t, y): ): solver.calculate_consistent_state(Model()) + def test_convert_to_casadi_format(self): + # Make sure model is converted to casadi format + model = pybamm.BaseModel() + v = pybamm.Variable("v") + model.rhs = {v: -1} + model.initial_conditions = {v: 1} + model.convert_to_format = "python" + + disc = pybamm.Discretisation() + disc.process_model(model) + + solver = pybamm.BaseSolver() + pybamm.set_logging_level("ERROR") + solver.set_up(model, {}) + self.assertEqual(model.convert_to_format, "casadi") + pybamm.set_logging_level("WARNING") + if __name__ == "__main__": print("Add -v for more debug output")