diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py index 4df4cbf0fb..9689d80f6b 100644 --- a/pybamm/solvers/idaklu_solver.py +++ b/pybamm/solvers/idaklu_solver.py @@ -321,10 +321,11 @@ def sensfn(resvalS, t, y, yp, yS, ypS): number_of_states = y0.size y_out = sol.y.reshape((number_of_timesteps, number_of_states)) - # return solution, we need to tranpose y to match scipy's interface + # return sensitivity solution, we need to flatten yS to + # (#timesteps * #states,) to match format used by Solution if number_of_sensitivity_parameters != 0: yS_out = { - name: sol.yS[i].transpose() for i, name in enumerate(sens0.keys()) + name: sol.yS[i].reshape(-1, 1) for i, name in enumerate(sens0.keys()) } else: yS_out = False diff --git a/pybamm/solvers/solution.py b/pybamm/solvers/solution.py index abc1e52a05..40ef2fc635 100644 --- a/pybamm/solvers/solution.py +++ b/pybamm/solvers/solution.py @@ -81,23 +81,10 @@ def __init__( else: self.all_inputs = all_inputs - # sensitivities - if isinstance(sensitivities, bool): - self._sensitivities = {} - # if solution consists of explicit sensitivity equations, extract them - if ( - sensitivities is True - and all_models[0] is not None - and not isinstance(all_ys[0], casadi.Function) - and all_models[0].len_rhs_and_alg != all_ys[0].shape[0] - and all_models[0].len_rhs_and_alg != 0 # for the dummy solver - ): - self.extract_explicit_sensitivities() - - elif isinstance(sensitivities, dict): - self._sensitivities = sensitivities - else: + # sensitivities must be a dict or bool + if not isinstance(sensitivities, (bool, dict)): raise TypeError('sensitivities arg needs to be a bool or dict') + self._sensitivities = sensitivities self._t_event = t_event self._y_event = y_event @@ -148,23 +135,10 @@ def __init__( pybamm.citations.register("Andersson2019") def extract_explicit_sensitivities(self): - for index, (model, ys, ts, inputs) in enumerate( - zip(self.all_models, self.all_ys, self.all_ts, - self.all_inputs) - ): - # TODO: only support sensitivities for one solution atm - # but make sure that sensitivities are removed for all - # solutions - if index == 0: - self._all_ys[index], self._sensitivities = \ - self._extract_explicit_sensitivities( - model, ys, ts, inputs - ) - else: - self._all_ys[index], _ = \ - self._extract_explicit_sensitivities( - model, ys, ts, inputs - ) + self._y, self._sensitivities = \ + self._extract_explicit_sensitivities( + self.all_models[0], self.y, self.t, self.all_inputs[0] + ) def _extract_explicit_sensitivities(self, model, y, t_eval, inputs): """ @@ -277,11 +251,18 @@ def y(self): return self._y except AttributeError: self.set_y() + + # if y is evaluated before sensitivities then need to extract them + if isinstance(self._sensitivities, bool) and self._sensitivities: + self.extract_explicit_sensitivities() + return self._y @property def sensitivities(self): """Values of the sensitivities. Returns a dict of param_name: np_array""" + if isinstance(self._sensitivities, bool) and self._sensitivities: + self.extract_explicit_sensitivities() return self._sensitivities def set_y(self): diff --git a/tests/integration/test_models/standard_model_tests.py b/tests/integration/test_models/standard_model_tests.py index 839e60f1d0..fa798c94ad 100644 --- a/tests/integration/test_models/standard_model_tests.py +++ b/tests/integration/test_models/standard_model_tests.py @@ -83,7 +83,6 @@ def test_solving(self, solver=None, t_eval=None, inputs=None, self.solution = self.solver.solve( self.model, t_eval, inputs=inputs, - calculate_sensitivities=calculate_sensitivities ) def test_outputs(self): @@ -93,40 +92,45 @@ def test_outputs(self): ) std_out_test.test_all() - def test_sensitivities(self): - param_name = "Negative electrode conductivity [S.m-1]" - neg_electrode_cond = 100.0 + def test_sensitivities(self, param_name, param_value): self.parameter_values.update({param_name: "[input]"}) - inputs = {param_name: neg_electrode_cond} + inputs = {param_name: param_value} self.test_processing_parameters() self.test_processing_disc() - self.test_solving(inputs=inputs, calculate_sensitivities=True) + # Use tighter default tolerances for testing + self.solver.rtol = 1e-8 + self.solver.atol = 1e-8 + + Crate = abs( + self.parameter_values["Current function [A]"] + / self.parameter_values["Nominal cell capacity [A.h]"] + ) + t_eval = np.linspace(0, 3600 / Crate, 100) + + self.solution = self.solver.solve( + self.model, t_eval, inputs=inputs, + calculate_sensitivities=True + ) # check via finite differencing - h = 1e-6 - inputs_plus = {param_name: neg_electrode_cond + 0.5 * h} - inputs_neg = {param_name: neg_electrode_cond - 0.5 * h} + h = 1e-6 * param_value + inputs_plus = {param_name: (param_value + 0.5 * h)} + inputs_neg = {param_name: (param_value - 0.5 * h)} sol_plus = self.solver.solve( - self.model, self.solution.all_ts[0], inputs=inputs_plus + self.model, t_eval, inputs=inputs_plus, ) sol_neg = self.solver.solve( - self.model, self.solution.all_ts[0], inputs=inputs_neg + self.model, t_eval, inputs=inputs_neg ) - n = self.solution.sensitivities[param_name].shape[0] - np.testing.assert_array_almost_equal( - self.solution.sensitivities[param_name], - ((sol_plus.y - sol_neg.y) / h).reshape((n, 1)) + fd = ((np.array(sol_plus.y) - np.array(sol_neg.y)) / h) + fd = fd.transpose().reshape(-1, 1) + np.testing.assert_allclose( + self.solution.sensitivities[param_name], fd, + rtol=1e-1, atol=1e-5, ) - if ( - isinstance( - self.model, (pybamm.lithium_ion.BaseModel, pybamm.lead_acid.BaseModel) - ) - ): - self.test_outputs() - def test_all( self, param=None, disc=None, solver=None, t_eval=None, skip_output_tests=False ): diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py index 48eca49167..732273f27d 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py @@ -31,7 +31,10 @@ def test_sensitivities(self): modeltest = tests.StandardModelTest( model, parameter_values=param, var_pts=var_pts ) - modeltest.test_sensitivities() + modeltest.test_sensitivities( + #'Separator thickness [m]', 2e-05, + 'Typical current [A]', 0.15652, + ) def test_basic_processing_1plus1D(self): options = {"current collector": "potential pair", "dimensionality": 1} diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py index 81bbe466d4..0bd8d1f5e3 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py @@ -20,10 +20,18 @@ def test_basic_processing(self): def test_sensitivities(self): options = {"thermal": "isothermal"} model = pybamm.lithium_ion.SPM(options) - # use Ecker parameters for nonlinear diffusion param = pybamm.ParameterValues(chemistry=pybamm.parameter_sets.Ecker2015) modeltest = tests.StandardModelTest(model, parameter_values=param) - modeltest.test_sensitivities() + modeltest.test_sensitivities( + #"Negative electrode conductivity [S.m-1]", 14.0 + 'Typical current [A]', 0.15652, + #"Typical electrolyte concentration [mol.m-3]", 1000.0 + #''Negative electrode diffusivity [m2.s-1]', 1e-3, + #'Negative electrode active material volume fraction', 0.372403, + #'Separator thickness [m]', 2e-05, + #'Negative electrode electrons in reaction', 1.0, + #'Outer SEI open-circuit potential [V]', 0.8, + ) def test_basic_processing_1plus1D(self): options = {"current collector": "potential pair", "dimensionality": 1} diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py index 42db8d9c79..c5e651a00e 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py @@ -24,7 +24,10 @@ def test_sensitivities(self): # use Ecker parameters for nonlinear diffusion param = pybamm.ParameterValues(chemistry=pybamm.parameter_sets.Ecker2015) modeltest = tests.StandardModelTest(model, parameter_values=param) - modeltest.test_sensitivities() + modeltest.test_sensitivities( + #'Separator thickness [m]', 2e-05, + 'Typical current [A]', 0.15652, + ) def test_basic_processing_python(self): options = {"thermal": "isothermal"}