#858 add check for time derivs in model.rhs, fix some evaluate arg bugs

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"""

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):

pybamm/expression_tree/concatenations.py

@@ -54,22 +54,22 @@ 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:
             if self.id not in known_evals:
                 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):

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):

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):

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"]