Merge pull request #4020 from pybamm-team/change-remove-independent-v…

…ars-default

change remove independent variables setting to False by default

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))

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,

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]

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"

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)
 

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")