Problems with the gradient operator

[isaacbasil]

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 :) 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])

[github-actions[bot]]

Hi! We have now moved our discussions to Discourse. Please post your question there.

[brosaplanella]

Hi! Do you mind posting the issue in Discourse and I will reply there? This is an issue that has come up a few times, so it would be good to have it archived there for future reference (please tag me here when you have opened the issue so I don't forget).

[isaacbasil]

Hi! No problem, I've posted it on Discourse. Thanks.

[brosaplanella]

I can't find it, can you post the link here?

[isaacbasil]

here's the link: https://pybamm.discourse.group/t/problems-with-the-gradient-operator/115