pybamm-team
[Bug]: `BoundaryGradient` is only first-order accurate
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
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.
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 https://github.com/pybamm-team/PyBaMM/issues/4433 as completed via #4614 https://github.com/pybamm-team/PyBaMM/pull/4614.
— Reply to this email directly, view it on GitHub https://github.com/pybamm-team/PyBaMM/issues/4433#event-15457434441, or unsubscribe https://github.com/notifications/unsubscribe-auth/BLHDY3UDVTY6KKCTA7IK3ZD2CYUM3AVCNFSM6AAAAABOB3W64GVHI2DSMVQWIX3LMV45UABCJFZXG5LFIV3GK3TUJZXXI2LGNFRWC5DJN5XDWMJVGQ2TONBTGQ2DIMI . You are receiving this because you authored the thread.Message ID: @.***>
Hi @mwverbrugge, the fix hasn't been released yet: https://github.com/pybamm-team/PyBaMM/pull/4614. It should be available with PyBaMM version 25.1.
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 https://github.com/pybamm-team/PyBaMM/pull/4614. It should be available with PyBaMM version 25.1.
— Reply to this email directly, view it on GitHub https://github.com/pybamm-team/PyBaMM/issues/4433#issuecomment-2508760448, or unsubscribe https://github.com/notifications/unsubscribe-auth/BLHDY3QWIKJEJ55GOHVXYBD2DEDLRAVCNFSM6AAAAABOB3W64GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKMBYG43DANBUHA . You are receiving this because you were mentioned.Message ID: @.***>