Repository: pybamm-team/PyBaMM

Issue:
"Throughput capacity [A.h]" is different for sol.cycles[0] and sol.cycles[0].steps[0]
### Discussed in https://github.com/pybamm-team/PyBaMM/discussions/2523

<div type='discussions-op-text'>

<sup>Originally posted by **RuiheLi** December  5, 2022</sup>
Dear PyBaMM,

I am using PyBaMM 22.9 and Python 3.9.6.

I find when using that set_initial_conditions_from, "Throughput capacity [A.h]" is different for sol.cycles[0] and sol.cycles[0].steps[0] for a later solution. Figure below:
![image](https://user-images.githubusercontent.com/84007676/205644495-558eabbe-1f36-4dd5-9935-44f131704896.png)


Script below:

```
import pybamm 
import matplotlib as mpl
fs=17;
font = {'family' : 'DejaVu Sans','size'   : fs}
mpl.rc('font', **font)
V_max = 4.0;  V_min = 2.7; 
model_0 = pybamm.lithium_ion.DFN(options={
    "calculate discharge energy":"true", },)
exp_AGE_text = [(f"Discharge at 1 C until {V_min} V", 
    f"Charge at 1 C until {V_max} V", ),  ]
Experiment    = pybamm.Experiment( exp_AGE_text * 1)
Para_0 = pybamm.ParameterValues("OKane2022" )

Sim_0= pybamm.Simulation(
    model_0,
    experiment = Experiment, 
    parameter_values=Para_0, 
    solver = pybamm.CasadiSolver(),)
Sol_0 = Sim_0.solve(calc_esoh=False,) 

model_1 = model_0.set_initial_conditions_from(Sol_0, inplace=False)
Sim_1= pybamm.Simulation(
    model_1,
    experiment = Experiment, 
    parameter_values=Para_0, 
    solver = pybamm.CasadiSolver(),)
Sol_1 = Sim_1.solve(calc_esoh=False,) 

model_2 = model_0.set_initial_conditions_from(Sol_1, inplace=False)
Sim_2= pybamm.Simulation(
    model_2,
    experiment = Experiment, 
    parameter_values=Para_0, 
    solver = pybamm.CasadiSolver(),)
Sol_2 = Sim_2.solve(calc_esoh=False,) 

label = ["Sol_0","Sol_1","Sol_2",] 
output_variables3 = [
    "Throughput energy [W.h]",
    "Throughput capacity [A.h]", 
    "Terminal voltage [V]", 
]
quick_plot = pybamm.QuickPlot(
    [Sol_0,Sol_1,Sol_2,], output_variables3,label,variable_limits='fixed',
    time_unit='hours',n_rows=1) #figsize = (18,12),
quick_plot.dynamic_plot()
label = ["Sol_2 whole","Sol_2 cycles","Sol_2 steps",] 
quick_plot = pybamm.QuickPlot(
    [Sol_2,Sol_2.cycles[0],Sol_2.cycles[0].steps[0]], output_variables3,label,variable_limits='fixed',
    time_unit='hours',n_rows=1) #figsize = (18,12),
quick_plot.dynamic_plot()
```

Do you have any ideas on this? Big thanks!

</div>

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ea458b4b6..a8d7d155d 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 1503c3048..e854c94e0 100644
--- a/pybamm/batch_study.py
+++ b/pybamm/batch_study.py
@@ -102,7 +102,6 @@ class BatchStudy:
         self,
         t_eval=None,
         solver=None,
-        check_model=True,
         save_at_cycles=None,
         calc_esoh=True,
         starting_solution=None,
@@ -158,7 +157,6 @@ class BatchStudy:
                 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 92e9e4e24..d61d8a4d6 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 @@ class Discretisation:
         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 @@ class Discretisation:
         # 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 @@ class Discretisation:
         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 @@ class Discretisation:
         # 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 @@ class Discretisation:
         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 @@ class Discretisation:
                     # 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 060d1cdd8..a81f7b62b 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 @@ class Simulation:
         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 @@ class Simulation:
         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 @@ class Simulation:
         # 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 @@ class Simulation:
 
         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 @@ class Simulation:
         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 @@ class Simulation:
             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 @@ class Simulation:
                 # 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 @@ class Simulation:
         self,
         t_eval=None,
         solver=None,
-        check_model=True,
         save_at_cycles=None,
         calc_esoh=True,
         starting_solution=None,
@@ -377,9 +384,6 @@ class Simulation:
             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 @@ class Simulation:
         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 @@ class Simulation:
 
         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"


Test patch:
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index e9da61e00..61e3c2bc6 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -1074,23 +1074,6 @@ class TestDiscretise(TestCase):
         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 @@ class TestDiscretise(TestCase):
     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 @@ class TestDiscretise(TestCase):
             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 6565a6d93..d617dc661 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -7,6 +7,7 @@ import sys
 import unittest
 import uuid
 from tempfile import TemporaryDirectory
+from scipy.integrate import cumulative_trapezoid
 
 
 class TestSimulation(TestCase):
@@ -69,8 +70,10 @@ class TestSimulation(TestCase):
                 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 @@ class TestSimulation(TestCase):
         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 b1311001a..6d891010d 100644
--- a/tests/unit/test_solvers/test_idaklu_jax.py
+++ b/tests/unit/test_solvers/test_idaklu_jax.py
@@ -30,7 +30,7 @@ if pybamm.have_idaklu() and pybamm.have_jax():
     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 8423c16b6..169762348 100644
--- a/tests/unit/test_solvers/test_idaklu_solver.py
+++ b/tests/unit/test_solvers/test_idaklu_solver.py
@@ -160,7 +160,7 @@ class TestIDAKLUSolver(TestCase):
             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 @@ class TestIDAKLUSolver(TestCase):
             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