Repository: pybamm-team/PyBaMM (python)

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))
```



Pull request description:
#793 fix simulation stepping
# Description

Fixes solution appending in the simulation

Fixes #793 

## 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. 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)
- [x] Bug fix (non-breaking change which fixes an issue)


# Key checklist:

- [x] No style issues: `$ flake8`
- [x] All tests pass: `$ python run-tests.py --unit`
- [x] The documentation builds: `$ cd docs` and then `$ make clean; make html`

You can run all three at once, using `$ python run-tests.py --quick`.

## Further checks:

- [x] Code is commented, particularly in hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works


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