Problems with the gradient operator

PyBaMM Q&A
1

Hi, I am trying to build a custom model with a bit of an unusual form, but I don’t see why it shouldn’t be solvable. However, I am having some problem with the grad operator, and getting an “off by one” shape error, which I believe is caused by a problem with the boundary conditions not being implemented properly.

Essentially, I am trying to add another grad term to my PDE for c_e, which will be made much clearer in my minimal working example below. In the example, the key line is the following:

a_eff = a + scale * pybamm.grad(c_e)

When I do this, I get the following error:

pybamm.expression_tree.exceptions.ModelError: rhs and initial conditions must have the same shape after discretisation but rhs.shape = (41, 1) and initial_conditions.shape = (40, 1) for variable ‘c_e’.

However, if I instead write,

a_eff = a

then everything works fine. I think this is because grad returns an object on the edges, but I need an object defined on the nodes - perhaps there’s a way I can map between these two types of objects?

If anyone has any ideas on how to tackle this that’d be greatly appreciated :slight_smile: Here is the example (the physics in this case is meaningless, I just want to get it to run without errors):

import pybamm

model = pybamm.BaseModel()

x_p = pybamm.SpatialVariable("x", domain="positive electrode", coord_sys="cartesian")
r_p = pybamm.SpatialVariable("r", domain="positive particle", coord_sys="cylindrical polar")

L_p = 40e-6
R_p = 10e-6

geometry = {"positive electrode": {x_p: {"min": 0, "max": L_p }},
    "positive particle": {r_p: {"min": 0, "max": R_p }},}

var_pts = {x_p: 40, r_p: 10}

submesh_types = {"positive electrode": pybamm.Uniform1DSubMesh, "positive particle": pybamm.Uniform1DSubMesh}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
spatial_methods = {"positive electrode": pybamm.FiniteVolume(),"positive particle": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)

c_e = pybamm.Variable("c_e", domain="positive electrode")
c_s = pybamm.Variable("c_s", domain="positive particle", auxiliary_domains={"secondary": "positive electrode"})

D_e = pybamm.Scalar(1)
D_s = pybamm.Scalar(1)
a = pybamm.PrimaryBroadcast(10000, "positive electrode")
scale = pybamm.Scalar(1e-5)

j = pybamm.PrimaryBroadcast(1, "positive electrode")

a_eff = a + scale * pybamm.grad(c_e)

dc_e_dt = pybamm.div(D_e * pybamm.grad(c_e)) + a_eff * j
dc_s_dt = pybamm.div(D_s * pybamm.grad(c_s))
model.rhs[c_e] = dc_e_dt
model.rhs[c_s] = dc_s_dt

model.boundary_conditions = {
    c_e: {"left": (0, "Neumann"), "right": (1, "Neumann")},
    c_s: {"left": (0, "Neumann"), "right": (- j, "Neumann")},
    }
model.initial_conditions = {c_e: 1000, c_s: 20000}

model.variables = {"c_e": c_e, "c_s": c_s}

disc.process_model(model)
sim = pybamm.Simulation(model)
sim.solve([0, 1])
2

The issue is that the gradient takes a variable that’s defined on the nodes and returns something that’s defined on the edges, that’s why you get the dimension mismatch. I am not sure I fully understand the equation you are implementing, but a way to get it working is to use the divergence instead of the gradient and then broadcast your variable to the edges.

I tried and I can’t broadcast c_e directly (don’t know why) but this seems to work:

u = pybamm.PrimaryBroadcastToEdges(pybamm.Scalar(1), "positive electrode")
a_eff = a + scale * pybamm.div(u * c_e)
1 Like
3

thanks very much! that works great for me :slight_smile: