Repository: tinosulzer/PyBaMM

Issue:
Model well-posedness
When checking that a PDE model is well-posed, make sure that all variables that appear in the rhs dict values expression trees (i.e. the model equations) also appear in the rhs keys. 
For DAE models, give some "free passes" corresponding to the number of algebraic equations (for which we don't explicitly define a variable)

Most submodels won't be well-posed, but they don't need to be

Discussion:

Patch:
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 95eefdba..b7488a6d 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -195,32 +195,57 @@ class BaseModel(object):
     def check_well_posedness(self):
         """
         Check that the model is well-posed by executing the following tests:
+        - All the variables that appear in the rhs and algebraic equations appear in
+        the rhs keys, with leeway for exactly n variables where n is the number of
+        algebraic equations. Overdetermined if more equations than variables,
+        underdetermined if more variables than equations.
         - There is an initial condition in self.initial_conditions for each
         variable/equation pair in self.rhs
         - There are appropriate boundary conditions in self.boundary_conditions for each
         variable/equation pair in self.rhs
         """
+        # Equations (differential and algebraic)
+        # Get all the variables from differential and algebraic equations
+        variable_ids_in_eqns = set()
+        for eqn in list(self.rhs.values()) + self.algebraic:
+            variable_ids_in_eqns.update(
+                [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
+            )
+        # Get all variables ids from rhs keys
+        variable_ids_in_keys = set([var.id for var in self.rhs.keys()])
+        # Compare eqns and keys
+        unaccounted_variables = variable_ids_in_eqns.difference(variable_ids_in_keys)
+        # Count how many variables are not accounted for by algebraic equations
+        n_extra_variables = len(unaccounted_variables) - len(self.algebraic)
+        # Make sure n_extra_variables == 0
+        if n_extra_variables > 0:
+            raise pybamm.ModelError("model is underdetermined")
+        elif n_extra_variables < 0:
+            raise pybamm.ModelError("model is overdetermined")
+
         # Initial conditions
         for var in self.rhs.keys():
-            assert var in self.initial_conditions.keys(), pybamm.ModelError(
-                """no initial condition given for variable '{}'""".format(var)
-            )
+            if var not in self.initial_conditions.keys():
+                raise pybamm.ModelError(
+                    """no initial condition given for variable '{}'""".format(var)
+                )
 
         # Boundary conditions
         for var, eqn in self.rhs.items():
             if eqn.has_spatial_derivatives():
                 # Variable must be in at least one expression in the boundary condition
                 # keys (to account for both Dirichlet and Neumann boundary conditions)
-                assert any(
+                if not any(
                     [
                         any([var.id == symbol.id for symbol in key.pre_order()])
                         for key in self.boundary_conditions.keys()
                     ]
-                ), pybamm.ModelError(
-                    """
-                    no boundary condition given for variable '{}'
-                    with equation '{}'
-                    """.format(
-                        var, eqn
+                ):
+                    raise pybamm.ModelError(
+                        """
+                        no boundary condition given for variable '{}'
+                        with equation '{}'
+                        """.format(
+                            var, eqn
+                        )
                     )
-                )


Test patch:
diff --git a/tests/test_models/test_base_model.py b/tests/test_models/test_base_model.py
index e7aa7555..c175476f 100644
--- a/tests/test_models/test_base_model.py
+++ b/tests/test_models/test_base_model.py
@@ -157,7 +157,45 @@ class TestBaseModel(unittest.TestCase):
         self.assertEqual(model.initial_conditions[e], submodel2.initial_conditions[e])
         self.assertEqual(model.boundary_conditions[e], submodel2.boundary_conditions[e])
 
-    def test_check_well_posedness(self):
+    def test_check_well_posedness_variables(self):
+        # Well-posed ODE model
+        model = pybamm.BaseModel()
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        c = pybamm.Variable("c", domain=whole_cell)
+        d = pybamm.Variable("d", domain=whole_cell)
+        model.rhs = {c: 5 * pybamm.div(pybamm.grad(d)) - 1, d: -c}
+        model.initial_conditions = {c: 1, d: 2}
+        model.boundary_conditions = {
+            c: {"left": 0, "right": 0},
+            d: {"left": 0, "right": 0},
+        }
+        model.check_well_posedness()
+
+        # Well-posed DAE model
+        e = pybamm.Variable("e", domain=whole_cell)
+        model.algebraic = [e - c - d]
+        model.check_well_posedness()
+
+        # Underdetermined model - not enough differential equations
+        model.rhs = {c: 5 * pybamm.div(pybamm.grad(d)) - 1}
+        model.algebraic = [e - c - d]
+        with self.assertRaisesRegex(pybamm.ModelError, "underdetermined"):
+            model.check_well_posedness()
+
+        # Underdetermined model - not enough algebraic equations
+        model.rhs = {c: 5 * pybamm.div(pybamm.grad(d)) - 1, d: -c - e}
+        model.algebraic = []
+        with self.assertRaisesRegex(pybamm.ModelError, "underdetermined"):
+            model.check_well_posedness()
+
+        # Overdetermined model - too many algebraic equations
+        # Model cannot be overdetermined by too many differential equations, as rhs
+        # keys are all unique variables
+        model.algebraic = [c - d, e + d]
+        with self.assertRaisesRegex(pybamm.ModelError, "overdetermined"):
+            model.check_well_posedness()
+
+    def test_check_well_posedness_initial_boundary_conditions(self):
         # Well-posed model - Dirichlet
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         model = pybamm.BaseModel()
@@ -174,28 +212,22 @@ class TestBaseModel(unittest.TestCase):
         # Model with bad initial conditions (expect assertion error)
         d = pybamm.Variable("d", domain=whole_cell)
         model.initial_conditions = {d: 3}
-        with self.assertRaises(AssertionError) as error:
+        with self.assertRaisesRegex(pybamm.ModelError, "initial condition"):
             model.check_well_posedness()
-        self.assertIsInstance(error.exception.args[0], pybamm.ModelError)
-        self.assertIn("initial condition", error.exception.args[0].args[0])
 
         # Model with bad boundary conditions - Dirichlet (expect assertion error)
         d = pybamm.Variable("d", domain=whole_cell)
         model.initial_conditions = {c: 3}
         model.boundary_conditions = {d: {"left": 0, "right": 0}}
-        with self.assertRaises(AssertionError) as error:
+        with self.assertRaisesRegex(pybamm.ModelError, "boundary condition"):
             model.check_well_posedness()
-        self.assertIsInstance(error.exception.args[0], pybamm.ModelError)
-        self.assertIn("boundary condition", error.exception.args[0].args[0])
 
         # Model with bad boundary conditions - Neumann (expect assertion error)
         d = pybamm.Variable("d", domain=whole_cell)
         model.initial_conditions = {c: 3}
         model.boundary_conditions = {4 * pybamm.grad(d): {"left": 0, "right": 0}}
-        with self.assertRaises(AssertionError) as error:
+        with self.assertRaisesRegex(pybamm.ModelError, "boundary condition"):
             model.check_well_posedness()
-        self.assertIsInstance(error.exception.args[0], pybamm.ModelError)
-        self.assertIn("boundary condition", error.exception.args[0].args[0])
 
 
 if __name__ == "__main__":