From 3e11efe5c587d17c7f59986d48be5d09fa763b05 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Tue, 5 Nov 2019 17:13:45 +0000 Subject: [PATCH 1/4] #705 fix sign error --- CHANGELOG.md | 1 + .../spatial_methods/scikit_finite_element.py | 4 +- .../test_scikit_finite_element.py | 37 +++++++++++++++++++ 3 files changed, 40 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index dde65eed1e..fbd3c1b491 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ ## Bug fixes +- Correct a sign error in Dirichlet boundary conditions in the Finite Element Methid ([]()) - Pass the correct dimensional temperature to open circuit potential ([#702](https://github.com/pybamm-team/PyBaMM/pull/702)) - Adds missing temperature dependence in electrolyte and interface submodels ([#698](https://github.com/pybamm-team/PyBaMM/pull/698)) - Fix differentiation of functions that have more than one argument ([#687](https://github.com/pybamm-team/PyBaMM/pull/687)) diff --git a/pybamm/spatial_methods/scikit_finite_element.py b/pybamm/spatial_methods/scikit_finite_element.py index ddf501049e..5d8331e489 100644 --- a/pybamm/spatial_methods/scikit_finite_element.py +++ b/pybamm/spatial_methods/scikit_finite_element.py @@ -121,7 +121,7 @@ def unit_bc_load_form(v, dv, w): # set Dirichlet value at facets corresponding to tab neg_bc_load = np.zeros(mesh.npts) neg_bc_load[mesh.negative_tab_dofs] = 1 - boundary_load = boundary_load - neg_bc_value * pybamm.Vector(neg_bc_load) + boundary_load = boundary_load + neg_bc_value * pybamm.Vector(neg_bc_load) else: raise ValueError( "boundary condition must be Dirichlet or Neumann, not '{}'".format( @@ -138,7 +138,7 @@ def unit_bc_load_form(v, dv, w): # set Dirichlet value at facets corresponding to tab pos_bc_load = np.zeros(mesh.npts) pos_bc_load[mesh.positive_tab_dofs] = 1 - boundary_load = boundary_load - pos_bc_value * pybamm.Vector(pos_bc_load) + boundary_load = boundary_load + pos_bc_value * pybamm.Vector(pos_bc_load) else: raise ValueError( "boundary condition must be Dirichlet or Neumann, not '{}'".format( 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 5a108e275d..bf9a93a535 100644 --- a/tests/unit/test_spatial_methods/test_scikit_finite_element.py +++ b/tests/unit/test_spatial_methods/test_scikit_finite_element.py @@ -461,6 +461,43 @@ def test_pure_neumann_poisson(self): u_exact = z ** 2 / 2 - 1 / 6 np.testing.assert_array_almost_equal(solution.y[:-1], u_exact, decimal=1) + def test_dirichlet_bcs(self): + # manufactured solution u = a*z^2 + b*z + c + model = pybamm.BaseModel() + a = 3 + b = 4 + c = 5 + u = pybamm.Variable("variable", domain="current collector") + model.algebraic = {u: -pybamm.laplacian(u) + pybamm.source(2 * a, u)} + # set boundary conditions ("negative tab" = bottom of unit square, + # "positive tab" = top of unit square, elsewhere normal derivative is zero) + model.boundary_conditions = { + u: { + "negative tab": (pybamm.Scalar(c), "Dirichlet"), + "positive tab": (pybamm.Scalar(a + b + c), "Dirichlet"), + } + } + # bad initial guess (on purpose) + model.initial_conditions = {u: pybamm.Scalar(1)} + model.variables = {"u": u} + # create discretisation + mesh = get_unit_2p1D_mesh_for_testing(ypts=8, zpts=32) + spatial_methods = { + "macroscale": pybamm.FiniteVolume, + "current collector": pybamm.ScikitFiniteElement, + } + disc = pybamm.Discretisation(mesh, spatial_methods) + disc.process_model(model) + + # solve model + solver = pybamm.AlgebraicSolver() + solution = solver.solve(model) + + # indepdent of y, so just check values for one y + z = mesh["current collector"][0].edges["z"][:, np.newaxis] + u_exact = a * z ** 2 + b * z + c + np.testing.assert_array_almost_equal(solution.y[0 : len(z)], u_exact) + if __name__ == "__main__": print("Add -v for more debug output") From e4cf50e32af6490a8da9d227beeca38c573d8960 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Tue, 5 Nov 2019 17:17:05 +0000 Subject: [PATCH 2/4] #705 changelog --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fbd3c1b491..ffd77754d4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,7 +21,7 @@ ## Bug fixes -- Correct a sign error in Dirichlet boundary conditions in the Finite Element Methid ([]()) +- Correct a sign error in Dirichlet boundary conditions in the Finite Element Method ([#706](https://github.com/pybamm-team/PyBaMM/pull/706)) - Pass the correct dimensional temperature to open circuit potential ([#702](https://github.com/pybamm-team/PyBaMM/pull/702)) - Adds missing temperature dependence in electrolyte and interface submodels ([#698](https://github.com/pybamm-team/PyBaMM/pull/698)) - Fix differentiation of functions that have more than one argument ([#687](https://github.com/pybamm-team/PyBaMM/pull/687)) From 3a2c03700a31fc8ccdf6a002a52865f3e328cea6 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Wed, 6 Nov 2019 10:45:02 +0000 Subject: [PATCH 3/4] #705 fix disc of spatial var --- .../spatial_methods/scikit_finite_element.py | 6 ++-- tests/unit/test_processed_variable.py | 28 +++++++------------ .../test_scikit_finite_element.py | 27 +++++++++++++++++- 3 files changed, 40 insertions(+), 21 deletions(-) diff --git a/pybamm/spatial_methods/scikit_finite_element.py b/pybamm/spatial_methods/scikit_finite_element.py index 5d8331e489..f6e0832cec 100644 --- a/pybamm/spatial_methods/scikit_finite_element.py +++ b/pybamm/spatial_methods/scikit_finite_element.py @@ -51,11 +51,13 @@ def spatial_variable(self, symbol): symbol_mesh = self.mesh if symbol.name == "y": vector = pybamm.Vector( - symbol_mesh["current collector"][0].edges["y"], domain=symbol.domain + symbol_mesh["current collector"][0].coordinates[0, :][:, np.newaxis] + # symbol_mesh["current collector"][0].edges["y"], domain=symbol.domain ) elif symbol.name == "z": vector = pybamm.Vector( - symbol_mesh["current collector"][0].edges["z"], domain=symbol.domain + symbol_mesh["current collector"][0].coordinates[1, :][:, np.newaxis] + # symbol_mesh["current collector"][0].edges["z"], domain=symbol.domain ) else: raise pybamm.GeometryError( diff --git a/tests/unit/test_processed_variable.py b/tests/unit/test_processed_variable.py index a891fdbc6b..87c29a84a0 100644 --- a/tests/unit/test_processed_variable.py +++ b/tests/unit/test_processed_variable.py @@ -140,13 +140,11 @@ def test_processed_variable_3D_x_z(self): def test_processed_variable_3D_scikit(self): var = pybamm.Variable("var", domain=["current collector"]) - y = pybamm.SpatialVariable("y", domain=["current collector"]) - z = pybamm.SpatialVariable("z", domain=["current collector"]) disc = tests.get_2p1d_discretisation_for_testing() disc.set_variable_slices([var]) - y_sol = disc.process_symbol(y).entries[:, 0] - z_sol = disc.process_symbol(z).entries[:, 0] + y = disc.mesh["current collector"][0].edges["y"] + z = disc.mesh["current collector"][0].edges["z"] var_sol = disc.process_symbol(var) t_sol = np.linspace(0, 1) u_sol = np.ones(var_sol.shape[0])[:, np.newaxis] * np.linspace(0, 5) @@ -154,25 +152,23 @@ def test_processed_variable_3D_scikit(self): processed_var = pybamm.ProcessedVariable(var_sol, t_sol, u_sol, mesh=disc.mesh) np.testing.assert_array_equal( processed_var.entries, - np.reshape(u_sol, [len(y_sol), len(z_sol), len(t_sol)]), + np.reshape(u_sol, [len(y), len(z), len(t_sol)]), ) def test_processed_variable_2Dspace_scikit(self): var = pybamm.Variable("var", domain=["current collector"]) - y = pybamm.SpatialVariable("y", domain=["current collector"]) - z = pybamm.SpatialVariable("z", domain=["current collector"]) disc = tests.get_2p1d_discretisation_for_testing() disc.set_variable_slices([var]) - y_sol = disc.process_symbol(y).entries[:, 0] - z_sol = disc.process_symbol(z).entries[:, 0] + y = disc.mesh["current collector"][0].edges["y"] + z = disc.mesh["current collector"][0].edges["z"] var_sol = disc.process_symbol(var) t_sol = np.array([0]) u_sol = np.ones(var_sol.shape[0])[:, np.newaxis] processed_var = pybamm.ProcessedVariable(var_sol, t_sol, u_sol, mesh=disc.mesh) np.testing.assert_array_equal( - processed_var.entries, np.reshape(u_sol, [len(y_sol), len(z_sol)]) + processed_var.entries, np.reshape(u_sol, [len(y), len(z)]) ) def test_processed_var_1D_interpolation(self): @@ -367,13 +363,11 @@ def test_processed_var_3D_r_first_dimension(self): def test_processed_var_3D_scikit_interpolation(self): var = pybamm.Variable("var", domain=["current collector"]) - y = pybamm.SpatialVariable("y", domain=["current collector"]) - z = pybamm.SpatialVariable("z", domain=["current collector"]) disc = tests.get_2p1d_discretisation_for_testing() disc.set_variable_slices([var]) - y_sol = disc.process_symbol(y).entries[:, 0] - z_sol = disc.process_symbol(z).entries[:, 0] + y_sol = disc.mesh["current collector"][0].edges["y"] + z_sol = disc.mesh["current collector"][0].edges["z"] var_sol = disc.process_symbol(var) t_sol = np.linspace(0, 1) u_sol = np.ones(var_sol.shape[0])[:, np.newaxis] * np.linspace(0, 5) @@ -406,13 +400,11 @@ def test_processed_var_3D_scikit_interpolation(self): def test_processed_var_2Dspace_scikit_interpolation(self): var = pybamm.Variable("var", domain=["current collector"]) - y = pybamm.SpatialVariable("y", domain=["current collector"]) - z = pybamm.SpatialVariable("z", domain=["current collector"]) disc = tests.get_2p1d_discretisation_for_testing() disc.set_variable_slices([var]) - y_sol = disc.process_symbol(y).entries[:, 0] - z_sol = disc.process_symbol(z).entries[:, 0] + y_sol = disc.mesh["current collector"][0].edges["y"] + z_sol = disc.mesh["current collector"][0].edges["z"] var_sol = disc.process_symbol(var) t_sol = np.array([0]) u_sol = np.ones(var_sol.shape[0])[:, np.newaxis] 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 bf9a93a535..21bd4b8799 100644 --- a/tests/unit/test_spatial_methods/test_scikit_finite_element.py +++ b/tests/unit/test_spatial_methods/test_scikit_finite_element.py @@ -493,11 +493,36 @@ def test_dirichlet_bcs(self): solver = pybamm.AlgebraicSolver() solution = solver.solve(model) - # indepdent of y, so just check values for one y + # indepedent of y, so just check values for one y z = mesh["current collector"][0].edges["z"][:, np.newaxis] u_exact = a * z ** 2 + b * z + c np.testing.assert_array_almost_equal(solution.y[0 : len(z)], u_exact) + def test_disc_spatial_var(self): + mesh = get_unit_2p1D_mesh_for_testing(ypts=4, zpts=5) + spatial_methods = { + "macroscale": pybamm.FiniteVolume, + "current collector": pybamm.ScikitFiniteElement, + } + disc = pybamm.Discretisation(mesh, spatial_methods) + + # discretise y and z + y = pybamm.SpatialVariable("y", ["current collector"]) + z = pybamm.SpatialVariable("z", ["current collector"]) + y_disc = disc.process_symbol(y) + z_disc = disc.process_symbol(z) + + # create expected meshgrid + y_vec = np.linspace(0, 1, 4) + z_vec = np.linspace(0, 1, 5) + Y, Z = np.meshgrid(y_vec, z_vec) + y_actual = np.transpose(Y).flatten()[:, np.newaxis] + z_actual = np.transpose(Z).flatten()[:, np.newaxis] + + # spatial vars should discretise to the flattend meshgrid + np.testing.assert_array_equal(y_disc.evaluate(), y_actual) + np.testing.assert_array_equal(z_disc.evaluate(), z_actual) + if __name__ == "__main__": print("Add -v for more debug output") From cb19c95f4ac1626a35bad70496b5eabf6b4290ee Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Wed, 6 Nov 2019 12:52:29 +0000 Subject: [PATCH 4/4] #705 remove old spatial var disc --- pybamm/spatial_methods/scikit_finite_element.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pybamm/spatial_methods/scikit_finite_element.py b/pybamm/spatial_methods/scikit_finite_element.py index f6e0832cec..464a1d68a7 100644 --- a/pybamm/spatial_methods/scikit_finite_element.py +++ b/pybamm/spatial_methods/scikit_finite_element.py @@ -52,12 +52,10 @@ def spatial_variable(self, symbol): if symbol.name == "y": vector = pybamm.Vector( symbol_mesh["current collector"][0].coordinates[0, :][:, np.newaxis] - # symbol_mesh["current collector"][0].edges["y"], domain=symbol.domain ) elif symbol.name == "z": vector = pybamm.Vector( symbol_mesh["current collector"][0].coordinates[1, :][:, np.newaxis] - # symbol_mesh["current collector"][0].edges["z"], domain=symbol.domain ) else: raise pybamm.GeometryError(