Issue 609 broadcast to cell edges

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added functionality to broadcast to edges ([#891](https://github.com/pybamm-team/PyBaMM/pull/891))
 -   Added additional notebooks showing how to create and compare models ([#877](https://github.com/pybamm-team/PyBaMM/pull/877))
 -   Added `Minimum`, `Maximum` and `Sign` operators ([#876](https://github.com/pybamm-team/PyBaMM/pull/876))
 -   Added a search feature to `FuzzyDict` ([#875](https://github.com/pybamm-team/PyBaMM/pull/875))

docs/source/expression_tree/broadcasts.rst

@@ -12,3 +12,14 @@ Broadcasting Operators
 
 .. autoclass:: pybamm.SecondaryBroadcast
   :members:
+
+.. autoclass:: pybamm.FullBroadcastToEdges
+  :members:
+
+.. autoclass:: pybamm.PrimaryBroadcastToEdges
+  :members:
+
+.. autoclass:: pybamm.SecondaryBroadcastToEdges
+  :members:
+
+.. autofunction:: pybamm.ones_like

docs/source/models/submodels/particle/fast_many_particles.rst

@@ -0,0 +1,5 @@
+Fast Many Particles
+===================
+
+.. autoclass:: pybamm.particle.FastManyParticles
+    :members:

docs/source/models/submodels/particle/fast_single_particle.rst

@@ -0,0 +1,5 @@
+Fast Single Particle
+====================
+
+.. autoclass:: pybamm.particle.FastSingleParticle
+    :members:

docs/source/models/submodels/particle/fickian_many_particles.rst

@@ -0,0 +1,7 @@
+Fickian Many Particles
+======================
+
+.. autoclass:: pybamm.particle.FickianManyParticles
+    :members:
+
+

docs/source/models/submodels/particle/fickian_single_particle.rst

@@ -0,0 +1,6 @@
+Fickian Single Particle
+=======================
+
+.. autoclass:: pybamm.particle.FickianSingleParticle
+    :members:
+

docs/source/models/submodels/particle/index.rst

@@ -5,5 +5,7 @@ Particle
   :maxdepth: 1
 
   base_particle
-  fickian/index
-  fast/index
+  fickian_single_particle
+  fickian_many_particles
+  fast_single_particle
+  fast_many_particles

examples/notebooks/using-submodels.ipynb

@@ -113,7 +113,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model.submodels[\"negative particle\"] = pybamm.particle.fast.SingleParticle(model.param, \"Negative\")"
+    "model.submodels[\"negative particle\"] = pybamm.particle.FastSingleParticle(model.param, \"Negative\")"
    ]
   },
   {
@@ -370,10 +370,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model.submodels[\"negative particle\"] = pybamm.particle.fast.SingleParticle(\n",
+    "model.submodels[\"negative particle\"] = pybamm.particle.FastSingleParticle(\n",
     "    model.param, \"Negative\"\n",
     ")\n",
-    "model.submodels[\"positive particle\"] = pybamm.particle.fast.SingleParticle(\n",
+    "model.submodels[\"positive particle\"] = pybamm.particle.FastSingleParticle(\n",
     "    model.param, \"Positive\"\n",
     ")"
    ]

examples/scripts/compare_lithium_ion.py

@@ -16,7 +16,7 @@
     pybamm.set_logging_level("INFO")
 
 # load models
-options = {"thermal": "isothermal"}
+options = {"thermal": "x-lumped"}
 models = [
     pybamm.lithium_ion.SPM(options),
     pybamm.lithium_ion.SPMe(options),

examples/scripts/custom_model.py

@@ -22,10 +22,10 @@
 model.submodels["positive electrode"] = pybamm.electrode.ohm.LeadingOrder(
     model.param, "Positive"
 )
-model.submodels["negative particle"] = pybamm.particle.fast.SingleParticle(
+model.submodels["negative particle"] = pybamm.particle.FastSingleParticle(
     model.param, "Negative"
 )
-model.submodels["positive particle"] = pybamm.particle.fast.SingleParticle(
+model.submodels["positive particle"] = pybamm.particle.FastSingleParticle(
     model.param, "Positive"
 )
 model.submodels["negative interface"] = pybamm.interface.InverseButlerVolmer(

pybamm/discretisations/discretisation.py

@@ -226,7 +226,7 @@ def set_variable_slices(self, variables):
                     for child, mesh in meshes.items():
                         for domain_mesh in mesh:
                             submesh = domain_mesh[i]
-                            end += submesh.npts_for_broadcast
+                            end += submesh.npts_for_broadcast_to_nodes
                         y_slices[child.id].append(slice(start, end))
                         start = end
             else:
@@ -249,7 +249,7 @@ def _get_variable_size(self, variable):
             size = 0
             for dom in variable.domain:
                 for submesh in self.spatial_methods[dom].mesh[dom]:
-                    size += submesh.npts_for_broadcast
+                    size += submesh.npts_for_broadcast_to_nodes
             return size
 
     def _preprocess_external_variables(self, model):

pybamm/expression_tree/broadcasts.py

@@ -37,7 +37,7 @@ def __init__(
         child,
         broadcast_domain,
         broadcast_auxiliary_domains=None,
-        broadcast_type="full",
+        broadcast_type="full to nodes",
         name=None,
     ):
         # Convert child to scalar if it is a number
@@ -84,7 +84,9 @@ class PrimaryBroadcast(Broadcast):
     """
 
     def __init__(self, child, broadcast_domain, name=None):
-        super().__init__(child, broadcast_domain, broadcast_type="primary", name=name)
+        super().__init__(
+            child, broadcast_domain, broadcast_type="primary to nodes", name=name
+        )
 
     def check_and_set_domains(
         self, child, broadcast_type, broadcast_domain, broadcast_auxiliary_domains
@@ -127,8 +129,8 @@ def check_and_set_domains(
         return domain, auxiliary_domains
 
     def _unary_new_copy(self, child):
-        """ See :meth:`pybamm.UnaryOperator.simplify()`. """
-        return PrimaryBroadcast(child, self.broadcast_domain)
+        """ See :meth:`pybamm.UnaryOperator._unary_new_copy()`. """
+        return self.__class__(child, self.broadcast_domain)
 
     def _evaluate_for_shape(self):
         """
@@ -140,6 +142,18 @@ def _evaluate_for_shape(self):
         return np.outer(child_eval, vec).reshape(-1, 1)
 
 
+class PrimaryBroadcastToEdges(PrimaryBroadcast):
+    "A primary broadcast onto the edges of the domain"
+
+    def __init__(self, child, broadcast_domain, name=None):
+        name = name or "broadcast to edges"
+        super().__init__(child, broadcast_domain, name)
+        self.broadcast_type = "primary to edges"
+
+    def evaluates_on_edges(self):
+        return True
+
+
 class SecondaryBroadcast(Broadcast):
     """A node in the expression tree representing a primary broadcasting operator.
     Broadcasts in a `secondary` dimension only. That is, makes explicit copies of the
@@ -162,7 +176,9 @@ class SecondaryBroadcast(Broadcast):
     """
 
     def __init__(self, child, broadcast_domain, name=None):
-        super().__init__(child, broadcast_domain, broadcast_type="secondary", name=name)
+        super().__init__(
+            child, broadcast_domain, broadcast_type="secondary to nodes", name=name
+        )
 
     def check_and_set_domains(
         self, child, broadcast_type, broadcast_domain, broadcast_auxiliary_domains
@@ -207,7 +223,7 @@ def check_and_set_domains(
         return domain, auxiliary_domains
 
     def _unary_new_copy(self, child):
-        """ See :meth:`pybamm.UnaryOperator.simplify()`. """
+        """ See :meth:`pybamm.UnaryOperator._unary_new_copy()`. """
         return SecondaryBroadcast(child, self.broadcast_domain)
 
     def _evaluate_for_shape(self):
@@ -220,6 +236,18 @@ def _evaluate_for_shape(self):
         return np.outer(vec, child_eval).reshape(-1, 1)
 
 
+class SecondaryBroadcastToEdges(SecondaryBroadcast):
+    "A secondary broadcast onto the edges of a domain"
+
+    def __init__(self, child, broadcast_domain, name=None):
+        name = name or "broadcast to edges"
+        super().__init__(child, broadcast_domain, name)
+        self.broadcast_type = "secondary to edges"
+
+    def evaluates_on_edges(self):
+        return True
+
+
 class FullBroadcast(Broadcast):
     "A class for full broadcasts"
 
@@ -230,7 +258,7 @@ def __init__(self, child, broadcast_domain, auxiliary_domains, name=None):
             child,
             broadcast_domain,
             broadcast_auxiliary_domains=auxiliary_domains,
-            broadcast_type="full",
+            broadcast_type="full to nodes",
             name=name,
         )
 
@@ -250,7 +278,7 @@ def check_and_set_domains(
         return domain, auxiliary_domains
 
     def _unary_new_copy(self, child):
-        """ See :meth:`pybamm.UnaryOperator.simplify()`. """
+        """ See :meth:`pybamm.UnaryOperator._unary_new_copy()`. """
         return FullBroadcast(child, self.broadcast_domain, self.auxiliary_domains)
 
     def _evaluate_for_shape(self):
@@ -266,6 +294,21 @@ def _evaluate_for_shape(self):
         return child_eval * vec
 
 
+class FullBroadcastToEdges(FullBroadcast):
+    """
+    A full broadcast onto the edges of a domain (edges of primary dimension, nodes of
+    other dimensions)
+    """
+
+    def __init__(self, child, broadcast_domain, auxiliary_domains, name=None):
+        name = name or "broadcast to edges"
+        super().__init__(child, broadcast_domain, auxiliary_domains, name)
+        self.broadcast_type = "full to edges"
+
+    def evaluates_on_edges(self):
+        return True
+
+
 def ones_like(*symbols):
     """
     Create a symbol with the same shape as the input symbol and with constant value '1',

pybamm/expression_tree/independent_variable.py

@@ -104,10 +104,32 @@ def __init__(self, name, domain=None, auxiliary_domains=None, coord_sys=None):
 
     def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
-        return SpatialVariable(
+        return self.__class__(
             self.name, self.domain, self.auxiliary_domains, self.coord_sys
         )
 
 
+class SpatialVariableEdge(SpatialVariable):
+    """A node in the expression tree representing a spatial variable, which evaluates
+    on the edges
+
+    Parameters
+    ----------
+    name : str
+        name of the node (e.g. "x", "y", "z", "r", "x_n", "x_s", "x_p", "r_n", "r_p")
+    domain : iterable of str
+        list of domains that this variable is valid over (e.g. "cartesian", "spherical
+        polar")
+
+    *Extends:* :class:`Symbol`
+    """
+
+    def __init__(self, name, domain=None, auxiliary_domains=None, coord_sys=None):
+        super().__init__(name, domain, auxiliary_domains, coord_sys)
+
+    def evaluates_on_edges(self):
+        return True
+
+
 # the independent variable time
 t = Time()