Repository: pybamm-team/PyBaMM (python)

Issue:
[Bug]: Upwind/downwind schemes are not working properly
### PyBaMM Version

develop

### Python Version

all

### Describe the bug

Checking the problem solved in `finite-volumes.ipynb` it does not correspond with what one would expect from an upwind scheme.

See example below.

### Steps to Reproduce
The example below compares a toy problem with its analytical solution. Note that for coarse meshes there is some disagreement but it converges if you refine the mesh. Commenting the lines that fix the system retrieves the previous behaviour which disagrees with the analytical solution.

```python3
import pybamm


model = pybamm.BaseModel()
x = pybamm.SpatialVariable("x", domain="domain", coord_sys="cartesian")
u = pybamm.Variable("u", domain="domain")
u_an = x + (pybamm.t - x) * ((x - pybamm.t) > 0)
v = pybamm.PrimaryBroadcastToEdges(1, ["domain"])
rhs = - pybamm.div(pybamm.upwind(u) * v) + 1
model.boundary_conditions = {
    u: {
        "left": (pybamm.Scalar(0), "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}
N = 5
var_pts = {x: N}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
spatial_methods = {"domain": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)
model_disc = disc.process_model(model, inplace=False)

# Fix the entries of the discretised system (uncomment to get original behaviour)
for i in range(N - 1):
    model_disc.rhs[u].children[1].children[0].entries[i + 1, i] = 0
    model_disc.rhs[u].children[1].children[0].entries[i + 1, i + 1] = -N
    model_disc.rhs[u].children[1].children[0].entries[i + 1, i + 2] = N
model_disc.rhs[u].children[1].children[0].entries[0, 0] = - 2 * N
model_disc.rhs[u].children[1].children[0].entries[0, 1] = 2 * N


print(model_disc.rhs[u].children[1].children[0].entries.toarray())
solver = pybamm.CasadiSolver()
solution =solver.solve(model_disc, [0, 10])

plot = pybamm.QuickPlot(solution,[["u", "analytical"]])
plot.dynamic_plot()
```

### Relevant log output

_No response_

Pull request:
Fix upwind/downwind finite volume operators
# Description

Fixed the issue with the upwind and downwind operators. I am now leveraging the `add_ghost_nodes` method, to which I had to do some minor changes.

Fixes #3526 

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [ ] New feature (non-breaking change which adds functionality)
- [ ] Optimization (back-end change that speeds up the code)
- [x] Bug fix (non-breaking change which fixes an issue)

# Key checklist:

- [x] No style issues: `$ pre-commit run` (or `$ nox -s pre-commit`) (see [CONTRIBUTING.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CONTRIBUTING.md#installing-and-using-pre-commit) for how to set this up to run automatically when committing locally, in just two lines of code)
- [ ] All tests pass: `$ python run-tests.py --all` (or `$ nox -s tests`)
- [ ] The documentation builds: `$ python run-tests.py --doctest` (or `$ nox -s doctests`)

You can run integration tests, unit tests, and doctests together at once, using `$ python run-tests.py --quick` (or `$ nox -s quick`).

## Further checks:

- [ ] Code is commented, particularly in hard-to-understand areas
- [ ] Tests added that prove fix is effective or that feature works


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index b474549bc..22d89ca9b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,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 2a2bb3664..f99159973 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 SpatialOperator(UnaryOperator):
         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 11313a145..d02a2b073 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -613,8 +613,13 @@ class FiniteVolume(pybamm.SpatialMethod):
         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 @@ class FiniteVolume(pybamm.SpatialMethod):
             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 @@ class FiniteVolume(pybamm.SpatialMethod):
             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 @@ class FiniteVolume(pybamm.SpatialMethod):
         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 @@ class FiniteVolume(pybamm.SpatialMethod):
             )
 
         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


Test patch:
diff --git a/tests/integration/test_spatial_methods/test_finite_volume.py b/tests/integration/test_spatial_methods/test_finite_volume.py
index 2f9f6ef1b..ce2d07c2d 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 @@ class TestFiniteVolumeLaplacian(TestCase):
         )
 
 
+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 16a3bbde2..52fa379ef 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 @@ class TestFiniteVolume(TestCase):
         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 5651a9dd1..47e1a2fda 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 @@ class TestGhostNodes(TestCase):
         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)