Disallow modifications to pybamm.Simulation object attributes after initialisation

CHANGELOG.md

@@ -10,14 +10,15 @@
 ## Bug fixes
 
 - Fixed a bug with `_Heaviside._evaluate_for_shape` which meant some expressions involving heaviside function and subtractions did not work ([#3306](https://github.com/pybamm-team/PyBaMM/pull/3306))
-- The `OneDimensionalX` thermal model has been updated to account for edge/tab cooling and account for the current collector volumetric heat capacity. It now gives the correct behaviour compared with a lumped model with the correct total heat transfer coefficient and surface area for cooling. ([#3042](https://github.com/pybamm-team/PyBaMM/pull/3042))
+- Attributes of `pybamm.Simulation` objects (models, parameter values, geometries, choice of solver, and output variables) are now private and as such cannot be edited in-place after the simulation has been created ([#3267](https://github.com/pybamm-team/PyBaMM/pull/3267)
 - Fixed a bug where the "basic" lithium-ion models gave incorrect results when using nonlinear particle diffusivity ([#3207](https://github.com/pybamm-team/PyBaMM/pull/3207))
 - Particle size distributions now work with SPMe and NewmanTobias models ([#3207](https://github.com/pybamm-team/PyBaMM/pull/3207))
 - Fix to simulate c_rate steps with drive cycles ([#3186](https://github.com/pybamm-team/PyBaMM/pull/3186))
-- Parameters in `Prada2013` have been updated to better match those given in the paper, which is a 2.3 Ah cell, instead of the mix-and-match with the 1.1 Ah cell from Lain2019.
-- Error generated when invalid parameter values are passed. ([#3132](https://github.com/pybamm-team/PyBaMM/pull/3132))
-- Thevenin() model is now constructed with standard variables: `Time [s], Time [min], Time [h]` ([#3143](https://github.com/pybamm-team/PyBaMM/pull/3143))
 - Fix SEI Example Notebook ([#3166](https://github.com/pybamm-team/PyBaMM/pull/3166))
+- Thevenin() model is now constructed with standard variables: `Time [s], Time [min], Time [h]` ([#3143](https://github.com/pybamm-team/PyBaMM/pull/3143))
+- Error generated when invalid parameter values are passed. ([#3132](https://github.com/pybamm-team/PyBaMM/pull/3132))
+- Parameters in `Prada2013` have been updated to better match those given in the paper, which is a 2.3 Ah cell, instead of the mix-and-match with the 1.1 Ah cell from Lain2019. ([#3096](https://https://github.com/pybamm-team/PyBaMM/pull/3096))
+- The `OneDimensionalX` thermal model has been updated to account for edge/tab cooling and account for the current collector volumetric heat capacity. It now gives the correct behaviour compared with a lumped model with the correct total heat transfer coefficient and surface area for cooling. ([#3042](https://github.com/pybamm-team/PyBaMM/pull/3042))
 
 ## Breaking changes
 

pybamm/simulation.py

@@ -4,7 +4,6 @@
 import pickle
 import pybamm
 import numpy as np
-import copy
 import warnings
 import sys
 from functools import lru_cache
@@ -74,8 +73,8 @@ def __init__(
         output_variables=None,
         C_rate=None,
     ):
-        self.parameter_values = parameter_values or model.default_parameter_values
-        self._unprocessed_parameter_values = self.parameter_values
+        self._parameter_values = parameter_values or model.default_parameter_values
+        self._unprocessed_parameter_values = self._parameter_values
 
         if isinstance(model, pybamm.lithium_ion.BasicDFNHalfCell):
             if experiment is not None:
@@ -114,14 +113,14 @@ def __init__(
             self.experiment = experiment.copy()
 
         self._unprocessed_model = model
-        self.model = model
+        self._model = model
 
-        self.geometry = geometry or self.model.default_geometry
-        self.submesh_types = submesh_types or self.model.default_submesh_types
-        self.var_pts = var_pts or self.model.default_var_pts
-        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._geometry = geometry or self._model.default_geometry
+        self._submesh_types = submesh_types or self._model.default_submesh_types
+        self._var_pts = var_pts or self._model.default_var_pts
+        self._spatial_methods = spatial_methods or self._model.default_spatial_methods
+        self._solver = solver or self._model.default_solver
+        self._output_variables = output_variables
 
         # Initialize empty built states
         self._model_with_set_params = None
@@ -201,16 +200,16 @@ def set_up_and_parameterise_experiment(self):
 
     def set_up_and_parameterise_model_for_experiment(self):
         """
-        Set up self.model to be able to run the experiment (new version).
+        Set up self._model to be able to run the experiment (new version).
         In this version, a new model is created for each step.
 
         This increases set-up time since several models to be processed, but
         reduces simulation time since the model formulation is efficient.
         """
         self.experiment_unique_steps_to_model = {}
         for op_number, op in enumerate(self.experiment.unique_steps):
-            new_model = self.model.new_copy()
-            new_parameter_values = self.parameter_values.copy()
+            new_model = self._model.new_copy()
+            new_parameter_values = self._parameter_values.copy()
 
             if op.type != "current":
                 # Voltage or power control
@@ -259,9 +258,9 @@ def set_up_and_parameterise_model_for_experiment(self):
 
         # Set up rest model if experiment has start times
         if self.experiment.initial_start_time:
-            new_model = self.model.new_copy()
+            new_model = self._model.new_copy()
             # Update parameter values
-            new_parameter_values = self.parameter_values.copy()
+            new_parameter_values = self._parameter_values.copy()
             self._original_temperature = new_parameter_values["Ambient temperature [K]"]
             new_parameter_values.update(
                 {"Current function [A]": 0, "Ambient temperature [K]": "[input]"},
@@ -360,8 +359,8 @@ def set_parameters(self):
         self._model_with_set_params = self._parameter_values.process_model(
             self._unprocessed_model, inplace=False
         )
-        self._parameter_values.process_geometry(self.geometry)
-        self.model = self._model_with_set_params
+        self._parameter_values.process_geometry(self._geometry)
+        self._model = self._model_with_set_params
 
     def set_initial_soc(self, initial_soc):
         if self._built_initial_soc != initial_soc:
@@ -371,8 +370,8 @@ def set_initial_soc(self, initial_soc):
             self.op_conds_to_built_models = None
             self.op_conds_to_built_solvers = None
 
-        param = self.model.param
-        self.parameter_values = (
+        param = self._model.param
+        self._parameter_values = (
             self._unprocessed_parameter_values.set_initial_stoichiometries(
                 initial_soc, param=param, inplace=False
             )
@@ -403,9 +402,9 @@ def build(self, check_model=True, initial_soc=None):
 
         if self.built_model:
             return
-        elif self.model.is_discretised:
-            self._model_with_set_params = self.model
-            self._built_model = self.model
+        elif self._model.is_discretised:
+            self._model_with_set_params = self._model
+            self._built_model = self._model
         else:
             self.set_parameters()
             self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts)
@@ -447,7 +446,7 @@ def build_for_experiment(self, check_model=True, initial_soc=None):
                 built_model = self._disc.process_model(
                     model_with_set_params, inplace=True, check_model=check_model
                 )
-                solver = self.solver.copy()
+                solver = self._solver.copy()
                 self.op_conds_to_built_solvers[op_cond] = solver
                 self.op_conds_to_built_models[op_cond] = built_model
 
@@ -519,7 +518,7 @@ def solve(
         """
         # Setup
         if solver is None:
-            solver = self.solver
+            solver = self._solver
 
         callbacks = pybamm.callbacks.setup_callbacks(callbacks)
         logs = {}
@@ -537,7 +536,7 @@ def solve(
                 )
             if (
                 self.operating_mode == "without experiment"
-                or self.model.name == "ElectrodeSOH model"
+                or self._model.name == "ElectrodeSOH model"
             ):
                 if t_eval is None:
                     raise pybamm.SolverError(
@@ -857,7 +856,7 @@ def solve(
                 # Calculate capacity_start using the first cycle
                 if cycle_num == 1:
                     # Note capacity_start could be defined as
-                    # self.parameter_values["Nominal cell capacity [A.h]"] instead
+                    # self._parameter_values["Nominal cell capacity [A.h]"] instead
                     if "capacity" in self.experiment.termination:
                         capacity_start = all_summary_variables[0]["Capacity [A.h]"]
                         logs["start capacity"] = capacity_start
@@ -927,7 +926,7 @@ def step(
             self.build()
 
         if solver is None:
-            solver = self.solver
+            solver = self._solver
 
         if starting_solution is None:
             starting_solution = self._solution
@@ -941,14 +940,14 @@ def step(
     def _get_esoh_solver(self, calc_esoh):
         if (
             calc_esoh is False
-            or isinstance(self.model, pybamm.lead_acid.BaseModel)
-            or isinstance(self.model, pybamm.equivalent_circuit.Thevenin)
-            or self.model.options["working electrode"] != "both"
+            or isinstance(self._model, pybamm.lead_acid.BaseModel)
+            or isinstance(self._model, pybamm.equivalent_circuit.Thevenin)
+            or self._model.options["working electrode"] != "both"
         ):
             return None
 
         return pybamm.lithium_ion.ElectrodeSOHSolver(
-            self.parameter_values, self.model.param
+            self._parameter_values, self._model.param
         )
 
     def plot(self, output_variables=None, **kwargs):
@@ -973,7 +972,7 @@ def plot(self, output_variables=None, **kwargs):
             )
 
         if output_variables is None:
-            output_variables = self.output_variables
+            output_variables = self._output_variables
 
         self.quick_plot = pybamm.dynamic_plot(
             self._solution, output_variables=output_variables, **kwargs
@@ -1009,10 +1008,6 @@ def create_gif(self, number_of_images=80, duration=0.1, output_filename="plot.gi
     def model(self):
         return self._model
 
-    @model.setter
-    def model(self, model):
-        self._model = copy.copy(model)
-
     @property
     def model_with_set_params(self):
         return self._model_with_set_params
@@ -1025,26 +1020,14 @@ def built_model(self):
     def geometry(self):
         return self._geometry
 
-    @geometry.setter
-    def geometry(self, geometry):
-        self._geometry = geometry
-
     @property
     def parameter_values(self):
         return self._parameter_values
 
-    @parameter_values.setter
-    def parameter_values(self, parameter_values):
-        self._parameter_values = parameter_values.copy()
-
     @property
     def submesh_types(self):
         return self._submesh_types
 
-    @submesh_types.setter
-    def submesh_types(self, submesh_types):
-        self._submesh_types = submesh_types.copy()
-
     @property
     def mesh(self):
         return self._mesh
@@ -1053,41 +1036,25 @@ def mesh(self):
     def var_pts(self):
         return self._var_pts
 
-    @var_pts.setter
-    def var_pts(self, var_pts):
-        self._var_pts = var_pts.copy()
-
     @property
     def spatial_methods(self):
         return self._spatial_methods
 
-    @spatial_methods.setter
-    def spatial_methods(self, spatial_methods):
-        self._spatial_methods = spatial_methods.copy()
-
     @property
     def solver(self):
         return self._solver
 
-    @solver.setter
-    def solver(self, solver):
-        self._solver = solver.copy()
-
     @property
     def output_variables(self):
         return self._output_variables
 
-    @output_variables.setter
-    def output_variables(self, output_variables):
-        self._output_variables = copy.copy(output_variables)
-
     @property
     def solution(self):
         return self._solution
 
     def save(self, filename):
         """Save simulation using pickle"""
-        if self.model.convert_to_format == "python":
+        if self._model.convert_to_format == "python":
             # We currently cannot save models in the 'python' format
             raise NotImplementedError(
                 """