Merge pull request #706 from pybamm-team/issue-705-fem-bc-bug

#705 fix sign error

CHANGELOG.md

@@ -21,6 +21,7 @@
 
 ## Bug fixes
 
+-   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))

pybamm/spatial_methods/scikit_finite_element.py

@@ -51,11 +51,11 @@ 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]
             )
         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]
             )
         else:
             raise pybamm.GeometryError(
@@ -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(

tests/unit/test_processed_variable.py

@@ -140,39 +140,35 @@ 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)
 
         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]

tests/unit/test_spatial_methods/test_scikit_finite_element.py

@@ -461,6 +461,68 @@ 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)
+
+        # 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")