diff --git a/CHANGELOG.md b/CHANGELOG.md index ea458b4b60..a8d7d155de 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- Added `plot_thermal_components` to plot the contributions to the total heat generation in a battery ([#4021](https://github.com/pybamm-team/PyBaMM/pull/4021)) - Added functions for normal probability density function (`pybamm.normal_pdf`) and cumulative distribution function (`pybamm.normal_cdf`) ([#3999](https://github.com/pybamm-team/PyBaMM/pull/3999)) - Updates multiprocess `Pool` in `BaseSolver.solve()` to be constructed with context `fork`. Adds small example for multiprocess inputs. ([#3974](https://github.com/pybamm-team/PyBaMM/pull/3974)) - Added custom experiment steps ([#3835](https://github.com/pybamm-team/PyBaMM/pull/3835)) @@ -14,6 +15,7 @@ ## Bug Fixes +- Set the `remove_independent_variables_from_rhs` to `False` by default, and moved the option from `Discretisation.process_model` to `Discretisation.__init__`. This fixes a bug related to the discharge capacity, but may make the simulation slower in some cases. To set the option to `True`, use `Simulation(..., discretisation_kwargs={"remove_independent_variables_from_rhs": True})`. ([#4020](https://github.com/pybamm-team/PyBaMM/pull/4020)) - Fixed a bug where independent variables were removed from models even if they appeared in events ([#4019](https://github.com/pybamm-team/PyBaMM/pull/4019)) - Fix bug with upwind and downwind schemes producing the wrong discretised system ([#3979](https://github.com/pybamm-team/PyBaMM/pull/3979)) - Allow evaluation of an `Interpolant` object with a number ([#3932](https://github.com/pybamm-team/PyBaMM/pull/3932)) @@ -33,6 +35,7 @@ ## Breaking changes +- Removed `check_model` argument from `Simulation.solve`. To change the `check_model` option, use `Simulation(..., discretisation_kwargs={"check_model": False})`. ([#4020](https://github.com/pybamm-team/PyBaMM/pull/4020)) - Removed multiple Docker images. Here on, a single Docker image tagged `pybamm/pybamm:latest` will be provided with both solvers (`IDAKLU` and `JAX`) pre-installed. ([#3992](https://github.com/pybamm-team/PyBaMM/pull/3992)) - Removed support for Python 3.8 ([#3961](https://github.com/pybamm-team/PyBaMM/pull/3961)) - Renamed "ocp_soc_0_dimensional" to "ocp_soc_0" and "ocp_soc_100_dimensional" to "ocp_soc_100" ([#3942](https://github.com/pybamm-team/PyBaMM/pull/3942)) diff --git a/pybamm/batch_study.py b/pybamm/batch_study.py index 1503c3048c..e854c94e00 100644 --- a/pybamm/batch_study.py +++ b/pybamm/batch_study.py @@ -102,7 +102,6 @@ def solve( self, t_eval=None, solver=None, - check_model=True, save_at_cycles=None, calc_esoh=True, starting_solution=None, @@ -158,7 +157,6 @@ def solve( sol = sim.solve( t_eval, solver, - check_model, save_at_cycles, calc_esoh, starting_solution, diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py index 92e9e4e24f..d61d8a4d60 100644 --- a/pybamm/discretisations/discretisation.py +++ b/pybamm/discretisations/discretisation.py @@ -26,14 +26,32 @@ class Discretisation: Parameters ---------- mesh : pybamm.Mesh - contains all submeshes to be used on each domain + contains all submeshes to be used on each domain spatial_methods : dict - a dictionary of the spatial methods to be used on each - domain. The keys correspond to the model domains and the - values to the spatial method. + a dictionary of the spatial methods to be used on each + domain. The keys correspond to the model domains and the + values to the spatial method. + check_model : bool, optional + If True, model checks are performed after discretisation. For large + systems these checks can be slow, so can be skipped by setting this + option to False. When developing, testing or debugging it is recommended + to leave this option as True as it may help to identify any errors. + Default is True. + remove_independent_variables_from_rhs : bool, optional + If True, model checks to see whether any variables from the RHS are used + in any other equation. If a variable meets all of the following criteria + (not used anywhere in the model, len(rhs)>1), then the variable + is moved to be explicitly integrated when called by the solution object. + Default is False. """ - def __init__(self, mesh=None, spatial_methods=None): + def __init__( + self, + mesh=None, + spatial_methods=None, + check_model=True, + remove_independent_variables_from_rhs=False, + ): self._mesh = mesh if mesh is None: self._spatial_methods = {} @@ -60,6 +78,10 @@ def __init__(self, mesh=None, spatial_methods=None): self._bcs = {} self.y_slices = {} self._discretised_symbols = {} + self._check_model_flag = check_model + self._remove_independent_variables_from_rhs_flag = ( + remove_independent_variables_from_rhs + ) @property def mesh(self): @@ -90,13 +112,7 @@ def bcs(self, value): # reset discretised_symbols self._discretised_symbols = {} - def process_model( - self, - model, - inplace=True, - check_model=True, - remove_independent_variables_from_rhs=True, - ): + def process_model(self, model, inplace=True): """ Discretise a model. Currently inplace, could be changed to return a new model. @@ -108,18 +124,6 @@ def process_model( inplace : bool, optional If True, discretise the model in place. Otherwise, return a new discretised model. Default is True. - check_model : bool, optional - If True, model checks are performed after discretisation. For large - systems these checks can be slow, so can be skipped by setting this - option to False. When developing, testing or debugging it is recommended - to leave this option as True as it may help to identify any errors. - Default is True. - remove_independent_variables_from_rhs : bool, optional - If True, model checks to see whether any variables from the RHS are used - in any other equation. If a variable meets all of the following criteria - (not used anywhere in the model, len(rhs)>1), then the variable - is moved to be explicitly integrated when called by the solution object. - Default is True. Returns ------- @@ -158,7 +162,7 @@ def process_model( # set variables (we require the full variable not just id) # Search Equations for Independence - if remove_independent_variables_from_rhs: + if self._remove_independent_variables_from_rhs_flag: model = self.remove_independent_variables_from_rhs(model) variables = list(model.rhs.keys()) + list(model.algebraic.keys()) # Find those RHS's that are constant @@ -240,7 +244,7 @@ def process_model( model_disc._geometry = getattr(self.mesh, "_geometry", None) # Check that resulting model makes sense - if check_model: + if self._check_model_flag: pybamm.logger.verbose(f"Performing model checks for {model.name}") self.check_model(model_disc) @@ -1123,6 +1127,7 @@ def remove_independent_variables_from_rhs(self, model): # in variables twice under different names for key in model.variables: if model.variables[key] == var: + print("here") model.variables[key] = model.variables[var.name] del model.rhs[var] del model.initial_conditions[var] diff --git a/pybamm/simulation.py b/pybamm/simulation.py index 060d1cdd8c..a81f7b62b2 100644 --- a/pybamm/simulation.py +++ b/pybamm/simulation.py @@ -63,6 +63,9 @@ class Simulation: A list of variables to plot automatically C_rate: float (optional) The C-rate at which you would like to run a constant current (dis)charge. + discretisation_kwargs: dict (optional) + Any keyword arguments to pass to the Discretisation class. + See :class:`pybamm.Discretisation` for details. """ def __init__( @@ -77,6 +80,7 @@ def __init__( solver=None, output_variables=None, C_rate=None, + discretisation_kwargs=None, ): self._parameter_values = parameter_values or model.default_parameter_values self._unprocessed_parameter_values = self._parameter_values @@ -126,6 +130,7 @@ def __init__( self._spatial_methods = spatial_methods or self._model.default_spatial_methods self._solver = solver or self._model.default_solver self._output_variables = output_variables + self._discretisation_kwargs = discretisation_kwargs or {} # Initialize empty built states self._model_with_set_params = None @@ -268,7 +273,7 @@ def set_initial_soc(self, initial_soc, inputs=None): # Save solved initial SOC in case we need to re-build the model self._built_initial_soc = initial_soc - def build(self, check_model=True, initial_soc=None, inputs=None): + def build(self, initial_soc=None, inputs=None): """ A method to build the model into a system of matrices and vectors suitable for performing numerical computations. If the model has already been built or @@ -278,13 +283,12 @@ def build(self, check_model=True, initial_soc=None, inputs=None): Parameters ---------- - check_model : bool, optional - If True, model checks are performed after discretisation (see - :meth:`pybamm.Discretisation.process_model`). Default is True. initial_soc : float, optional Initial State of Charge (SOC) for the simulation. Must be between 0 and 1. If given, overwrites the initial concentrations provided in the parameter set. + inputs : dict, optional + A dictionary of input parameters to pass to the model when solving. """ if initial_soc is not None: self.set_initial_soc(initial_soc, inputs=inputs) @@ -297,14 +301,16 @@ def build(self, check_model=True, initial_soc=None, inputs=None): else: self.set_parameters() self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts) - self._disc = pybamm.Discretisation(self._mesh, self._spatial_methods) + self._disc = pybamm.Discretisation( + self._mesh, self._spatial_methods, **self._discretisation_kwargs + ) self._built_model = self._disc.process_model( - self._model_with_set_params, inplace=False, check_model=check_model + self._model_with_set_params, inplace=False ) # rebuilt model so clear solver setup self._solver._model_set_up = {} - def build_for_experiment(self, check_model=True, initial_soc=None, inputs=None): + def build_for_experiment(self, initial_soc=None, inputs=None): """ Similar to :meth:`Simulation.build`, but for the case of simulating an experiment, where there may be several models and solvers to build. @@ -322,7 +328,9 @@ def build_for_experiment(self, check_model=True, initial_soc=None, inputs=None): self._parameter_values.process_geometry(self._geometry) # Only needs to set up mesh and discretisation once self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts) - self._disc = pybamm.Discretisation(self._mesh, self._spatial_methods) + self._disc = pybamm.Discretisation( + self._mesh, self._spatial_methods, **self._discretisation_kwargs + ) # Process all the different models self.steps_to_built_models = {} self.steps_to_built_solvers = {} @@ -333,7 +341,7 @@ def build_for_experiment(self, check_model=True, initial_soc=None, inputs=None): # It's ok to modify the model with set parameters in place as it's # not returned anywhere built_model = self._disc.process_model( - model_with_set_params, inplace=True, check_model=check_model + model_with_set_params, inplace=True ) solver = self._solver.copy() self.steps_to_built_solvers[step] = solver @@ -343,7 +351,6 @@ def solve( self, t_eval=None, solver=None, - check_model=True, save_at_cycles=None, calc_esoh=True, starting_solution=None, @@ -377,9 +384,6 @@ def solve( provided in the data. solver : :class:`pybamm.BaseSolver`, optional The solver to use to solve the model. If None, Simulation.solver is used - check_model : bool, optional - If True, model checks are performed after discretisation (see - :meth:`pybamm.Discretisation.process_model`). Default is True. save_at_cycles : int or list of ints, optional Which cycles to save the full sub-solutions for. If None, all cycles are saved. If int, every multiple of save_at_cycles is saved. If list, every @@ -416,7 +420,7 @@ def solve( inputs = inputs or {} if self.operating_mode in ["without experiment", "drive cycle"]: - self.build(check_model=check_model, initial_soc=initial_soc, inputs=inputs) + self.build(initial_soc=initial_soc, inputs=inputs) if save_at_cycles is not None: raise ValueError( "'save_at_cycles' option can only be used if simulating an " @@ -493,9 +497,7 @@ def solve( elif self.operating_mode == "with experiment": callbacks.on_experiment_start(logs) - self.build_for_experiment( - check_model=check_model, initial_soc=initial_soc, inputs=inputs - ) + self.build_for_experiment(initial_soc=initial_soc, inputs=inputs) if t_eval is not None: pybamm.logger.warning( "Ignoring t_eval as solution times are specified by the experiment" diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py index e9da61e001..61e3c2bc6e 100644 --- a/tests/unit/test_discretisations/test_discretisation.py +++ b/tests/unit/test_discretisations/test_discretisation.py @@ -1074,23 +1074,6 @@ def test_check_tab_bcs_error(self): with self.assertRaisesRegex(pybamm.ModelError, "Boundary conditions"): disc.check_tab_conditions(b, bcs) - def test_process_with_no_check(self): - # create model - whole_cell = ["negative electrode", "separator", "positive electrode"] - c = pybamm.Variable("c", domain=whole_cell) - N = pybamm.grad(c) - model = pybamm.BaseModel() - model.rhs = {c: pybamm.div(N)} - model.initial_conditions = {c: pybamm.Scalar(3)} - model.boundary_conditions = { - c: {"left": (0, "Neumann"), "right": (0, "Neumann")} - } - model.variables = {"c": c, "N": N} - - # create discretisation - disc = get_discretisation_for_testing() - disc.process_model(model, check_model=False) - def test_mass_matrix_inverse(self): # get mesh mesh = get_2p1d_mesh_for_testing(ypts=5, zpts=5) @@ -1229,17 +1212,34 @@ def test_length_scale_errors(self): def test_independent_rhs(self): a = pybamm.Variable("a") b = pybamm.Variable("b") - c = pybamm.Variable("c") + # Include a concatenation for the test + c_n = pybamm.Variable("c_n", ["negative electrode"]) + c_s = pybamm.Variable("c_s", ["separator"]) + c = pybamm.concatenation(c_n, c_s) + model = pybamm.BaseModel() - model.rhs = {a: b, b: c, c: -c} - model.initial_conditions = { - a: pybamm.Scalar(0), - b: pybamm.Scalar(1), - c: pybamm.Scalar(1), - } - disc = pybamm.Discretisation() + model.rhs = {a: b, b: 0, c: -c} + model.initial_conditions = {a: 0, b: 1, c: 1} + # test edge case where variable appears twice with different names + model.variables = {"a": a, "a again": a} + mesh = get_mesh_for_testing() + spatial_methods = {"macroscale": SpatialMethodForTesting()} + disc = pybamm.Discretisation( + mesh, spatial_methods, remove_independent_variables_from_rhs=True + ) disc.process_model(model) self.assertEqual(len(model.rhs), 2) + self.assertEqual(model.variables["a"], model.variables["a again"]) + + def test_independent_rhs_one_equation(self): + # Test that if there is only one equation, it is not removed + a = pybamm.Variable("a") + model = pybamm.BaseModel() + model.rhs = {a: 0} + model.initial_conditions = {a: 0} + disc = pybamm.Discretisation(remove_independent_variables_from_rhs=True) + disc.process_model(model) + self.assertEqual(len(model.rhs), 1) def test_independent_rhs_with_event(self): a = pybamm.Variable("a") @@ -1253,7 +1253,7 @@ def test_independent_rhs_with_event(self): c: pybamm.Scalar(1), } model.events = [pybamm.Event("a=1", a - 1)] - disc = pybamm.Discretisation() + disc = pybamm.Discretisation(remove_independent_variables_from_rhs=True) disc.process_model(model) self.assertEqual(len(model.rhs), 3) diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py index 6565a6d930..d617dc6619 100644 --- a/tests/unit/test_simulation.py +++ b/tests/unit/test_simulation.py @@ -7,6 +7,7 @@ import unittest import uuid from tempfile import TemporaryDirectory +from scipy.integrate import cumulative_trapezoid class TestSimulation(TestCase): @@ -69,8 +70,10 @@ def test_solve(self): self.assertTrue(val.has_symbol_of_classes(pybamm.Matrix)) # test solve without check - sim = pybamm.Simulation(pybamm.lithium_ion.SPM()) - sol = sim.solve(t_eval=[0, 600], check_model=False) + sim = pybamm.Simulation( + pybamm.lithium_ion.SPM(), discretisation_kwargs={"check_model": False} + ) + sol = sim.solve(t_eval=[0, 600]) for val in list(sim.built_model.rhs.values()): self.assertFalse(val.has_symbol_of_classes(pybamm.Parameter)) # skip test for scalar variables (e.g. discharge capacity) @@ -83,6 +86,19 @@ def test_solve(self): with self.assertRaisesRegex(ValueError, "starting_solution"): sim.solve(starting_solution=sol) + def test_solve_remove_independent_variables_from_rhs(self): + sim = pybamm.Simulation( + pybamm.lithium_ion.SPM(), + discretisation_kwargs={"remove_independent_variables_from_rhs": True}, + ) + sol = sim.solve([0, 600]) + t = sol["Time [s]"].data + I = sol["Current [A]"].data + q = sol["Discharge capacity [A.h]"].data + np.testing.assert_array_almost_equal( + q, cumulative_trapezoid(I, t, initial=0) / 3600 + ) + def test_solve_non_battery_model(self): model = pybamm.BaseModel() v = pybamm.Variable("v") diff --git a/tests/unit/test_solvers/test_idaklu_jax.py b/tests/unit/test_solvers/test_idaklu_jax.py index b1311001aa..6d891010d6 100644 --- a/tests/unit/test_solvers/test_idaklu_jax.py +++ b/tests/unit/test_solvers/test_idaklu_jax.py @@ -30,7 +30,7 @@ model.initial_conditions = {u1: 0, u2: 0, v: 1} model.variables = {"v": v, "u1": u1, "u2": u2} disc = pybamm.Discretisation() - disc.process_model(model, remove_independent_variables_from_rhs=False) + disc.process_model(model) t_eval = np.linspace(0, 1, 100) idaklu_solver = pybamm.IDAKLUSolver(rtol=1e-6, atol=1e-6) diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py index 8423c16b68..1697623486 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -160,7 +160,7 @@ def test_input_params(self): model.initial_conditions = {u1: 0, u2: 0, u3: 0, v: 1} disc = pybamm.Discretisation() - disc.process_model(model, remove_independent_variables_from_rhs=False) + disc.process_model(model) solver = pybamm.IDAKLUSolver(root_method=root_method) @@ -563,7 +563,11 @@ def construct_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 = pybamm.Discretisation( + mesh, + model.default_spatial_methods, + remove_independent_variables_from_rhs=True, + ) disc.process_model(model) return model