Repository: pybamm-team/PyBaMM

Issue:
Solution[<variable>] doesn't work when some parameters are inputs
**Describe the bug**
The inputs are not making their way to all the things that need to be evaluated somehow. The example below will work if current is the only input but fails when electrode height is one. Something to do with mesh

**To Reproduce**
```
import pybamm
import numpy as np

pybamm.set_logging_level('DEBUG')

I_typical = 1.0

def current_function(t):
    return pybamm.InputParameter("Current")

model = pybamm.lithium_ion.SPM()
geometry = model.default_geometry
param = model.default_parameter_values
param.update(
    {
        "Typical current [A]": I_typical,
        "Current function [A]": current_function,
        "Current": "[input]",
        "Electrode height [m]": "[input]",
    }
)
param.process_model(model)
param.process_geometry(geometry)
var = pybamm.standard_spatial_vars
var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 10, var.r_p: 10}
spatial_methods = model.default_spatial_methods

solver = model.default_solver
sim = pybamm.Simulation(
    model=model,
    geometry=geometry,
    parameter_values=param,
    var_pts=var_pts,
    spatial_methods=spatial_methods,
    solver=solver,
)

inputs = {
        'Current': I_typical,
        'Electrode height [m]': 0.1
        }

sim.solve(t_eval=np.linspace(0, 0.01, 100), inputs=inputs)
print(sim.solution['Time [h]'](sim.solution.t))
```



Discussion:

Patch:
diff --git a/pybamm/processed_variable.py b/pybamm/processed_variable.py
index 847fdef80..96f97f18e 100644
--- a/pybamm/processed_variable.py
+++ b/pybamm/processed_variable.py
@@ -93,7 +93,7 @@ class ProcessedVariable(object):
         for idx in range(len(self.t_sol)):
             t = self.t_sol[idx]
             u = self.u_sol[:, idx]
-            inputs = {name: inp[0] 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, known_evals=self.known_evals[t]
@@ -117,7 +117,7 @@ class ProcessedVariable(object):
         for idx in range(len(self.t_sol)):
             t = self.t_sol[idx]
             u = self.u_sol[:, idx]
-            inputs = {name: inp[0] 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, known_evals=self.known_evals[t]
@@ -221,7 +221,7 @@ class ProcessedVariable(object):
         for idx in range(len(self.t_sol)):
             t = self.t_sol[idx]
             u = self.u_sol[:, idx]
-            inputs = {name: inp[0] 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, known_evals=self.known_evals[t]
@@ -258,7 +258,11 @@ class ProcessedVariable(object):
         len_z = len(z_sol)
 
         # Evaluate the base_variable
-        entries = np.reshape(self.base_variable.evaluate(0, self.u_sol), [len_y, len_z])
+        inputs = {name: inp[0] for name, inp in self.inputs.items()}
+
+        entries = np.reshape(
+            self.base_variable.evaluate(0, self.u_sol, inputs), [len_y, len_z]
+        )
 
         # assign attributes for reference
         self.entries = entries
@@ -284,15 +288,17 @@ class ProcessedVariable(object):
         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()}
+
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, self.inputs, known_evals=self.known_evals[t]
+                    t, u, inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(eval_and_known_evals[0], [len_y, len_z])
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
                 entries[:, :, idx] = np.reshape(
-                    self.base_variable.evaluate(t, u, self.inputs), [len_y, len_z]
+                    self.base_variable.evaluate(t, u, inputs), [len_y, len_z]
                 )
 
         # assign attributes for reference
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 7eafa8fb2..f7080bf11 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -225,22 +225,10 @@ class Simulation:
         elif self._solution.t[-1] == solution.t[-1]:
             pass
         else:
-            self._update_solution(solution)
+            self._solution.append(solution)
 
         self._made_first_step = True
 
-    def _update_solution(self, solution):
-
-        self._solution.set_up_time += solution.set_up_time
-        self._solution.solve_time += solution.solve_time
-        self._solution.t = np.append(self._solution.t, solution.t[-1])
-        self._solution.t_event = solution.t_event
-        self._solution.termination = solution.termination
-        self._solution.y = np.concatenate(
-            [self._solution.y, solution.y[:, -1][:, np.newaxis]], axis=1
-        )
-        self._solution.y_event = solution.y_event
-
     def get_variable_array(self, *variables):
         """
         A helper function to easily obtain a dictionary of arrays of values
diff --git a/pybamm/solvers/solution.py b/pybamm/solvers/solution.py
index 64ab97f35..b3f798544 100644
--- a/pybamm/solvers/solution.py
+++ b/pybamm/solvers/solution.py
@@ -142,6 +142,10 @@ class Solution(object):
             self.inputs[name] = np.concatenate((inp, solution_inp[1:]))
         # Update solution time
         self.solve_time += solution.solve_time
+        # Update events
+        self.termination = solution.termination
+        self.t_event = solution.t_event
+        self.y_event = solution.y_event
 
         # Update known_evals
         for t, evals in solution.known_evals.items():


Test patch:
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index fcfaf2831..ee7fa188d 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -223,6 +223,31 @@ class TestSimulation(unittest.TestCase):
         self.assertEqual(sim.solution.t[0], 2 * dt)
         self.assertEqual(sim.solution.t[1], 3 * dt)
 
+    def test_step_with_inputs(self):
+        def current_function(t):
+            return pybamm.InputParameter("Current")
+
+        dt = 0.001
+        model = pybamm.lithium_ion.SPM()
+        param = model.default_parameter_values
+        param.update({"Current function [A]": current_function, "Current": "[input]"})
+        sim = pybamm.Simulation(model, parameter_values=param)
+        sim.step(dt, inputs={"Current": 1})  # 1 step stores first two points
+        self.assertEqual(sim.solution.t.size, 2)
+        self.assertEqual(sim.solution.y[0, :].size, 2)
+        self.assertEqual(sim.solution.t[0], 0)
+        self.assertEqual(sim.solution.t[1], dt)
+        np.testing.assert_array_equal(sim.solution.inputs["Current"], 1)
+        sim.step(dt, inputs={"Current": 2})  # automatically append the next step
+        self.assertEqual(sim.solution.t.size, 3)
+        self.assertEqual(sim.solution.y[0, :].size, 3)
+        self.assertEqual(sim.solution.t[0], 0)
+        self.assertEqual(sim.solution.t[1], dt)
+        self.assertEqual(sim.solution.t[2], 2 * dt)
+        np.testing.assert_array_equal(
+            sim.solution.inputs["Current"], np.array([1, 1, 2])
+        )
+
     def test_save_load(self):
         model = pybamm.lead_acid.LOQS()
         model.use_jacobian = True