Skip to content

Commit b112fd7

Browse files
martinjrobinsmartinjrobins
authored
Merge pull request #1897 from pybamm-team/i1887-sens-subset
fix for when only a subset of input params are used for sensitivities
2 parents 524839d + b8e097e commit b112fd7

File tree

4 files changed

+131
-24
lines changed

4 files changed

+131
-24
lines changed

‎CHANGELOG.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
12
# [Unreleased](https://github.com/pybamm-team/PyBaMM/)
23

34
## Features
@@ -7,6 +8,7 @@
78

89
## Bug fixes
910

11+
- Fix bug where sensitivity calculation failed if len of `calculate_sensitivities` was less than `inputs` ([#1897](https://github.com/pybamm-team/PyBaMM/pull/1897))
1012
- Fixed a bug in the eSOH variable calculation when OCV is given as data ([#1975](https://github.com/pybamm-team/PyBaMM/pull/1975))
1113
- Fixed a bug where isothermal models did not compute any heat source terms ([#1958](https://github.com/pybamm-team/PyBaMM/pull/1958))
1214

@@ -50,6 +52,7 @@
5052
- Parameters can now be imported from any given path in `Windows` ([#1900](https://github.com/pybamm-team/PyBaMM/pull/1900))
5153
- Fixed initial conditions for the EC SEI model ([#1895](https://github.com/pybamm-team/PyBaMM/pull/1895))
5254
- Fixed issue in extraction of sensitivites ([#1894](https://github.com/pybamm-team/PyBaMM/pull/1894))
55+
>>>>>>> develop
5356
5457
# [v21.12](https://github.com/pybamm-team/PyBaMM/tree/v21.11) - 2021-12-29
5558

‎pybamm/solvers/base_solver.py‎

Lines changed: 26 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -248,7 +248,32 @@ def set_up(self, model, inputs=None, t_eval=None, ics_only=False):
248248
y_and_S = y_casadi
249249

250250
# if we will change the equations to include the explicit sensitivity
251-
# equations, then we also need to update the mass matrix and bounds
251+
# equations, then we also need to update the mass matrix and bounds.
252+
# First, we reset the mass matrix and bounds back to their original form
253+
# if they have been extended
254+
if model.bounds[0].shape[0] > model.len_rhs_and_alg:
255+
model.bounds = (
256+
model.bounds[0][: model.len_rhs_and_alg],
257+
model.bounds[1][: model.len_rhs_and_alg],
258+
)
259+
if (
260+
model.mass_matrix is not None
261+
and model.mass_matrix.shape[0] > model.len_rhs_and_alg
262+
):
263+
if model.mass_matrix_inv is not None:
264+
model.mass_matrix_inv = pybamm.Matrix(
265+
model.mass_matrix_inv.entries[
266+
: model.len_rhs, : model.len_rhs
267+
]
268+
)
269+
model.mass_matrix = pybamm.Matrix(
270+
model.mass_matrix.entries[
271+
: model.len_rhs_and_alg, : model.len_rhs_and_alg
272+
]
273+
)
274+
275+
# now we can extend them by the number of sensitivity parameters
276+
# if needed
252277
if calculate_sensitivities_explicit:
253278
if model.len_rhs != 0:
254279
n_inputs = model.len_rhs_sens // model.len_rhs
@@ -276,28 +301,6 @@ def set_up(self, model, inputs=None, t_eval=None, ics_only=False):
276301
[model.mass_matrix.entries] * (n_inputs + 1), format="csr"
277302
)
278303
)
279-
else:
280-
# take care if calculate_sensitivites used then not used
281-
if model.bounds[0].shape[0] > model.len_rhs_and_alg:
282-
model.bounds = (
283-
model.bounds[0][: model.len_rhs_and_alg],
284-
model.bounds[1][: model.len_rhs_and_alg],
285-
)
286-
if (
287-
model.mass_matrix is not None
288-
and model.mass_matrix.shape[0] > model.len_rhs_and_alg
289-
):
290-
if model.mass_matrix_inv is not None:
291-
model.mass_matrix_inv = pybamm.Matrix(
292-
model.mass_matrix_inv.entries[
293-
: model.len_rhs, : model.len_rhs
294-
]
295-
)
296-
model.mass_matrix = pybamm.Matrix(
297-
model.mass_matrix.entries[
298-
: model.len_rhs_and_alg, : model.len_rhs_and_alg
299-
]
300-
)
301304

302305
def process(symbol, name, use_jacobian=None):
303306
def report(string):

‎pybamm/solvers/casadi_algebraic_solver.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -83,7 +83,7 @@ def _integrate(self, model, t_eval, inputs_dict=None):
8383
if model.len_rhs_and_alg == y0.shape[0]:
8484
len_rhs = model.len_rhs
8585
else:
86-
len_rhs = model.len_rhs * (inputs.shape[0] + 1)
86+
len_rhs = model.len_rhs + model.len_rhs_sens
8787
y0_diff = y0[:len_rhs]
8888
y0_alg = y0[len_rhs:]
8989

‎tests/unit/test_solvers/test_casadi_solver.py‎

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -865,6 +865,64 @@ def test_solve_sensitivity_then_no_sensitivity(self):
865865
solution["var squared"].data, np.exp(0.1 * solution.t) ** 2
866866
)
867867

868+
def test_solve_sensitivity_subset(self):
869+
# Create model
870+
model = pybamm.BaseModel()
871+
var = pybamm.Variable("var")
872+
p = pybamm.InputParameter("p")
873+
q = pybamm.InputParameter("q")
874+
r = pybamm.InputParameter("r")
875+
model.rhs = {var: p * q}
876+
model.initial_conditions = {var: r}
877+
878+
# only calculate the sensitivities of a subset of parameters
879+
solver = pybamm.CasadiSolver(rtol=1e-10, atol=1e-10)
880+
t_eval = np.linspace(0, 1, 80)
881+
solution = solver.solve(
882+
model,
883+
t_eval,
884+
inputs={"q": 2, "r": -1, "p": 0.1},
885+
calculate_sensitivities=["q", "p"],
886+
)
887+
np.testing.assert_allclose(solution.y[0], -1 + 0.2 * solution.t)
888+
np.testing.assert_allclose(
889+
solution.sensitivities["p"],
890+
(2 * solution.t)[:, np.newaxis],
891+
)
892+
np.testing.assert_allclose(
893+
solution.sensitivities["q"],
894+
(0.1 * solution.t)[:, np.newaxis],
895+
)
896+
self.assertTrue("r" not in solution.sensitivities)
897+
np.testing.assert_allclose(
898+
solution.sensitivities["all"],
899+
np.hstack(
900+
[
901+
solution.sensitivities["p"],
902+
solution.sensitivities["q"],
903+
]
904+
),
905+
)
906+
907+
solution = solver.solve(
908+
model,
909+
t_eval,
910+
inputs={"q": 2, "r": -1, "p": 0.1},
911+
calculate_sensitivities=["r"],
912+
)
913+
np.testing.assert_allclose(solution.y[0], -1 + 0.2 * solution.t)
914+
self.assertTrue("p" not in solution.sensitivities)
915+
self.assertTrue("q" not in solution.sensitivities)
916+
np.testing.assert_allclose(solution.sensitivities["r"], 1)
917+
np.testing.assert_allclose(
918+
solution.sensitivities["all"],
919+
np.hstack(
920+
[
921+
solution.sensitivities["r"],
922+
]
923+
),
924+
)
925+
868926

869927
class TestCasadiSolverDAEsWithForwardSensitivityEquations(unittest.TestCase):
870928
def test_solve_sensitivity_scalar_var_scalar_input(self):
@@ -939,6 +997,49 @@ def test_solve_sensitivity_algebraic(self):
939997
atol=1e-7,
940998
)
941999

1000+
def test_solve_sensitivity_subset(self):
1001+
# Create model
1002+
model = pybamm.BaseModel()
1003+
var = pybamm.Variable("var")
1004+
var2 = pybamm.Variable("var2")
1005+
p = pybamm.InputParameter("p")
1006+
q = pybamm.InputParameter("q")
1007+
r = pybamm.InputParameter("r")
1008+
model.rhs = {var: p * q}
1009+
model.algebraic = {var2: 2 * var - var2}
1010+
model.initial_conditions = {var: r, var2: 2 * r}
1011+
1012+
# only calculate the sensitivities of a subset of parameters
1013+
solver = pybamm.CasadiSolver(rtol=1e-10, atol=1e-10)
1014+
t_eval = np.linspace(0, 1, 80)
1015+
solution = solver.solve(
1016+
model,
1017+
t_eval,
1018+
inputs={"p": 0.1, "q": 2, "r": -1, "s": 0.5},
1019+
calculate_sensitivities=["p", "q"],
1020+
)
1021+
np.testing.assert_allclose(solution.y[0], -1 + 0.2 * solution.t)
1022+
np.testing.assert_allclose(solution.y[-1], 2 * (-1 + 0.2 * solution.t))
1023+
np.testing.assert_allclose(
1024+
solution.sensitivities["p"][::2],
1025+
(2 * solution.t)[:, np.newaxis],
1026+
)
1027+
np.testing.assert_allclose(
1028+
solution.sensitivities["q"][::2],
1029+
(0.1 * solution.t)[:, np.newaxis],
1030+
)
1031+
self.assertTrue("r" not in solution.sensitivities)
1032+
self.assertTrue("s" not in solution.sensitivities)
1033+
np.testing.assert_allclose(
1034+
solution.sensitivities["all"],
1035+
np.hstack(
1036+
[
1037+
solution.sensitivities["p"],
1038+
solution.sensitivities["q"],
1039+
]
1040+
),
1041+
)
1042+
9421043

9431044
if __name__ == "__main__":
9441045
print("Add -v for more debug output")

0 commit comments

Comments
 (0)