Repository: pybamm-team/PyBaMM

Issue:
Error while simulating cycling
### PyBaMM Version

22.5

### Python Version

3.8.10

### Describe the bug

I'm struggling with fast-charge (2C) aging (SEI+irreversible plating) simulation for long cycling like 50 cycles or above. I tried different things like:
- Different models: DFN, SPMe
- Different ambient temperatures: 0-25 C
- Different i0s for SEI and plating
- Different solver modes and dt_max: Only Casadi "safe" works

The problem has two parts: first, I could run it in COMSOL with more aggressive parameters e.g. higher anode tortuosity, lower temperatures, higher charge C-rates, etc.; however, in Pybamm, I am getting the error concentration reached zero with much smaller tortuosity. 
Second, when I use smaller anode tortuosity it works without the irreversible Li plating submodel only if I include CV discharge at 2.5V at the end, otherwise, it raises the error "Experiment is infeasible" in the middle of cycling and stops at cycle 26 (reversible plating also works fine in this case). So, the problem comes from the irreversible Li plating model. I use the O'Kane model and parameters but it seems that Pybamm has some issues with it although LLI is not big at the end. The script and error for 25C are below:

### Steps to Reproduce


```
# Solving for 2C fast charge

parameter_values["Initial concentration in negative electrode [mol.m-3]"] = 300
parameter_values["Initial concentration in positive electrode [mol.m-3]"] = 49000

model = pybamm.lithium_ion.DFN({"thermal": "x-full", "SEI": "solvent-diffusion limited","lithium plating":"irreversible"})

solver= pybamm.CasadiSolver("safe")

experiment = pybamm.Experiment(
    [
        
            "Charge at 2 C until 4.2 V",
            "Hold at 4.2V until C/20",
            "Rest for 15 minutes",
            "Discharge at 0.5 C until 2.5 V",
            "Hold at 2.5 V until C/20",
            "Rest for 15 minutes"
        
    ]
    * 50)

var_pts = {
    "x_n": 40,  # negative electrode
    "x_s": 5,  # separator 
    "x_p": 10,  # positive electrode
    "r_n": 20,  # negative particle
    "r_p": 10,  # positive particle
}

sim = pybamm.Simulation(model, var_pts=var_pts, experiment=experiment, parameter_values=parameter_values, solver=solver)

sim.solve(calc_esoh=False)
sim.plot();

##################
At t = 0.00934147 and h = 1.58982e-018, the corrector convergence failed repeatedly or with |h| = hmin.
At t = 0.00949702, , mxstep steps taken before reaching tout.
At t = 0.00854179 and h = 6.02689e-026, the corrector convergence failed repeatedly or with |h| = hmin.
The linesearch algorithm failed with too small a step.
2022-06-26 22:07:10.050 - [WARNING] callbacks.on_experiment_infeasible(237): 

	Experiment is infeasible: 'final time' was triggered during 'Discharge at 0.5 C until 2.5 V'. The returned solution only contains the first 153 cycles. Try reducing the current, shortening the time interval, or reducing the period.
``` 

### Relevant log output

_No response_

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ad51621a4..96edd382c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,14 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Optimizations
+
+-   Added error for when solution vector gets too large, to help debug solver errors ([#2138](https://github.com/pybamm-team/PyBaMM/pull/2138))
+
+## Bug fixes
+
+-   Fixed bug causing cut-off voltage to change after setting up a simulation with a model ([#2138](https://github.com/pybamm-team/PyBaMM/pull/2138))
+-   A single solution cycle can now be used as a starting solution for a simulation ([#2138](https://github.com/pybamm-team/PyBaMM/pull/2138))
+
 ## Breaking changes
 
 -   Exchange-current density functions (and some other functions) now take an additional argument, the maximum particle concentration for that phase ([#2134](https://github.com/pybamm-team/PyBaMM/pull/2134))
diff --git a/pybamm/models/submodels/particle_mechanics/crack_propagation.py b/pybamm/models/submodels/particle_mechanics/crack_propagation.py
index 63a5dd819..c723a2e32 100644
--- a/pybamm/models/submodels/particle_mechanics/crack_propagation.py
+++ b/pybamm/models/submodels/particle_mechanics/crack_propagation.py
@@ -92,3 +92,17 @@ class CrackPropagation(BaseMechanics):
             l_cr = variables[self.domain + " particle crack length"]
             l0 = pybamm.PrimaryBroadcast(1, self.domain.lower() + " electrode")
         self.initial_conditions = {l_cr: l0}
+
+    def set_events(self, variables):
+        domain = self.domain.lower()
+        if self.x_average is True:
+            l_cr = variables[f"X-averaged {domain} particle crack length"]
+        else:
+            l_cr = variables[self.domain + " particle crack length"]
+        self.events.append(
+            pybamm.Event(
+                f"{domain} particle crack length larger than particle radius",
+                pybamm.max(l_cr) - self.domain_param.R_typ / self.domain_param.l_cr_0,
+                pybamm.EventType.TERMINATION,
+            )
+        )
diff --git a/pybamm/settings.py b/pybamm/settings.py
index da2f80e41..c3fd20bf7 100644
--- a/pybamm/settings.py
+++ b/pybamm/settings.py
@@ -11,6 +11,7 @@ class Settings(object):
     _heaviside_smoothing = "exact"
     _abs_smoothing = "exact"
     max_words_in_line = 4
+    max_y_value = 1e5
 
     @property
     def debug_mode(self):
@@ -37,16 +38,19 @@ class Settings(object):
         self.heaviside_smoothing = k
         self.abs_smoothing = k
 
+    def check_k(self, k):
+        if k != "exact" and k <= 0:
+            raise ValueError(
+                "smoothing parameter must be 'exact' or a strictly positive number"
+            )
+
     @property
     def min_smoothing(self):
         return self._min_smoothing
 
     @min_smoothing.setter
     def min_smoothing(self, k):
-        if k != "exact" and k <= 0:
-            raise ValueError(
-                "smoothing parameter must be 'exact' or a strictly positive number"
-            )
+        self.check_k(k)
         self._min_smoothing = k
 
     @property
@@ -55,10 +59,7 @@ class Settings(object):
 
     @max_smoothing.setter
     def max_smoothing(self, k):
-        if k != "exact" and k <= 0:
-            raise ValueError(
-                "smoothing parameter must be 'exact' or a strictly positive number"
-            )
+        self.check_k(k)
         self._max_smoothing = k
 
     @property
@@ -67,10 +68,7 @@ class Settings(object):
 
     @heaviside_smoothing.setter
     def heaviside_smoothing(self, k):
-        if k != "exact" and k <= 0:
-            raise ValueError(
-                "smoothing parameter must be 'exact' or a strictly positive number"
-            )
+        self.check_k(k)
         self._heaviside_smoothing = k
 
     @property
@@ -79,10 +77,7 @@ class Settings(object):
 
     @abs_smoothing.setter
     def abs_smoothing(self, k):
-        if k != "exact" and k <= 0:
-            raise ValueError(
-                "smoothing parameter must be 'exact' or a strictly positive number"
-            )
+        self.check_k(k)
         self._abs_smoothing = k
 
 
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index fc9c7f0f6..9b80574cf 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -407,11 +407,15 @@ class Simulation:
                 # Keep the min and max voltages as safeguards but add some tolerances
                 # so that they are not triggered before the voltage limits in the
                 # experiment
-                for event in new_model.events:
+                for i, event in enumerate(new_model.events):
                     if event.name == "Minimum voltage":
-                        event._expression += 1
+                        new_model.events[i] = pybamm.Event(
+                            event.name, event.expression + 1, event.event_type
+                        )
                     elif event.name == "Maximum voltage":
-                        event._expression -= 1
+                        new_model.events[i] = pybamm.Event(
+                            event.name, event.expression - 1, event.event_type
+                        )
 
                 # Update parameter values
                 new_parameter_values = self.parameter_values.copy()
@@ -723,10 +727,28 @@ class Simulation:
             inputs = kwargs.get("inputs", {})
             timer = pybamm.Timer()
 
+            # Set up eSOH model (for summary variables)
+            if calc_esoh is True:
+                esoh_model = pybamm.lithium_ion.ElectrodeSOH()
+                esoh_sim = pybamm.Simulation(
+                    esoh_model, parameter_values=self.parameter_values
+                )
+            else:
+                esoh_sim = None
+
             if starting_solution is None:
                 starting_solution_cycles = []
                 starting_solution_summary_variables = []
                 starting_solution_first_states = []
+            elif not hasattr(starting_solution, "all_summary_variables"):
+                (
+                    cycle_solution,
+                    cycle_sum_vars,
+                    cycle_first_state,
+                ) = pybamm.make_cycle_solution(starting_solution.steps, esoh_sim, True)
+                starting_solution_cycles = [cycle_solution]
+                starting_solution_summary_variables = [cycle_sum_vars]
+                starting_solution_first_states = [cycle_first_state]
             else:
                 starting_solution_cycles = starting_solution.cycles.copy()
                 starting_solution_summary_variables = (
@@ -742,15 +764,6 @@ class Simulation:
             all_first_states = starting_solution_first_states
             current_solution = starting_solution
 
-            # Set up eSOH model (for summary variables)
-            if calc_esoh is True:
-                esoh_model = pybamm.lithium_ion.ElectrodeSOH()
-                esoh_sim = pybamm.Simulation(
-                    esoh_model, parameter_values=self.parameter_values
-                )
-            else:
-                esoh_sim = None
-
             voltage_stop = self.experiment.termination.get("voltage")
             logs["stopping conditions"] = {"voltage": voltage_stop}
 
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index e10508848..a78d89070 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -180,6 +180,7 @@ class CasadiSolver(pybamm.BaseSolver):
             # Check if the sign of an event changes, if so find an accurate
             # termination point and exit
             solution = self._solve_for_event(solution, init_event_signs)
+            solution.check_ys_are_not_too_large()
             return solution
         elif self.mode in ["safe", "safe without grid"]:
             y0 = model.y0
@@ -289,6 +290,7 @@ class CasadiSolver(pybamm.BaseSolver):
             if bool(model.calculate_sensitivities):
                 solution.sensitivities = True
 
+            solution.check_ys_are_not_too_large()
             return solution
 
     def _solve_for_event(self, coarse_solution, init_event_signs):
@@ -704,6 +706,7 @@ class CasadiSolver(pybamm.BaseSolver):
                     model,
                     inputs_dict,
                     sensitivities=extract_sensitivities_in_solution,
+                    check_solution=False,
                 )
                 sol.integration_time = integration_time
                 return sol
@@ -738,6 +741,7 @@ class CasadiSolver(pybamm.BaseSolver):
                 model,
                 inputs_dict,
                 sensitivities=extract_sensitivities_in_solution,
+                check_solution=False,
             )
             sol.integration_time = integration_time
             return sol
diff --git a/pybamm/solvers/solution.py b/pybamm/solvers/solution.py
index 01d5acb3c..4d816a53f 100644
--- a/pybamm/solvers/solution.py
+++ b/pybamm/solvers/solution.py
@@ -73,6 +73,7 @@ class Solution(object):
         y_event=None,
         termination="final time",
         sensitivities=False,
+        check_solution=True,
     ):
         if not isinstance(all_ts, list):
             all_ts = [all_ts]
@@ -105,6 +106,10 @@ class Solution(object):
             isinstance(v, casadi.MX) for v in self.all_inputs[0].values()
         )
 
+        # Check no ys are too large
+        if check_solution and not self.has_symbolic_inputs:
+            self.check_ys_are_not_too_large()
+
         # Copy the timescale_eval and lengthscale_evals if they exist
         if hasattr(all_models[0], "timescale_eval"):
             self.timescale_eval = all_models[0].timescale_eval
@@ -312,6 +317,24 @@ class Solution(object):
                 "computed explicitly."
             )
 
+    def check_ys_are_not_too_large(self):
+        # Only check last one so that it doesn't take too long
+        # We only care about the cases where y is growing too large without any
+        # restraint, so if y gets large in the middle then comes back down that is ok
+        y, model = self.all_ys[-1], self.all_models[-1]
+        y = y[:, -1]
+        if np.any(y > pybamm.settings.max_y_value):
+            for var in [*model.rhs.keys(), *model.algebraic.keys()]:
+                y_var = y[model.variables[var.name].y_slices[0]]
+                if np.any(y_var > pybamm.settings.max_y_value):
+                    pybamm.logger.error(
+                        f"Solution for '{var}' exceeds the maximum allowed value "
+                        f"of `{pybamm.settings.max_y_value}. This could be due to "
+                        "incorrect nondimensionalisation, model formulation, or "
+                        "parameter values. The maximum allowed value is set by "
+                        "'pybammm.settings.max_y_value'."
+                    )
+
     @property
     def all_ts(self):
         return self._all_ts


Test patch:
diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py
index dd6a01cfe..f0c7bdf2d 100644
--- a/tests/unit/test_experiments/test_simulation_with_experiment.py
+++ b/tests/unit/test_experiments/test_simulation_with_experiment.py
@@ -109,6 +109,12 @@ class TestSimulationExperiment(unittest.TestCase):
         self.assertGreater(sol2.t[-1], sol.t[-1])
         self.assertEqual(sol2.cycles[0], sol.cycles[0])
         self.assertEqual(len(sol2.cycles), 2)
+        # Solve again starting from solution but only inputting the cycle
+        sol2 = sim.solve(starting_solution=sol.cycles[-1])
+        self.assertEqual(sol2.termination, "final time")
+        self.assertGreater(sol2.t[-1], sol.t[-1])
+        self.assertEqual(len(sol2.cycles), 2)
+
         # Check starting solution is unchanged
         self.assertEqual(len(sol.cycles), 1)
 
diff --git a/tests/unit/test_plotting/test_quick_plot.py b/tests/unit/test_plotting/test_quick_plot.py
index d876d4035..5f3217dc0 100644
--- a/tests/unit/test_plotting/test_quick_plot.py
+++ b/tests/unit/test_plotting/test_quick_plot.py
@@ -16,7 +16,7 @@ class TestQuickPlot(unittest.TestCase):
         c = pybamm.Variable("c", domain=[])
 
         # Simple ODEs
-        model.rhs = {a: pybamm.Scalar(2), b: pybamm.Scalar(0), c: -c}
+        model.rhs = {a: pybamm.Scalar(0.2), b: pybamm.Scalar(0), c: -c}
 
         # Simple initial conditions
         model.initial_conditions = {
@@ -28,6 +28,7 @@ class TestQuickPlot(unittest.TestCase):
         # Broadcast some of the variables
         model.variables = {
             "a": a,
+            "b": b,
             "b broadcasted": pybamm.FullBroadcast(b, whole_cell, "current collector"),
             "c broadcasted": pybamm.FullBroadcast(
                 c, ["negative electrode", "separator"], "current collector"
@@ -146,7 +147,7 @@ class TestQuickPlot(unittest.TestCase):
             quick_plot.plots[("a",)][0][0].get_xdata(), t_eval
         )
         np.testing.assert_array_almost_equal(
-            quick_plot.plots[("a",)][0][0].get_ydata(), 2 * t_eval
+            quick_plot.plots[("a",)][0][0].get_ydata(), 0.2 * t_eval
         )
         quick_plot = pybamm.QuickPlot(solution, ["a"], time_unit="minutes")
         quick_plot.plot(0)
@@ -155,7 +156,7 @@ class TestQuickPlot(unittest.TestCase):
             quick_plot.plots[("a",)][0][0].get_xdata(), t_eval / 60
         )
         np.testing.assert_array_almost_equal(
-            quick_plot.plots[("a",)][0][0].get_ydata(), 2 * t_eval
+            quick_plot.plots[("a",)][0][0].get_ydata(), 0.2 * t_eval
         )
         quick_plot = pybamm.QuickPlot(solution, ["a"], time_unit="hours")
         quick_plot.plot(0)
@@ -164,7 +165,7 @@ class TestQuickPlot(unittest.TestCase):
             quick_plot.plots[("a",)][0][0].get_xdata(), t_eval / 3600
         )
         np.testing.assert_array_almost_equal(
-            quick_plot.plots[("a",)][0][0].get_ydata(), 2 * t_eval
+            quick_plot.plots[("a",)][0][0].get_ydata(), 0.2 * t_eval
         )
         with self.assertRaisesRegex(ValueError, "time unit"):
             pybamm.QuickPlot(solution, ["a"], time_unit="bad unit")
diff --git a/tests/unit/test_solvers/test_solution.py b/tests/unit/test_solvers/test_solution.py
index 244c70bd3..9fea56d51 100644
--- a/tests/unit/test_solvers/test_solution.py
+++ b/tests/unit/test_solvers/test_solution.py
@@ -39,6 +39,16 @@ class TestSolution(unittest.TestCase):
         ):
             sol.set_t()
 
+        ts = [np.array([1, 2, 3])]
+        bad_ys = [(pybamm.settings.max_y_value + 1) * np.ones((1, 3))]
+        model = pybamm.BaseModel()
+        var = pybamm.StateVector(slice(0, 1))
+        model.rhs = {var: 0}
+        model.variables = {var.name: var}
+        with self.assertLogs() as captured:
+            pybamm.Solution(ts, bad_ys, model, {})
+        self.assertIn("exceeds the maximum", captured.records[0].getMessage())
+
     def test_add_solutions(self):
         # Set up first solution
         t1 = np.linspace(0, 1)
@@ -279,11 +289,9 @@ class TestSolution(unittest.TestCase):
         csv_str = solution.save_data(variables=["c", "2c"], to_format="csv")
 
         # check string is the same as the file
-        with open('test.csv') as f:
+        with open("test.csv") as f:
             # need to strip \r chars for windows
-            self.assertEqual(
-                csv_str.replace('\r', ''), f.read()
-            )
+            self.assertEqual(csv_str.replace("\r", ""), f.read())
 
         # read csv
         df = pd.read_csv("test.csv")
@@ -295,11 +303,9 @@ class TestSolution(unittest.TestCase):
         json_str = solution.save_data(to_format="json")
 
         # check string is the same as the file
-        with open('test.json') as f:
+        with open("test.json") as f:
             # need to strip \r chars for windows
-            self.assertEqual(
-                json_str.replace('\r', ''), f.read()
-            )
+            self.assertEqual(json_str.replace("\r", ""), f.read())
 
         # check if string has the right values
         json_data = json.loads(json_str)