diff --git a/examples/scripts/3E_cell.py b/examples/scripts/3E_cell.py new file mode 100644 index 0000000000..625e25a0a7 --- /dev/null +++ b/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]", + ] +) diff --git a/pybamm/__init__.py b/pybamm/__init__.py index cab7914cd9..f9b809d121 100644 --- a/pybamm/__init__.py +++ b/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 # diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py index bb6e678f4c..09f0e37496 100644 --- a/pybamm/discretisations/discretisation.py +++ b/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( diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py index 2705685814..1832181754 100644 --- a/pybamm/expression_tree/unary_operators.py +++ b/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. diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py index 5c32e5a2c0..0e25a7b3fb 100644 --- a/pybamm/spatial_methods/finite_volume.py +++ b/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 diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py index acb0227bc2..4945c7e1bb 100644 --- a/pybamm/spatial_methods/spatial_method.py +++ b/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. diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py index 3c34db7dcd..f39dba335e 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 @@ 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") 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 37b4eb6a0b..d48ea69a7b 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 @@ 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 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 91a5b70044..b98cfa2abe 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,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")