#546 make discretised delta function dense

pybamm/discretisations/discretisation.py

@@ -539,7 +539,11 @@ def process_symbol(self, symbol):
         except KeyError:
             discretised_symbol = self._process_symbol(symbol)
             self._discretised_symbols[symbol.id] = discretised_symbol
-            discretised_symbol.test_shape()
+            try:
+                discretised_symbol.test_shape()
+            except:
+                import ipdb; ipdb.set_trace()
+                discretised_symbol.test_shape()
             return discretised_symbol
 
     def _process_symbol(self, symbol):

pybamm/expression_tree/broadcasts.py

@@ -124,9 +124,6 @@ class PrimaryBroadcast(Broadcast):
         child node
     broadcast_domain : iterable of str
         Primary domain for broadcast. This will become the domain of the symbol
-    broadcast_type : str, optional
-        Whether to broadcast to the full domain (primary and secondary) or only in the
-        primary direction. Default is "full".
     name : str
         name of the node
 

pybamm/expression_tree/concatenations.py

@@ -4,7 +4,7 @@
 import copy
 import numpy as np
 import pybamm
-from scipy.sparse import vstack
+from scipy.sparse import vstack, issparse
 from collections import defaultdict
 
 

pybamm/expression_tree/unary_operators.py

@@ -632,9 +632,13 @@ class DeltaFunction(SpatialOperator):
     **Extends:** :class:`SpatialOperator`
     """
 
-    def __init__(self, child, side, domain, auxiliary_domains=None):
+    def __init__(self, child, side, domain):
         self.side = side
-        super().__init__("delta function", child, domain, auxiliary_domains)
+        if child.domain != []:
+            auxiliary_domains = {"secondary": child.domain}
+        else:
+            auxiliary_domains = {}
+        super().__init__("delta_function", child, domain, auxiliary_domains)
 
     def set_id(self):
         """ See :meth:`pybamm.Symbol.set_id()` """
@@ -650,14 +654,20 @@ def evaluates_on_edges(self):
 
     def _unary_simplify(self, simplified_child):
         """ See :meth:`UnaryOperator._unary_simplify()`. """
-        return self.__class__(
-            simplified_child, self.side, self.domain, self.auxiliary_domains
-        )
+        return self.__class__(simplified_child, self.side, self.domain)
 
     def _unary_new_copy(self, child):
         """ See :meth:`UnaryOperator._unary_new_copy()`. """
-        return self.__class__(child, self.side, self.domain, self.auxiliary_domains)
+        return self.__class__(child, self.side, self.domain)
+
+    def evaluate_for_shape(self):
+        """
+        See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()`
+        """
+        child_eval = self.children[0].evaluate_for_shape()
+        vec = pybamm.evaluate_for_shape_using_domain(self.domain)
 
+        return np.outer(child_eval, vec).reshape(-1, 1)
 
 class BoundaryOperator(SpatialOperator):
     """A node in the expression tree which gets the boundary value of a variable.

pybamm/models/submodels/interface/diffusion_limited/full_diffusion_limited.py

@@ -33,9 +33,6 @@ def _get_diffusion_limited_current_density(self, variables):
                 * param.curlyD_ox
                 * pybamm.BoundaryGradient(c_ox_s, "left")
             )
-            import ipdb
-
-            ipdb.set_trace()
             j = -pybamm.DeltaFunction(
                 N_ox_neg_sep_interface / param.C_e / param.s_ox_Ox / param.l_n,
                 "right",

pybamm/spatial_methods/finite_volume.py

@@ -364,7 +364,7 @@ def delta_function(self, symbol, discretised_symbol):
         # not supported by the default kron format
         # Note that this makes column-slicing inefficient, but this should not be an
         # issue
-        matrix = csr_matrix(kron(eye(sec_pts), sub_matrix))
+        matrix = kron(eye(sec_pts), sub_matrix).toarray()
 
         # Return delta function, keep domains
         delta_fn = pybamm.Matrix(1 / dx * matrix) * discretised_symbol

tests/unit/test_expression_tree/test_simplify.py

@@ -82,7 +82,7 @@ def myfunction(x, y):
 
         # Delta function
         self.assertIsInstance(
-            (pybamm.DeltaFunction(v_neg, "right", None, None)).simplify(),
+            (pybamm.DeltaFunction(v_neg, "right", None)).simplify(),
             pybamm.DeltaFunction,
         )
 

tests/unit/test_expression_tree/test_unary_operators.py

@@ -186,7 +186,7 @@ def test_id(self):
 
     def test_delta_function(self):
         a = pybamm.Symbol("a")
-        delta_a = pybamm.DeltaFunction(a, "right", "some domain", {})
+        delta_a = pybamm.DeltaFunction(a, "right", "some domain")
         self.assertEqual(delta_a.side, "right")
         self.assertEqual(delta_a.child.id, a.id)
         self.assertFalse(delta_a.evaluates_on_edges())

tests/unit/test_parameters/test_parameter_values.py

@@ -93,7 +93,7 @@ def test_process_symbol(self):
 
         # process delta function
         aa = pybamm.Parameter("a")
-        delta_aa = pybamm.DeltaFunction(aa, "left", "some domain", {})
+        delta_aa = pybamm.DeltaFunction(aa, "left", "some domain")
         processed_delta_aa = parameter_values.process_symbol(delta_aa)
         self.assertIsInstance(processed_delta_aa, pybamm.DeltaFunction)
         self.assertEqual(processed_delta_aa.side, "left")

tests/unit/test_spatial_methods/test_finite_volume.py

@@ -1362,8 +1362,8 @@ def test_delta_function(self):
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
         var = pybamm.Variable("var")
-        delta_fn_left = pybamm.DeltaFunction(var, "left", "negative electrode", {})
-        delta_fn_right = pybamm.DeltaFunction(var, "right", "negative electrode", {})
+        delta_fn_left = pybamm.DeltaFunction(var, "left", "negative electrode")
+        delta_fn_right = pybamm.DeltaFunction(var, "right", "negative electrode")
         disc.set_variable_slices([var])
         delta_fn_left_disc = disc.process_symbol(delta_fn_left)
         delta_fn_right_disc = disc.process_symbol(delta_fn_right)