diff --git a/CHANGELOG.md b/CHANGELOG.md index 4e3e90f035..e1bb4214af 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,7 @@ ## Bug Fixes +- Fix bug with upwind and downwind schemes producing the wrong discretised system ([#3979](https://github.com/pybamm-team/PyBaMM/pull/3979)) - Allow evaluation of an `Interpolant` object with a number ([#3932](https://github.com/pybamm-team/PyBaMM/pull/3932)) - `plot_voltage_components` now works even if the time does not start at 0 ([#3915](https://github.com/pybamm-team/PyBaMM/pull/3915)) - Fixed bug where separator porosity was used in calculation instead of transport efficiency ([#3905](https://github.com/pybamm-team/PyBaMM/pull/3905)) diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py index 2a2bb36640..f991599735 100644 --- a/pybamm/expression_tree/unary_operators.py +++ b/pybamm/expression_tree/unary_operators.py @@ -24,6 +24,8 @@ class UnaryOperator(pybamm.Symbol): name of the node child : :class:`Symbol` child node + domains : dict + A dictionary equivalent to {'primary': domain, auxiliary_domains}. """ def __init__( @@ -398,6 +400,8 @@ class with a :class:`Matrix` name of the node child : :class:`Symbol` child node + domains : dict + A dictionary equivalent to {'primary': domain, auxiliary_domains}. """ def __init__( diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py index 11313a1450..d02a2b0734 100644 --- a/pybamm/spatial_methods/finite_volume.py +++ b/pybamm/spatial_methods/finite_volume.py @@ -613,8 +613,13 @@ def add_ghost_nodes(self, symbol, discretised_symbol, bcs): n = submesh.npts second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains) - lbc_value, lbc_type = bcs["left"] - rbc_value, rbc_type = bcs["right"] + # Catch if no boundary conditions are defined + if "left" not in bcs.keys() and "right" not in bcs.keys(): + raise ValueError(f"No boundary conditions have been provided for {symbol}") + + # Allow to only pass one boundary condition (for upwind/downwind) + lbc_value, lbc_type = bcs.get("left", (None, None)) + rbc_value, rbc_type = bcs.get("right", (None, None)) # Add ghost node(s) to domain where necessary and count number of # Dirichlet boundary conditions @@ -637,7 +642,7 @@ def add_ghost_nodes(self, symbol, discretised_symbol, bcs): else: left_ghost_constant = 2 * lbc_value lbc_vector = pybamm.Matrix(lbc_matrix) @ left_ghost_constant - elif lbc_type == "Neumann": + elif lbc_type in ["Neumann", None]: lbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats)) else: raise ValueError( @@ -656,7 +661,7 @@ def add_ghost_nodes(self, symbol, discretised_symbol, bcs): else: right_ghost_constant = 2 * rbc_value rbc_vector = pybamm.Matrix(rbc_matrix) @ right_ghost_constant - elif rbc_type == "Neumann": + elif rbc_type in ["Neumann", None]: rbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats)) else: raise ValueError( @@ -1390,8 +1395,6 @@ def upwind_or_downwind(self, symbol, discretised_symbol, bcs, direction): direction : str Direction in which to apply the operator (upwind or downwind) """ - submesh = self.mesh[symbol.domain] - n = submesh.npts if symbol not in bcs: raise pybamm.ModelError( @@ -1399,36 +1402,17 @@ def upwind_or_downwind(self, symbol, discretised_symbol, bcs, direction): ) if direction == "upwind": - bc, typ = bcs[symbol]["left"] - if typ != "Dirichlet": - raise pybamm.ModelError( - "Dirichlet boundary conditions must be provided for " - f"upwinding '{symbol}'" - ) - - concat_bc = pybamm.NumpyConcatenation(bc, discretised_symbol) - - upwind_mat = vstack( - [ - csr_matrix(([1], ([0], [0])), shape=(1, n + 1)), - diags([-0.5, 1.5], [0, 1], shape=(n, n + 1)), - ] - ) - symbol_out = pybamm.Matrix(upwind_mat) @ concat_bc + bc_side = "left" elif direction == "downwind": - bc, typ = bcs[symbol]["right"] - if typ != "Dirichlet": - raise pybamm.ModelError( - "Dirichlet boundary conditions must be provided for " - f"downwinding '{symbol}'" - ) + bc_side = "right" - concat_bc = pybamm.NumpyConcatenation(discretised_symbol, bc) - downwind_mat = vstack( - [ - diags([1.5, -0.5], [0, 1], shape=(n, n + 1)), - csr_matrix(([1], ([0], [n])), shape=(1, n + 1)), - ] + if bcs[symbol][bc_side][1] != "Dirichlet": + raise pybamm.ModelError( + "Dirichlet boundary conditions must be provided for " + f"{direction}ing '{symbol}'" ) - symbol_out = pybamm.Matrix(downwind_mat) @ concat_bc + + # Extract only the relevant boundary condition as the model might have both + bc_subset = {bc_side: bcs[symbol][bc_side]} + symbol_out, _ = self.add_ghost_nodes(symbol, discretised_symbol, bc_subset) return symbol_out diff --git a/tests/integration/test_spatial_methods/test_finite_volume.py b/tests/integration/test_spatial_methods/test_finite_volume.py index 2f9f6ef1b0..ce2d07c2de 100644 --- a/tests/integration/test_spatial_methods/test_finite_volume.py +++ b/tests/integration/test_spatial_methods/test_finite_volume.py @@ -339,6 +339,55 @@ def test_laplacian_spherical(self): ) +def solve_advection_equation(direction="upwind", source=1, bc=0): + model = pybamm.BaseModel() + x = pybamm.SpatialVariable("x", domain="domain", coord_sys="cartesian") + u = pybamm.Variable("u", domain="domain") + if direction == "upwind": + bc_side = "left" + y = x + v = pybamm.PrimaryBroadcastToEdges(1, ["domain"]) + rhs = -pybamm.div(pybamm.upwind(u) * v) + source + elif direction == "downwind": + bc_side = "right" + y = 1 - x + v = pybamm.PrimaryBroadcastToEdges(-1, ["domain"]) + rhs = -pybamm.div(pybamm.downwind(u) * v) + source + + u_an = (bc + source * y) - (bc + source * (y - pybamm.t)) * ((y - pybamm.t) > 0) + model.boundary_conditions = { + u: { + bc_side: (pybamm.Scalar(bc), "Dirichlet"), + } + } + model.rhs = {u: rhs} + model.initial_conditions = {u: pybamm.Scalar(0)} + model.variables = {"u": u, "x": x, "analytical": u_an} + geometry = {"domain": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}} + submesh_types = {"domain": pybamm.Uniform1DSubMesh} + var_pts = {x: 1000} + mesh = pybamm.Mesh(geometry, submesh_types, var_pts) + spatial_methods = {"domain": pybamm.FiniteVolume()} + disc = pybamm.Discretisation(mesh, spatial_methods) + disc.process_model(model) + solver = pybamm.CasadiSolver() + return solver.solve(model, [0, 1]) + + +class TestUpwindDownwind(TestCase): + def test_upwind(self): + solution = solve_advection_equation("upwind") + np.testing.assert_array_almost_equal( + solution["u"].entries, solution["analytical"].entries, decimal=2 + ) + + def test_downwind(self): + solution = solve_advection_equation("downwind") + np.testing.assert_array_almost_equal( + solution["u"].entries, solution["analytical"].entries, decimal=2 + ) + + if __name__ == "__main__": print("Add -v for more debug output") import sys diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py index 16a3bbde2c..52fa379ef6 100644 --- a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py +++ b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py @@ -410,11 +410,11 @@ def test_upwind_downwind(self): y_test = 2 * np.ones_like(nodes) np.testing.assert_array_equal( disc_upwind.evaluate(y=y_test), - np.concatenate([np.array([5, 0.5]), 2 * np.ones(n - 1)])[:, np.newaxis], + np.concatenate([np.array([8]), 2 * np.ones(n)])[:, np.newaxis], ) np.testing.assert_array_equal( disc_downwind.evaluate(y=y_test), - np.concatenate([2 * np.ones(n - 1), np.array([1.5, 3])])[:, np.newaxis], + np.concatenate([2 * np.ones(n), np.array([4])])[:, np.newaxis], ) # Remove boundary conditions and check error is raised diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_ghost_nodes_and_neumann.py b/tests/unit/test_spatial_methods/test_finite_volume/test_ghost_nodes_and_neumann.py index 5651a9dd1e..47e1a2fda1 100644 --- a/tests/unit/test_spatial_methods/test_finite_volume/test_ghost_nodes_and_neumann.py +++ b/tests/unit/test_spatial_methods/test_finite_volume/test_ghost_nodes_and_neumann.py @@ -52,6 +52,8 @@ def test_add_ghost_nodes(self): bcs = {"left": (pybamm.Scalar(0), "Neumann"), "right": (pybamm.Scalar(3), "x")} with self.assertRaisesRegex(ValueError, "boundary condition must be"): sp_meth.add_ghost_nodes(var, discretised_symbol, bcs) + with self.assertRaisesRegex(ValueError, "No boundary conditions"): + sp_meth.add_ghost_nodes(var, discretised_symbol, {}) with self.assertRaisesRegex(ValueError, "boundary condition must be"): sp_meth.add_neumann_values(var, discretised_symbol, bcs, var.domain)