Issue 2193 starting solution

CHANGELOG.md

@@ -21,6 +21,7 @@
 
 ## Bug fixes
 
+-   Fixed a bug in `Simulation` that caused initial conditions to change when solving an experiment multiple times ([#2204](https://github.com/pybamm-team/PyBaMM/pull/2204))
 -   Fixed labels and ylims in `plot_voltage_components`([#2183](https://github.com/pybamm-team/PyBaMM/pull/2183))
 -   Fixed 2D interpolant ([#2180](https://github.com/pybamm-team/PyBaMM/pull/2180))
 -   Fixes a bug where the SPMe always builds even when `build=False` ([#2169](https://github.com/pybamm-team/PyBaMM/pull/2169))

pybamm/models/base_model.py

@@ -579,7 +579,7 @@ def _build_model(self):
 
         self.build_model_equations()
 
-    def set_initial_conditions_from(self, solution, inplace=True):
+    def set_initial_conditions_from(self, solution, inplace=True, return_type="model"):
         """
         Update initial conditions with the final states from a Solution object or from
         a dictionary.
@@ -590,14 +590,12 @@ def set_initial_conditions_from(self, solution, inplace=True):
         ----------
         solution : :class:`pybamm.Solution`, or dict
             The solution to use to initialize the model
-        inplace : bool
-            Whether to modify the model inplace or create a new model
+        inplace : bool, optional
+            Whether to modify the model inplace or create a new model (default True)
+        return_type : str, optional
+            Whether to return the model (default) or initial conditions ("ics")
         """
-        if inplace is True:
-            model = self
-        else:
-            model = self.new_copy()
-
+        initial_conditions = {}
         if isinstance(solution, pybamm.Solution):
             solution = solution.last_state
         for var in self.initial_conditions:
@@ -623,7 +621,7 @@ def set_initial_conditions_from(self, solution, inplace=True):
                     final_state_eval = final_state[:, :, -1].flatten(order="F")
                 else:
                     raise NotImplementedError("Variable must be 0D, 1D, or 2D")
-                model.initial_conditions[var] = pybamm.Vector(final_state_eval)
+                initial_conditions[var] = pybamm.Vector(final_state_eval)
             elif isinstance(var, pybamm.Concatenation):
                 children = []
                 for child in var.orphans:
@@ -645,7 +643,7 @@ def set_initial_conditions_from(self, solution, inplace=True):
                             "Variable in concatenation must be 1D"
                         )
                     children.append(final_state_eval)
-                model.initial_conditions[var] = pybamm.Vector(np.concatenate(children))
+                initial_conditions[var] = pybamm.Vector(np.concatenate(children))
 
             else:
                 raise NotImplementedError(
@@ -670,14 +668,26 @@ def set_initial_conditions_from(self, solution, inplace=True):
                     )
                 else:
                     slices.append(y_slices[symbol][0])
-            equations = list(model.initial_conditions.values())
+            equations = list(initial_conditions.values())
             # sort equations according to slices
             sorted_equations = [eq for _, eq in sorted(zip(slices, equations))]
-            model.concatenated_initial_conditions = pybamm.NumpyConcatenation(
+            concatenated_initial_conditions = pybamm.NumpyConcatenation(
                 *sorted_equations
             )
+        else:
+            concatenated_initial_conditions = None
+
+        if return_type == "model":
+            if inplace is True:
+                model = self
+            else:
+                model = self.new_copy()
 
-        return model
+            model.initial_conditions = initial_conditions
+            model.concatenated_initial_conditions = concatenated_initial_conditions
+            return model
+        elif return_type == "ics":
+            return initial_conditions, concatenated_initial_conditions
 
     def check_and_combine_dict(self, dict1, dict2):
         # check that the key ids are distinct

pybamm/solvers/base_solver.py

@@ -1296,9 +1296,12 @@ def step(
                 # initialize with old solution
                 model.y0 = old_solution.all_ys[-1][:, -1]
             else:
-                model.y0 = model.set_initial_conditions_from(
-                    old_solution
-                ).concatenated_initial_conditions.evaluate(0, inputs=ext_and_inputs)
+                _, concatenated_initial_conditions = model.set_initial_conditions_from(
+                    old_solution, return_type="ics"
+                )
+                model.y0 = concatenated_initial_conditions.evaluate(
+                    0, inputs=ext_and_inputs
+                )
         set_up_time = timer.time()
 
         # (Re-)calculate consistent initial conditions

tests/unit/test_experiments/test_simulation_with_experiment.py

@@ -124,6 +124,26 @@ def test_run_experiment(self):
         self.assertEqual(len(sol3.cycles), 2)
         os.remove("test_experiment.sav")
 
+    def test_run_experiment_multiple_times(self):
+        experiment = pybamm.Experiment(
+            [
+                (
+                    "Discharge at C/20 for 1 hour",
+                    "Charge at C/20 until 4.1 V",
+                )
+            ]
+            * 3
+        )
+        model = pybamm.lithium_ion.DFN()
+        sim = pybamm.Simulation(model, experiment=experiment)
+
+        # Test that solving twice gives the same solution (see #2193)
+        sol1 = sim.solve()
+        sol2 = sim.solve()
+        np.testing.assert_array_equal(
+            sol1["Terminal voltage [V]"].data, sol2["Terminal voltage [V]"].data
+        )
+
     def test_run_experiment_cccv_ode(self):
         experiment_2step = pybamm.Experiment(
             [