From e44775fb8d5df85ad0e8e88eaaa020628bdb4bd4 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Thu, 19 Nov 2020 16:57:41 +0000
Subject: [PATCH 01/31] Make inputs argument of base_solver.solve a list

---
 pybamm/solvers/base_solver.py | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 99ce106c08..7d31d01a34 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -495,7 +495,7 @@ def calculate_consistent_state(self, model, time=0, inputs=None):
             y0 = y0.flatten()
         return y0
 
-    def solve(self, model, t_eval=None, external_variables=None, inputs=None):
+    def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
         """
         Execute the solver setup and calculate the solution of the model at
         specified times.
@@ -555,14 +555,17 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
             raise pybamm.SolverError("t_eval must increase monotonically")
 
         # Set up external variables and inputs
-        ext_and_inputs = self._set_up_ext_and_inputs(model, external_variables, inputs)
+        ext_and_inputs_list = [
+            self._set_up_ext_and_inputs(model, external_variables, inputs)
+            for inputs in inputs_list
+        ]
 
         # Set up
         timer = pybamm.Timer()
 
         # Set up (if not done already)
         if model not in self.models_set_up:
-            self.set_up(model, ext_and_inputs, t_eval)
+            self.set_up(model, ext_and_inputs_list[0], t_eval)
             self.models_set_up.update(
                 {model: {"initial conditions": model.concatenated_initial_conditions}}
             )
@@ -573,7 +576,7 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
                 # If the new initial conditions are different, set up again
                 # Doing the whole setup again might be slow, but no need to prematurely
                 # optimize this
-                self.set_up(model, ext_and_inputs, t_eval)
+                self.set_up(model, ext_and_inputs_list[0], t_eval)
                 self.models_set_up[model][
                     "initial conditions"
                 ] = model.concatenated_initial_conditions
@@ -581,14 +584,14 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         timer.reset()
 
         # (Re-)calculate consistent initial conditions
-        self._set_initial_conditions(model, ext_and_inputs, update_rhs=True)
+        self._set_initial_conditions(model, ext_and_inputs_list[0], update_rhs=True)
 
         # Non-dimensionalise time
         t_eval_dimensionless = t_eval / model.timescale_eval
 
         # Calculate discontinuities
         discontinuities = [
-            event.expression.evaluate(inputs=inputs)
+            event.expression.evaluate(inputs=inputs_list[0])
             for event in model.discontinuity_events_eval
         ]
 

From 412f3038a7df2c61ea4bc821f7ae4ad3def2ac7c Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Thu, 19 Nov 2020 17:02:09 +0000
Subject: [PATCH 02/31] Make base_solver.solve() return a list of solutions
 objects

Call to integrate is now made inside of a map() wich results in a list of
solutions corresponfing to the list of inputs given to solve().
---
 pybamm/solvers/base_solver.py | 63 ++++++++++++++++++++---------------
 1 file changed, 36 insertions(+), 27 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 7d31d01a34..455be800e5 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -643,7 +643,7 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
         # object, restarting the solver at each discontinuity (and recalculating a
         # consistent state afterwards if a dae)
         old_y0 = model.y0
-        solution = None
+        solutions = None
         for start_index, end_index in zip(start_indices, end_indices):
             pybamm.logger.info(
                 "Calling solver for {} < t < {}".format(
@@ -651,46 +651,55 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
                     t_eval_dimensionless[end_index - 1] * model.timescale_eval,
                 )
             )
-            new_solution = self._integrate(
-                model, t_eval_dimensionless[start_index:end_index], ext_and_inputs
+            ninputs = len(ext_and_inputs_list)
+            new_solutions = list(
+                map(
+                    self._integrate,
+                    [model] * ninputs,
+                    [t_eval_dimensionless[start_index:end_index]] * ninputs,
+                    ext_and_inputs_list
+                )
             )
-            new_solution.solve_time = timer.time()
-            if solution is None:
-                solution = new_solution
+            new_solutions.solve_time = timer.time()
+            if start_index == start_indices[0]:
+                solutions = [sol for sol in new_solutions]
             else:
-                solution.append(new_solution, start_index=0)
+                for i, new_solution in enumerate(new_solutions):
+                    solutions[i].append(new_solution, start_index=0)
 
-            if solution.termination != "final time":
+            if solutions[0].termination != "final time":
                 break
 
             if end_index != len(t_eval_dimensionless):
                 # setup for next integration subsection
-                last_state = solution.y[:, -1]
+                last_states = [sol.y[:, -1] for sol in solutions]
                 # update y0 (for DAE solvers, this updates the initial guess for the
                 # rootfinder)
-                model.y0 = last_state
                 if len(model.algebraic) > 0:
-                    model.y0 = self.calculate_consistent_state(
-                        model, t_eval_dimensionless[end_index], ext_and_inputs
-                    )
-
-        # Assign times
-        solution.set_up_time = set_up_time
-        solution.solve_time = timer.time()
-
-        # restore old y0
-        model.y0 = old_y0
+                    for i, last_state in enumerate(last_states):
+                        model.y0 = last_state
+                        last_states[i] = self.calculate_consistent_state(
+                            model, t_eval_dimensionless[end_index], ext_and_inputs_list[0]
+                        )
+                model.y0 = last_states
 
-        # Add model and inputs to solution
-        solution.model = model
-        solution.inputs = ext_and_inputs
+        for solution in solutions:
+            # Assign times
+            solution.set_up_time = set_up_time
+            solution.solve_time = timer.time()
+            # Add model and inputs to solution
+            solution.model = model
+            solution.inputs = ext_and_inputs_list
 
-        # Copy the timescale_eval and lengthscale_evals
-        solution.timescale_eval = model.timescale_eval
-        solution.length_scales_eval = model.length_scales_eval
+            # Copy the timescale_eval and lengthscale_evals
+            solution.timescale_eval = model.timescale_eval
+            solution.length_scales_eval = model.length_scales_eval
 
         # Identify the event that caused termination
-        termination = self.get_termination_reason(solution, model.events)
+        termination = self.get_termination_reason(solution[0], model.events)
+
+        # restore old y0
+        model.y0 = old_y0
 
         pybamm.logger.info("Finish solving {} ({})".format(model.name, termination))
         pybamm.logger.info(

From 2218b05e2ca69f05e89c2639cb4a5732cfcf14aa Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Fri, 20 Nov 2020 10:32:49 +0000
Subject: [PATCH 03/31] Do not record time in between discontinuity events

---
 pybamm/solvers/base_solver.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 455be800e5..10c555be3a 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -660,7 +660,6 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
                     ext_and_inputs_list
                 )
             )
-            new_solutions.solve_time = timer.time()
             if start_index == start_indices[0]:
                 solutions = [sol for sol in new_solutions]
             else:

From 89fb481daeec39686fabe83950741471ab44503b Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Fri, 20 Nov 2020 10:35:27 +0000
Subject: [PATCH 04/31] Assign individual inputs to solutions

---
 pybamm/solvers/base_solver.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 10c555be3a..3cdf7db8e0 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -682,13 +682,13 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
                         )
                 model.y0 = last_states
 
-        for solution in solutions:
+        for i, solution in enumerate(solutions):
             # Assign times
             solution.set_up_time = set_up_time
             solution.solve_time = timer.time()
             # Add model and inputs to solution
             solution.model = model
-            solution.inputs = ext_and_inputs_list
+            solution.inputs = ext_and_inputs_list[i]
 
             # Copy the timescale_eval and lengthscale_evals
             solution.timescale_eval = model.timescale_eval

From 359934a603cf14ec5883045d892b7d286fbf7cf9 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Fri, 20 Nov 2020 10:35:48 +0000
Subject: [PATCH 05/31] solution -> solutions

---
 pybamm/solvers/base_solver.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 3cdf7db8e0..657528e736 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -695,7 +695,7 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
             solution.length_scales_eval = model.length_scales_eval
 
         # Identify the event that caused termination
-        termination = self.get_termination_reason(solution[0], model.events)
+        termination = self.get_termination_reason(solutions[0], model.events)
 
         # restore old y0
         model.y0 = old_y0
@@ -717,7 +717,7 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
                 "Check whether simulation terminated too early."
             )
 
-        return solution
+        return solutions
 
     def step(
         self,

From f8666a416cb91ff0bf422a8d2323017789417975 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Fri, 20 Nov 2020 18:20:57 +0000
Subject: [PATCH 06/31] Replace map by multiprocessing.Pool.starmap

For now max number of processes is set to 2 because I only have four cores on my machine.
---
 pybamm/solvers/base_solver.py | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 657528e736..b8f1b7e4be 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -8,6 +8,7 @@
 import numpy as np
 import sys
 import itertools
+import multiprocessing as mp
 
 
 class BaseSolver(object):
@@ -652,14 +653,16 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
                 )
             )
             ninputs = len(ext_and_inputs_list)
-            new_solutions = list(
-                map(
+            with mp.Pool(processes=2) as p:
+                new_solutions = p.starmap(
                     self._integrate,
-                    [model] * ninputs,
-                    [t_eval_dimensionless[start_index:end_index]] * ninputs,
-                    ext_and_inputs_list
+                    zip(
+                        [model] * ninputs,
+                        [t_eval_dimensionless[start_index:end_index]] * ninputs,
+                        ext_and_inputs_list
+                    )
                 )
-            )
+
             if start_index == start_indices[0]:
                 solutions = [sol for sol in new_solutions]
             else:

From 8276ffaaf4049ec91a0e14f294a4dd9d5939d26e Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Wed, 25 Nov 2020 13:49:09 +0000
Subject: [PATCH 07/31] Make it possible to supply argument inputs_list as a
 single dictionary

---
 pybamm/solvers/base_solver.py | 31 +++++++++++++++++++++----------
 1 file changed, 21 insertions(+), 10 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index b8f1b7e4be..ae5d403141 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -511,8 +511,9 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
         external_variables : dict
             A dictionary of external variables and their corresponding
             values at the current time
-        inputs : dict, optional
-            Any input parameters to pass to the model when solving
+        inputs : dict or list, optional
+            A dictionary or list of dictionaries describing any input parameters to
+            pass to the model when solving
 
         Raises
         ------
@@ -556,6 +557,8 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
             raise pybamm.SolverError("t_eval must increase monotonically")
 
         # Set up external variables and inputs
+        if not isinstance(inputs_list, list):
+            inputs_list = [inputs_list]
         ext_and_inputs_list = [
             self._set_up_ext_and_inputs(model, external_variables, inputs)
             for inputs in inputs_list
@@ -653,15 +656,23 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
                 )
             )
             ninputs = len(ext_and_inputs_list)
-            with mp.Pool(processes=2) as p:
-                new_solutions = p.starmap(
-                    self._integrate,
-                    zip(
-                        [model] * ninputs,
-                        [t_eval_dimensionless[start_index:end_index]] * ninputs,
-                        ext_and_inputs_list
-                    )
+            if ninputs == 1:
+                new_solution = self._integrate(
+                    model,
+                    t_eval_dimensionless[start_index:end_index],
+                    ext_and_inputs_list[0],
                 )
+                new_solutions = [new_solution]
+            else:
+                with mp.Pool(processes=2) as p:
+                    new_solutions = p.starmap(
+                        self._integrate,
+                        zip(
+                            [model] * ninputs,
+                            [t_eval_dimensionless[start_index:end_index]] * ninputs,
+                            ext_and_inputs_list,
+                        ),
+                    )
 
             if start_index == start_indices[0]:
                 solutions = [sol for sol in new_solutions]

From fc2f2af5b1635eb79cb67f64a73c646327d7b246 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Wed, 25 Nov 2020 13:50:08 +0000
Subject: [PATCH 08/31] [style] black

---
 pybamm/solvers/base_solver.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index ae5d403141..0a9f49b7e2 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -692,7 +692,9 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
                     for i, last_state in enumerate(last_states):
                         model.y0 = last_state
                         last_states[i] = self.calculate_consistent_state(
-                            model, t_eval_dimensionless[end_index], ext_and_inputs_list[0]
+                            model,
+                            t_eval_dimensionless[end_index],
+                            ext_and_inputs_list[0],
                         )
                 model.y0 = last_states
 

From a30479633aec72bd0592aee688adea66eff67884 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Wed, 25 Nov 2020 14:06:30 +0000
Subject: [PATCH 09/31] Fix solve() API

---
 pybamm/solvers/base_solver.py | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 0a9f49b7e2..19f9fa8664 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -496,7 +496,7 @@ def calculate_consistent_state(self, model, time=0, inputs=None):
             y0 = y0.flatten()
         return y0
 
-    def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
+    def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         """
         Execute the solver setup and calculate the solution of the model at
         specified times.
@@ -557,8 +557,12 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
             raise pybamm.SolverError("t_eval must increase monotonically")
 
         # Set up external variables and inputs
-        if not isinstance(inputs_list, list):
-            inputs_list = [inputs_list]
+        #
+        # Argument "inputs" can be either a list of input dicts or
+        # a single dict. The remaining of this function is only working
+        # with variable "input_list", which is a list of dictionaries.
+        # If "inputs" is a single dict, "inputs_list" is a list of only one dict.
+        inputs_list = inputs if isinstance(inputs, list) else [inputs]
         ext_and_inputs_list = [
             self._set_up_ext_and_inputs(model, external_variables, inputs)
             for inputs in inputs_list
@@ -733,7 +737,10 @@ def solve(self, model, t_eval=None, external_variables=None, inputs_list=None):
                 "Check whether simulation terminated too early."
             )
 
-        return solutions
+        if ninputs == 1:
+            return solutions[0]
+        else:
+            return solutions
 
     def step(
         self,

From 8229c69e97f543c04b204694ff21b70a3b030463 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Wed, 25 Nov 2020 14:11:11 +0000
Subject: [PATCH 10/31] Fix logging of solve and setup time

---
 pybamm/solvers/base_solver.py | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 19f9fa8664..0d2eaf7e2a 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -702,10 +702,11 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
                         )
                 model.y0 = last_states
 
+        solve_time = timer.time()
         for i, solution in enumerate(solutions):
             # Assign times
             solution.set_up_time = set_up_time
-            solution.solve_time = timer.time()
+            solution.solve_time = solve_time
             # Add model and inputs to solution
             solution.model = model
             solution.inputs = ext_and_inputs_list[i]
@@ -723,15 +724,15 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         pybamm.logger.info("Finish solving {} ({})".format(model.name, termination))
         pybamm.logger.info(
             "Set-up time: {}, Solve time: {}, Total time: {}".format(
-                timer.format(solution.set_up_time),
-                timer.format(solution.solve_time),
-                timer.format(solution.total_time),
+                timer.format(solutions[0].set_up_time),
+                timer.format(solutions[0].solve_time),
+                timer.format(solutions[0].total_time),
             )
         )
 
-        # Raise error if solution only contains one timestep (except for algebraic
+        # Raise error if solutions[0] only contains one timestep (except for algebraic
         # solvers, where we may only expect one time in the solution)
-        if self.algebraic_solver is False and len(solution.t) == 1:
+        if self.algebraic_solver is False and len(solutions[0].t) == 1:
             raise pybamm.SolverError(
                 "Solution time vector has length 1. "
                 "Check whether simulation terminated too early."

From 60850e780c5e011a7d11d48dbbaf513bfb1dc9fb Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Wed, 25 Nov 2020 16:00:06 +0000
Subject: [PATCH 11/31] Add argument so solve() to set number of processes

---
 pybamm/solvers/base_solver.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 0d2eaf7e2a..da480c1d56 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -496,7 +496,9 @@ def calculate_consistent_state(self, model, time=0, inputs=None):
             y0 = y0.flatten()
         return y0
 
-    def solve(self, model, t_eval=None, external_variables=None, inputs=None):
+    def solve(
+        self, model, t_eval=None, external_variables=None, inputs=None, nproc=None
+    ):
         """
         Execute the solver setup and calculate the solution of the model at
         specified times.
@@ -514,6 +516,9 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         inputs : dict or list, optional
             A dictionary or list of dictionaries describing any input parameters to
             pass to the model when solving
+        nproc : int, optional
+            Number of processes to use when solving for more than one set of input
+            parameters. Defaults to value returned by "os.cpu_count()".
 
         Raises
         ------
@@ -668,7 +673,7 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
                 )
                 new_solutions = [new_solution]
             else:
-                with mp.Pool(processes=2) as p:
+                with mp.Pool(processes=nproc) as p:
                     new_solutions = p.starmap(
                         self._integrate,
                         zip(

From 6dd96a3625e1c2e80fc9f9a9d90c40c17d75d3df Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Fri, 27 Nov 2020 18:22:43 +0000
Subject: [PATCH 12/31] Not parallel processing of ensemble of inputs when
 having discontinuities

---
 pybamm/solvers/base_solver.py | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index da480c1d56..25adb8d72b 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -629,6 +629,13 @@ def solve(
             pybamm.logger.info(
                 "Discontinuity events found at t = {}".format(discontinuities)
             )
+            if isinstance(inputs, list):
+                with RuntimeError as e:
+                    e.message = (
+                        "Cannot solve for a list of input parameters"
+                        " sets with discontinuities"
+                    )
+                    raise e
         else:
             pybamm.logger.info("No discontinuity events found")
 
@@ -694,18 +701,14 @@ def solve(
 
             if end_index != len(t_eval_dimensionless):
                 # setup for next integration subsection
-                last_states = [sol.y[:, -1] for sol in solutions]
+                last_state = solutions[0].y[:, -1]
                 # update y0 (for DAE solvers, this updates the initial guess for the
                 # rootfinder)
+                model.y0 = last_state
                 if len(model.algebraic) > 0:
-                    for i, last_state in enumerate(last_states):
-                        model.y0 = last_state
-                        last_states[i] = self.calculate_consistent_state(
-                            model,
-                            t_eval_dimensionless[end_index],
-                            ext_and_inputs_list[0],
-                        )
-                model.y0 = last_states
+                    model.y0 = self.calculate_consistent_state(
+                        model, t_eval_dimensionless[end_index], ext_and_inputs_list[0]
+                    )
 
         solve_time = timer.time()
         for i, solution in enumerate(solutions):

From 1560d61d53e193eeee037c559de9e7dc15aeebb2 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Fri, 27 Nov 2020 18:24:51 +0000
Subject: [PATCH 13/31] Record solve time for each segment between
 discontinuities

this is required because pybamm.solution.append adds solve time.
---
 pybamm/solvers/base_solver.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 25adb8d72b..4e26295cc6 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -689,7 +689,12 @@ def solve(
                             ext_and_inputs_list,
                         ),
                     )
-
+            # Setting the solve time for each segment
+            # Not sure why required here since overall solve time
+            # is set further down
+            solve_time = timer.time()
+            for sol in new_solutions:
+                sol.solve_time = solve_time
             if start_index == start_indices[0]:
                 solutions = [sol for sol in new_solutions]
             else:

From 8577fea3d35c6814b6a4708819c9f30179ab008f Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Wed, 9 Dec 2020 10:43:03 +0000
Subject: [PATCH 14/31] Add test for ScipySolver.solve() with multiple inputs

---
 tests/unit/test_solvers/test_scipy_solver.py | 31 ++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index fae2775b3f..27fccbdd95 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -217,6 +217,37 @@ def test_model_solver_with_inputs(self):
         np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
         np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
 
+    def test_model_solver_with_multiple_inputs(self):
+        # Create model
+        model = pybamm.BaseModel()
+        model.convert_to_format = "python"
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=domain)
+        model.rhs = {var: -pybamm.InputParameter("rate") * var}
+        model.initial_conditions = {var: 1}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+        t_eval = np.linspace(0, 10, 100)
+        ninputs = 8
+        inputs_list = [{"rate": 0.01 * (i + 1)} for i in range(ninputs)]
+        solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
+        for i in range(ninputs):
+            with self.subTest(i=i):
+                solution = solutions[i]
+                self.assertLess(len(solution.t), len(t_eval))
+                np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+                np.testing.assert_allclose(
+                    solution.y[0], np.exp(-0.01 * (i + 1) * solution.t)
+                )
+
     def test_model_solver_with_external(self):
         # Create model
         model = pybamm.BaseModel()

From c176badaedd6313b3f4c77aee3f97c076cb0086f Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Wed, 9 Dec 2020 10:45:45 +0000
Subject: [PATCH 15/31] solution.t should be the same length as t_eval since
 there are no events

---
 tests/unit/test_solvers/test_scipy_solver.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 27fccbdd95..3ebb2e0230 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -242,8 +242,7 @@ def test_model_solver_with_multiple_inputs(self):
         for i in range(ninputs):
             with self.subTest(i=i):
                 solution = solutions[i]
-                self.assertLess(len(solution.t), len(t_eval))
-                np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+                np.testing.assert_array_equal(solution.t, t_eval)
                 np.testing.assert_allclose(
                     solution.y[0], np.exp(-0.01 * (i + 1) * solution.t)
                 )

From 8bfe384da4299a191d4e917f1dbbe1f69f6c5c1e Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Wed, 9 Dec 2020 11:28:48 +0000
Subject: [PATCH 16/31] Set format to casadi when testing solve with multiple
 outputs. See #1283.

---
 tests/unit/test_solvers/test_scipy_solver.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 3ebb2e0230..46fe974db6 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -220,7 +220,9 @@ def test_model_solver_with_inputs(self):
     def test_model_solver_with_multiple_inputs(self):
         # Create model
         model = pybamm.BaseModel()
-        model.convert_to_format = "python"
+        # Covert to casadi instead of python to avoid pickling of
+        # "EvaluatorPython" objects.
+        model.convert_to_format = "casadi"
         domain = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=domain)
         model.rhs = {var: -pybamm.InputParameter("rate") * var}

From 2962389608c58f2c94c1425cda6bfdf8868dc649 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Thu, 10 Dec 2020 12:16:07 +0000
Subject: [PATCH 17/31] Test for SolverError if discontinuity

---
 pybamm/solvers/base_solver.py                |  4 +---
 tests/unit/test_solvers/test_scipy_solver.py | 22 +++++++++++++++++++-
 2 files changed, 22 insertions(+), 4 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index abb168d4fb..fa8576f984 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -630,12 +630,10 @@ def solve(
                 "Discontinuity events found at t = {}".format(discontinuities)
             )
             if isinstance(inputs, list):
-                with RuntimeError as e:
-                    e.message = (
+                raise pybamm.SolverError(
                         "Cannot solve for a list of input parameters"
                         " sets with discontinuities"
                     )
-                    raise e
         else:
             pybamm.logger.info("No discontinuity events found")
 
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 3a39b455ce..151d212a0a 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -235,11 +235,31 @@ def test_model_solver_with_multiple_inputs(self):
         disc = pybamm.Discretisation(mesh, spatial_methods)
         disc.process_model(model)
 
-        # Solve
         solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
         t_eval = np.linspace(0, 10, 100)
         ninputs = 8
         inputs_list = [{"rate": 0.01 * (i + 1)} for i in range(ninputs)]
+
+        model.events = [
+            pybamm.Event(
+                "discontinuity",
+                pybamm.Scalar(t_eval[-1] / 2),
+                event_type=pybamm.EventType.DISCONTINUITY,
+            )
+        ]
+        with self.assertRaisesRegex(
+            pybamm.SolverError,
+            (
+                "Cannot solve for a list of input parameters"
+                " sets with discontinuities"
+            ),
+        ):
+            solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
+
+        model.events = []
+        # Must set up model again to update attribute
+        # "discontinuity_events_eval" of object "model".
+        solver.set_up(model, inputs_list[0], t_eval)
         solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
         for i in range(ninputs):
             with self.subTest(i=i):

From 45a2e5ec795d57cb6e9b637794d546c0fe1b58ab Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Thu, 10 Dec 2020 14:20:07 +0000
Subject: [PATCH 18/31] Add description of returned object in docstring of
 solve

---
 pybamm/solvers/base_solver.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index fa8576f984..ce93f6dc9f 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -499,8 +499,7 @@ def calculate_consistent_state(self, model, time=0, inputs=None):
     def solve(
         self, model, t_eval=None, external_variables=None, inputs=None, nproc=None
     ):
-        """
-        Execute the solver setup and calculate the solution of the model at
+        """Execute the solver setup and calculate the solution of the model at
         specified times.
 
         Parameters
@@ -520,6 +519,13 @@ def solve(
             Number of processes to use when solving for more than one set of input
             parameters. Defaults to value returned by "os.cpu_count()".
 
+        Returns
+        ------- 
+        :class:`pybamm.Solution` or list of :class:`pybamm.Solution` objects.
+             If type of `inputs` is `list`, return a list of corresponding 
+             :class:`pybamm.Solution` objects.
+            
+
         Raises
         ------
         :class:`pybamm.ModelError`

From 34af020cf7e5bd934ac1e769113e4a4133921055 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Thu, 10 Dec 2020 14:52:07 +0000
Subject: [PATCH 19/31] Add comments explaining why only first element of list
 of inputs is passed.

---
 pybamm/solvers/base_solver.py | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index ce93f6dc9f..77c946282b 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -584,6 +584,11 @@ def solve(
 
         # Set up (if not done already)
         if model not in self.models_set_up:
+            # It is assumed that when len(inputs_list) > 1, model set
+            # up (initial condition, time-scale and length-scale) does
+            # not depend on input parameters. Thefore only `ext_and_inputs[0]`
+            # is passed to `set_up`.
+            # See https://github.com/pybamm-team/PyBaMM/pull/1261
             self.set_up(model, ext_and_inputs_list[0], t_eval)
             self.models_set_up.update(
                 {model: {"initial conditions": model.concatenated_initial_conditions}}
@@ -603,6 +608,10 @@ def solve(
         timer.reset()
 
         # (Re-)calculate consistent initial conditions
+        # Assuming initial conditions do not depend on input parameters
+        # when len(inputs_list) > 1, only `ext_and_inputs_list[0]`
+        # is passed to `_set_initial_conditions`.
+        # See https://github.com/pybamm-team/PyBaMM/pull/1261
         self._set_initial_conditions(model, ext_and_inputs_list[0], update_rhs=True)
 
         # Non-dimensionalise time
@@ -610,6 +619,10 @@ def solve(
 
         # Calculate discontinuities
         discontinuities = [
+            # Assuming that discontinuities do not depend on
+            # input parameters when len(input_list) > 1, only
+            # `input_list[0]` is passed to `evaluate`.
+            # See https://github.com/pybamm-team/PyBaMM/pull/1261
             event.expression.evaluate(inputs=inputs_list[0])
             for event in model.discontinuity_events_eval
         ]

From 110d81eb8cc2118cfb43034bbae1c2dee44609e9 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Thu, 10 Dec 2020 15:37:48 +0000
Subject: [PATCH 20/31] Add comment justifying recording of solve_time for each
 integration segment

---
 pybamm/solvers/base_solver.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 77c946282b..d2c73db078 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -520,11 +520,11 @@ def solve(
             parameters. Defaults to value returned by "os.cpu_count()".
 
         Returns
-        ------- 
+        -------
         :class:`pybamm.Solution` or list of :class:`pybamm.Solution` objects.
-             If type of `inputs` is `list`, return a list of corresponding 
+             If type of `inputs` is `list`, return a list of corresponding
              :class:`pybamm.Solution` objects.
-            
+
 
         Raises
         ------
@@ -650,9 +650,9 @@ def solve(
             )
             if isinstance(inputs, list):
                 raise pybamm.SolverError(
-                        "Cannot solve for a list of input parameters"
-                        " sets with discontinuities"
-                    )
+                    "Cannot solve for a list of input parameters"
+                    " sets with discontinuities"
+                )
         else:
             pybamm.logger.info("No discontinuity events found")
 
@@ -706,9 +706,9 @@ def solve(
                             ext_and_inputs_list,
                         ),
                     )
-            # Setting the solve time for each segment
-            # Not sure why required here since overall solve time
-            # is set further down
+            # Setting the solve time for each segment.
+            # pybamm.Solution.append assumes attribute
+            # solve_time.
             solve_time = timer.time()
             for sol in new_solutions:
                 sol.solve_time = solve_time

From 0837b78b5924a7949268d0a8771a5573076ea517 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Fri, 11 Dec 2020 17:38:41 +0000
Subject: [PATCH 21/31] Split test for solve with multiple inputs in two.

---
 tests/unit/test_solvers/test_scipy_solver.py | 46 ++++++++++++++------
 1 file changed, 33 insertions(+), 13 deletions(-)

diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 151d212a0a..5c2821647a 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -217,7 +217,39 @@ def test_model_solver_with_inputs(self):
         np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
         np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
 
-    def test_model_solver_with_multiple_inputs(self):
+    def test_model_solver_multiple_inputs_happy_path(self):
+        # Create model
+        model = pybamm.BaseModel()
+        # Covert to casadi instead of python to avoid pickling of
+        # "EvaluatorPython" objects.
+        model.convert_to_format = "casadi"
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=domain)
+        model.rhs = {var: -pybamm.InputParameter("rate") * var}
+        model.initial_conditions = {var: 1}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+        t_eval = np.linspace(0, 10, 100)
+        ninputs = 8
+        inputs_list = [{"rate": 0.01 * (i + 1)} for i in range(ninputs)]
+
+        solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
+        for i in range(ninputs):
+            with self.subTest(i=i):
+                solution = solutions[i]
+                np.testing.assert_array_equal(solution.t, t_eval)
+                np.testing.assert_allclose(
+                    solution.y[0], np.exp(-0.01 * (i + 1) * solution.t)
+                )
+
+    def test_model_solver_multiple_inputs_discontinuity_error(self):
         # Create model
         model = pybamm.BaseModel()
         # Covert to casadi instead of python to avoid pickling of
@@ -256,18 +288,6 @@ def test_model_solver_with_multiple_inputs(self):
         ):
             solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
 
-        model.events = []
-        # Must set up model again to update attribute
-        # "discontinuity_events_eval" of object "model".
-        solver.set_up(model, inputs_list[0], t_eval)
-        solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
-        for i in range(ninputs):
-            with self.subTest(i=i):
-                solution = solutions[i]
-                np.testing.assert_array_equal(solution.t, t_eval)
-                np.testing.assert_allclose(
-                    solution.y[0], np.exp(-0.01 * (i + 1) * solution.t)
-                )
 
     def test_model_solver_with_external(self):
         # Create model

From b47c30b0f2e9526911ab84e68fed4576e6a69340 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Fri, 11 Dec 2020 17:40:41 +0000
Subject: [PATCH 22/31] Raise SolverError if input parameters appear in initial
 conditions

---
 pybamm/solvers/base_solver.py                | 15 ++++++++++
 tests/unit/test_solvers/test_scipy_solver.py | 31 ++++++++++++++++++++
 2 files changed, 46 insertions(+)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 7caa8b114c..927c3569c9 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -622,6 +622,21 @@ def solve(
         # when len(inputs_list) > 1, only `ext_and_inputs_list[0]`
         # is passed to `_set_initial_conditions`.
         # See https://github.com/pybamm-team/PyBaMM/pull/1261
+        if len(inputs_list) > 1:
+            all_inputs_names = set(
+                itertools.chain.from_iterable(
+                    [ext_and_inputs.keys() for ext_and_inputs in ext_and_inputs_list]
+                )
+            )
+            initial_conditions_node_names = set(
+                [it.name for it in model.concatenated_initial_conditions.pre_order()]
+            )
+            if all_inputs_names.issubset(initial_conditions_node_names):
+                raise pybamm.SolverError(
+                    "Input parameters cannot appear in expression "
+                    "for initial conditions."
+                )
+
         self._set_initial_conditions(model, ext_and_inputs_list[0], update_rhs=True)
 
         # Non-dimensionalise time
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 5c2821647a..7d676cad47 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -288,6 +288,37 @@ def test_model_solver_multiple_inputs_discontinuity_error(self):
         ):
             solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
 
+    def test_model_solver_multiple_inputs_initial_conditions_error(self):
+        # Create model
+        model = pybamm.BaseModel()
+        # Covert to casadi instead of python to avoid pickling of
+        # "EvaluatorPython" objects.
+        model.convert_to_format = "casadi"
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=domain)
+        model.rhs = {var: -pybamm.InputParameter("rate") * var}
+        model.initial_conditions = {var: 2 * pybamm.InputParameter("rate")}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+        t_eval = np.linspace(0, 10, 100)
+        ninputs = 8
+        inputs_list = [{"rate": 0.01 * (i + 1)} for i in range(ninputs)]
+
+        with self.assertRaisesRegex(
+            pybamm.SolverError,
+            (
+                "Input parameters cannot appear in expression "
+                "for initial conditions."
+            ),
+        ):
+            solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
 
     def test_model_solver_with_external(self):
         # Create model

From 48e002e26568bb6488142264505735fde2271ebb Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Fri, 11 Dec 2020 18:10:45 +0000
Subject: [PATCH 23/31] Do not return solutions list when testing for
 exceptions

---
 tests/unit/test_solvers/test_scipy_solver.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 7d676cad47..4f9d2a0a47 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -286,7 +286,7 @@ def test_model_solver_multiple_inputs_discontinuity_error(self):
                 " sets with discontinuities"
             ),
         ):
-            solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
+            solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
 
     def test_model_solver_multiple_inputs_initial_conditions_error(self):
         # Create model
@@ -318,7 +318,7 @@ def test_model_solver_multiple_inputs_initial_conditions_error(self):
                 "for initial conditions."
             ),
         ):
-            solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
+            solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
 
     def test_model_solver_with_external(self):
         # Create model

From 35491ba0cedd6ab2da9e25da57710de2e41cb96f Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Mon, 14 Dec 2020 10:43:53 +0000
Subject: [PATCH 24/31] Update changelog

---
 CHANGELOG.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index ad6f3848bd..1b3079ee13 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Solver method `solve()` can take a list of inputs as the `inputs` keyword argument, in the form of a list of dictionaries. In this case the model is solved for each element in the list, and `solve()` returns a list of solutions mapping the set of inputs to the solutions of the model. Note that `solve()` can still take a single dictionary as the `inputs` keyword argument. In this case the behaviour is unchanged compared to previous versions.
 -   Added option to make the total interfacial current density a state ([#1280](https://github.com/pybamm-team/PyBaMM/pull/1280))
 -   Added functionality to initialize a model using the solution from another model ([#1278](https://github.com/pybamm-team/PyBaMM/pull/1278))
 -   Added submodels for active material ([#1262](https://github.com/pybamm-team/PyBaMM/pull/1262))
@@ -9,6 +10,7 @@
 
 ## Optimizations
 
+-   If solver method `solve()` is passed a list of inputs as the `inputs` keyword argument, the resolution of the model for each input set is spread acrosss several Python processes, usually running in parallel on different processors. The default number of processes is the number of processors available. `solve()` takes a new keyword argument `nproc` which can be used to set this number a manually.
 -   Operations such as `1*x` and `0+x` now directly return `x` ([#1252](https://github.com/pybamm-team/PyBaMM/pull/1252))
 
 ## Bug fixes

From dab096d31a669a153bce0b8a119023345f4a7482 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Fri, 11 Dec 2020 16:11:43 +0000
Subject: [PATCH 25/31] Add methods getstate and setstate to avoid including
 method _evaluate to pickle

---
 pybamm/expression_tree/operations/evaluate.py | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/pybamm/expression_tree/operations/evaluate.py b/pybamm/expression_tree/operations/evaluate.py
index a18010ae71..08b40eb605 100644
--- a/pybamm/expression_tree/operations/evaluate.py
+++ b/pybamm/expression_tree/operations/evaluate.py
@@ -485,6 +485,7 @@ def __init__(self, symbol):
         python_str = python_str + "\nself._evaluate = evaluate"
 
         self._python_str = python_str
+        self._result_var = result_var
         self._symbol = symbol
 
         # compile and run the generated python code,
@@ -507,6 +508,23 @@ def evaluate(self, t=None, y=None, y_dot=None, inputs=None, known_evals=None):
         else:
             return result
 
+    def __getstate__(self):
+        # Control the state of instances of EvaluatorPython
+        # before pickling. Method "_evaluate" cannot be pickled.
+        # See https://github.com/pybamm-team/PyBaMM/issues/1283
+        state = self.__dict__.copy()
+        del state["_evaluate"]
+        return state
+
+    def __setstate__(self, state):
+        # Restore pickled attributes and
+        # compile code from "python_str"
+        # Execution of bytecode (re)adds attribute
+        # "_method"
+        self.__dict__.update(state)
+        compiled_function = compile(self._python_str, self._result_var, "exec")
+        exec(compiled_function)
+
 
 class EvaluatorJax:
     """

From 35b26ddb56a8e9b7ba1c071de8ac09599593ba30 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Mon, 14 Dec 2020 11:47:28 +0000
Subject: [PATCH 26/31] Test solving multiple inputs using multiprocessing with
 model in both python and casadi mode

---
 tests/unit/test_solvers/test_scipy_solver.py | 64 ++++++++++----------
 1 file changed, 31 insertions(+), 33 deletions(-)

diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 4f9d2a0a47..02dd106274 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -218,36 +218,37 @@ def test_model_solver_with_inputs(self):
         np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
 
     def test_model_solver_multiple_inputs_happy_path(self):
-        # Create model
-        model = pybamm.BaseModel()
-        # Covert to casadi instead of python to avoid pickling of
-        # "EvaluatorPython" objects.
-        model.convert_to_format = "casadi"
-        domain = ["negative electrode", "separator", "positive electrode"]
-        var = pybamm.Variable("var", domain=domain)
-        model.rhs = {var: -pybamm.InputParameter("rate") * var}
-        model.initial_conditions = {var: 1}
-        # No need to set parameters; can use base discretisation (no spatial
-        # operators)
-        # create discretisation
-        mesh = get_mesh_for_testing()
-        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
-        disc = pybamm.Discretisation(mesh, spatial_methods)
-        disc.process_model(model)
-
-        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
-        t_eval = np.linspace(0, 10, 100)
-        ninputs = 8
-        inputs_list = [{"rate": 0.01 * (i + 1)} for i in range(ninputs)]
+        for convert_to_format in ["python", "casadi"]:
+            # Create model
+            model = pybamm.BaseModel()
+            # Covert to casadi instead of python to avoid pickling of
+            # "EvaluatorPython" objects.
+            model.convert_to_format = convert_to_format
+            domain = ["negative electrode", "separator", "positive electrode"]
+            var = pybamm.Variable("var", domain=domain)
+            model.rhs = {var: -pybamm.InputParameter("rate") * var}
+            model.initial_conditions = {var: 1}
+            # No need to set parameters; can use base discretisation (no spatial
+            # operators)
+            # create discretisation
+            mesh = get_mesh_for_testing()
+            spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+            disc = pybamm.Discretisation(mesh, spatial_methods)
+            disc.process_model(model)
 
-        solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
-        for i in range(ninputs):
-            with self.subTest(i=i):
-                solution = solutions[i]
-                np.testing.assert_array_equal(solution.t, t_eval)
-                np.testing.assert_allclose(
-                    solution.y[0], np.exp(-0.01 * (i + 1) * solution.t)
-                )
+            solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+            t_eval = np.linspace(0, 10, 100)
+            ninputs = 8
+            inputs_list = [{"rate": 0.01 * (i + 1)} for i in range(ninputs)]
+
+            solutions = solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
+            for i in range(ninputs):
+                with self.subTest(i=i):
+                    solution = solutions[i]
+                    np.testing.assert_array_equal(solution.t, t_eval)
+                    np.testing.assert_allclose(
+                        solution.y[0], np.exp(-0.01 * (i + 1) * solution.t)
+                    )
 
     def test_model_solver_multiple_inputs_discontinuity_error(self):
         # Create model
@@ -313,10 +314,7 @@ def test_model_solver_multiple_inputs_initial_conditions_error(self):
 
         with self.assertRaisesRegex(
             pybamm.SolverError,
-            (
-                "Input parameters cannot appear in expression "
-                "for initial conditions."
-            ),
+            ("Input parameters cannot appear in expression " "for initial conditions."),
         ):
             solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
 

From 35db29dddc66200fb1d31e98dd9dc88ff42f994c Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Mon, 14 Dec 2020 11:59:23 +0000
Subject: [PATCH 27/31] Raise SolverError if attempt to solve for list of
 inputs with model in jax format.

---
 pybamm/solvers/base_solver.py                |  6 ++++
 tests/unit/test_solvers/test_scipy_solver.py | 32 ++++++++++++++++++++
 2 files changed, 38 insertions(+)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 927c3569c9..01402ee1fc 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -545,6 +545,12 @@ def solve(
         """
         pybamm.logger.info("Start solving {} with {}".format(model.name, self.name))
 
+        # Cannot use multiprocessing with model in "jax" format
+        if(len(inputs) > 1) and model.convert_to_format == "jax":
+            raise pybamm.SolverError(
+                "Cannot solve list of inputs with multiprocessing "
+                "when model in format \"jax\"."
+            )
         # Make sure model isn't empty
         if len(model.rhs) == 0 and len(model.algebraic) == 0:
             if not isinstance(self, pybamm.DummySolver):
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 02dd106274..2af2aa803b 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -318,6 +318,38 @@ def test_model_solver_multiple_inputs_initial_conditions_error(self):
         ):
             solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
 
+    def test_model_solver_multiple_inputs_jax_format_error(self):
+        # Create model
+        model = pybamm.BaseModel()
+        # Covert to casadi instead of python to avoid pickling of
+        # "EvaluatorPython" objects.
+        model.convert_to_format = "jax"
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=domain)
+        model.rhs = {var: -pybamm.InputParameter("rate") * var}
+        model.initial_conditions = {var: 2 * pybamm.InputParameter("rate")}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+        t_eval = np.linspace(0, 10, 100)
+        ninputs = 8
+        inputs_list = [{"rate": 0.01 * (i + 1)} for i in range(ninputs)]
+
+        with self.assertRaisesRegex(
+            pybamm.SolverError,
+            (
+                "Cannot solve list of inputs with multiprocessing "
+                'when model in format "jax".'
+            ),
+        ):
+            solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
+
     def test_model_solver_with_external(self):
         # Create model
         model = pybamm.BaseModel()

From d3e132ae9da1e50dc2b6558c8ecb1a324f293937 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Mon, 14 Dec 2020 15:31:54 +0000
Subject: [PATCH 28/31] relocate check for jax format

So that `inputs_list` is defined
---
 pybamm/solvers/base_solver.py | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 01402ee1fc..f38a972eb7 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -545,12 +545,6 @@ def solve(
         """
         pybamm.logger.info("Start solving {} with {}".format(model.name, self.name))
 
-        # Cannot use multiprocessing with model in "jax" format
-        if(len(inputs) > 1) and model.convert_to_format == "jax":
-            raise pybamm.SolverError(
-                "Cannot solve list of inputs with multiprocessing "
-                "when model in format \"jax\"."
-            )
         # Make sure model isn't empty
         if len(model.rhs) == 0 and len(model.algebraic) == 0:
             if not isinstance(self, pybamm.DummySolver):
@@ -595,6 +589,13 @@ def solve(
             for inputs in inputs_list
         ]
 
+        # Cannot use multiprocessing with model in "jax" format
+        if(len(inputs_list) > 1) and model.convert_to_format == "jax":
+            raise pybamm.SolverError(
+                "Cannot solve list of inputs with multiprocessing "
+                "when model in format \"jax\"."
+            )
+
         # Set up
         timer = pybamm.Timer()
 

From 43da47a336511347e1e29e1cc2fc4321eac6aaf3 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Tue, 5 Jan 2021 20:30:29 +0000
Subject: [PATCH 29/31] Add call to close() and join() for coverage to pick up
 processes forkred by Pool

---
 pybamm/solvers/base_solver.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index f38a972eb7..bb3418922e 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -738,6 +738,8 @@ def solve(
                             ext_and_inputs_list,
                         ),
                     )
+                    p.close()
+                    p.join()
             # Setting the solve time for each segment.
             # pybamm.Solution.append assumes attribute
             # solve_time.

From 7729e87d32546a1ce4623580c25de29b89d2bda7 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Wed, 6 Jan 2021 07:38:28 +0000
Subject: [PATCH 30/31] Enable multiprocessing mode for coverage

---
 tox.ini | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/tox.ini b/tox.ini
index 4e90a958c3..b7f094fcee 100644
--- a/tox.ini
+++ b/tox.ini
@@ -50,6 +50,11 @@ deps =
      scikits.odes
 commands = 
      coverage run run-tests.py --nosub
+     # Some tests make use of multiple processes through
+     # multiprocessing. Coverage data is then generated for each
+     # process separately and data must then be combined into one
+     # single coverage data file.
+     coverage combine
      coverage xml
 
 [testenv:docs]
@@ -114,4 +119,7 @@ ignore=
     W605,
 
 [coverage:run]
-source = pybamm
\ No newline at end of file
+source = pybamm
+# By default coverage data isn't collected in forked processes, see
+# https://coverage.readthedocs.io/en/coverage-5.3.1/subprocess.html
+concurrency = multiprocessing

From e5ae0a1c2e7b6e7ac50f33571695fe6fbe48de59 Mon Sep 17 00:00:00 2001
From: Thibault Lestang <thibault.lestang@cs.ox.ac.uk>
Date: Wed, 6 Jan 2021 09:34:32 +0000
Subject: [PATCH 31/31] Remove out of date comments

---
 tests/unit/test_solvers/test_scipy_solver.py | 17 +----------------
 1 file changed, 1 insertion(+), 16 deletions(-)

diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 2af2aa803b..b3851c6f0f 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -200,10 +200,9 @@ def test_model_solver_with_inputs(self):
         var = pybamm.Variable("var", domain=domain)
         model.rhs = {var: -pybamm.InputParameter("rate") * var}
         model.initial_conditions = {var: 1}
-        model.events = [pybamm.Event("var=0.5", pybamm.min(var - 0.5))]
         # No need to set parameters; can use base discretisation (no spatial
         # operators)
-
+        model.events = [pybamm.Event("var=0.5", pybamm.min(var - 0.5))]
         # create discretisation
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
@@ -221,15 +220,11 @@ def test_model_solver_multiple_inputs_happy_path(self):
         for convert_to_format in ["python", "casadi"]:
             # Create model
             model = pybamm.BaseModel()
-            # Covert to casadi instead of python to avoid pickling of
-            # "EvaluatorPython" objects.
             model.convert_to_format = convert_to_format
             domain = ["negative electrode", "separator", "positive electrode"]
             var = pybamm.Variable("var", domain=domain)
             model.rhs = {var: -pybamm.InputParameter("rate") * var}
             model.initial_conditions = {var: 1}
-            # No need to set parameters; can use base discretisation (no spatial
-            # operators)
             # create discretisation
             mesh = get_mesh_for_testing()
             spatial_methods = {"macroscale": pybamm.FiniteVolume()}
@@ -253,15 +248,11 @@ def test_model_solver_multiple_inputs_happy_path(self):
     def test_model_solver_multiple_inputs_discontinuity_error(self):
         # Create model
         model = pybamm.BaseModel()
-        # Covert to casadi instead of python to avoid pickling of
-        # "EvaluatorPython" objects.
         model.convert_to_format = "casadi"
         domain = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=domain)
         model.rhs = {var: -pybamm.InputParameter("rate") * var}
         model.initial_conditions = {var: 1}
-        # No need to set parameters; can use base discretisation (no spatial
-        # operators)
         # create discretisation
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
@@ -292,15 +283,11 @@ def test_model_solver_multiple_inputs_discontinuity_error(self):
     def test_model_solver_multiple_inputs_initial_conditions_error(self):
         # Create model
         model = pybamm.BaseModel()
-        # Covert to casadi instead of python to avoid pickling of
-        # "EvaluatorPython" objects.
         model.convert_to_format = "casadi"
         domain = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=domain)
         model.rhs = {var: -pybamm.InputParameter("rate") * var}
         model.initial_conditions = {var: 2 * pybamm.InputParameter("rate")}
-        # No need to set parameters; can use base discretisation (no spatial
-        # operators)
         # create discretisation
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
@@ -321,8 +308,6 @@ def test_model_solver_multiple_inputs_initial_conditions_error(self):
     def test_model_solver_multiple_inputs_jax_format_error(self):
         # Create model
         model = pybamm.BaseModel()
-        # Covert to casadi instead of python to avoid pickling of
-        # "EvaluatorPython" objects.
         model.convert_to_format = "jax"
         domain = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=domain)