rtimms
commented
2 months ago 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.
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