From 3a1ad703778b41b72342845c34490e868130c5d2 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Thu, 26 Sep 2019 10:21:17 -0400 Subject: [PATCH] #618 fix discretisation --- pybamm/spatial_methods/finite_volume.py | 12 +++++++++--- .../unit/test_spatial_methods/test_finite_volume.py | 12 +++++------- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py index 81675eedea..9698b4ffd2 100644 --- a/pybamm/spatial_methods/finite_volume.py +++ b/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,15 +362,18 @@ 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 = 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 + delta_fn = pybamm.Matrix(domain_width / dx * matrix) * discretised_symbol delta_fn.domain = symbol.domain delta_fn.auxiliary_domains = symbol.auxiliary_domains diff --git a/tests/unit/test_spatial_methods/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume.py index 0746e4e0f7..c8215bd963 100644 --- a/tests/unit/test_spatial_methods/test_finite_volume.py +++ b/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):