Repository: pybamm-team/PyBaMM

Issue:
Introduction of 3 electrode setup
### Description

The feature involves additional potential variables. The first one equates to a potential difference between the current collector of the anode and the middle of separator called "Negative electrode potential vs separator reference". The second one is the potential difference between the middle of the separator and the current collector of the positive electrode called "Positive electrode potential vs sparator reference".


### Motivation

Motivation for this comes from the 3-electrode setup in the lab. 

In a 2-electrode setup in the lab one can not resolve the individual SOCs of each electrode. You only get the potential difference between both current collectors. In a 3-electrode setup a reference electrode is introduced in the middle of the separator for example a little place of lithium or a whole lithium reference ring around the separator. With the reference electrode, the potentials of both electrodes can be monitored separately. These potentials can for example be used for plating detection on lithium ion graphite anodes.

Modelling this feature directly in PyBaMM is not yet possible. There are already a lot of fancy ways in pybamm to resolve processes on the electrodes directly for example different overpotential contributions. These though can not be compared directly to measured data in the lab. With the mentioned feature there is a way to describe these electrode potentials in PyBaMM and relate them directly to measured potentials in the lab.  The additional potential variables could help a great deal for fitting models better to measurements and increase model validity.


### Possible Implementation

As a rough draft for my own simulations i implemented the following lines in "base_electrolyte_conductivity". It could be placed in the current-collector submodel or maybe also "base electrode" but since the separator electrolyte potential is needed this would involve including a new submodel at these places. For now these lines worked for me:

```
        phi_e_s_mid = pybamm.x_average(phi_e_s)  # CAVEAT
        phi_e_s_mid_dim = -param.U_n_ref + phi_e_s_mid * param.potential_scale
        v_n = pybamm.boundary_value(variables["Negative electrode potential [V]"], 'left') \
              - phi_e_s_mid_dim
        v_p = pybamm.boundary_value(variables["Positive electrode potential [V]"], 'right') \
              - phi_e_s_mid_dim
        variables.update(
            {
                "Negative electrode potential vs reference [V]": v_n,
                "Positive electrode potential vs reference [V]": v_p,
            }
        )
```
There is one caveat though. For a point of reference we need a single value (here termed "phi_e_s_mid"). Typically this would just be the evaluation of the electrolyte in the middle of the separator but until now i did not manage to get it working. I tried to do it with `phi_e_s.evaluate(param.L_n+param.L_s/2)' but it didn't want to be evaluated as it said the method "evaluate is not impleted for "Current function"'. I did not try it with another approach as this wouldn't be working for my system.




### Additional context

Provided there is feedback and discussion and it's not too much work to put the feature inside while making sure that everything else in the code still works i would implement this and create a pull request. 

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index bd9c751af..a7bd875f6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,8 @@
 
 ## Features
 
+- Added a new unary operator, `EvaluateAt`, that evaluates a spatial variable at a given position ([#3573](https://github.com/pybamm-team/PyBaMM/pull/3573))
+- Added a method, `insert_reference_electrode`, to `pybamm.lithium_ion.BaseModel` that insert a reference electrode to measure the electrolyte potential at a given position in space and adds new variables that mimic a 3E cell setup. ([#3573](https://github.com/pybamm-team/PyBaMM/pull/3573))
 - Serialisation added so models can be written to/read from JSON ([#3397](https://github.com/pybamm-team/PyBaMM/pull/3397))
 
 ## Bug fixes
diff --git a/docs/source/api/expression_tree/unary_operator.rst b/docs/source/api/expression_tree/unary_operator.rst
index ad5bb0a48..e6a3cbe55 100644
--- a/docs/source/api/expression_tree/unary_operator.rst
+++ b/docs/source/api/expression_tree/unary_operator.rst
@@ -34,6 +34,9 @@ Unary Operators
 .. autoclass:: pybamm.Mass
   :members:
 
+.. autoclass:: pybamm.BoundaryMass
+  :members:
+
 .. autoclass:: pybamm.Integral
   :members:
 
@@ -58,6 +61,9 @@ Unary Operators
 .. autoclass:: pybamm.BoundaryGradient
   :members:
 
+.. autoclass:: pybamm.EvaluateAt
+  :members:
+
 .. autoclass:: pybamm.UpwindDownwind
   :members:
 
diff --git a/docs/source/examples/index.rst b/docs/source/examples/index.rst
index 4afaa6eee..e025ea71b 100644
--- a/docs/source/examples/index.rst
+++ b/docs/source/examples/index.rst
@@ -65,6 +65,7 @@ The notebooks are organised into subfolders, and can be viewed in the galleries
     notebooks/models/rate-capability.ipynb
     notebooks/models/saving_models.ipynb
     notebooks/models/SEI-on-cracks.ipynb
+    notebooks/models/simulate-3E-cell.ipynb
     notebooks/models/simulating-ORegan-2022-parameter-set.ipynb
     notebooks/models/SPM.ipynb
     notebooks/models/SPMe.ipynb
diff --git a/docs/source/examples/notebooks/models/simulate-3E-cell.ipynb b/docs/source/examples/notebooks/models/simulate-3E-cell.ipynb
new file mode 100644
index 000000000..c92cb5346
--- /dev/null
+++ b/docs/source/examples/notebooks/models/simulate-3E-cell.ipynb
@@ -0,0 +1,210 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Simulating a 3E cell\n",
+    "\n",
+    "In this notebook we show how to insert a reference electrode to mimic a 3E cell."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Note: you may need to restart the kernel to use updated packages.\n"
+     ]
+    }
+   ],
+   "source": [
+    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
+    "import pybamm"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We first load a model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model = pybamm.lithium_ion.DFN()"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we use the helper function `insert_reference_electrode` to insert a reference electrode into the model. This function takes the position of the reference electrode as an optional argument. If no position is given, the reference electrode is inserted at the midpoint of the separator. The helper function adds the new variables \"Reference electrode potential [V]\", \"Negative electrode 3E potential [V]\" and \"Positive electrode 3E potential [V]\" to the model.\n",
+    "\n",
+    "In this example we will explicitly pass a position to show how it is done"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "L_n = model.param.n.L  # Negative electrode thickness [m]\n",
+    "L_s = model.param.s.L  # Separator thickness [m]\n",
+    "L_ref = L_n + L_s / 2  # Reference electrode position [m]\n",
+    "\n",
+    "model.insert_reference_electrode(L_ref)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we can set up a simulation and solve the model as usual"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<pybamm.solvers.solution.Solution at 0x169dc5950>"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "sim = pybamm.Simulation(model)\n",
+    "sim.solve([0, 3600])"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's plot a comparison of the 3E potentials and the potential difference between the solid and electrolyte phases at the electrode/separator interfaces"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "d1f4e1ed03764660b87ee56b135b24a7",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=1.0, step=0.01), Output()), _dom_classes=('w…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<pybamm.plotting.quick_plot.QuickPlot at 0x169816a90>"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "sim.plot(\n",
+    "    [\n",
+    "        [\n",
+    "            \"Negative electrode surface potential difference at separator interface [V]\",\n",
+    "            \"Negative electrode 3E potential [V]\",\n",
+    "        ],\n",
+    "        [\n",
+    "            \"Positive electrode surface potential difference at separator interface [V]\",\n",
+    "            \"Positive electrode 3E potential [V]\",\n",
+    "        ],\n",
+    "        \"Voltage [V]\",\n",
+    "    ]\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n",
+      "[2] Marc Doyle, Thomas F. Fuller, and John Newman. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. Journal of the Electrochemical society, 140(6):1526–1533, 1993. doi:10.1149/1.2221597.\n",
+      "[3] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n",
+      "[4] Scott G. Marquis, Valentin Sulzer, Robert Timms, Colin P. Please, and S. Jon Chapman. An asymptotic derivation of a single particle model with electrolyte. Journal of The Electrochemical Society, 166(15):A3693–A3706, 2019. doi:10.1149/2.0341915jes.\n",
+      "[5] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "pybamm.print_citations()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "venv",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.6"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "9ff3d0c7e37de5f5aa47f4f719e4c84fc6cba7b39c571a05173422444e82fa58"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/source/examples/notebooks/models/simulating-ORegan-2022-parameter-set.ipynb b/docs/source/examples/notebooks/models/simulating-ORegan-2022-parameter-set.ipynb
index 7eb647fc9..f20f38560 100644
--- a/docs/source/examples/notebooks/models/simulating-ORegan-2022-parameter-set.ipynb
+++ b/docs/source/examples/notebooks/models/simulating-ORegan-2022-parameter-set.ipynb
@@ -163,7 +163,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3.7.4 ('dev': venv)",
+   "display_name": "dev",
    "language": "python",
    "name": "python3"
   },
@@ -177,7 +177,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.9.16"
   },
   "toc": {
    "base_numbering": 1,
@@ -194,7 +194,7 @@
   },
   "vscode": {
    "interpreter": {
-    "hash": "0f0e5a277ebcf03e91e138edc3d4774b5dee64e7d6640c0d876f99a9f6b2a4dc"
+    "hash": "bca2b99bfac80e18288b793d52fa0653ab9b5fe5d22e7b211c44eb982a41c00c"
    }
   }
  },
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index cab7914cd..f9b809d12 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -54,6 +54,7 @@ from .util import (
 from .logger import logger, set_logging_level, get_new_logger
 from .settings import settings
 from .citations import Citations, citations, print_citations
+
 #
 # Classes for the Expression Tree
 #
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index bb6e678f4..62110b167 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -865,6 +865,10 @@ class Discretisation(object):
                 return child_spatial_method.boundary_value_or_flux(
                     symbol, disc_child, self.bcs
                 )
+            elif isinstance(symbol, pybamm.EvaluateAt):
+                return child_spatial_method.evaluate_at(
+                    symbol, disc_child, symbol.position
+                )
             elif isinstance(symbol, pybamm.UpwindDownwind):
                 direction = symbol.name  # upwind or downwind
                 return spatial_method.upwind_or_downwind(
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 270568581..319429183 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -325,7 +325,14 @@ class Index(UnaryOperator):
     def set_id(self):
         """See :meth:`pybamm.Symbol.set_id()`"""
         self._id = hash(
-            (self.__class__, self.name, self.slice.start, self.slice.stop, self.children[0].id, *tuple(self.domain))
+            (
+                self.__class__,
+                self.name,
+                self.slice.start,
+                self.slice.stop,
+                self.children[0].id,
+                *tuple(self.domain),
+            )
         )
 
     def _unary_evaluate(self, child):
@@ -465,7 +472,9 @@ class Divergence(SpatialOperator):
 
     def _sympy_operator(self, child):
         """Override :meth:`pybamm.UnaryOperator._sympy_operator`"""
-        sympy_Divergence = have_optional_dependency("sympy.vector.operators", "Divergence")
+        sympy_Divergence = have_optional_dependency(
+            "sympy.vector.operators", "Divergence"
+        )
         return sympy_Divergence(child)
 
 
@@ -616,7 +625,18 @@ class Integral(SpatialOperator):
     def set_id(self):
         """See :meth:`pybamm.Symbol.set_id()`"""
         self._id = hash(
-            (self.__class__, self.name, *tuple([integration_variable.id for integration_variable in self.integration_variable]), self.children[0].id, *tuple(self.domain))
+            (
+                self.__class__,
+                self.name,
+                *tuple(
+                    [
+                        integration_variable.id
+                        for integration_variable in self.integration_variable
+                    ]
+                ),
+                self.children[0].id,
+                *tuple(self.domain),
+            )
         )
 
     def _unary_new_copy(self, child):
@@ -757,7 +777,13 @@ class DefiniteIntegralVector(SpatialOperator):
     def set_id(self):
         """See :meth:`pybamm.Symbol.set_id()`"""
         self._id = hash(
-            (self.__class__, self.name, self.vector_type, self.children[0].id, *tuple(self.domain))
+            (
+                self.__class__,
+                self.name,
+                self.vector_type,
+                self.children[0].id,
+                *tuple(self.domain),
+            )
         )
 
     def _unary_new_copy(self, child):
@@ -851,7 +877,13 @@ class DeltaFunction(SpatialOperator):
     def set_id(self):
         """See :meth:`pybamm.Symbol.set_id()`"""
         self._id = hash(
-            (self.__class__, self.name, self.side, self.children[0].id, *tuple([(k, tuple(v)) for k, v in self.domains.items()]))
+            (
+                self.__class__,
+                self.name,
+                self.side,
+                self.children[0].id,
+                *tuple([(k, tuple(v)) for k, v in self.domains.items()]),
+            )
         )
 
     def _evaluates_on_edges(self, dimension):
@@ -872,7 +904,8 @@ class DeltaFunction(SpatialOperator):
 
 class BoundaryOperator(SpatialOperator):
     """
-    A node in the expression tree which gets the boundary value of a variable.
+    A node in the expression tree which gets the boundary value of a variable on its
+    primary domain.
 
     Parameters
     ----------
@@ -909,7 +942,13 @@ class BoundaryOperator(SpatialOperator):
     def set_id(self):
         """See :meth:`pybamm.Symbol.set_id()`"""
         self._id = hash(
-            (self.__class__, self.name, self.side, self.children[0].id, *tuple([(k, tuple(v)) for k, v in self.domains.items()]))
+            (
+                self.__class__,
+                self.name,
+                self.side,
+                self.children[0].id,
+                *tuple([(k, tuple(v)) for k, v in self.domains.items()]),
+            )
         )
 
     def _unary_new_copy(self, child):
@@ -923,7 +962,8 @@ class BoundaryOperator(SpatialOperator):
 
 class BoundaryValue(BoundaryOperator):
     """
-    A node in the expression tree which gets the boundary value of a variable.
+    A node in the expression tree which gets the boundary value of a variable on its
+    primary domain.
 
     Parameters
     ----------
@@ -998,7 +1038,8 @@ class ExplicitTimeIntegral(UnaryOperator):
 
 class BoundaryGradient(BoundaryOperator):
     """
-    A node in the expression tree which gets the boundary flux of a variable.
+    A node in the expression tree which gets the boundary flux of a variable on its
+    primary domain.
 
     Parameters
     ----------
@@ -1012,6 +1053,59 @@ class BoundaryGradient(BoundaryOperator):
         super().__init__("boundary flux", child, side)
 
 
+class EvaluateAt(SpatialOperator):
+    """
+    A node in the expression tree which evaluates a symbol at a given position in space
+    in its primary domain. Currently this is only implemented for 1D primary domains.
+
+    Parameters
+    ----------
+    child : :class:`pybamm.Symbol`
+        The variable to evaluate
+    position : :class:`pybamm.Symbol`
+        The position in space on the symbol's primary domain at which to evaluate
+        the symbol.
+    """
+
+    def __init__(self, child, position):
+        self.position = position
+
+        # "evaluate at" of a child takes the primary domain from secondary domain
+        # of the child
+        # tertiary auxiliary domain shift down to secondary, quarternary to tertiary
+        domains = {
+            "primary": child.domains["secondary"],
+            "secondary": child.domains["tertiary"],
+            "tertiary": child.domains["quaternary"],
+        }
+
+        super().__init__("evaluate", child, domains)
+
+    def set_id(self):
+        """See :meth:`pybamm.Symbol.set_id()`"""
+        self._id = hash(
+            (
+                self.__class__,
+                self.name,
+                self.position,
+                self.children[0].id,
+                *tuple([(k, tuple(v)) for k, v in self.domains.items()]),
+            )
+        )
+
+    def _unary_jac(self, child_jac):
+        """See :meth:`pybamm.UnaryOperator._unary_jac()`."""
+        return pybamm.Scalar(0)
+
+    def _unary_new_copy(self, child):
+        """See :meth:`UnaryOperator._unary_new_copy()`."""
+        return self.__class__(child, self.position)
+
+    def _evaluate_for_shape(self):
+        """See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()`"""
+        return pybamm.evaluate_for_shape_using_domain(self.domains)
+
+
 class UpwindDownwind(SpatialOperator):
     """
     A node in the expression tree representing an upwinding or downwinding operator.
diff --git a/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py b/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
index cc736d6d0..fbe19b0d4 100644
--- a/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
+++ b/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
@@ -462,3 +462,51 @@ class BaseModel(pybamm.BaseBatteryModel):
         self.submodels[
             "through-cell convection"
         ] = pybamm.convection.through_cell.NoConvection(self.param, self.options)
+
+    def insert_reference_electrode(self, position=None):
+        """
+        Insert a reference electrode to measure the electrolyte potential at a given
+        position in space. Adds model variables for the electrolyte potential at the
+        reference electrode and for the potential difference between the electrode
+        potentials measured at the electrode/current collector interface and the
+        reference electrode. Only implemented for 1D models (i.e. where the
+        'dimensionality' option is 0).
+
+        Parameters
+        ----------
+        position : :class:`pybamm.Symbol`, optional
+            The position in space at which to measure the electrolyte potential. If
+            None, defaults to the mid-point of the separator.
+        """
+        if self.options["dimensionality"] != 0:
+            raise NotImplementedError(
+                "Reference electrode can only be inserted for models where "
+                "'dimensionality' is 0. For other models, please add a reference "
+                "electrode manually."
+            )
+
+        param = self.param
+        if position is None:
+            position = param.n.L + param.s.L / 2
+
+        phi_e_ref = pybamm.EvaluateAt(
+            self.variables["Electrolyte potential [V]"], position
+        )
+        phi_p = pybamm.boundary_value(
+            self.variables["Positive electrode potential [V]"], "right"
+        )
+        variables = {
+            "Positive electrode 3E potential [V]": phi_p - phi_e_ref,
+            "Reference electrode potential [V]": phi_e_ref,
+        }
+
+        if self.options["working electrode"] == "both":
+            phi_n = pybamm.boundary_value(
+                self.variables["Negative electrode potential [V]"], "left"
+            )
+            variables.update(
+                {
+                    "Negative electrode 3E potential [V]": phi_n - phi_e_ref,
+                }
+            )
+        self.variables.update(variables)
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index d5f12f362..049910ae9 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -721,6 +721,15 @@ class ParameterValues:
             # f_a_dist in the size average needs to be processed
             if isinstance(new_symbol, pybamm.SizeAverage):
                 new_symbol.f_a_dist = self.process_symbol(new_symbol.f_a_dist)
+            # position in evaluate at needs to be processed, and should be a Scalar
+            if isinstance(new_symbol, pybamm.EvaluateAt):
+                new_symbol_position = self.process_symbol(new_symbol.position)
+                if not isinstance(new_symbol_position, pybamm.Scalar):
+                    raise ValueError(
+                        "'position' in 'EvaluateAt' must evaluate to a scalar"
+                    )
+                else:
+                    new_symbol.position = new_symbol_position
             return new_symbol
 
         # Functions
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 5c32e5a2c..636243f82 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -1023,6 +1023,46 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         return boundary_value
 
+    def evaluate_at(self, symbol, discretised_child, position):
+        """
+        Returns the symbol evaluated at a given position in space.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The boundary value or flux symbol
+        discretised_child : :class:`pybamm.StateVector`
+            The discretised variable from which to calculate the boundary value
+        position : :class:`pybamm.Scalar`
+            The point in one-dimensional space at which to evaluate the symbol.
+
+        Returns
+        -------
+        :class:`pybamm.MatrixMultiplication`
+            The variable representing the value at the given point.
+        """
+        # Get mesh nodes
+        domain = discretised_child.domain
+        mesh = self.mesh[domain]
+        nodes = mesh.nodes
+        repeats = self._get_auxiliary_domain_repeats(discretised_child.domains)
+
+        # Find the index of the node closest to the value
+        index = np.argmin(np.abs(nodes - position.value))
+
+        # Create a sparse matrix with a 1 at the index
+        sub_matrix = csr_matrix(([1], ([0], [index])), shape=(1, mesh.npts))
+        # repeat across auxiliary domains
+        matrix = csr_matrix(kron(eye(repeats), sub_matrix))
+
+        # Index into the discretised child
+        out = pybamm.Matrix(matrix) @ discretised_child
+
+        # `EvaluateAt` removes domain
+        out.clear_domains()
+
+        return out
+
     def process_binary_operators(self, bin_op, left, right, disc_left, disc_right):
         """Discretise binary operators in model equations.  Performs appropriate
         averaging of diffusivities if one of the children is a gradient operator, so
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index acb0227bc..a461d6c15 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -331,7 +331,7 @@ class SpatialMethod:
 
     def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
         """
-        Returns the boundary value or flux using the approriate expression for the
+        Returns the boundary value or flux using the appropriate expression for the
         spatial method. To do this, we create a sparse vector 'bv_vector' that extracts
         either the first (for side="left") or last (for side="right") point from
         'discretised_child'.
@@ -377,6 +377,26 @@ class SpatialMethod:
         out.clear_domains()
         return out
 
+    def evaluate_at(self, symbol, discretised_child, position):
+        """
+        Returns the symbol evaluated at a given position in space.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The boundary value or flux symbol
+        discretised_child : :class:`pybamm.StateVector`
+            The discretised variable from which to calculate the boundary value
+        position : :class:`pybamm.Scalar`
+            The point in one-dimensional space at which to evaluate the symbol.
+
+        Returns
+        -------
+        :class:`pybamm.MatrixMultiplication`
+            The variable representing the value at the given point.
+        """
+        raise NotImplementedError
+
     def mass_matrix(self, symbol, boundary_conditions):
         """
         Calculates the mass matrix for a spatial method.


Test patch:
diff --git a/tests/unit/test_expression_tree/test_operations/test_copy.py b/tests/unit/test_expression_tree/test_operations/test_copy.py
index 0340e56bb..6800f9092 100644
--- a/tests/unit/test_expression_tree/test_operations/test_copy.py
+++ b/tests/unit/test_expression_tree/test_operations/test_copy.py
@@ -60,6 +60,7 @@ class TestCopy(TestCase):
             pybamm.maximum(a, b),
             pybamm.SparseStack(mat, mat),
             pybamm.Equality(a, b),
+            pybamm.EvaluateAt(a, 0),
         ]:
             self.assertEqual(symbol, symbol.new_copy())
 
diff --git a/tests/unit/test_expression_tree/test_operations/test_jac.py b/tests/unit/test_expression_tree/test_operations/test_jac.py
index c6e04d331..503e7321e 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac.py
@@ -236,6 +236,12 @@ class TestJacobian(TestCase):
         jac = ind.jac(vec).evaluate(y=np.linspace(0, 2, 5)).toarray()
         np.testing.assert_array_equal(jac, np.array([[0, 0, 0, 0, 0]]))
 
+    def test_evaluate_at(self):
+        y = pybamm.StateVector(slice(0, 4))
+        expr = pybamm.EvaluateAt(y, 2)
+        jac = expr.jac(y).evaluate(y=np.linspace(0, 2, 4))
+        np.testing.assert_array_equal(jac, 0)
+
     def test_jac_of_number(self):
         """Jacobian of a number should be zero"""
         a = pybamm.Scalar(1)
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 3c34db7dc..6ae6b62d0 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -454,6 +454,11 @@ class TestUnaryOperators(TestCase):
             pybamm.Index(vec, 5)
         pybamm.settings.debug_mode = debug_mode
 
+    def test_evaluate_at(self):
+        a = pybamm.Symbol("a", domain=["negative electrode"])
+        f = pybamm.EvaluateAt(a, 1)
+        self.assertEqual(f.position, 1)
+
     def test_upwind_downwind(self):
         # upwind of scalar symbol should fail
         a = pybamm.Symbol("a")
@@ -673,9 +678,10 @@ class TestUnaryOperators(TestCase):
         self.assertFalse((2 * a).is_constant())
 
     def test_to_equation(self):
-
         sympy = have_optional_dependency("sympy")
-        sympy_Divergence = have_optional_dependency("sympy.vector.operators", "Divergence")
+        sympy_Divergence = have_optional_dependency(
+            "sympy.vector.operators", "Divergence"
+        )
         sympy_Gradient = have_optional_dependency("sympy.vector.operators", "Gradient")
 
         a = pybamm.Symbol("a", domain="negative particle")
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
index 315896b29..fbc916d4a 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
@@ -29,6 +29,25 @@ class TestBaseLithiumIonModel(TestCase):
         )
         os.chdir(cwd)
 
+    def test_insert_reference_electrode(self):
+        model = pybamm.lithium_ion.SPM()
+        model.insert_reference_electrode()
+        self.assertIn("Negative electrode 3E potential [V]", model.variables)
+        self.assertIn("Positive electrode 3E potential [V]", model.variables)
+        self.assertIn("Reference electrode potential [V]", model.variables)
+
+        model = pybamm.lithium_ion.SPM({"working electrode": "positive"})
+        model.insert_reference_electrode()
+        self.assertNotIn("Negative electrode potential [V]", model.variables)
+        self.assertIn("Positive electrode 3E potential [V]", model.variables)
+        self.assertIn("Reference electrode potential [V]", model.variables)
+
+        model = pybamm.lithium_ion.SPM({"dimensionality": 2})
+        with self.assertRaisesRegex(
+            NotImplementedError, "Reference electrode can only be"
+        ):
+            model.insert_reference_electrode()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index 37ec89068..40964e8d6 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -267,6 +267,15 @@ class TestParameterValues(TestCase):
         self.assertEqual(processed_a.value, 4)
         self.assertEqual(processed_x, x)
 
+        # process EvaluateAt
+        evaluate_at = pybamm.EvaluateAt(x, aa)
+        processed_evaluate_at = parameter_values.process_symbol(evaluate_at)
+        self.assertIsInstance(processed_evaluate_at, pybamm.EvaluateAt)
+        self.assertEqual(processed_evaluate_at.children[0], x)
+        self.assertEqual(processed_evaluate_at.position, 4)
+        with self.assertRaisesRegex(ValueError, "'position' in 'EvaluateAt'"):
+            parameter_values.process_symbol(pybamm.EvaluateAt(x, x))
+
         # process broadcast
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         broad = pybamm.PrimaryBroadcast(a, whole_cell)
diff --git a/tests/unit/test_spatial_methods/test_base_spatial_method.py b/tests/unit/test_spatial_methods/test_base_spatial_method.py
index 37b4eb6a0..d48ea69a7 100644
--- a/tests/unit/test_spatial_methods/test_base_spatial_method.py
+++ b/tests/unit/test_spatial_methods/test_base_spatial_method.py
@@ -36,6 +36,8 @@ class TestSpatialMethod(TestCase):
             spatial_method.delta_function(None, None)
         with self.assertRaises(NotImplementedError):
             spatial_method.internal_neumann_condition(None, None, None, None)
+        with self.assertRaises(NotImplementedError):
+            spatial_method.evaluate_at(None, None, None)
 
     def test_get_auxiliary_domain_repeats(self):
         # Test the method to read number of repeats from auxiliary domains
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
index 91a5b7004..16a3bbde2 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
@@ -551,6 +551,30 @@ class TestFiniteVolume(TestCase):
         disc = pybamm.Discretisation(mesh, spatial_methods)
         disc.process_model(model)
 
+    def test_evaluate_at(self):
+        mesh = get_p2d_mesh_for_testing()
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume(),
+            "negative particle": pybamm.FiniteVolume(),
+            "positive particle": pybamm.FiniteVolume(),
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        n = mesh["negative electrode"].npts
+        var = pybamm.StateVector(slice(0, n), domain="negative electrode")
+
+        idx = 3
+        position = pybamm.Scalar(mesh["negative electrode"].nodes[idx])
+        evaluate_at = pybamm.EvaluateAt(var, position)
+        evaluate_at_disc = disc.process_symbol(evaluate_at)
+
+        self.assertIsInstance(evaluate_at_disc, pybamm.MatrixMultiplication)
+        self.assertIsInstance(evaluate_at_disc.left, pybamm.Matrix)
+        self.assertIsInstance(evaluate_at_disc.right, pybamm.StateVector)
+
+        y = np.arange(n)[:, np.newaxis]
+        self.assertEqual(evaluate_at_disc.evaluate(y=y), y[idx])
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")