[Bug]: Upwind/downwind schemes are not working properly

[brosaplanella]

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.

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

[brosaplanella]

The current approach creates an "upwind matrix" that converts the standard divergence matrix to the upwind divergence. However, the standard divergence matrix is singular (or at least for the problem I tested above) so I don't think such transformation exists.

I think a better approach would be to allow a keyword argument when calling pybamm.div that allows specifying upwind/downwind (or other methods in the future) so different matrices are assembled (similar on how different matrices are assembled for different coordinate systems).

[brosaplanella]

I have been playing with it and I can think of 2 possible ways of implementing this, and it boils down to how we conceptualise the model. The entry point for both would be to allow an extra keyword method.

1. Define UpwindDivergence and DownwindDivergence operator
   The divergence function would contain an if statement that defines a different symbol depending on the value of method. Then the spatial method should have separate methods for upwind and downwind. This would enable defining these methods as NotImplemented in the base class so it throws an error if the operator is used with a spatial method that doesn't support it.

2. Keep Divergence operator and have the if statement in the spatial method
   It has the advantage of being more generalisable and applicable to other operators, but not sure if there are any other cases we might want to implement that would benefit from this. I am thinking, for example, of a TVD scheme, but not sure if that's something we want to implement to a given operator or is a whole new spatial method. It would also require a lot of checks in the spatial methods to ensure if a wrong method is passed we throw an error.

I am not sure which way to go. I am tempted to say that, unless we can think of a good near-future feature that would benefit from approach 2, let's use 1 for now as it is cleaner and less prone to errors.

What do you think @tinosulzer @rtimms @martinjrobins?

[valentinsulzer]

It's been too long since I thought hard about spatial methods so just do whatever you think is simplest and safest. Does adding method require a user input or is there a sensible default?