[Bug]: BoundaryGradient is only first-order accurate

[mwverbrugge]

PyBaMM Version

24.1 and 24.9

Python Version

3.12.4

Describe the bug

The following program illustrates a bug in the routine pybamm.BoundaryGradient. We have tested this routine on the most recent pybamm version 24.9, but the problem exists in earlier versions as well.

In the program, a differential equation for c(y) is solved, such that the exact solution is
c(y)=y**2/2+y-3/2

The function
grad(c)=y+1
BoundaryGradient is set to give the gradient of the variable c(y) at y=0, which should be equal to 1, but BoundaryGradient is instead calculating the gradient of c(y) between the first two nodes for c(y), which, in this program, is at y=0.2.

As a result, one sees that
pybamm.BoundaryGradient(c,"left")=1.2. In a similar manner, BoundaryGradient is set to give the gradient of the variable c(y) at y=1, which should be equal to 2, but BoundaryGradient is instead calculating the gradient of c(y) between the last two nodes for c(y), which, in this program, is at y=0.8. As a result, one sees that
pybamm.BoundaryGradient(c,"right")=1.8.

The program makes plots of the solution (attached Figure_1) and prints out the values of pybamm.BoundaryGradient(c,"left") and pybamm.BoundaryGradient(c,"right") to verify the above assertions.

Steps to Reproduce

import matplotlib.pyplot as plt
import numpy as np
import pybamm as pb

start use of PyBaMM

model = pb.BaseModel()

domains and associated spatial coordinates

x_p = pb.SpatialVariable("y", domain="positive electrode", coord_sys="cartesian")

spatial coordinate ranges

geometry = { "positive electrode": {x_p: {"min": 0, "max": 1 }}}

elements for each domain

var_pts = {x_p: 5}

t_eval = np.linspace(0,1,3)

c = pb.Variable("c", domain="positive electrode")
f = pb.Variable("f")

dcdy = pb.grad(c)
dcdy_BGl = pb.BoundaryGradient(c, "left")
dcdy_BGr = pb.BoundaryGradient(c, "right")

"""
The function f below is not relevant to the BoundaryGradient problem. It is only
there because the CasadiSolver used below requires a transient component to any
problem it solves.
"""
model.rhs[f]=f-1
model.algebraic[c]=pb.div(pb.grad(c))-1

model.boundary_conditions = {
c: {"left": (1, "Neumann"),
"right": (0, "Dirichlet")}, }

model.initial_conditions = {c: 1, f: 0}

model.variables = {"c":c, "dcdy":dcdy, "dcdy_BGl":dcdy_BGl,"dcdy_BGr":dcdy_BGr}

submesh_types = {
"positive electrode": pb.Uniform1DSubMesh, }

mesh = pb.Mesh(geometry, submesh_types, var_pts)

spatial_methods = {
"positive electrode": pb.FiniteVolume(), }

disc = pb.Discretisation(mesh, spatial_methods)
disc.process_model(model)

solver = pb.CasadiSolver(mode="safe", atol=1e-6, rtol=1e-6,)

solver = pb.CasadiSolver()
solution = solver.solve(model, t_eval)

c_out = solution["c"]
dcdy_out = solution["dcdy"]
dcdy_BGl_out = solution["dcdy_BGl"]
dcdy_BGr_out = solution["dcdy_BGr"]
#the solution for c is independent of time, so below we only give it at time t=0.5

plt.subplot(1, 2, 1)
y_n=np.array([0,.1,.3,.5,.7,.9,1])
c_exact=y_n**2/2+y_n-3/2
plt.figure(1)
plt.plot(y_n[1:-1],c_out(t=0.5,y=y_n[1:-1]),label="c numerical",
marker="o",linestyle='None',)
plt.plot(y_n,c_exact,label="exact solution")
plt.xlabel("y")
plt.grid(axis="both")
plt.legend()

plt.subplot(1, 2, 2)
y_g=np.array([0,.2,.4,.6,.8,1])
plt.plot(y_g[1:-1],dcdy_out(t=0.5,y=y_g[1:-1]), label="dc/dy numerical", marker="o")
dcdy_exact=y_g+1
plt.plot(y_g,dcdy_exact,label="exact dcdy")
plt.plot(0,dcdy_BGl_out(0),marker='s',linestyle='None',
label='BoundaryGradient(c, "left")')
plt.plot(1,dcdy_BGr_out(0),marker='^',linestyle='None',
label='BoundaryGradient(c, "right")')
plt.xlabel("y")
plt.grid(axis="both")
plt.legend()

print("nodes :",y_n[1:-1])
print("numerical c at nodes:",c_out(t=0.5,y=y_n[1:-1]))
print("pb.BoundaryGradient(c, ""left"")=",dcdy_BGl_out(t=0.5))
print("numerical (c(1)-c(0))/dy=(-1.16+1.4)/.2=",(-1.16+1.4)/.2)
print("numerical grad(c) at y=0.2=",dcdy_out(t=0.5,y=0.2))
print("analytic dcdy(y=0)=",y_g[0]+1)

print("pb.BoundaryGradient(c, ""right"")=",dcdy_BGr_out(t=0.5))
print("numerical (c(4)-c(3))/dy=(-0.2+0.56)/.2=",(-0.2+0.56)/.2)
print("numerical grad(c) at y=0.8=",dcdy_out(t=0.5,y=0.8))
print("analytic dcdy(y=1)=",y_g[-1]+1)

Relevant log output

No response

[rtimms]

It looks like the current implementation of BoundaryValue is only first-order accurate, see https://github.com/pybamm-team/PyBaMM/blob/develop/src/pybamm/spatial_methods/finite_volume.py#L954.

The finite volume code could do with some cleaning up as it is pretty difficult to read and should be updated to ensure second-order accuracy everywhere.

There's an argument to be made that we should just rely on a dependency instead of having our own method, but since we have it already it's probably simplest just to update our won.