diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py index 80ef5dc8fe..efbfe46772 100644 --- a/pybamm/discretisations/discretisation.py +++ b/pybamm/discretisations/discretisation.py @@ -514,6 +514,7 @@ def process_rhs_and_algebraic(self, model): equations) and processed_concatenated_algebraic """ + # Discretise right-hand sides, passing domain from variable processed_rhs = self.process_dict(model.rhs) @@ -996,6 +997,16 @@ def check_model(self, model): self.check_initial_conditions(model) self.check_initial_conditions_rhs(model) self.check_variables(model) + self.check_for_time_derivatives_in_rhs(model) + + def check_for_time_derivatives_in_rhs(self, model): + # Check that no variable time derivatives exist in the rhs equations + for eq in model.rhs.values(): + for node in eq.pre_order(): + if isinstance(node, pybamm.VariableDot): + raise pybamm.ModelError( + "time derivative of variable found ({}) in rhs equation" + ) def check_initial_conditions(self, model): """Check initial conditions are a numpy array""" diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py index beb322faa8..90ea3c41c0 100644 --- a/pybamm/expression_tree/binary_operators.py +++ b/pybamm/expression_tree/binary_operators.py @@ -162,21 +162,21 @@ def _binary_new_copy(self, left, right): "Default behaviour for new_copy" return self.__class__(left, right) - def evaluate(self, t=None, y=None, u=None, known_evals=None): + def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None): """ See :meth:`pybamm.Symbol.evaluate()`. """ if known_evals is not None: id = self.id try: return known_evals[id], known_evals except KeyError: - left, known_evals = self.left.evaluate(t, y, u, known_evals) - right, known_evals = self.right.evaluate(t, y, u, known_evals) + left, known_evals = self.left.evaluate(t, y, y_dot, u, known_evals) + right, known_evals = self.right.evaluate(t, y, y_dot, u, known_evals) value = self._binary_evaluate(left, right) known_evals[id] = value return value, known_evals else: - left = self.left.evaluate(t, y, u) - right = self.right.evaluate(t, y, u) + left = self.left.evaluate(t, y, y_dot, u) + right = self.right.evaluate(t, y, y_dot, u) return self._binary_evaluate(left, right) def _evaluate_for_shape(self): diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py index e1225bf03a..4a74de84cf 100644 --- a/pybamm/expression_tree/concatenations.py +++ b/pybamm/expression_tree/concatenations.py @@ -54,7 +54,7 @@ def _concatenation_evaluate(self, children_eval): else: return self.concatenation_function(children_eval) - def evaluate(self, t=None, y=None, u=None, known_evals=None): + def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None): """ See :meth:`pybamm.Symbol.evaluate()`. """ children = self.cached_children if known_evals is not None: @@ -62,14 +62,14 @@ def evaluate(self, t=None, y=None, u=None, known_evals=None): children_eval = [None] * len(children) for idx, child in enumerate(children): children_eval[idx], known_evals = child.evaluate( - t, y, u, known_evals + t, y, y_dot, u, known_evals ) known_evals[self.id] = self._concatenation_evaluate(children_eval) return known_evals[self.id], known_evals else: children_eval = [None] * len(children) for idx, child in enumerate(children): - children_eval[idx] = child.evaluate(t, y, u) + children_eval[idx] = child.evaluate(t, y, y_dot, u) return self._concatenation_evaluate(children_eval) def new_copy(self): diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py index 2601034b98..c21b473ef4 100644 --- a/pybamm/expression_tree/functions.py +++ b/pybamm/expression_tree/functions.py @@ -152,19 +152,19 @@ def _function_jac(self, children_jacs): return jacobian - def evaluate(self, t=None, y=None, u=None, known_evals=None): + def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None): """ See :meth:`pybamm.Symbol.evaluate()`. """ if known_evals is not None: if self.id not in known_evals: evaluated_children = [None] * len(self.children) for i, child in enumerate(self.children): evaluated_children[i], known_evals = child.evaluate( - t, y, u, known_evals=known_evals + t, y, y_dot, u, known_evals=known_evals ) known_evals[self.id] = self._function_evaluate(evaluated_children) return known_evals[self.id], known_evals else: - evaluated_children = [child.evaluate(t, y, u) for child in self.children] + evaluated_children = [child.evaluate(t, y, y_dot, u) for child in self.children] return self._function_evaluate(evaluated_children) def _evaluate_for_shape(self): diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py index ae3a1af9a7..a4c7f8574b 100644 --- a/pybamm/expression_tree/unary_operators.py +++ b/pybamm/expression_tree/unary_operators.py @@ -63,15 +63,15 @@ def _unary_evaluate(self, child): """Perform unary operation on a child. """ raise NotImplementedError - def evaluate(self, t=None, y=None, u=None, known_evals=None): + def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None): """ See :meth:`pybamm.Symbol.evaluate()`. """ if known_evals is not None: if self.id not in known_evals: - child, known_evals = self.child.evaluate(t, y, u, known_evals) + child, known_evals = self.child.evaluate(t, y, y_dot, u, known_evals) known_evals[self.id] = self._unary_evaluate(child) return known_evals[self.id], known_evals else: - child = self.child.evaluate(t, y, u) + child = self.child.evaluate(t, y, y_dot, u) return self._unary_evaluate(child) def _evaluate_for_shape(self): diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py index 17122a3c74..28df731ec9 100644 --- a/tests/unit/test_discretisations/test_discretisation.py +++ b/tests/unit/test_discretisations/test_discretisation.py @@ -701,6 +701,24 @@ def test_process_model_ode(self): with self.assertRaises(pybamm.ModelError): disc.process_model(model) + # test that ill possed model with time derivatives of variables in rhs raises an + # error + model = pybamm.BaseModel() + model.rhs = {c: pybamm.div(N) + pybamm.d_dt(c), T: pybamm.div(q), S: pybamm.div(p)} + model.initial_conditions = { + c: pybamm.Scalar(2), + T: pybamm.Scalar(5), + S: pybamm.Scalar(8), + } + model.boundary_conditions = { + c: {"left": (0, "Neumann"), "right": (0, "Neumann")}, + T: {"left": (0, "Neumann"), "right": (0, "Neumann")}, + S: {"left": (0, "Neumann"), "right": (0, "Neumann")}, + } + model.variables = {"ST": S * T} + with self.assertRaises(pybamm.ModelError): + disc.process_model(model) + def test_process_model_dae(self): # one rhs equation and one algebraic whole_cell = ["negative electrode", "separator", "positive electrode"]