Repository: pybamm-team/PyBaMM

Issue:
Order of sensitivites
- PyBaMM version - 21.12:
- Python version - 3.8.10:

**Describe the bug**
The order of sensitivities seems to not align correctly when there are multiple sensitivities. The example below shows the sensitivity of terminal voltage SPM to D_n in two cases. Case 1 where only D_n is an input and Case 2 where D_n and b_n (negative electrolyte bruggeman coefficient) are inputs. In Case 2, it appears that the sensitivites of the terminal voltage w.r.t D_n and b_n have been switched. Seems likely a switching of order in either the building of the sensitivity rhs or when extracting from casadi. 

**To Reproduce**
```
import pybamm
import numpy as np
import matplotlib.pyplot as plt

# Case 1
chemistry = pybamm.parameter_sets.Ecker2015
parameter_values = pybamm.ParameterValues(chemistry=chemistry)
parameter_values["Negative electrode diffusivity [m2.s-1]"] = "[input]"

model = pybamm.lithium_ion.SPM()
sim = pybamm.Simulation(model, parameter_values=parameter_values)

inputs = {
    "Negative electrode diffusivity [m2.s-1]": 1e-13,
}

t_eval = np.linspace(0, 30, 1000)
sim.solve(t_eval=t_eval, inputs=inputs, calculate_sensitivities=True)
sol1 = sim.solution


# Case 2
chemistry = pybamm.parameter_sets.Ecker2015
parameter_values = pybamm.ParameterValues(chemistry=chemistry)
parameter_values["Negative electrode diffusivity [m2.s-1]"] = "[input]"
parameter_values["Negative electrode Bruggeman coefficient (electrolyte)"] = "[input]"

model = pybamm.lithium_ion.SPM()
sim = pybamm.Simulation(model, parameter_values=parameter_values)

inputs = {
    "Negative electrode diffusivity [m2.s-1]": 1e-13,
    "Negative electrode Bruggeman coefficient (electrolyte)": 1.5,
}

sim.solve(t_eval=t_eval, inputs=inputs, calculate_sensitivities=True)
sol2 = sim.solution

# Plot
fig, ax = plt.subplots(2, 1, figsize=(16, 9), sharex=True)
ax[0].plot(
    t_eval,
    sol1["Terminal voltage [V]"].sensitivities[
        "Negative electrode diffusivity [m2.s-1]"
    ],
    label="sol1 D_n",
)
ax[0].plot(
    t_eval,
    sol2["Terminal voltage [V]"].sensitivities[
        "Negative electrode diffusivity [m2.s-1]"
    ],
    linestyle="--",
    label="sol2 D_n",
)

ax[1].plot(
    t_eval,
    sol1["Terminal voltage [V]"].sensitivities[
        "Negative electrode diffusivity [m2.s-1]"
    ],
    label="sol1 D_n",
)
ax[1].plot(
    t_eval,
    sol2["Terminal voltage [V]"].sensitivities[
        "Negative electrode Bruggeman coefficient (electrolyte)"
    ],
    linestyle="--",
    label="sol2 b",
)

ax[0].grid()
ax[1].grid()

ax[0].set_xlabel("Time [s]")
ax[1].set_xlabel("Time [s]")

ax[0].legend()
ax[1].legend()

plt.show()
```


**Expected behaviour**
The sensitivity of terminal voltage w.r.t D_n should be independent of the number of other inputs. 

**Screenshots**
Plot generated by above example. 

![image](https://user-images.githubusercontent.com/22661308/148528381-ddfa43d6-2bd7-42a6-a329-c6a38944e13c.png)




Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ef7d7c5a2..29d1f3f99 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@
 
 -   Parameters can now be imported from any given path in `Windows` ([#1900](https://github.com/pybamm-team/PyBaMM/pull/1900))
 -   Fixed initial conditions for the EC SEI model ([#1895](https://github.com/pybamm-team/PyBaMM/pull/1895))
+-   Fixed issue in extraction of sensitivites ([#1894](https://github.com/pybamm-team/PyBaMM/pull/1894))
 
 # [v21.12](https://github.com/pybamm-team/PyBaMM/tree/v21.11) - 2021-12-29
 
@@ -14,7 +15,6 @@
 -   Added cylindrical geometry and finite volume method ([#1824](https://github.com/pybamm-team/PyBaMM/pull/1824))
 
 ## Bug fixes
-
 -   `PyBaMM` is now importable in `Linux` systems where `jax` is already installed ([#1874](https://github.com/pybamm-team/PyBaMM/pull/1874))
 -   Simulations with drive cycles now support `initial_soc` ([#1842](https://github.com/pybamm-team/PyBaMM/pull/1842))
 -   Fixed bug in expression tree simplification ([#1831](https://github.com/pybamm-team/PyBaMM/pull/1831))
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 98fa02dd5..1ed978516 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -1459,7 +1459,12 @@ class BaseSolver(object):
                 inputs[name] = casadi.MX.sym(name, input_param._expected_size)
 
         external_variables = external_variables or {}
-        ext_and_inputs = {**external_variables, **inputs}
+
+        ordered_inputs_names = list(inputs.keys())
+        ordered_inputs_names.sort()
+        ordered_inputs = {name: inputs[name] for name in ordered_inputs_names}
+
+        ext_and_inputs = {**external_variables, **ordered_inputs}
         return ext_and_inputs
 
 
diff --git a/pybamm/solvers/processed_variable.py b/pybamm/solvers/processed_variable.py
index 7261a9aaa..20a615020 100644
--- a/pybamm/solvers/processed_variable.py
+++ b/pybamm/solvers/processed_variable.py
@@ -351,10 +351,14 @@ class ProcessedVariable(object):
             self.second_dimension = "x"
             self.R_sol = first_dim_pts
             self.x_sol = second_dim_pts
-        elif self.domain[0] in [
-            "negative particle size",
-            "positive particle size",
-        ] and self.auxiliary_domains["secondary"] == ["current collector"]:
+        elif (
+            self.domain[0]
+            in [
+                "negative particle size",
+                "positive particle size",
+            ]
+            and self.auxiliary_domains["secondary"] == ["current collector"]
+        ):
             self.first_dimension = "R"
             self.second_dimension = "z"
             self.R_sol = first_dim_pts
@@ -563,6 +567,7 @@ class ProcessedVariable(object):
             name: casadi.MX.sym(name, value.shape[0])
             for name, value in self.all_inputs[0].items()
         }
+
         p_casadi_stacked = casadi.vertcat(*[p for p in p_casadi.values()])
 
         # Convert variable to casadi format for differentiating


Test patch:
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 4678e7fec..1b4444511 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -589,9 +589,10 @@ class TestCasadiSolverODEsWithForwardSensitivityEquations(unittest.TestCase):
         solution = solver.solve(
             model,
             t_eval,
-            inputs={"p": 0.1, "q": 2, "r": -1, "s": 0.5},
+            inputs={"r": -1, "s": 0.5, "q": 2, "p": 0.1},
             calculate_sensitivities=True,
         )
+
         np.testing.assert_allclose(solution.y[0], -1 + 0.2 * solution.t)
         np.testing.assert_allclose(
             solution.sensitivities["p"],