Repository: tinosulzer/PyBaMM

Issue:
refactor BaseModel to hold an expression tree
See #16. 

The base model will internally hold an expression tree representing the rhs of an ode model. e.g.

```python
class BaseModel:
   def __init__(self):
         self._rhs = ... # lhs is a dictionary mapping variables to expressions
         self._initial_conditions = ... # a dictionary mapping expression to expressions representing the initial conditions
         self._boundary_conditions = ... # a dictionary mapping expression to expressions representing the boundary conditions
   @property
   def rhs(self):
          return self._rhs
   ...
```

Note: rather than having `rhs` return a dict, it might be nicer to have `BaseModel` emulate a dictionary, so you could write `model[variable]` and get the rhs expression associated with that variable (see https://docs.python.org/3/reference/datamodel.html?emulating-container-types#emulating-container-types)

Discussion:

Patch:
diff --git a/docs/source/expression_tree/binary_operator.rst b/docs/source/expression_tree/binary_operator.rst
index f56a2f25..61d0bb1c 100644
--- a/docs/source/expression_tree/binary_operator.rst
+++ b/docs/source/expression_tree/binary_operator.rst
@@ -7,9 +7,6 @@ Binary Operators
 .. autoclass:: pybamm.Addition
   :members:
 
-.. autoclass:: pybamm.Power
-  :members:
-
 .. autoclass:: pybamm.Subtraction
   :members:
 
diff --git a/docs/source/expression_tree/unary_operator.rst b/docs/source/expression_tree/unary_operator.rst
index d4aaa793..b1d69b08 100644
--- a/docs/source/expression_tree/unary_operator.rst
+++ b/docs/source/expression_tree/unary_operator.rst
@@ -4,12 +4,6 @@ Unary Operators
 .. autoclass:: pybamm.UnaryOperator
   :members:
 
-.. autoclass:: pybamm.Negate
-  :members:
-
-.. autoclass:: pybamm.AbsoluteValue
-  :members:
-
 .. autoclass:: pybamm.SpatialOperator
   :members:
 
@@ -17,7 +11,6 @@ Unary Operators
   :members:
 
 .. autoclass:: pybamm.Divergence
-  :members:
 
 .. autofunction:: pybamm.grad
 
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 52442596..9f248d9a 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -53,7 +53,6 @@ from .expression_tree.symbol import Symbol
 from .expression_tree.binary_operators import (
     BinaryOperator,
     Addition,
-    Power,
     Subtraction,
     Multiplication,
     Division,
@@ -65,8 +64,6 @@ from .expression_tree.domain import Domain
 from .expression_tree.parameter import Parameter
 from .expression_tree.unary_operators import (
     UnaryOperator,
-    Negate,
-    AbsoluteValue,
     SpatialOperator,
     Gradient,
     Divergence,
diff --git a/pybamm/discretisations/base_discretisation.py b/pybamm/discretisations/base_discretisation.py
index de878c94..6803c2f5 100644
--- a/pybamm/discretisations/base_discretisation.py
+++ b/pybamm/discretisations/base_discretisation.py
@@ -37,18 +37,23 @@ class BaseDiscretisation(object):
             Model to dicretise. Must have attributes rhs, initial_conditions and
             boundary_conditions (all dicts of {variable: equation})
 
+        Returns
+        -------
+        y0 : :class:`numpy.array`
+            Vector of initial conditions
+
         """
         # Set the y split for variables
         y_slices = self.get_variable_slices(model.rhs.keys())
 
         # Discretise and concatenate initial conditions, passing domain from variable
-        model.initial_conditions = self.process_initial_conditions(
-            model.initial_conditions
-        )
+        y0 = self.process_initial_conditions(model.initial_conditions)
 
         # Discretise right-hand sides, passing domain from variable
         model.rhs = self.process_rhs(model.rhs, model.boundary_conditions, y_slices)
 
+        return y0
+
     def get_variable_slices(self, variables):
         """Set the slicing for variables.
 
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index 890b01e7..2262a186 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -33,21 +33,6 @@ class BinaryOperator(pybamm.Symbol):
         return "{!s} {} {!s}".format(self.children[0], self.name, self.children[1])
 
 
-class Power(BinaryOperator):
-    """A node in the expression tree representing a `**` power operator
-
-    **Extends:** :class:`BinaryOperator`
-    """
-
-    def __init__(self, left, right):
-        """ See :meth:`pybamm.BinaryOperator.__init__()`. """
-        super().__init__("**", left, right)
-
-    def evaluate(self, t=None, y=None):
-        """ See :meth:`pybamm.Symbol.evaluate()`. """
-        return self.children[0].evaluate(t, y)**self.children[1].evaluate(t, y)
-
-
 class Addition(BinaryOperator):
     """A node in the expression tree representing an addition operator
 
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 79491358..f27baec2 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -80,9 +80,8 @@ class Symbol(anytree.NodeMixin):
         return self._name
 
     def __repr__(self):
-        """returns the string `__class__(id, name, parent expression)`"""
-        return "{!s}({}, {!s}, {!s})".format(self.__class__.__name__,
-                                             hex(self.id), self._name, self.parent)
+        """returns the string `Symbol(name, parent expression)`"""
+        return "Symbol({!s}, {!s})".format(self._name, self.parent)
 
     def __add__(self, other):
         """return an :class:`Addition` object"""
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 926cdc90..86aaf2e4 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -4,7 +4,6 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 import pybamm
-import numpy as np
 
 
 class UnaryOperator(pybamm.Symbol):
@@ -28,44 +27,6 @@ class UnaryOperator(pybamm.Symbol):
     def __init__(self, name, child):
         super().__init__(name, children=[child])
 
-    def __str__(self):
-        """ See :meth:`pybamm.Symbol.__str__()`. """
-        return "{}({!s})".format(self.name, self.children[0])
-
-
-class Negate(UnaryOperator):
-    """A node in the expression tree representing a `-` negation operator
-
-    **Extends:** :class:`UnaryOperator`
-    """
-
-    def __init__(self, child):
-        """ See :meth:`pybamm.UnaryOperator.__init__()`. """
-        super().__init__("-", child)
-
-    def evaluate(self, t=None, y=None):
-        """ See :meth:`pybamm.Symbol.evaluate()`. """
-        return -self.children[0].evaluate(t, y)
-
-    def __str__(self):
-        """ See :meth:`pybamm.Symbol.__str__()`. """
-        return "{}{!s}".format(self.name, self.children[0])
-
-
-class AbsoluteValue(UnaryOperator):
-    """A node in the expression tree representing an `abs` operator
-
-    **Extends:** :class:`UnaryOperator`
-    """
-
-    def __init__(self, child):
-        """ See :meth:`pybamm.UnaryOperator.__init__()`. """
-        super().__init__("abs", child)
-
-    def evaluate(self, t=None, y=None):
-        """ See :meth:`pybamm.Symbol.evaluate()`. """
-        return np.abs(self.children[0].evaluate(t, y))
-
 
 class SpatialOperator(UnaryOperator):
     """A node in the expression tree representing a unary spatial operator
@@ -92,6 +53,10 @@ class SpatialOperator(UnaryOperator):
         super().__init__(name, child)
         # self.domain = child.domain
 
+    def __str__(self):
+        """ See :meth:`pybamm.Symbol.__str__()`. """
+        return "{}({!s})".format(self.name, self.children[0])
+
 
 class Gradient(SpatialOperator):
     """A node in the expression tree representing an grad operator
diff --git a/pybamm/models/core.py b/pybamm/models/core.py
index 6ef77e20..b023c0ac 100644
--- a/pybamm/models/core.py
+++ b/pybamm/models/core.py
@@ -4,131 +4,52 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 
-import numpy as np
-
 
 class BaseModel(object):
     """Base model class for other models to extend.
 
-    Parameters
+    Attributes
     ----------
-    tests : dict, optional
-        A dictionary for testing the convergence of the numerical solution:
-            * {} (default): We are not running in test mode, use built-ins.
-            * {'inits': dict of initial conditions,
-               'bcs': dict of boundary conditions,
-               'sources': dict of source terms
-               }: To be used for testing convergence to an exact solution.
-    """
 
-    def __init__(self, tests={}):
-        self.name = "Base Model"
-        # Assign tests as an attribute
-        if tests:
-            assert set(tests.keys()) == {
-                "inits",
-                "bcs",
-                "sources",
-            }, "tests.keys() must include, 'inits', 'bcs' and 'sources'"
-        self.tests = tests
+    _rhs: dict
+        A dictionary that maps expressions (variables) to expressions that represent
+        the rhs
+    _initial_conditions: dict
+        A dictionary that maps expressions (variables) to expressions that represent
+        the initial conditions
+    _boundary_conditions: dict
+        A dictionary that maps expressions (variables) to expressions that represent
+        the boundary conditions
+    """
 
-    def __str__(self):
-        return self.name
+    def __init__(self):
+        self._rhs = {}
+        self._initial_conditions = {}
+        self._boundary_conditions = {}
 
     @property
-    def pde_variables(self):
-        """The variables of the model, as defined by submodels."""
-        return [
-            (variable, submodel.submesh)
-            for variable, submodel in self.submodels["pdes"].items()
-        ]
-
-    # @property
-    # def domains(self):
-    #     """The domain(s) in which the model holds."""
-    #     return set([domain for self.domain in self.variables])
-
-    def set_simulation(self, param, operators, mesh):
-        """
-        Assign simulation-specific objects as attributes.
+    def rhs(self):
+        return self._rhs
 
-        Parameters
-        ----------
-        param : :class:`pybamm.BaseParameterValues` instance
-            The parameters of the simulation
-        operators : :class:`pybamm.Operators` instance
-            The spatial operators.
-        mesh : :class:`pybamm.Mesh` instance
-            The spatial and temporal discretisation.
-        """
-        self.param = param
-        self.operators = operators
-        self.mesh = mesh
-
-        self.simulation_set = True
+    @rhs.setter
+    def rhs(self, rhs):
+        self._rhs = rhs
 
+    @property
     def initial_conditions(self):
-        """
-        Calculate the initial conditions for the simulation.
-
-        Returns
-        -------
-        array_like
-            A concatenated vector of all the initial conditions.
-
-        """
-        return np.concatenate(
-            [
-                submodel.initial_conditions()
-                for submodel in self.submodels["pdes"].values()
-            ]
-        )
+        return self._initial_conditions
 
-    def pdes_rhs(self, vars):
-        """
-        Calculate the spatial derivates of the spatial terms in the PDEs and returns
-        the right-hand side to be used by the ODE solver (Method of Lines).
+    @initial_conditions.setter
+    def initial_conditions(self, initial_conditions):
+        self._initial_conditions = initial_conditions
 
-        Parameters
-        ----------
-        vars : pybamm.variables.Variables() instance
-            The variables of the model.
-
-        Returns
-        -------
-        dydt : array_like
-            A concatenated vector of all the derivatives.
-
-        """
-        j = self.reactions(vars)
-        vars.set_reaction_vars({"j": j})
-        dydt = np.concatenate(
-            [submodel.pdes_rhs(vars) for submodel in self.submodels["pdes"].values()]
-        )
-
-        return dydt
-
-    def reactions(self, vars):
-        """
-        Calculate the interfacial current density using the reactions defined in
-        self.submodels.
-
-        Parameters
-        ----------
-        vars : pybamm.variables.Variables() instance
-            The variables of the model.
-
-        Returns
-        -------
-        array_like
-            The interfacial current density.
-
-        """
-        jn = self.submodels["reactions"]["neg"].reaction(vars.neg)
-        jp = self.submodels["reactions"]["pos"].reaction(vars.pos)
+    @property
+    def boundary_conditions(self):
+        return self._boundary_conditions
 
-        return np.concatenate([jn, np.zeros_like(self.mesh.xs.nodes), jp])
+    @boundary_conditions.setter
+    def boundary_conditions(self, boundary_conditions):
+        self._boundary_conditions = boundary_conditions
 
-    @property
-    def submodels(self):
-        raise NotImplementedError
+    def __getitem__(self, key):
+        return self.rhs[key]


Test patch:
diff --git a/tests/test_base_model.py b/tests/test_base_model.py
new file mode 100644
index 00000000..9948dcbb
--- /dev/null
+++ b/tests/test_base_model.py
@@ -0,0 +1,34 @@
+#
+# Tests for the base model class
+#
+import pybamm
+
+import unittest
+
+
+class TestBaseModel(unittest.TestCase):
+
+    def test_rhs_set_get(self):
+        model = pybamm.BaseModel()
+        rhs = {"c": "alpha", "d": "beta"}
+        model.rhs = rhs
+        self.assertEqual(rhs, model.rhs)
+
+    def test_initial_conditions_set_get(self):
+        model = pybamm.BaseModel()
+        initial_conditions = {"c0": "gamma", "d0": "delta"}
+        model.initial_conditions = initial_conditions
+        self.assertEqual(initial_conditions, model.initial_conditions)
+
+    def test_boundary_conditions_set_get(self):
+        model = pybamm.BaseModel()
+        boundary_conditions = {"c_left": "epsilon", "c_right": "eta"}
+        model.boundary_conditions = boundary_conditions
+        self.assertEqual(boundary_conditions, model.boundary_conditions)
+
+    def test_model_dict_behaviour(self):
+        model = pybamm.BaseModel()
+        rhs = {"c": "alpha", "d": "beta"}
+        model.rhs = rhs
+        self.assertEqual(model["c"], rhs["c"])
+        self.assertEqual(model["c"], model.rhs["c"])
diff --git a/tests/test_discretisations/test_base_discretisations.py b/tests/test_discretisations/test_base_discretisations.py
index 3a3811af..e3fbe3f0 100644
--- a/tests/test_discretisations/test_base_discretisations.py
+++ b/tests/test_discretisations/test_base_discretisations.py
@@ -194,8 +194,7 @@ class TestDiscretise(unittest.TestCase):
         mesh = MeshForTesting()
         disc = DiscretisationForTesting(mesh)
 
-        disc.process_model(model)
-        y0 = model.initial_conditions
+        y0 = disc.process_model(model)
         np.testing.assert_array_equal(y0, 3 * np.ones_like(mesh["whole cell"].nodes))
         np.testing.assert_array_equal(y0, model.rhs.evaluate(None, y0))
 
@@ -207,8 +206,7 @@ class TestDiscretise(unittest.TestCase):
         boundary_conditions = {}
         model = ModelForTesting(rhs, initial_conditions, boundary_conditions)
 
-        disc.process_model(model)
-        y0 = model.initial_conditions
+        y0 = disc.process_model(model)
         np.testing.assert_array_equal(
             y0,
             np.concatenate(
diff --git a/tests/test_expression_tree/test_binary_operators.py b/tests/test_expression_tree/test_binary_operators.py
index 716a1786..0d554ac9 100644
--- a/tests/test_expression_tree/test_binary_operators.py
+++ b/tests/test_expression_tree/test_binary_operators.py
@@ -23,19 +23,6 @@ class TestBinaryOperators(unittest.TestCase):
         self.assertEqual(sum.children[0].name, a.name)
         self.assertEqual(sum.children[1].name, b.name)
 
-    def test_power(self):
-        a = pybamm.Symbol("a")
-        b = pybamm.Symbol("b")
-        pow = pybamm.Power(a, b)
-        self.assertEqual(pow.name, "**")
-        self.assertEqual(pow.children[0].name, a.name)
-        self.assertEqual(pow.children[1].name, b.name)
-
-        a = pybamm.Scalar(4)
-        b = pybamm.Scalar(2)
-        pow = pybamm.Power(a, b)
-        self.assertEqual(pow.evaluate(), 16)
-
     def test_addition_printing(self):
         a = pybamm.Symbol("a")
         b = pybamm.Symbol("b")
diff --git a/tests/test_expression_tree/test_symbol.py b/tests/test_expression_tree/test_symbol.py
index ac9df67d..06ecde11 100644
--- a/tests/test_expression_tree/test_symbol.py
+++ b/tests/test_expression_tree/test_symbol.py
@@ -61,42 +61,6 @@ class TestSymbol(unittest.TestCase):
         with self.assertRaises(NotImplementedError):
             a.evaluate()
 
-    def test_symbol_repr(self):
-        """
-        test that __repr___ returns the string
-        `__class__(id, name, parent expression)`
-        """
-        a = pybamm.Symbol("a")
-        b = pybamm.Symbol("b")
-        hex_regex = r"\-?0x[0-9,a-f]+"
-        self.assertRegex(a.__repr__(), r"Symbol\(" + hex_regex + r", a, None\)")
-        self.assertRegex(b.__repr__(), r"Symbol\(" + hex_regex + r", b, None\)")
-        self.assertRegex(
-            (a + b).__repr__(), r"Addition\(" + hex_regex + r", \+, None\)"
-        )
-        self.assertRegex(
-            (a + b).children[0].__repr__(), r"Symbol\(" + hex_regex + r", a, a \+ b\)"
-        )
-        self.assertRegex(
-            (a + b).children[1].__repr__(), r"Symbol\(" + hex_regex + r", b, a \+ b\)"
-        )
-        self.assertRegex(
-            (a * b).__repr__(), r"Multiplication\(" + hex_regex + r", \*, None\)"
-        )
-        self.assertRegex(
-            (a * b).children[0].__repr__(), r"Symbol\(" + hex_regex + r", a, a \* b\)"
-        )
-        self.assertRegex(
-            (a * b).children[1].__repr__(), r"Symbol\(" + hex_regex + ", b, a \* b\)"
-        )
-        self.assertRegex(
-            pybamm.grad(a).__repr__(), r"Gradient\(" + hex_regex + ", grad, None\)"
-        )
-        self.assertRegex(
-            pybamm.grad(a).children[0].__repr__(),
-            r"Symbol\(" + hex_regex + ", a, grad\(a\)\)",
-        )
-
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/test_expression_tree/test_unary_operators.py b/tests/test_expression_tree/test_unary_operators.py
index 2f92ae6c..dac4877e 100644
--- a/tests/test_expression_tree/test_unary_operators.py
+++ b/tests/test_expression_tree/test_unary_operators.py
@@ -14,26 +14,6 @@ class TestUnaryOperators(unittest.TestCase):
         un = pybamm.UnaryOperator("unary test", a)
         self.assertEqual(un.children[0].name, a.name)
 
-    def test_negation(self):
-        a = pybamm.Symbol("a")
-        nega = pybamm.Negate(a)
-        self.assertEqual(nega.name, "-")
-        self.assertEqual(nega.children[0].name, a.name)
-
-        b = pybamm.Scalar(4)
-        negb = pybamm.Negate(b)
-        self.assertEqual(negb.evaluate(), -4)
-
-    def test_absolute(self):
-        a = pybamm.Symbol("a")
-        nega = pybamm.AbsoluteValue(a)
-        self.assertEqual(nega.name, "abs")
-        self.assertEqual(nega.children[0].name, a.name)
-
-        b = pybamm.Scalar(-4)
-        negb = pybamm.AbsoluteValue(b)
-        self.assertEqual(negb.evaluate(), 4)
-
     def test_gradient(self):
         a = pybamm.Symbol("a")
         grad = pybamm.Gradient(a)
diff --git a/tests/test_reaction_diffusion.py b/tests/test_reaction_diffusion.py
new file mode 100644
index 00000000..42564e69
--- /dev/null
+++ b/tests/test_reaction_diffusion.py
@@ -0,0 +1,98 @@
+#
+# Test for the reaction-diffusion model
+#
+import pybamm
+from tests.shared import pdes_io
+
+import unittest
+import numpy as np
+from numpy.linalg import norm
+import scipy.integrate as it
+
+
+@unittest.skip("not yet implemented")
+class TestReactionDiffusion(unittest.TestCase):
+    def setUp(self):
+        self.model = pybamm.ReactionDiffusionModel()
+        self.param = pybamm.Parameters()
+        target_npts = 10
+        tsteps = 10
+        tend = 1
+        self.mesh = pybamm.Mesh(self.param, target_npts, tsteps=tsteps, tend=tend)
+        self.param.set_mesh(self.mesh)
+
+    def tearDown(self):
+        del self.model
+        del self.param
+        del self.mesh
+
+    def test_model_shape(self):
+        for spatial_discretisation in pybamm.KNOWN_SPATIAL_DISCRETISATIONS:
+            solver = pybamm.Solver(spatial_discretisation=spatial_discretisation)
+            sim = pybamm.Simulation(self.model, solver=solver)
+            y, dydt = pdes_io(sim.model)
+            self.assertEqual(y.shape, dydt.shape)
+
+    def test_model_physics(self):
+        """Check that the average concentration is as expected"""
+        sim = pybamm.Simulation(self.model)
+        sim.run()
+        sim.average()
+
+        c_avg_expected = self.param.electrolyte.c0 + (
+            self.param.electrolyte.sn - self.param.electrolyte.sp
+        ) * it.cumtrapz(self.param.icell(sim.vars.t), sim.vars.t, initial=0.0)
+
+        np.testing.assert_allclose(sim.vars.c_avg, c_avg_expected, atol=4e-16)
+
+    def test_model_convergence(self):
+        """
+        Exact solution: c = exp(-4*pi**2*t) * cos(2*pi*x)
+        Initial conditions: c0 = cos(2*pi*x)
+        Boundary conditions: Zero flux
+        Source terms: 0
+
+        Can achieve "convergence" in time by changing the integrator tolerance
+        Can't get h**2 convergence in space
+        """
+        param = pybamm.Parameters(tests="convergence")
+        param.set_mesh(self.mesh)
+
+        def c_exact(t):
+            return np.exp(-4 * np.pi ** 2 * t) * np.cos(2 * np.pi * self.mesh.x.centres)
+
+        inits = {"concentration": c_exact(0)}
+
+        def bcs(t):
+            return {"concentration": (np.array([0]), np.array([0]))}
+
+        def sources(t):
+            return {"concentration": 0}
+
+        tests = {"inits": inits, "bcs": bcs, "sources": sources}
+
+        model = pybamm.ReactionDiffusionModel(tests=tests)
+
+        ns = [1, 2, 3]
+        errs = [0] * len(ns)
+        for i, n in enumerate(ns):
+            solver = pybamm.Solver(
+                integrator="BDF",
+                spatial_discretisation="Finite Volumes",
+                tol=10 ** (-n),
+            )
+            sim = pybamm.Simulation(model, param=param, mesh=self.mesh, solver=solver)
+            sim.run()
+            errs[i] = norm(
+                sim.vars.c.T - c_exact(self.mesh.time[:, np.newaxis])
+            ) / norm(c_exact(self.mesh.time[:, np.newaxis]))
+        [self.assertLess(errs[i + 1] / errs[i], 0.14) for i in range(len(errs) - 1)]
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
new file mode 100644
index 00000000..0f745817
--- /dev/null
+++ b/tests/test_simulation.py
@@ -0,0 +1,40 @@
+#
+# Test the simulation class
+#
+import pybamm
+
+import numpy as np
+
+import unittest
+
+
+@unittest.skip("not yet implemented")
+class TestSimulation(unittest.TestCase):
+    """Test the simulation class."""
+
+    def test_simulation_name(self):
+        model = pybamm.ReactionDiffusionModel()
+        param = pybamm.Parameters()
+        mesh = pybamm.Mesh(param, target_npts=50)
+        solver = pybamm.Solver()
+        sim = pybamm.Simulation(
+            model, param=param, mesh=mesh, solver=solver, name="test name"
+        )
+
+        self.assertEqual(sim.param.electrolyte.s.shape, sim.mesh.x.centres.shape)
+        np.testing.assert_array_equal(
+            sim.operators.x.div(mesh.x.edges), np.ones_like(mesh.x.centres)
+        )
+        self.assertEqual(param, sim.model.param)
+        self.assertEqual(mesh, sim.model.mesh)
+
+        self.assertEqual(str(sim), "test name")
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()