I3790 idaklu outvars

CHANGELOG.md

@@ -10,6 +10,7 @@
 ## Bug Fixes
 
 - Fixed a bug where 1+1D and 2+1D models would not work with voltage or power controlled experiments([#3829](https://github.com/pybamm-team/PyBaMM/pull/3829))
+- Update IDAKLU solver to fail gracefully when a variable is requested that was not in the solves `output_variables` list ([#3803](https://github.com/pybamm-team/PyBaMM/pull/3803))
 - Updated `_steps_util.py` to throw a specific exception when drive cycle starts at t>0 ([#3756](https://github.com/pybamm-team/PyBaMM/pull/3756))
 - Updated `plot_voltage_components.py` to support both `Simulation` and `Solution` objects. Added new methods in both `Simulation` and `Solution` classes for allow the syntax `simulation.plot_voltage_components` and `solution.plot_voltage_components`. Updated `test_plot_voltage_components.py` to reflect these changes ([#3723](https://github.com/pybamm-team/PyBaMM/pull/3723)).
 - The SEI thickness decreased at some intervals when the 'electron-migration limited' model was used. It has been corrected ([#3622](https://github.com/pybamm-team/PyBaMM/pull/3622))

pybamm/solvers/solution.py

@@ -478,7 +478,15 @@ def update(self, variables):
             for i, (model, ys, inputs, var_pybamm) in enumerate(
                 zip(self.all_models, self.all_ys, self.all_inputs, vars_pybamm)
             ):
-                if isinstance(var_pybamm, pybamm.ExplicitTimeIntegral):
+                if ys.size == 0 and var_pybamm.has_symbol_of_classes(
+                    pybamm.expression_tree.state_vector.StateVector
+                ):
+                    raise KeyError(
+                        f"Cannot process variable '{key}' as it was not part of the "
+                        "solve. Please re-run the solve with `output_variables` set to "
+                        "include this variable."
+                    )
+                elif isinstance(var_pybamm, pybamm.ExplicitTimeIntegral):
                     cumtrapz_ic = var_pybamm.initial_condition
                     cumtrapz_ic = cumtrapz_ic.evaluate()
                     var_pybamm = var_pybamm.child

tests/unit/test_solvers/test_idaklu_solver.py

@@ -550,21 +550,23 @@ def test_with_output_variables(self):
         # Construct a model and solve for all vairables, then test
         # the 'output_variables' option for each variable in turn, confirming
         # equivalence
-
-        # construct model
-        model = pybamm.lithium_ion.DFN()
-        geometry = model.default_geometry
-        param = model.default_parameter_values
         input_parameters = {}  # Sensitivities dictionary
-        param.update({key: "[input]" for key in input_parameters})
-        param.process_model(model)
-        param.process_geometry(geometry)
-        var_pts = {"x_n": 50, "x_s": 50, "x_p": 50, "r_n": 5, "r_p": 5}
-        mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-        disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-        disc.process_model(model)
         t_eval = np.linspace(0, 3600, 100)
 
+        # construct model
+        def construct_model():
+            model = pybamm.lithium_ion.DFN()
+            geometry = model.default_geometry
+            param = model.default_parameter_values
+            param.update({key: "[input]" for key in input_parameters})
+            param.process_model(model)
+            param.process_geometry(geometry)
+            var_pts = {"x_n": 50, "x_s": 50, "x_p": 50, "r_n": 5, "r_p": 5}
+            mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
+            disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+            disc.process_model(model)
+            return model
+
         options = {
             "linear_solver": "SUNLinSol_KLU",
             "jacobian": "sparse",
@@ -585,14 +587,32 @@ def test_with_output_variables(self):
             "Throughput capacity [A.h]",  # ExplicitTimeIntegral
         ]
 
+        # vars that are not in the output_variables list, but are still accessible as
+        # they are either model parameters, or do not require access to the state vector
+        model_vars = [
+            "Time [s]",
+            "C-rate",
+            "Ambient temperature [K]",
+            "Porosity",
+        ]
+
+        # A list of variables that are not in the model and cannot be computed
+        inaccessible_vars = [
+            "Terminal voltage [V]",
+            "Negative particle surface stoichiometry",
+            "Electrode current density [A.m-2]",
+            "Power [W]",
+            "Resistance [Ohm]",
+        ]
+
         # Use the full model as comparison (tested separately)
         solver_all = pybamm.IDAKLUSolver(
             atol=1e-8,
             rtol=1e-8,
             options=options,
         )
         sol_all = solver_all.solve(
-            model,
+            construct_model(),
             t_eval,
             inputs=input_parameters,
             calculate_sensitivities=True,
@@ -606,15 +626,20 @@ def test_with_output_variables(self):
             output_variables=output_variables,
         )
         sol = solver.solve(
-            model,
+            construct_model(),
             t_eval,
             inputs=input_parameters,
         )
 
         # Compare output to sol_all
-        for varname in output_variables:
+        for varname in [*output_variables, *model_vars]:
             self.assertTrue(np.allclose(sol[varname].data, sol_all[varname].data))
 
+        # Check that the missing variables are not available in the solution
+        for varname in inaccessible_vars:
+            with self.assertRaises(KeyError):
+                sol[varname].data
+
         # Mock a 1D current collector and initialise (none in the model)
         sol["x_s [m]"].domain = ["current collector"]
         sol["x_s [m]"].initialise_1D()