#546 redo charge results with delta function

examples/scripts/compare_lead_acid.py

@@ -28,7 +28,7 @@
 
 # load parameter values and process models and geometry
 param = models[0].default_parameter_values
-param.update({"Typical current [A]": -20, "Initial State of Charge": 0.5})
+param.update({"Typical current [A]": -17 * 4, "Initial State of Charge": 0.5})
 for model in models:
     param.process_model(model)
 

pybamm/spatial_methods/finite_volume.py

@@ -341,7 +341,10 @@ def indefinite_integral_matrix_edges(self, domain):
     def delta_function(self, symbol, discretised_symbol):
         """
         Delta function. Implemented as a vector whose only non-zero element is the
-        first (if symbol.side = "left") or last (if symbol.side = "right").
+        first (if symbol.side = "left") or last (if symbol.side = "right"), with
+        appropriate value so that the integral of the delta function across the whole
+        domain is the same as the integral of the discretised symbol across the whole
+        domain.
 
         See :meth:`pybamm.SpatialMethod.delta_function`
         """
@@ -359,13 +362,15 @@ def delta_function(self, symbol, discretised_symbol):
             dx = submesh_list[0].d_nodes[-1]
             sub_matrix = csr_matrix(([1], ([prim_pts - 1], [0])), shape=(prim_pts, 1))
 
+        # Calculate domain width, to make sure that the integral of the delta function
+        # is the same as the integral of the child
+        domain_width = submesh_list[0].edges[-1] - submesh_list[0].edges[0]
         # Generate full matrix from the submatrix
         # Convert to csr_matrix so that we can take the index (row-slicing), which is
         # not supported by the default kron format
         # Note that this makes column-slicing inefficient, but this should not be an
         # issue
         matrix = kron(eye(sec_pts), sub_matrix).toarray()
-        domain_width = submesh_list[0].edges[-1] - submesh_list[0].edges[0]
 
         # Return delta function, keep domains
         delta_fn = pybamm.Matrix(domain_width / dx * matrix) * discretised_symbol

results/2019_09_sulzer_thesis/lead_acid_charge.py

@@ -16,7 +16,7 @@
 
 
 def plot_voltages(all_variables, t_eval):
-    Crates = [-0.1, -0.2, -0.5, -1, -2, -5]
+    Crates = [-0.1, -0.2, -0.5, -1, -2, -4]
     all_variables = {k: v for k, v in all_variables.items() if k in Crates}
     shared_plotting.plot_voltages(
         all_variables, t_eval, linestyles=["k-", "g--", "r-."]
@@ -27,7 +27,7 @@ def plot_voltages(all_variables, t_eval):
 
 
 def plot_interfacial_currents(all_variables, t_eval):
-    Crates = [-0.1, -2, -5]
+    Crates = [-0.1, -2, -4]
     all_variables = {Crate: v for Crate, v in all_variables.items() if Crate in Crates}
     file_name = "charge_interfacial_current_density_comparison.eps"
     output_vars = [
@@ -59,7 +59,7 @@ def plot_interfacial_currents(all_variables, t_eval):
 
 def plot_variables(all_variables, t_eval):
     # Set up
-    Crates = [-0.1, -2, -5]
+    Crates = [-0.1, -2, -4]
     times = np.linspace(0, 2, 4)
     var_file_names = {
         "Electrolyte concentration [Molar]"
@@ -86,7 +86,7 @@ def plot_variables(all_variables, t_eval):
 
 
 def plot_voltage_components(all_variables, t_eval):
-    Crates = [-0.1, -2, -5]
+    Crates = [-0.1, -2, -4]
     model = "Composite"
     shared_plotting.plot_voltage_components(all_variables, t_eval, model, Crates)
     file_name = "charge_voltage_components.eps"
@@ -97,9 +97,7 @@ def plot_voltage_components(all_variables, t_eval):
 def charge_states(compute):
     if compute:
         models = [
-            pybamm.lead_acid.Full(
-                {"side reactions": ["oxygen"]}, name="Full"
-            ),
+            pybamm.lead_acid.Full({"side reactions": ["oxygen"]}, name="Full"),
             pybamm.lead_acid.LOQS(
                 {"surface form": "algebraic", "side reactions": ["oxygen"]}, name="LOQS"
             ),
@@ -108,7 +106,7 @@ def charge_states(compute):
                 name="Composite",
             ),
         ]
-        Crates = [-0.1, -0.2, -0.5, -1, -2, -5]
+        Crates = [-0.1, -0.2, -0.5, -1, -2, -4, -5]
         t_eval = np.linspace(0, 3, 100)
         extra_parameter_values = {
             "Positive electrode"

tests/unit/test_spatial_methods/test_finite_volume.py

@@ -1377,9 +1377,7 @@ def test_delta_function(self):
         )
         self.assertIsInstance(delta_fn_left_disc, pybamm.Multiplication)
         self.assertIsInstance(delta_fn_left_disc.left, pybamm.Matrix)
-        np.testing.assert_array_equal(
-            delta_fn_left_disc.left.evaluate().toarray()[:, 1:], 0
-        )
+        np.testing.assert_array_equal(delta_fn_left_disc.left.evaluate()[:, 1:], 0)
         self.assertEqual(delta_fn_left_disc.shape, y.shape)
         # Right
         self.assertEqual(delta_fn_right_disc.domain, delta_fn_right.domain)
@@ -1388,17 +1386,17 @@ def test_delta_function(self):
         )
         self.assertIsInstance(delta_fn_right_disc, pybamm.Multiplication)
         self.assertIsInstance(delta_fn_right_disc.left, pybamm.Matrix)
-        np.testing.assert_array_equal(
-            delta_fn_right_disc.left.evaluate().toarray()[:, :-1], 0
-        )
+        np.testing.assert_array_equal(delta_fn_right_disc.left.evaluate()[:, :-1], 0)
         self.assertEqual(delta_fn_right_disc.shape, y.shape)
 
         # Value tests
+        # Delta function should integrate to the same thing as variable
         var_disc = disc.process_symbol(var)
         x = pybamm.standard_spatial_vars.x_n
         delta_fn_int_disc = disc.process_symbol(pybamm.Integral(delta_fn_left, x))
         np.testing.assert_array_equal(
-            var_disc.evaluate(y=y), np.sum(delta_fn_int_disc.evaluate(y=y))
+            var_disc.evaluate(y=y) * mesh["negative electrode"][0].edges[-1],
+            np.sum(delta_fn_int_disc.evaluate(y=y)),
         )
 
     def test_grad_div_with_bcs_on_tab(self):