Repository: pybamm-team/PyBaMM

Issue:
Add simulation class 
**Summary**
Create a simulation class which wraps model, discretization, parameters etc. 

**Motivation** 
Can be extended to allow linking of the experiment, model, discretization, and solver together and allow each to reference each other. Also simplifies the pybamm user interface so that users only need to see what the actually are setting


Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0c656afe1..7af4022ad 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added Simulation class ([#693](https://github.com/pybamm-team/PyBaMM/pull/693))
 -   Add interface to CasADi solver ([#687](https://github.com/pybamm-team/PyBaMM/pull/687))
 -   Add option to use CasADi's Algorithmic Differentiation framework to calculate Jacobians ([#687](https://github.com/pybamm-team/PyBaMM/pull/687))
 -   Add method to evaluate parameters more easily ([#669](https://github.com/pybamm-team/PyBaMM/pull/669))
diff --git a/README.md b/README.md
index 2946fc01d..8e4ca96df 100644
--- a/README.md
+++ b/README.md
@@ -11,7 +11,17 @@ Python Battery Mathematical Modelling solves continuum models for batteries, usi
 
 ## How do I use PyBaMM?
 
-PyBaMM comes with a number of [detailed examples](examples/notebooks/README.md), hosted here on
+The easiest way to use PyBaMM is to run a 1C constant-current discharge with a model of your choice with all the default settings:
+```python3
+import pybamm
+model = pybamm.lithium_ion.DFN() # Doyle-Fuller-Newman model
+sim = pybamm.Simulation(model)
+sim.solve()
+sim.plot()
+```
+However, much greater customisation is available. It is possible to change the physics, parameter values, geometry, submesh type,  number of submesh points, methods for spatial discretisation and solver for integration (see DFN [script](examples/scripts/DFN.py) or [notebook](examples/notebooks/models/dfn.ipynb)).
+
+Further details can be found in a number of [detailed examples](examples/notebooks/README.md), hosted here on
 github. In addition, there is a [full API documentation](http://pybamm.readthedocs.io/),
 hosted on [Read The Docs](readthedocs.io). A set of slides giving an overview of PyBaMM
 can be found
diff --git a/docs/index.rst b/docs/index.rst
index 49b6913f1..d991331d6 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -32,6 +32,7 @@ Contents
     source/solvers/index
     source/processed_variable
     source/util
+    source/simulation
 
 Examples
 ========
diff --git a/docs/source/simulation.rst b/docs/source/simulation.rst
new file mode 100644
index 000000000..a50a67eb2
--- /dev/null
+++ b/docs/source/simulation.rst
@@ -0,0 +1,5 @@
+Simulation
+==========
+
+.. autoclass:: pybamm.Simulation
+  :members:
\ No newline at end of file
diff --git a/examples/notebooks/models/compare-lithium-ion.ipynb b/examples/notebooks/models/compare-lithium-ion.ipynb
index 718009401..076b1eb46 100644
--- a/examples/notebooks/models/compare-lithium-ion.ipynb
+++ b/examples/notebooks/models/compare-lithium-ion.ipynb
@@ -463,7 +463,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,
diff --git a/examples/notebooks/parameter-values.ipynb b/examples/notebooks/parameter-values.ipynb
index 1a2e2b798..0ece13828 100644
--- a/examples/notebooks/parameter-values.ipynb
+++ b/examples/notebooks/parameter-values.ipynb
@@ -482,7 +482,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,
diff --git a/examples/scripts/run_simulation.py b/examples/scripts/run_simulation.py
new file mode 100644
index 000000000..09e2f3ed4
--- /dev/null
+++ b/examples/scripts/run_simulation.py
@@ -0,0 +1,7 @@
+import pybamm
+
+model = pybamm.lithium_ion.SPM()
+
+sim = pybamm.Simulation(model)
+sim.solve()
+sim.plot()
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 3a94c5dc9..69cda501f 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -268,7 +268,6 @@ from .solvers.scikits_ode_solver import ScikitsOdeSolver, have_scikits_odes
 from .solvers.scipy_solver import ScipySolver
 from .solvers.idaklu_solver import IDAKLUSolver, have_idaklu
 
-
 #
 # Current profiles
 #
@@ -285,6 +284,8 @@ from .parameters.standard_current_functions.get_current_data import GetCurrentDa
 from .processed_variable import post_process_variables, ProcessedVariable
 from .quick_plot import QuickPlot, ax_min, ax_max
 
+from .simulation import Simulation
+
 #
 # Remove any imported modules, so we don't expose them as part of pybamm
 #
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 284e375cb..0c1ab0c6a 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -122,8 +122,8 @@ class Discretisation(object):
             # since they point to the same object
             model_disc = model
         else:
-            # create a blank model so that original model is unchanged
-            model_disc = pybamm.BaseModel()
+            # create a model of the same class as the original model
+            model_disc = model.__class__(model.options)
             model_disc.name = model.name
             model_disc.options = model.options
             model_disc.use_jacobian = model.use_jacobian
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 5aff085d3..9880746c0 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -211,13 +211,13 @@ class ParameterValues(dict):
                 values["C-rate"] = float(values["Typical current [A]"]) / capacity
         return values
 
-    def process_model(self, model, processing="process"):
+    def process_model(self, unprocessed_model, processing="process", inplace=True):
         """Assign parameter values to a model.
         Currently inplace, could be changed to return a new model.
 
         Parameters
         ----------
-        model : :class:`pybamm.BaseModel`
+        unprocessed_model : :class:`pybamm.BaseModel`
             Model to assign parameter values for
         processing : str, optional
             Flag to indicate how to process model (default 'process')
@@ -226,6 +226,9 @@ class ParameterValues(dict):
             and replace any Parameter with a Value)
             * 'update': Calls :meth:`update_scalars()` for use on already-processed \
             model (update the value of any Scalars in the expression tree.)
+        inplace: bool, optional
+            If True, replace the parameters in the model in place. Otherwise, return a
+            new model with parameter values set. Default is True.
 
         Raises
         ------
@@ -233,9 +236,25 @@ class ParameterValues(dict):
             If an empty model is passed (`model.rhs = {}` and `model.algebraic={}`)
 
         """
-        pybamm.logger.info("Start setting parameters for {}".format(model.name))
-
-        if len(model.rhs) == 0 and len(model.algebraic) == 0:
+        pybamm.logger.info(
+            "Start setting parameters for {}".format(unprocessed_model.name)
+        )
+
+        # set up inplace vs not inplace
+        if inplace:
+            # any changes to model_disc attributes will change model attributes
+            # since they point to the same object
+            model = unprocessed_model
+        else:
+            # create a blank model of the same class
+            model = unprocessed_model.__class__(unprocessed_model.options)
+            model.name = unprocessed_model.name
+            model.options = unprocessed_model.options
+            model.use_jacobian = unprocessed_model.use_jacobian
+            model.use_simplify = unprocessed_model.use_simplify
+            model.convert_to_format = unprocessed_model.convert_to_format
+
+        if len(unprocessed_model.rhs) == 0 and len(unprocessed_model.algebraic) == 0:
             raise pybamm.ModelError("Cannot process parameters for empty model")
 
         if processing == "process":
@@ -243,13 +262,13 @@ class ParameterValues(dict):
         elif processing == "update":
             processing_function = self.update_scalars
 
-        for variable, equation in model.rhs.items():
+        for variable, equation in unprocessed_model.rhs.items():
             pybamm.logger.debug(
                 "{} parameters for {!r} (rhs)".format(processing.capitalize(), variable)
             )
             model.rhs[variable] = processing_function(equation)
 
-        for variable, equation in model.algebraic.items():
+        for variable, equation in unprocessed_model.algebraic.items():
             pybamm.logger.debug(
                 "{} parameters for {!r} (algebraic)".format(
                     processing.capitalize(), variable
@@ -257,7 +276,7 @@ class ParameterValues(dict):
             )
             model.algebraic[variable] = processing_function(equation)
 
-        for variable, equation in model.initial_conditions.items():
+        for variable, equation in unprocessed_model.initial_conditions.items():
             pybamm.logger.debug(
                 "{} parameters for {!r} (initial conditions)".format(
                     processing.capitalize(), variable
@@ -270,7 +289,7 @@ class ParameterValues(dict):
         # small number of variables, e.g. {"negative tab": neg. tab bc,
         # "positive tab": pos. tab bc "no tab": no tab bc}.
         new_boundary_conditions = {}
-        for variable, bcs in model.boundary_conditions.items():
+        for variable, bcs in unprocessed_model.boundary_conditions.items():
             processed_variable = processing_function(variable)
             new_boundary_conditions[processed_variable] = {}
             for side in ["left", "right", "negative tab", "positive tab", "no tab"]:
@@ -288,14 +307,14 @@ class ParameterValues(dict):
 
         model.boundary_conditions = new_boundary_conditions
 
-        for variable, equation in model.variables.items():
+        for variable, equation in unprocessed_model.variables.items():
             pybamm.logger.debug(
                 "{} parameters for {!r} (variables)".format(
                     processing.capitalize(), variable
                 )
             )
             model.variables[variable] = processing_function(equation)
-        for event, equation in model.events.items():
+        for event, equation in unprocessed_model.events.items():
             pybamm.logger.debug(
                 "{} parameters for event '{}''".format(processing.capitalize(), event)
             )
@@ -303,6 +322,8 @@ class ParameterValues(dict):
 
         pybamm.logger.info("Finish setting parameters for {}".format(model.name))
 
+        return model
+
     def update_model(self, model, disc):
         """Process a discretised model.
         Currently inplace, could be changed to return a new model.
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
new file mode 100644
index 000000000..774a5750b
--- /dev/null
+++ b/pybamm/simulation.py
@@ -0,0 +1,279 @@
+import pybamm
+import numpy as np
+import copy
+
+
+class Simulation:
+    """A Simulation class for easy building and running of PyBaMM simulations.
+
+    Parameters
+    ----------
+    model : :class:`pybamm.BaseModel`
+        The model to be simulated
+    """
+
+    def __init__(
+        self,
+        model,
+        geometry=None,
+        parameter_values=None,
+        submesh_types=None,
+        var_pts=None,
+        spatial_methods=None,
+        solver=None,
+        quick_plot_vars=None,
+    ):
+        self.model = model
+
+        self.geometry = geometry or model.default_geometry
+        self._parameter_values = parameter_values or model.default_parameter_values
+        self._submesh_types = submesh_types or model.default_submesh_types
+        self._var_pts = var_pts or model.default_var_pts
+        self._spatial_methods = spatial_methods or model.default_spatial_methods
+        self._solver = solver or self._model.default_solver
+        self._quick_plot_vars = quick_plot_vars
+
+        self.reset()
+
+    def set_defaults(self):
+        """
+        A method to set all the simulation specs to default values for the
+        supplied model.
+        """
+        self.geometry = self._model.default_geometry
+        self._parameter_values = self._model.default_parameter_values
+        self._submesh_types = self._model.default_submesh_types
+        self._var_pts = self._model.default_var_pts
+        self._spatial_methods = self._model.default_spatial_methods
+        self._solver = self._model.default_solver
+        self._quick_plot_vars = None
+
+    def reset(self):
+        """
+        A method to reset a simulation back to its unprocessed state.
+        """
+        self.model = self._model_class(self._model_options)
+        self.geometry = copy.deepcopy(self._unprocessed_geometry)
+        self._model_with_set_params = None
+        self._built_model = None
+        self._mesh = None
+        self._disc = None
+        self._solution = None
+
+    def set_parameters(self):
+        """
+        A method to set the parameters in the model and the associated geometry. If
+        the model has already been built or solved then this will first reset to the
+        unprocessed state and then set the parameter values.
+        """
+
+        if self.model_with_set_params:
+            return None
+
+        self._model_with_set_params = self._parameter_values.process_model(
+            self._model, inplace=True
+        )
+        self._parameter_values.process_geometry(self._geometry)
+
+    def build(self):
+        """
+        A method to build the model into a system of matrices and vectors suitable for
+        performing numerical computations. If the model has already been built or
+        solved then this function will have no effect. If you want to rebuild,
+        first use "reset()". This method will automatically set the parameters
+        if they have not already been set.
+        """
+
+        if self.built_model:
+            return None
+
+        self.set_parameters()
+        self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts)
+        self._disc = pybamm.Discretisation(self._mesh, self._spatial_methods)
+        self._built_model = self._disc.process_model(self._model, inplace=False)
+
+    def solve(self, t_eval=None, solver=None):
+        """
+        A method to solve the model. This method will automatically build
+        and set the model parameters if not already done so.
+
+        Parameters
+        ----------
+        t_eval : numeric type (optional)
+            The times at which to compute the solution
+        solver : :class:`pybamm.BaseSolver`
+            The solver to use to solve the model.
+        """
+        self.build()
+
+        if t_eval is None:
+            t_eval = np.linspace(0, 1, 100)
+
+        if solver is None:
+            solver = self.solver
+
+        self._solution = solver.solve(self.built_model, t_eval)
+
+    def plot(self, quick_plot_vars=None):
+        """
+        A method to quickly plot the outputs of the simulation.
+
+        Parameters
+        ----------
+        quick_plot_vars: list
+            A list of the variables to plot.
+        """
+
+        if self._solution is None:
+            raise ValueError(
+                "Model has not been solved, please solve the model before plotting."
+            )
+
+        if quick_plot_vars is None:
+            quick_plot_vars = self.quick_plot_vars
+
+        plot = pybamm.QuickPlot(
+            self.built_model,
+            self._mesh,
+            self._solution,
+            output_variables=quick_plot_vars,
+        )
+        plot.dynamic_plot()
+
+    @property
+    def model(self):
+        return self._model
+
+    @model.setter
+    def model(self, model):
+        self._model = model
+        self._model_class = model.__class__
+        self._model_options = model.options
+
+    @property
+    def model_with_set_params(self):
+        return self._model_with_set_params
+
+    @property
+    def built_model(self):
+        return self._built_model
+
+    @property
+    def model_options(self):
+        return self._model_options
+
+    @property
+    def geometry(self):
+        return self._geometry
+
+    @geometry.setter
+    def geometry(self, geometry):
+        self._geometry = geometry
+        self._unprocessed_geometry = copy.deepcopy(geometry)
+
+    @property
+    def unprocessed_geometry(self):
+        return self._unprocessed_geometry
+
+    @property
+    def parameter_values(self):
+        return self._parameter_values
+
+    @property
+    def submesh_types(self):
+        return self._submesh_types
+
+    @property
+    def var_pts(self):
+        return self._var_pts
+
+    @property
+    def spatial_methods(self):
+        return self._spatial_methods
+
+    @property
+    def solver(self):
+        return self._solver
+
+    @solver.setter
+    def solver(self, solver):
+        self._solver = solver
+
+    @property
+    def quick_plot_vars(self):
+        return self._quick_plot_vars
+
+    @quick_plot_vars.setter
+    def quick_plot_vars(self, quick_plot_vars):
+        self._quick_plot_vars = quick_plot_vars
+
+    @property
+    def solution(self):
+        return self._solution
+
+    def specs(
+        self,
+        model_options=None,
+        geometry=None,
+        parameter_values=None,
+        submesh_types=None,
+        var_pts=None,
+        spatial_methods=None,
+        solver=None,
+        quick_plot_vars=None,
+    ):
+        """
+        A method to set the various specs of the simulation. This method
+        automatically resets the model after the new specs have been set.
+
+        Parameters
+        ----------
+        model_options: dict (optional)
+            A dictionary of options to tweak the model you are using
+        geometry: :class:`pybamm.Geometry` (optional)
+            The geometry upon which to solve the model
+        parameter_values: dict (optional)
+            A dictionary of parameters and their corresponding numerical
+            values
+        submesh_types: dict (optional)
+            A dictionary of the types of submesh to use on each subdomain
+        var_pts: dict (optional)
+            A dictionary of the number of points used by each spatial
+            variable
+        spatial_methods: dict (optional)
+            A dictionary of the types of spatial method to use on each
+            domain (e.g. pybamm.FiniteVolume)
+        solver: :class:`pybamm.BaseSolver`
+            The solver to use to solve the model.
+        quick_plot_vars: list
+            A list of variables to plot automatically
+        """
+
+        if model_options:
+            self._model_options = model_options
+
+        if geometry:
+            self.geometry = geometry
+
+        if parameter_values:
+            self._parameter_values = parameter_values
+        if submesh_types:
+            self._submesh_types = submesh_types
+        if var_pts:
+            self._var_pts = var_pts
+        if spatial_methods:
+            self._spatial_methods = spatial_methods
+        if solver:
+            self._solver = solver
+        if quick_plot_vars:
+            self._quick_plot_vars = quick_plot_vars
+
+        if (
+            model_options
+            or geometry
+            or parameter_values
+            or submesh_types
+            or var_pts
+            or spatial_methods
+        ):
+            self.reset()


Test patch:
diff --git a/tests/unit/test_parameters/test_update_parameters.py b/tests/unit/test_parameters/test_update_parameters.py
index ccc13a3c7..bab24a8f4 100644
--- a/tests/unit/test_parameters/test_update_parameters.py
+++ b/tests/unit/test_parameters/test_update_parameters.py
@@ -94,6 +94,17 @@ class TestUpdateParameters(unittest.TestCase):
         # results should be different
         self.assertNotEqual(np.linalg.norm(Y1 - Y3), 0)
 
+    def test_inplace(self):
+        model = pybamm.lithium_ion.SPM()
+        param = model.default_parameter_values
+        new_model = param.process_model(model, inplace=False)
+
+        for val in list(model.rhs.values()):
+            self.assertTrue(val.has_symbol_of_classes(pybamm.Parameter))
+
+        for val in list(new_model.rhs.values()):
+            self.assertFalse(val.has_symbol_of_classes(pybamm.Parameter))
+
     def test_update_geometry(self):
         # test on simple lead-acid model
         model1 = pybamm.lead_acid.LOQS()
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
new file mode 100644
index 000000000..8dd3aa1b3
--- /dev/null
+++ b/tests/unit/test_simulation.py
@@ -0,0 +1,140 @@
+import pybamm
+import unittest
+
+
+class TestSimulation(unittest.TestCase):
+    def test_basic_ops(self):
+
+        model = pybamm.lithium_ion.SPM()
+        sim = pybamm.Simulation(model)
+
+        self.assertEqual(model.__class__, sim._model_class)
+        self.assertEqual(model.options, sim._model_options)
+
+        # check that the model is unprocessed
+        self.assertEqual(sim._mesh, None)
+        self.assertEqual(sim._disc, None)
+        for val in list(sim.model.rhs.values()):
+            self.assertTrue(val.has_symbol_of_classes(pybamm.Parameter))
+            self.assertFalse(val.has_symbol_of_classes(pybamm.Matrix))
+
+        sim.set_parameters()
+        self.assertEqual(sim._mesh, None)
+        self.assertEqual(sim._disc, None)
+        for val in list(sim.model_with_set_params.rhs.values()):
+            self.assertFalse(val.has_symbol_of_classes(pybamm.Parameter))
+            self.assertFalse(val.has_symbol_of_classes(pybamm.Matrix))
+
+        sim.build()
+        self.assertFalse(sim._mesh is None)
+        self.assertFalse(sim._disc is None)
+        for val in list(sim.built_model.rhs.values()):
+            self.assertFalse(val.has_symbol_of_classes(pybamm.Parameter))
+            self.assertTrue(val.has_symbol_of_classes(pybamm.Matrix))
+
+        sim.reset()
+        sim.set_parameters()
+        self.assertEqual(sim._mesh, None)
+        self.assertEqual(sim._disc, None)
+        self.assertEqual(sim.built_model, None)
+
+        for val in list(sim.model_with_set_params.rhs.values()):
+            self.assertFalse(val.has_symbol_of_classes(pybamm.Parameter))
+            self.assertFalse(val.has_symbol_of_classes(pybamm.Matrix))
+
+        sim.build()
+        sim.reset()
+        self.assertEqual(sim._mesh, None)
+        self.assertEqual(sim._disc, None)
+        self.assertEqual(sim.model_with_set_params, None)
+        self.assertEqual(sim.built_model, None)
+        for val in list(sim.model.rhs.values()):
+            self.assertTrue(val.has_symbol_of_classes(pybamm.Parameter))
+            self.assertFalse(val.has_symbol_of_classes(pybamm.Matrix))
+
+    def test_solve(self):
+
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
+        sim.solve()
+        self.assertFalse(sim._solution is None)
+        for val in list(sim.built_model.rhs.values()):
+            self.assertFalse(val.has_symbol_of_classes(pybamm.Parameter))
+            self.assertTrue(val.has_symbol_of_classes(pybamm.Matrix))
+
+        sim.reset()
+        self.assertEqual(sim.model_with_set_params, None)
+        self.assertEqual(sim.built_model, None)
+        for val in list(sim.model.rhs.values()):
+            self.assertTrue(val.has_symbol_of_classes(pybamm.Parameter))
+            self.assertFalse(val.has_symbol_of_classes(pybamm.Matrix))
+
+        self.assertEqual(sim._solution, None)
+
+    def test_reuse_commands(self):
+
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
+
+        sim.set_parameters()
+        sim.set_parameters()
+
+        sim.build()
+        sim.build()
+
+        sim.solve()
+        sim.solve()
+
+        sim.build()
+        sim.solve()
+        sim.set_parameters()
+
+    def test_specs(self):
+        # test can rebuild after setting specs
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
+        sim.build()
+
+        model_options = {"thermal": "lumped"}
+        sim.specs(model_options=model_options)
+        sim.build()
+        self.assertEqual(sim.model.options["thermal"], "lumped")
+
+        params = sim.parameter_values
+        # normally is 0.0001
+        params.update({"Negative electrode thickness [m]": 0.0002})
+        sim.specs(parameter_values=params)
+
+        self.assertEqual(
+            sim.parameter_values["Negative electrode thickness [m]"], 0.0002
+        )
+        sim.build()
+
+        geometry = sim.unprocessed_geometry
+        custom_geometry = {}
+        x_n = pybamm.standard_spatial_vars.x_n
+        custom_geometry["negative electrode"] = {
+            "primary": {
+                x_n: {"min": pybamm.Scalar(0), "max": pybamm.geometric_parameters.l_n}
+            }
+        }
+        geometry.update(custom_geometry)
+        sim.specs(geometry=geometry)
+        sim.build()
+
+        var_pts = sim.var_pts
+        var_pts[pybamm.standard_spatial_vars.x_n] = 5
+        sim.specs(var_pts=var_pts)
+        sim.build()
+
+        spatial_methods = sim.spatial_methods
+        # nothing to change this to at the moment but just reload in
+        sim.specs(spatial_methods=spatial_methods)
+        sim.build()
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
+