Add simulation save and load

.github/PULL_REQUEST_TEMPLATE.md

@@ -6,7 +6,7 @@ Fixes # (issue)
 
 ## Type of change
 
-Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/master/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s.
+Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/master/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.
 
 - [ ] New feature (non-breaking change which adds functionality)
 - [ ] Optimization (back-end change that speeds up the code)

CHANGELOG.md

@@ -7,7 +7,7 @@
 -   Simplified the interface for setting current functions ([#723](https://github.com/pybamm-team/PyBaMM/pull/723))
 -   Added Heaviside operator ([#723](https://github.com/pybamm-team/PyBaMM/pull/723))
 -   Added some "Getting Started" documentation ([#703](https://github.com/pybamm-team/PyBaMM/pull/703))
--   Added Simulation class ([#693](https://github.com/pybamm-team/PyBaMM/pull/693))
+-   Added Simulation class ([#693](https://github.com/pybamm-team/PyBaMM/pull/693)) with load/save functionality ([#732](https://github.com/pybamm-team/PyBaMM/pull/732))
 -   Added interface to CasADi solver ([#687](https://github.com/pybamm-team/PyBaMM/pull/687), [#691](https://github.com/pybamm-team/PyBaMM/pull/691), [#714](https://github.com/pybamm-team/PyBaMM/pull/714)). This makes the SUNDIALS DAE solvers (Scikits and KLU) truly optional (though IDA KLU is recommended for solving the DFN).
 -   Added option to use CasADi's Algorithmic Differentiation framework to calculate Jacobians ([#687](https://github.com/pybamm-team/PyBaMM/pull/687))
 -   Added method to evaluate parameters more easily ([#669](https://github.com/pybamm-team/PyBaMM/pull/669))
@@ -36,6 +36,13 @@
 -   Added warning if `ProcessedVariable` is called outside its interpolation range ([#681](https://github.com/pybamm-team/PyBaMM/pull/681))
 -   Improved the way `ProcessedVariable` objects are created in higher dimensions ([#581](https://github.com/pybamm-team/PyBaMM/pull/581))
 
+## Breaking changes
+
+-   The parameters "Bruggeman coefficient" must now be specified separately as "Bruggeman coefficient (electrolyte)" and "Bruggeman coefficient (electrode)"
+-   The current classes (`GetConstantCurrent`, `GetUserCurrent` and `GetUserData`) have now been removed. Please refer to the [`change-input-current` notebook](https://github.com/pybamm-team/PyBaMM/blob/master/examples/notebooks/change-input-current.ipynb) for information on how to specify an input current
+-   Parameter functions must now use pybamm functions instead of numpy functions (e.g. `pybamm.exp` instead of `numpy.exp`), as these are then used to construct the expression tree directly. Generally, pybamm syntax follows numpy syntax; please get in touch if a function you need is missing.
+
+
 # [v0.1.0](https://github.com/pybamm-team/PyBaMM/tree/v0.1.0) - 2019-10-08
 
 This is the first official version of PyBaMM.

pybamm/__init__.py

@@ -246,7 +246,7 @@ def version(formatted=False):
 from .processed_variable import post_process_variables, ProcessedVariable
 from .quick_plot import QuickPlot, ax_min, ax_max
 
-from .simulation import Simulation
+from .simulation import Simulation, load_sim
 
 #
 # Remove any imported modules, so we don't expose them as part of pybamm

pybamm/discretisations/discretisation.py

@@ -137,13 +137,8 @@ def process_model(self, model, inplace=True):
             # since they point to the same object
             model_disc = model
         else:
-            # 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
-            model_disc.use_simplify = model.use_simplify
-            model_disc.convert_to_format = model.convert_to_format
+            # create an empty copy of the original model
+            model_disc = model.new_copy()
 
         model_disc.bcs = self.bcs
 

pybamm/expression_tree/binary_operators.py

@@ -227,7 +227,9 @@ def _binary_jac(self, left_jac, right_jac):
 
     def _binary_evaluate(self, left, right):
         """ See :meth:`pybamm.BinaryOperator._binary_evaluate()`. """
-        return left ** right
+        # don't raise RuntimeWarning for NaNs
+        with np.errstate(invalid="ignore"):
+            return left ** right
 
     def _binary_simplify(self, left, right):
         """ See :meth:`pybamm.BinaryOperator._binary_simplify()`. """

pybamm/models/base_model.py

@@ -1,12 +1,22 @@
 #
 # Base model class
 #
-import pybamm
-
+import inspect
 import numbers
+import pybamm
 import warnings
 
 
+class ParamClass:
+    """Class for converting a module of parameters into a class. For pickling."""
+
+    def __init__(self, methods):
+        for k, v in methods.__dict__.items():
+            # don't save module attributes (e.g. pybamm, numpy)
+            if not (k.startswith("__") or inspect.ismodule(v)):
+                self.__dict__[k] = v
+
+
 class BaseModel(object):
     """Base model class for other models to extend.
 
@@ -263,8 +273,14 @@ def jacobian_algebraic(self, jacobian_algebraic):
         self._jacobian_algebraic = jacobian_algebraic
 
     @property
-    def set_of_parameters(self):
-        return self._set_of_parameters
+    def param(self):
+        return self._param
+
+    @param.setter
+    def param(self, values):
+        # convert module into a class
+        # (StackOverflow: https://tinyurl.com/yk3euon3)
+        self._param = ParamClass(values)
 
     @property
     def options(self):
@@ -277,6 +293,16 @@ def options(self, options):
     def __getitem__(self, key):
         return self.rhs[key]
 
+    def new_copy(self, options=None):
+        "Create an empty copy with identical options, or new options if specified"
+        options = options or self.options
+        new_model = self.__class__(options)
+        new_model.name = self.name
+        new_model.use_jacobian = self.use_jacobian
+        new_model.use_simplify = self.use_simplify
+        new_model.convert_to_format = self.convert_to_format
+        return new_model
+
     def update(self, *submodels):
         """
         Update model to add new physics from submodels

pybamm/parameters/parameter_values.py

@@ -271,12 +271,7 @@ def process_model(self, unprocessed_model, processing="process", inplace=True):
             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
+            model = unprocessed_model.new_copy()
 
         if len(unprocessed_model.rhs) == 0 and len(unprocessed_model.algebraic) == 0:
             raise pybamm.ModelError("Cannot process parameters for empty model")

pybamm/simulation.py

@@ -1,3 +1,7 @@
+#
+# Simulation class
+#
+import pickle
 import pybamm
 import numpy as np
 import copy
@@ -71,7 +75,7 @@ def reset(self):
         """
         A method to reset a simulation back to its unprocessed state.
         """
-        self.model = self._model_class(self._model_options)
+        self.model = self.model.new_copy(self._model_options)
         self.geometry = copy.deepcopy(self._unprocessed_geometry)
         self._model_with_set_params = None
         self._built_model = None
@@ -321,3 +325,23 @@ def specs(
             or spatial_methods
         ):
             self.reset()
+
+    def save(self, filename):
+        """Save simulation using pickle"""
+        if self.model.convert_to_format == "python":
+            # We currently cannot save models in the 'python'
+            raise NotImplementedError(
+                """
+                Cannot save simulation if model format is python.
+                Set model.convert_to_format = 'casadi' instead.
+                """
+            )
+        with open(filename, "wb") as f:
+            pickle.dump(self, f, pickle.HIGHEST_PROTOCOL)
+
+
+def load_sim(filename):
+    """Load a saved simulation"""
+    with open(filename, "rb") as f:
+        sim = pickle.load(f)
+    return sim

pybamm/solvers/casadi_solver.py

@@ -194,7 +194,11 @@ def compute_solution(self, model, t_eval):
         solve_start_time = timer.time()
         pybamm.logger.debug("Calling DAE solver")
         solution = self.integrate_casadi(
-            self.casadi_problem, self.y0, t_eval, mass_matrix=model.mass_matrix.entries
+            self.casadi_rhs,
+            self.casadi_algebraic,
+            self.y0,
+            t_eval,
+            mass_matrix=model.mass_matrix.entries,
         )
         solve_time = timer.time() - solve_start_time
 
@@ -203,7 +207,7 @@ def compute_solution(self, model, t_eval):
 
         return solution, solve_time, termination
 
-    def integrate_casadi(self, problem, y0, t_eval, mass_matrix=None):
+    def integrate_casadi(self, rhs, algebraic, y0, t_eval, mass_matrix=None):
         """
         Solve a DAE model defined by residuals with initial conditions y0.
 
@@ -233,11 +237,24 @@ def integrate_casadi(self, problem, y0, t_eval, mass_matrix=None):
             options["calc_ic"] = True
 
         # set up and solve
+        t = casadi.MX.sym("t")
+        y_diff = casadi.MX.sym("y_diff", rhs(0, y0).shape[0])
+        if algebraic is None:
+            problem = {"t": t, "x": y_diff, "ode": rhs(t, y_diff)}
+        else:
+            y_alg = casadi.MX.sym("y_alg", algebraic(0, y0).shape[0])
+            y = casadi.vertcat(y_diff, y_alg)
+            problem = {
+                "t": t,
+                "x": y_diff,
+                "z": y_alg,
+                "ode": rhs(t, y),
+                "alg": algebraic(t, y),
+            }
         integrator = casadi.integrator("F", self.method, problem, options)
         try:
             # Try solving
-            len_rhs = problem["x"].size()[0]
-            y0_diff, y0_alg = np.split(y0, [len_rhs])
+            y0_diff, y0_alg = np.split(y0, [y_diff.shape[0]])
             sol = integrator(x0=y0_diff, z0=y0_alg)
             y_values = np.concatenate([sol["xf"].full(), sol["zf"].full()])
             return pybamm.Solution(t_eval, y_values, None, None, "final time")