Merge pull request #1026 from pybamm-team/issue-1024-plot-external-variables

#1024 fix external variables

Author: valentinsulzer

CHANGELOG.md

@@ -27,6 +27,8 @@
 
 ## Bug fixes
 
+-   Fix storing and plotting external variables in the solution ([#1026](https://github.com/pybamm-team/PyBaMM/pull/1026))
+-   Fix running a simulation with a model that is already discretized ([#1025](https://github.com/pybamm-team/PyBaMM/pull/1025))
 -   Fix CI not triggering for PR. ([#1013](https://github.com/pybamm-team/PyBaMM/pull/1013))
 -   Fix schedule testing running too often. ([#1010](https://github.com/pybamm-team/PyBaMM/pull/1010))
 -   Fix doctests failing due to mismatch in unsorted output.([#990](https://github.com/pybamm-team/PyBaMM/pull/990))

pybamm/expression_tree/variable.py

@@ -197,6 +197,8 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, inputs=None):
                     )
                 )
             else:
+                if isinstance(out, np.ndarray) and out.ndim == 1:
+                    out = out[:, np.newaxis]
                 return out
         # raise more informative error if can't find name in dict
         except KeyError:

pybamm/simulation.py

@@ -418,7 +418,12 @@ def solve(
                 t_eval = np.linspace(0, t_end, 100)
 
             self.t_eval = t_eval
-            self._solution = solver.solve(self.built_model, t_eval, inputs=inputs)
+            self._solution = solver.solve(
+                self.built_model,
+                t_eval,
+                external_variables=external_variables,
+                inputs=inputs,
+            )
 
         elif self.operating_mode == "with experiment":
             if t_eval is not None:

pybamm/solvers/base_solver.py

@@ -775,7 +775,7 @@ def get_termination_reason(self, solution, events):
                         event.expression.evaluate(
                             solution.t_event,
                             solution.y_event,
-                            inputs={k: v[-1] for k, v in solution.inputs.items()},
+                            inputs={k: v[:, -1] for k, v in solution.inputs.items()},
                         )
                     )
             termination_event = min(final_event_values, key=final_event_values.get)

pybamm/solvers/processed_variable.py

@@ -59,14 +59,14 @@ def __init__(self, base_variable, solution, known_evals=None):
             self.base_eval, self.known_evals[solution.t[0]] = base_variable.evaluate(
                 solution.t[0],
                 solution.y[:, 0],
-                inputs={name: inp[0] for name, inp in solution.inputs.items()},
+                inputs={name: inp[:, 0] for name, inp in solution.inputs.items()},
                 known_evals=self.known_evals[solution.t[0]],
             )
         else:
             self.base_eval = base_variable.evaluate(
                 solution.t[0],
                 solution.y[:, 0],
-                inputs={name: inp[0] for name, inp in solution.inputs.items()},
+                inputs={name: inp[:, 0] for name, inp in solution.inputs.items()},
             )
 
         # handle 2D (in space) finite element variables differently
@@ -115,7 +115,7 @@ def initialise_0D(self):
         for idx in range(len(self.t_sol)):
             t = self.t_sol[idx]
             u = self.u_sol[:, idx]
-            inputs = {name: inp[idx] for name, inp in self.inputs.items()}
+            inputs = {name: inp[:, idx] for name, inp in self.inputs.items()}
             if self.known_evals:
                 entries[idx], self.known_evals[t] = self.base_variable.evaluate(
                     t, u, inputs=inputs, known_evals=self.known_evals[t]
@@ -151,7 +151,7 @@ def initialise_1D(self, fixed_t=False):
         for idx in range(len(self.t_sol)):
             t = self.t_sol[idx]
             u = self.u_sol[:, idx]
-            inputs = {name: inp[idx] for name, inp in self.inputs.items()}
+            inputs = {name: inp[:, idx] for name, inp in self.inputs.items()}
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
                     t, u, inputs=inputs, known_evals=self.known_evals[t]
@@ -270,7 +270,7 @@ def initialise_2D(self):
         for idx in range(len(self.t_sol)):
             t = self.t_sol[idx]
             u = self.u_sol[:, idx]
-            inputs = {name: inp[idx] for name, inp in self.inputs.items()}
+            inputs = {name: inp[:, idx] for name, inp in self.inputs.items()}
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
                     t, u, inputs=inputs, known_evals=self.known_evals[t]
@@ -330,7 +330,7 @@ def initialise_2D_scikit_fem(self):
         for idx in range(len(self.t_sol)):
             t = self.t_sol[idx]
             u = self.u_sol[:, idx]
-            inputs = {name: inp[idx] for name, inp in self.inputs.items()}
+            inputs = {name: inp[:, idx] for name, inp in self.inputs.items()}
 
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(

pybamm/solvers/solution.py

@@ -108,8 +108,12 @@ def inputs(self, inputs):
             self.has_symbolic_inputs = False
             self._inputs = {}
             for name, inp in inputs.items():
+                # Convert number to vector of the right shape
                 if isinstance(inp, numbers.Number):
-                    inp = inp * np.ones_like(self.t)
+                    inp = inp * np.ones((1, len(self.t)))
+                # Tile a vector
+                else:
+                    inp = np.tile(inp, len(self.t))
                 self._inputs[name] = inp
 
     @property
@@ -316,7 +320,7 @@ def append(self, solution, start_index=1, create_sub_solutions=False):
         self._y = np.concatenate((self._y, solution.y[:, start_index:]), axis=1)
         for name, inp in self.inputs.items():
             solution_inp = solution.inputs[name]
-            self.inputs[name] = np.concatenate((inp, solution_inp[start_index:]))
+            self.inputs[name] = np.c_[inp, solution_inp[:, start_index:]]
         # Update solution time
         self.solve_time += solution.solve_time
         # Update termination

tests/integration/test_solvers/test_solution.py

@@ -59,6 +59,65 @@ def test_append(self):
             decimal=4,
         )
 
+    def test_append_external_variables(self):
+        model = pybamm.lithium_ion.SPM(
+            {
+                "thermal": "lumped",
+                "external submodels": ["thermal", "negative particle"],
+            }
+        )
+        # create geometry
+        geometry = model.default_geometry
+
+        # load parameter values and process model and geometry
+        param = model.default_parameter_values
+        param.process_model(model)
+        param.process_geometry(geometry)
+
+        # set mesh
+        mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
+
+        # discretise model
+        disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+        disc.process_model(model)
+
+        # solve model
+        solver = model.default_solver
+
+        T_av = 0
+        c_s_n_av = np.ones((10, 1)) * 0.6
+        external_variables = {
+            "Volume-averaged cell temperature": T_av,
+            "X-averaged negative particle concentration": c_s_n_av,
+        }
+
+        # Step
+        dt = 0.1
+        sol_step = None
+        for _ in range(5):
+            sol_step = solver.step(
+                sol_step, model, dt, external_variables=external_variables
+            )
+        np.testing.assert_array_equal(
+            sol_step.inputs["Volume-averaged cell temperature"],
+            np.zeros((1, len(sol_step.t))),
+        )
+        np.testing.assert_array_equal(
+            sol_step.inputs["X-averaged negative particle concentration"],
+            np.ones((mesh["negative particle"][0].npts, len(sol_step.t))) * 0.6,
+        )
+
+        # Solve
+        t_eval = np.linspace(0, 3600)
+        sol = solver.solve(model, t_eval, external_variables=external_variables)
+        np.testing.assert_array_equal(
+            sol.inputs["Volume-averaged cell temperature"], np.zeros((1, len(sol.t))),
+        )
+        np.testing.assert_array_equal(
+            sol.inputs["X-averaged negative particle concentration"],
+            np.ones((mesh["negative particle"][0].npts, len(sol.t))) * 0.6,
+        )
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

tests/unit/test_simulation.py

@@ -235,16 +235,30 @@ def test_get_variable_array(self):
         self.assertIsInstance(c_e, np.ndarray)
 
     def test_set_external_variable(self):
-        model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
+        model_options = {
+            "thermal": "lumped",
+            "external submodels": ["thermal", "negative particle"],
+        }
         model = pybamm.lithium_ion.SPMe(model_options)
         sim = pybamm.Simulation(model)
 
         T_av = 0
+        c_s_n_av = np.ones((10, 1)) * 0.5
+        external_variables = {
+            "Volume-averaged cell temperature": T_av,
+            "X-averaged negative particle concentration": c_s_n_av,
+        }
 
-        dt = 0.001
+        # Step
+        dt = 0.1
+        for _ in range(5):
+            sim.step(dt, external_variables=external_variables)
+        sim.plot(testing=True)
 
-        external_variables = {"Volume-averaged cell temperature": T_av}
-        sim.step(dt, external_variables=external_variables)
+        # Solve
+        t_eval = np.linspace(0, 3600)
+        sim.solve(t_eval, external_variables=external_variables)
+        sim.plot(testing=True)
 
     def test_step(self):
 
@@ -301,7 +315,7 @@ def test_step_with_inputs(self):
         self.assertEqual(sim.solution.t[1], dt / tau)
         self.assertEqual(sim.solution.t[2], 2 * dt / tau)
         np.testing.assert_array_equal(
-            sim.solution.inputs["Current function [A]"], np.array([1, 1, 2])
+            sim.solution.inputs["Current function [A]"], np.array([[1, 1, 2]])
         )
 
     def test_save_load(self):

tests/unit/test_solvers/test_solution.py

@@ -48,7 +48,9 @@ def test_append(self):
         np.testing.assert_array_equal(sol1.y, np.concatenate([y1, y2[:, 1:]], axis=1))
         np.testing.assert_array_equal(
             sol1.inputs["a"],
-            np.concatenate([1 * np.ones_like(t1), 2 * np.ones_like(t2[1:])]),
+            np.concatenate([1 * np.ones_like(t1), 2 * np.ones_like(t2[1:])])[
+                np.newaxis, :
+            ],
         )
 
         # Test sub-solutions
@@ -57,11 +59,11 @@ def test_append(self):
         np.testing.assert_array_equal(sol1.sub_solutions[1].t, t2)
         self.assertEqual(sol1.sub_solutions[0].model, sol1.model)
         np.testing.assert_array_equal(
-            sol1.sub_solutions[0].inputs["a"], 1 * np.ones_like(t1)
+            sol1.sub_solutions[0].inputs["a"], 1 * np.ones_like(t1)[np.newaxis, :]
         )
         self.assertEqual(sol1.sub_solutions[1].model, sol2.model)
         np.testing.assert_array_equal(
-            sol1.sub_solutions[1].inputs["a"], 2 * np.ones_like(t2)
+            sol1.sub_solutions[1].inputs["a"], 2 * np.ones_like(t2)[np.newaxis, :]
         )
 
     def test_total_time(self):