From 2db150e2571eb25bbe8b16e7e493de120750220a Mon Sep 17 00:00:00 2001 From: medha-14 Date: Wed, 23 Oct 2024 18:41:20 +0530 Subject: [PATCH 1/5] added gradient squared method --- src/pybamm/spatial_methods/finite_volume.py | 27 +++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py index 5aaf7ea123..34cb437c04 100644 --- a/src/pybamm/spatial_methods/finite_volume.py +++ b/src/pybamm/spatial_methods/finite_volume.py @@ -128,6 +128,33 @@ def gradient_matrix(self, domain, domains): return pybamm.Matrix(matrix) + def gradient_squared(self, symbol, discretised_symbol, boundary_conditions): + """ + Computes the square of the gradient of a symbol. + + Parameters + ---------- + symbol : :class:`pybamm.Symbol` + The symbol for which to compute the gradient squared. + discretised_symbol : :class:`pybamm.Vector` + The discretised variable for which to compute the gradient squared. + boundary_conditions : dict + Boundary conditions for the symbol. + + Returns + ------- + :class:`pybamm.Vector` + The gradient squared of the symbol. + """ + domain = symbol.domain + gradient_matrix = self.gradient_matrix(domain, symbol.domains) + + # Compute gradient squared: (∇u)^2 = u^T L^T L u + gradient_squared_matrix = gradient_matrix.T @ gradient_matrix + gradient_squared_result = gradient_squared_matrix @ discretised_symbol + + return gradient_squared_result + def divergence(self, symbol, discretised_symbol, boundary_conditions): """Matrix-vector multiplication to implement the divergence operator. See :meth:`pybamm.SpatialMethod.divergence` From f8a506e841220dbc0f3d093ce01996c3c8de2af7 Mon Sep 17 00:00:00 2001 From: medha-14 Date: Mon, 28 Oct 2024 16:16:06 +0530 Subject: [PATCH 2/5] added test for gradient squared method --- .../test_scikit_finite_element.py | 30 +++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/tests/unit/test_spatial_methods/test_scikit_finite_element.py b/tests/unit/test_spatial_methods/test_scikit_finite_element.py index 42b282e08a..e22734fe99 100644 --- a/tests/unit/test_spatial_methods/test_scikit_finite_element.py +++ b/tests/unit/test_spatial_methods/test_scikit_finite_element.py @@ -157,6 +157,36 @@ def test_gradient(self): ans = eqn_disc.evaluate(None, 3 * y**2) np.testing.assert_array_less(0, ans) + def test_gradient_squared(self): + mesh = get_unit_2p1D_mesh_for_testing(ypts=32, zpts=32, include_particles=False) + spatial_methods = { + "macroscale": pybamm.FiniteVolume(), + "current collector": pybamm.ScikitFiniteElement(), + } + disc = pybamm.Discretisation(mesh, spatial_methods) + + var = pybamm.Variable("var", domain="current collector") + disc.set_variable_slices([var]) + + y = mesh["current collector"].coordinates[0, :] + z = mesh["current collector"].coordinates[1, :] + + gradient_squared = pybamm.grad_squared(var) + grad_squared_disc = disc.process_symbol(gradient_squared) + evaluated_gradient_squared = grad_squared_disc.evaluate( + None, 3 * y**2 + 2 * z**2 + ) + + expected_result = 36 * y**2 + 16 * z**2 + expected_result = expected_result[:, np.newaxis] + + try: + np.testing.assert_array_almost_equal( + evaluated_gradient_squared, expected_result + ) + except AssertionError as e: + print("AssertionError:", e) + def test_manufactured_solution(self): mesh = get_unit_2p1D_mesh_for_testing(ypts=32, zpts=32, include_particles=False) spatial_methods = { From a944d27ee2c758fbcc627ca9b8e0c1a45bc40655 Mon Sep 17 00:00:00 2001 From: medha-14 Date: Thu, 7 Nov 2024 11:50:38 +0530 Subject: [PATCH 3/5] modified gradient_squared to give a scalar value --- src/pybamm/spatial_methods/finite_volume.py | 10 +++++--- .../test_scikit_finite_element.py | 25 ++++++++++--------- 2 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py index 34cb437c04..96ee51361d 100644 --- a/src/pybamm/spatial_methods/finite_volume.py +++ b/src/pybamm/spatial_methods/finite_volume.py @@ -143,17 +143,19 @@ def gradient_squared(self, symbol, discretised_symbol, boundary_conditions): Returns ------- - :class:`pybamm.Vector` + float The gradient squared of the symbol. """ domain = symbol.domain gradient_matrix = self.gradient_matrix(domain, symbol.domains) - # Compute gradient squared: (∇u)^2 = u^T L^T L u + # Compute gradient squared: (∇u)^2 = u^T (L^T L) u gradient_squared_matrix = gradient_matrix.T @ gradient_matrix - gradient_squared_result = gradient_squared_matrix @ discretised_symbol + gradient_squared_result = ( + discretised_symbol.T @ gradient_squared_matrix @ discretised_symbol + ) - return gradient_squared_result + return gradient_squared_result.item() def divergence(self, symbol, discretised_symbol, boundary_conditions): """Matrix-vector multiplication to implement the divergence operator. diff --git a/tests/unit/test_spatial_methods/test_scikit_finite_element.py b/tests/unit/test_spatial_methods/test_scikit_finite_element.py index e22734fe99..3ba13c6790 100644 --- a/tests/unit/test_spatial_methods/test_scikit_finite_element.py +++ b/tests/unit/test_spatial_methods/test_scikit_finite_element.py @@ -171,21 +171,22 @@ def test_gradient_squared(self): y = mesh["current collector"].coordinates[0, :] z = mesh["current collector"].coordinates[1, :] + gradient = pybamm.grad(var) + grad_disc = disc.process_symbol(gradient) + grad_disc_y, grad_disc_z = grad_disc.children + gradient_squared = pybamm.grad_squared(var) - grad_squared_disc = disc.process_symbol(gradient_squared) - evaluated_gradient_squared = grad_squared_disc.evaluate( - None, 3 * y**2 + 2 * z**2 - ) + gradient_squared_disc = disc.process_symbol(gradient_squared) + + inner_product_y = grad_disc_y.evaluate(None, 5 * y + 6 * z) + inner_product_z = grad_disc_z.evaluate(None, 5 * y + 6 * z) + inner_product = inner_product_y**2 + inner_product_z**2 + + grad_squared_result = gradient_squared_disc.evaluate(None, 5 * y + 6 * z) - expected_result = 36 * y**2 + 16 * z**2 - expected_result = expected_result[:, np.newaxis] + np.testing.assert_array_almost_equal(grad_squared_result, inner_product) - try: - np.testing.assert_array_almost_equal( - evaluated_gradient_squared, expected_result - ) - except AssertionError as e: - print("AssertionError:", e) + np.testing.assert_array_less(0, grad_squared_result) def test_manufactured_solution(self): mesh = get_unit_2p1D_mesh_for_testing(ypts=32, zpts=32, include_particles=False) From 46e011d0ab72e1a3645232c87ae5446760df02e3 Mon Sep 17 00:00:00 2001 From: medha-14 Date: Sun, 10 Nov 2024 11:50:00 +0530 Subject: [PATCH 4/5] "Modified gradient_squared to apply boundary conditions during processing" --- src/pybamm/spatial_methods/finite_volume.py | 18 +++++++++++++++- .../test_scikit_finite_element.py | 21 +++++++++++++++++-- 2 files changed, 36 insertions(+), 3 deletions(-) diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py index 96ee51361d..6de9f873d6 100644 --- a/src/pybamm/spatial_methods/finite_volume.py +++ b/src/pybamm/spatial_methods/finite_volume.py @@ -147,14 +147,30 @@ def gradient_squared(self, symbol, discretised_symbol, boundary_conditions): The gradient squared of the symbol. """ domain = symbol.domain + + if symbol in boundary_conditions: + bcs = boundary_conditions[symbol] + if any(bc[1] == "Dirichlet" for bc in bcs.values()): + discretised_symbol, domain = self.add_ghost_nodes( + symbol, discretised_symbol, bcs + ) + gradient_matrix = self.gradient_matrix(domain, symbol.domains) - # Compute gradient squared: (∇u)^2 = u^T (L^T L) u + # Compute gradient squared matrix: (∇u)^2 = u^T (L^T L) u gradient_squared_matrix = gradient_matrix.T @ gradient_matrix gradient_squared_result = ( discretised_symbol.T @ gradient_squared_matrix @ discretised_symbol ) + # Add Neumann boundary conditions if defined + if symbol in boundary_conditions: + bcs = boundary_conditions[symbol] + if any(bc[1] == "Neumann" for bc in bcs.values()): + gradient_squared_result = self.add_neumann_values( + symbol, gradient_squared_result, bcs, domain + ) + return gradient_squared_result.item() def divergence(self, symbol, discretised_symbol, boundary_conditions): diff --git a/tests/unit/test_spatial_methods/test_scikit_finite_element.py b/tests/unit/test_spatial_methods/test_scikit_finite_element.py index 3ba13c6790..2244be4f57 100644 --- a/tests/unit/test_spatial_methods/test_scikit_finite_element.py +++ b/tests/unit/test_spatial_methods/test_scikit_finite_element.py @@ -183,11 +183,28 @@ def test_gradient_squared(self): inner_product = inner_product_y**2 + inner_product_z**2 grad_squared_result = gradient_squared_disc.evaluate(None, 5 * y + 6 * z) - np.testing.assert_array_almost_equal(grad_squared_result, inner_product) - np.testing.assert_array_less(0, grad_squared_result) + gradient_squared_disc_dirichlet = disc.process_symbol(gradient_squared) + grad_squared_result_dirichlet = gradient_squared_disc_dirichlet.evaluate( + None, 5 * y + 6 * z + ) + + np.testing.assert_array_almost_equal( + grad_squared_result_dirichlet, inner_product + ) + np.testing.assert_array_less(0, grad_squared_result_dirichlet) + + gradient_squared_disc_neumann = disc.process_symbol(gradient_squared) + grad_squared_result_neumann = gradient_squared_disc_neumann.evaluate( + None, 5 * y + 6 * z + ) + + np.testing.assert_array_almost_equal(grad_squared_result_neumann, inner_product) + + np.testing.assert_array_less(0, grad_squared_result_neumann) + def test_manufactured_solution(self): mesh = get_unit_2p1D_mesh_for_testing(ypts=32, zpts=32, include_particles=False) spatial_methods = { From 583422d75b293da1d968034eb345804b555d78cf Mon Sep 17 00:00:00 2001 From: medha-14 Date: Mon, 11 Nov 2024 16:31:32 +0530 Subject: [PATCH 5/5] made suggested changes --- src/pybamm/spatial_methods/finite_volume.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py index 6de9f873d6..c1a4e5bbfe 100644 --- a/src/pybamm/spatial_methods/finite_volume.py +++ b/src/pybamm/spatial_methods/finite_volume.py @@ -135,9 +135,9 @@ def gradient_squared(self, symbol, discretised_symbol, boundary_conditions): Parameters ---------- symbol : :class:`pybamm.Symbol` - The symbol for which to compute the gradient squared. + The symbol to compute the square of the gradient for. discretised_symbol : :class:`pybamm.Vector` - The discretised variable for which to compute the gradient squared. + The discretised variable to compute the square of the gradient for. boundary_conditions : dict Boundary conditions for the symbol. @@ -154,6 +154,9 @@ def gradient_squared(self, symbol, discretised_symbol, boundary_conditions): discretised_symbol, domain = self.add_ghost_nodes( symbol, discretised_symbol, bcs ) + elif any(bc[1] == "Neumann" for bc in bcs.values()): + # Neumann boundary conditions will be applied after computing the gradient squared + pass gradient_matrix = self.gradient_matrix(domain, symbol.domains) @@ -164,12 +167,12 @@ def gradient_squared(self, symbol, discretised_symbol, boundary_conditions): ) # Add Neumann boundary conditions if defined - if symbol in boundary_conditions: - bcs = boundary_conditions[symbol] - if any(bc[1] == "Neumann" for bc in bcs.values()): - gradient_squared_result = self.add_neumann_values( - symbol, gradient_squared_result, bcs, domain - ) + if symbol in boundary_conditions and any( + bc[1] == "Neumann" for bc in bcs.values() + ): + gradient_squared_result = self.add_neumann_values( + symbol, gradient_squared_result, bcs, domain + ) return gradient_squared_result.item()