[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.

[mwverbrugge]

Thanks Eric/PyBaMM team. Can you let me know when the solution will be pushed to a released PyBaMM version? I just installed 24.11.1 to see if the fix might show up in this recent PB release, but the same issue prevails. Mark

…

On Wed, Nov 27, 2024 at 12:36 PM Eric G. Kratz ***@***.***> wrote: Closed #4433 <#4433> as completed via #4614 <#4614>. — Reply to this email directly, view it on GitHub <#4433 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BLHDY3UDVTY6KKCTA7IK3ZD2CYUM3AVCNFSM6AAAAABOB3W64GVHI2DSMVQWIX3LMV45UABCJFZXG5LFIV3GK3TUJZXXI2LGNFRWC5DJN5XDWMJVGQ2TONBTGQ2DIMI> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[agriyakhetarpal]

Hi @mwverbrugge, the fix hasn't been released yet: #4614. It should be available with PyBaMM version 25.1.

[mwverbrugge]

Got it...thanks

…

On Fri, Nov 29, 2024, 4:47 PM Agriya Khetarpal ***@***.***> wrote: Hi @mwverbrugge <https://github.com/mwverbrugge>, the fix hasn't been released yet: #4614 <#4614>. It should be available with PyBaMM version 25.1. — Reply to this email directly, view it on GitHub <#4433 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BLHDY3QWIKJEJ55GOHVXYBD2DEDLRAVCNFSM6AAAAABOB3W64GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKMBYG43DANBUHA> . You are receiving this because you were mentioned.Message ID: ***@***.***>