Discretisation problems when using both algebraic and time-derivative variables

[CSchwet]

Hi,

I am trying to implement a model in pybamm, which relies on three base variables, two governed by time-derivative equations of the form dX_a/dt = div(N), where N = sum_a A_a grad(X_a) and the final variable is determined by a poisson equation div(grad(Y)) = X_a. I have attached a simplified version below with only two variables, but reflecting the shape I need.

However, I need line 11 to be of the form N = -pybamm.grad(c - d / k) or N = -pybamm.grad(c) - pybamm.grad(d / k). This produces an error "pybamm.expression_tree.exceptions.ShapeError: Cannot find shape (original error: operands could not be broadcast together with shapes (21,1) (19,1) )", throwing two previous errors "AttributeError: 'Subtraction' object has no attribute '_saved_evaluate_for_shape'. Did you mean: '_evaluate_for_shape'?" and "ValueError: operands could not be broadcast together with shapes (21,1) (19,1)".

I suspect the solver uses different discretisations for the variables c and d. Is there a way to force the discretisations to be equal, or another way to handle this?

Thank you!

code

import pybamm

model = pybamm.BaseModel()

c = pybamm.Variable("Concentration", domain="negative particle") # the two variables I am solving for
d = pybamm.Variable("Algebraic", domain="negative particle")

k = pybamm.Parameter("Constant k")

N = -pybamm.grad(c) - pybamm.grad(d) / k # combine the gradients to get the flux, this line leads to issues

dcdt = -pybamm.div(N) # define the rhs equation

poisson = pybamm.div(pybamm.grad(d)) + c # define the algebraic equation

model.rhs = {c: dcdt} # add the equation to rhs dictionary with the variable as the key
model.algebraic = {d: poisson}

#boundary conditions
model.boundary_conditions = {
c: {"left": (0, "Neumann"), "right": (-1, "Neumann")},
d: {"left": (0, "Dirichlet"), "right": (0, "Dirichlet")}
}

#initial conditions
model.initial_conditions = {c: 0.9, d: 0}

model.variables = {
"Concentration": c,
"Algebraic": d,
}

r = pybamm.SpatialVariable("r", domain=["negative particle"], coord_sys="cartesian")

geometry = {"negative particle": {r: {"min": 0, "max": 1}}}

spatial_methods = {"negative particle": pybamm.FiniteVolume()}
submesh_types = {"negative particle": pybamm.Uniform1DSubMesh}
var_pts = {r: 20}
#create a mesh of our geometry, using a uniform grid with 20 volumes
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

parameter_values = pybamm.ParameterValues({
"Constant k": 5,
})

solver = pybamm.CasadiSolver()

sim = pybamm.Simulation(
model,
geometry=geometry,
parameter_values=parameter_values,
submesh_types=submesh_types,
var_pts=var_pts,
spatial_methods=spatial_methods,
solver=solver,
)

solution = sim.solve([0, 1])

sim.plot(["Concentration", "Algebraic"])

[valentinsulzer]

Your code works fine for me on the latest version of pybamm

[CSchwet]

That's the working version, if I change line 11 to N = -pybamm.grad(c) - pybamm.grad(d / k), it breaks

[valentinsulzer]

Ah I see. The boundary conditions are applied to the variable inside the grad, so you have to write your model in a way that the variable in the grad is the one that you have boundary conditions for (e.g. by applying the chain/product rule). Off-by-two shape error means that the boundary conditions have not been correctly applied.

[CSchwet]

Thank you! So I should do that using the diff operator (for Neumann b.c.)?

[valentinsulzer]

What are you trying to do exactly?

[CSchwet]

I would like to automate the evaluation of the b.c., so if I define a new expression

y = d / k

how can I evaluate this for a given boundary value of d to put into the model.boundary_conditions? Let's say

d0 = pybamm.Parameter("Initial algebraic")

model.boundary_conditions = {
c: {"left": (0, "Neumann"), "right": (-1, "Neumann")},
d: {"left": (0, "Dirichlet"), "right": (d0, "Dirichlet")},
y: {} # y(d=0) and y(d=d0)?
}