pybamm-team#2188 implement evaluate at in 1D

examples/scripts/3E_cell.py

@@ -0,0 +1,53 @@
+#
+# Simulate insertion of a reference electrode in the middle of the cell
+#
+import pybamm
+
+# load model
+model = pybamm.lithium_ion.SPM()
+
+# load parameters and evaluate the mid-point of the cell
+parameter_values = pybamm.ParameterValues("Chen2020")
+L_n = model.param.n.L
+L_s = model.param.s.L
+L_mid = parameter_values.evaluate(L_n + L_s / 2)
+
+# extract the potential in the negative and positive electrode at the electrode/current
+# collector interfaces
+phi_n = pybamm.boundary_value(
+    model.variables["Negative electrode potential [V]"], "left"
+)
+phi_p = pybamm.boundary_value(
+    model.variables["Positive electrode potential [V]"], "right"
+)
+
+# evaluate the electrolyte potential at the mid-point of the cell
+phi_e_mid = pybamm.EvaluateAt(model.variables["Electrolyte potential [V]"], L_mid)
+
+# add the new variables to the model
+model.variables.update(
+    {
+        "Negative electrode 3E potential [V]": phi_n - phi_e_mid,
+        "Positive electrode 3E potential [V]": phi_p - phi_e_mid,
+    }
+)
+
+# solve
+sim = pybamm.Simulation(model)
+sim.solve([0, 3600])
+
+# plot a comparison of the 3E potential and the potential difference between the solid
+# and electrolyte phases at the electrode/separator interfaces
+sim.plot(
+    [
+        [
+            "Negative electrode surface potential difference at separator interface [V]",
+            "Negative electrode 3E potential [V]",
+        ],
+        [
+            "Positive electrode surface potential difference at separator interface [V]",
+            "Positive electrode 3E potential [V]",
+        ],
+        "Voltage [V]",
+    ]
+)

pybamm/__init__.py

@@ -54,6 +54,7 @@
 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
 #

pybamm/discretisations/discretisation.py

@@ -865,6 +865,10 @@ def _process_symbol(self, symbol):
                 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.value
+                )
             elif isinstance(symbol, pybamm.UpwindDownwind):
                 direction = symbol.name  # upwind or downwind
                 return spatial_method.upwind_or_downwind(

pybamm/expression_tree/unary_operators.py

@@ -325,7 +325,14 @@ def _unary_jac(self, child_jac):
     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 @@ def _unary_new_copy(self, child):
 
     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 @@ def integration_variable(self):
     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 @@ def __init__(self, child, vector_type="row"):
     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):
@@ -827,6 +853,18 @@ def _evaluates_on_edges(self, dimension):
         return False
 
 
+class ExplicitTimeIntegral(UnaryOperator):
+    def __init__(self, children, initial_condition):
+        super().__init__("explicit time integral", children)
+        self.initial_condition = initial_condition
+
+    def _unary_new_copy(self, child):
+        return self.__class__(child, self.initial_condition)
+
+    def is_constant(self):
+        return False
+
+
 class DeltaFunction(SpatialOperator):
     """
     Delta function. Currently can only be implemented at the edge of a domain.
@@ -851,7 +889,13 @@ def __init__(self, child, side, domain):
     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):
@@ -909,7 +953,13 @@ def __init__(self, name, child, side):
     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):
@@ -1012,6 +1062,51 @@ def __init__(self, child, side):
         super().__init__("boundary flux", child, side)
 
 
+class EvaluateAt(SpatialOperator):
+    """
+    A node in the expression tree which evaluates a symbol at a given position. Only
+    implemented for variables that depend on a single spatial dimension.
+
+    Parameters
+    ----------
+    child : :class:`pybamm.Symbol`
+        The variable whose boundary value to take
+    value : float
+        The point in one-dimensional space at which to evaluate the symbol.
+    """
+
+    def __init__(self, child, value):
+        self.value = value
+
+        super().__init__("evaluate", child)
+
+        # evaluating removes the domain
+        self.clear_domains()
+
+    def set_id(self):
+        """See :meth:`pybamm.Symbol.set_id()`"""
+        self._id = hash(
+            (
+                self.__class__,
+                self.name,
+                self.value,
+                self.children[0].id,
+            )
+        )
+
+    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.value)
+
+    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.

pybamm/spatial_methods/finite_volume.py

@@ -1023,6 +1023,50 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
 
         return boundary_value
 
+    def evaluate_at(self, symbol, discretised_child, value):
+        """
+        Returns the symbol evaluated at a given position in space. In the Finite
+        Volume method, the symbol is evaluated at the nearest node to the given value.
+
+        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
+        value : float
+            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.
+        """
+        # Check dimension
+        if self._get_auxiliary_domain_repeats(discretised_child.domains) > 1:
+            raise NotImplementedError(
+                "'EvaluateAt' is only implemented for 1D variables."
+            )
+
+        # Get mesh nodes
+        domain = discretised_child.domain
+        mesh = self.mesh[domain]
+        nodes = mesh.nodes
+
+        # Find the index of the node closest to the value
+        index = np.argmin(np.abs(nodes - value))
+
+        # Create a sparse matrix with a 1 at the index
+        matrix = csr_matrix(([1], ([0], [index])), shape=(1, mesh.npts))
+
+        # 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

pybamm/spatial_methods/spatial_method.py

@@ -331,7 +331,7 @@ def internal_neumann_condition(
 
     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 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
         out.clear_domains()
         return out
 
+    def evaluate_at(self, symbol, discretised_child, value):
+        """
+        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
+        value : float
+            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.

tests/unit/test_expression_tree/test_unary_operators.py

@@ -454,6 +454,11 @@ def test_index(self):
             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.value, 1)
+
     def test_upwind_downwind(self):
         # upwind of scalar symbol should fail
         a = pybamm.Symbol("a")
@@ -673,9 +678,10 @@ def test_not_constant(self):
         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")

tests/unit/test_spatial_methods/test_base_spatial_method.py

@@ -36,6 +36,8 @@ def test_basics(self):
             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

tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py

@@ -551,6 +551,41 @@ def test_full_broadcast_domains(self):
         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
+        value = mesh["negative electrode"].nodes[idx]
+        evaluate_at = pybamm.EvaluateAt(var, value)
+        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])
+
+        # test fail if not 1D
+        var = pybamm.Variable(
+            "var",
+            domain=["negative particle"],
+            auxiliary_domains={"secondary": "negative electrode"},
+        )
+        disc.set_variable_slices([var])
+        evaluate_at = pybamm.EvaluateAt(var, value)
+        with self.assertRaisesRegex(NotImplementedError, "'EvaluateAt' is only"):
+            disc.process_symbol(evaluate_at)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")