Repository: pybamm-team/PyBaMM (python)
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'
```
Pull request description: bug: make sure symbol always has required mesh attribute
# Description

Fixes #4637

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [x] Bug fix (non-breaking change which fixes an issue)

# Key checklist:

- [x] No style issues: `$ pre-commit run` (or `$ nox -s pre-commit`) (see [CONTRIBUTING.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CONTRIBUTING.md#installing-and-using-pre-commit) for how to set this up to run automatically when committing locally, in just two lines of code)
- [x] All tests pass: `$ python run-tests.py --all` (or `$ nox -s tests`)
- [x] The documentation builds: `$ python run-tests.py --doctest` (or `$ nox -s doctests`)

You can run integration tests, unit tests, and doctests together at once, using `$ python run-tests.py --quick` (or `$ nox -s quick`).

## Further checks:

- [x] Code is commented, particularly in hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works


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])}