Repository: pybamm-team/PyBaMM

Issue:
[Bug]: post-processing variables require pybamm.Symbol to have a mesh
### PyBaMM Version

develop

### Python Version

3.10

### Describe the bug

Adding a variable to a pre-built pybamm model causes a failure in `pybamm.solvers.processed_variable.py::process_variable()` because it doesn't have a `mesh` attribute

### Steps to Reproduce

1. create a pybamm model (e.g. SPM)
2. manually add a new variable (model.variables.update({"name": expression}))
3. solve the model then try and access the variable (sol["name"])

### Relevant log output

```shell
  File "/home/mrobins/git/PyBOP/pybop/problems/fitting_problem.py", line 133, in evaluate
    return {
  File "/home/mrobins/git/PyBOP/pybop/problems/fitting_problem.py", line 134, in <dictcomp>
    signal: sol[signal].data
  File "/home/mrobins/git/PyBOP/env/lib/python3.10/site-packages/pybamm/solvers/solution.py", line 710, in __getitem__
    self.update(key)
  File "/home/mrobins/git/PyBOP/env/lib/python3.10/site-packages/pybamm/solvers/solution.py", line 590, in update
    self._update_variable(variable)
  File "/home/mrobins/git/PyBOP/env/lib/python3.10/site-packages/pybamm/solvers/solution.py", line 640, in _update_variable
    var = pybamm.process_variable(
  File "/home/mrobins/git/PyBOP/env/lib/python3.10/site-packages/pybamm/solvers/processed_variable.py", line 881, in process_variable
    mesh = base_variables[0].mesh
AttributeError: 'Power' object has no attribute 'mesh'
```

Discussion:
review-notebook-app[bot]: Check out this pull request on&nbsp; <a href="https://app.reviewnb.com/pybamm-team/PyBaMM/pull/4644"><img align="absmiddle"  alt="ReviewNB" height="28" class="BotMessageButtonImage" src="https://raw.githubusercontent.com/ReviewNB/support/master/images/button_reviewnb.png"/></a> 

 See visual diffs & provide feedback on Jupyter Notebooks. 

---

 <i>Powered by <a href='https://www.reviewnb.com/?utm_source=gh'>ReviewNB</a></i>
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4644?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 99.21%. Comparing base [(`72c23ea`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/72c23ea9ecfb81057b7d149e4988000b669b4e6b?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) to head [(`ea15ec9`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/ea15ec9d7767083a8e0d10e44e790dcef873461f?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).
> Report is 8 commits behind head on develop.

<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff            @@
##           develop    #4644   +/-   ##
========================================
  Coverage    99.21%   99.21%           
========================================
  Files          302      302           
  Lines        22858    22866    +8     
========================================
+ Hits         22679    22687    +8     
  Misses         179      179           
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4644?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7954585f0..508c5f305 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,7 @@
 ## Bug fixes
 
 - Remove internal use of deprecated `set_parameters` function in the `Simulation` class which caused warnings. ([#4638](https://github.com/pybamm-team/PyBaMM/pull/4638))
+- Provide default value for `Symbol.mesh` attribute to avoid errors when adding variables after discretisation. ([#4644](https://github.com/pybamm-team/PyBaMM/pull/4644))
 
 # [v24.11.2](https://github.com/pybamm-team/PyBaMM/tree/v24.11.2) - 2024-11-27
 
diff --git a/docs/source/examples/notebooks/models/pouch-cell-model.ipynb b/docs/source/examples/notebooks/models/pouch-cell-model.ipynb
index fa394a341..6b71f763e 100644
--- a/docs/source/examples/notebooks/models/pouch-cell-model.ipynb
+++ b/docs/source/examples/notebooks/models/pouch-cell-model.ipynb
@@ -301,7 +301,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -311,8 +311,6 @@
     "    pybamm.t,\n",
     "    name=\"voltage_comsol\",\n",
     ")\n",
-    "comsol_voltage.mesh = None\n",
-    "comsol_voltage.secondary_mesh = None\n",
     "comsol_phi_s_cn = get_interp_fun_curr_coll(\"phi_s_cn\")\n",
     "comsol_phi_s_cp = get_interp_fun_curr_coll(\"phi_s_cp\")\n",
     "comsol_current = get_interp_fun_curr_coll(\"current\")\n",
diff --git a/examples/scripts/compare_comsol/compare_comsol_DFN.py b/examples/scripts/compare_comsol/compare_comsol_DFN.py
index c2730e71e..a59b867e1 100644
--- a/examples/scripts/compare_comsol/compare_comsol_DFN.py
+++ b/examples/scripts/compare_comsol/compare_comsol_DFN.py
@@ -102,9 +102,6 @@ comsol_voltage = pybamm.Interpolant(
     comsol_t, np.array(comsol_variables["voltage"]), pybamm.t
 )
 
-comsol_voltage.mesh = None
-comsol_voltage.secondary_mesh = None
-
 # Create comsol model with dictionary of Matrix variables
 comsol_model = pybamm.lithium_ion.BaseModel()
 comsol_model._geometry = pybamm_model.default_geometry
diff --git a/src/pybamm/expression_tree/symbol.py b/src/pybamm/expression_tree/symbol.py
index aa9ebe66d..c4e0043f1 100644
--- a/src/pybamm/expression_tree/symbol.py
+++ b/src/pybamm/expression_tree/symbol.py
@@ -236,6 +236,9 @@ class Symbol:
         # Set domains (and hence id)
         self.domains = self.read_domain_or_domains(domain, auxiliary_domains, domains)
 
+        # mesh required for solution and processed variables classes
+        self.mesh = None
+
         self._saved_evaluates_on_edges: dict = {}
         self._print_name = None
 
diff --git a/src/pybamm/models/base_model.py b/src/pybamm/models/base_model.py
index 6cb3af2fd..d489b51f3 100644
--- a/src/pybamm/models/base_model.py
+++ b/src/pybamm/models/base_model.py
@@ -116,8 +116,6 @@ class BaseModel:
             for var in instance._variables.values():
                 if var.domain != []:
                     var.mesh = properties["mesh"][var.domain]
-                else:
-                    var.mesh = None
 
                 if var.domains["secondary"] != []:
                     var.secondary_mesh = properties["mesh"][var.domains["secondary"]]


Test patch:
diff --git a/tests/unit/test_solvers/test_processed_variable.py b/tests/unit/test_solvers/test_processed_variable.py
index b3ab85060..55fc87833 100644
--- a/tests/unit/test_solvers/test_processed_variable.py
+++ b/tests/unit/test_solvers/test_processed_variable.py
@@ -109,7 +109,6 @@ class TestProcessedVariable:
         t = pybamm.t
         y = pybamm.StateVector(slice(0, 1))
         var = t * y
-        var.mesh = None
         model = pybamm.BaseModel()
         t_sol = np.linspace(0, 1)
         y_sol = np.array([np.linspace(0, 5)])
@@ -124,7 +123,6 @@ class TestProcessedVariable:
 
         # scalar value
         var = y
-        var.mesh = None
         t_sol = np.array([0])
         y_sol = np.array([1])[:, np.newaxis]
         yp_sol = np.array([1])[:, np.newaxis]
@@ -180,7 +178,6 @@ class TestProcessedVariable:
         data = pybamm.DiscreteTimeData(data_t, data_v, "test_data")
         var = (y - data) ** order
         expected_entries = (y_sol - data_v) ** order
-        var.mesh = None
         model = pybamm.BaseModel()
         var_casadi = to_casadi(var, y_sol)
         processed_var = pybamm.process_variable(
@@ -206,7 +203,6 @@ class TestProcessedVariable:
         var = (y - data) ** order
         expected = (y_sol_interp - data_v) ** order
         expected_entries = (y_sol - data_v_interp) ** order
-        var.mesh = None
         model = pybamm.BaseModel()
         var_casadi = to_casadi(var, y_sol)
         processed_var = pybamm.process_variable(
@@ -227,7 +223,6 @@ class TestProcessedVariable:
         t = pybamm.t
         y = pybamm.StateVector(slice(0, 1))
         var = t * y
-        var.mesh = None
         t_sol = np.linspace(0, 1)
         y_sol = np.array([np.linspace(0, 5)])
         yp_sol = self._get_yps(y_sol, hermite_interp)
@@ -246,7 +241,6 @@ class TestProcessedVariable:
         y = pybamm.StateVector(slice(0, 1))
         a = pybamm.InputParameter("a")
         var = t * y * a
-        var.mesh = None
         t_sol = np.linspace(0, 1)
         y_sol = np.array([np.linspace(0, 5)])
         inputs = {"a": np.array([1.0])}
@@ -593,8 +587,6 @@ class TestProcessedVariable:
         y = pybamm.StateVector(slice(0, 1))
         var = y
         eqn = t * y
-        var.mesh = None
-        eqn.mesh = None
 
         t_sol = np.linspace(0, 1, 1000)
         y_sol = np.array([5 * t_sol])
@@ -631,10 +623,7 @@ class TestProcessedVariable:
     @pytest.mark.parametrize("hermite_interp", _hermite_args)
     def test_processed_var_0D_fixed_t_interpolation(self, hermite_interp):
         y = pybamm.StateVector(slice(0, 1))
-        var = y
         eqn = 2 * y
-        var.mesh = None
-        eqn.mesh = None
 
         t_sol = np.array([10])
         y_sol = np.array([[100]])
@@ -1191,7 +1180,6 @@ class TestProcessedVariable:
         t = pybamm.t
         y = pybamm.StateVector(slice(0, 1))
         var = t * y
-        var.mesh = None
         t_sol = np.linspace(0, 1)
         y_sol = np.array([np.linspace(0, 5)])
         var_casadi = to_casadi(var, y_sol)
@@ -1235,12 +1223,8 @@ class TestProcessedVariable:
             return sol
 
         # without spatial dependence
-        t = pybamm.t
         y = pybamm.StateVector(slice(0, 1))
         var = y
-        eqn = t * y
-        var.mesh = None
-        eqn.mesh = None
 
         sign1 = +1
         t_sol1 = np.linspace(0, 1, 100)
diff --git a/tests/unit/test_solvers/test_processed_variable_computed.py b/tests/unit/test_solvers/test_processed_variable_computed.py
index f28c053fd..a94817143 100644
--- a/tests/unit/test_solvers/test_processed_variable_computed.py
+++ b/tests/unit/test_solvers/test_processed_variable_computed.py
@@ -72,7 +72,6 @@ class TestProcessedVariableComputed:
         # without space
         y = pybamm.StateVector(slice(0, 1))
         var = y
-        var.mesh = None
         t_sol = np.array([0])
         y_sol = np.array([1])[:, np.newaxis]
         var_casadi = to_casadi(var, y_sol)
@@ -117,7 +116,6 @@ class TestProcessedVariableComputed:
         t = pybamm.t
         y = pybamm.StateVector(slice(0, 1))
         var = t * y
-        var.mesh = None
         t_sol = np.linspace(0, 1)
         y_sol = np.array([np.linspace(0, 5)])
         var_casadi = to_casadi(var, y_sol)
@@ -136,7 +134,6 @@ class TestProcessedVariableComputed:
         y = pybamm.StateVector(slice(0, 1))
         a = pybamm.InputParameter("a")
         var = t * y * a
-        var.mesh = None
         t_sol = np.linspace(0, 1)
         y_sol = np.array([np.linspace(0, 5)])
         inputs = {"a": np.array([1.0])}