From 716f8ea6226b754eefff664e0d11d98d487f1ee0 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Fri, 21 Feb 2020 17:02:03 -0500 Subject: [PATCH] #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()