[Bug]: rhs and initial_conditions have shape mismatch after discretisation

[pvprajwal]

PyBaMM Version

25.2

Python Version

3.11.0

Describe the bug

When I implement spherical particle diffusion equation with elastic-stress coupling and set model.rhs and model.initial_conditions, I get an error that rhs and initial conditions have shape mismatch discretisation. They differ by one row for the both the particle variables (rhs (21,1) vs initial (20,1)).

I think this is due to a primary/secondary broadcasting node mismatch introduced by pybamm.grad/pybamm.div, but I am not able to find a fix

Steps to Reproduce

I have written a sample code that uses the diffusion equation in expanded form. It combines pybamm.grad (dc/dr) and pybamm.div (d2c/dr2) expressions. This is throwing the mismatch error.

import pybamm

model = pybamm.BaseModel()
c = pybamm.Variable("Concentration", domain="negative particle")
r_n = pybamm.SpatialVariable("r_n", domain="negative particle")

D = pybamm.Scalar(1)
theta = pybamm.Scalar(1)
dc_dr = pybamm.grad(c)
d2c_dr2 = pybamm.div(pybamm.grad(c))

dc_dt = D * (d2c_dr2 + 2 / r_n * dc_dr + theta * dc_dr**2 + theta * c * (d2c_dr2 + 2 / r_n * dc_dr))

model.rhs = {c: dc_dt}
model.initial_conditions = {c: pybamm.Scalar(1)}
model.boundary_conditions = {
    c: {
        "left": (pybamm.Scalar(0), "Neumann"),
        "right": (pybamm.Scalar(0), "Neumann"),
    }
}

geometry = {"negative particle": {"r_n": {"min": 0, "max": 1}}}
var_pts = {"r_n": 20}
submesh_types = {"negative particle": pybamm.Uniform1DSubMesh}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
spatial_methods = {"negative particle": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)

disc.process_model(model)

pybamm.expression_tree.exceptions.ModelError: rhs and initial conditions must have the same shape after discretisation but rhs.shape = (21, 1) and initial_conditions.shape = (20, 1) for variable 'Concentration'.

Relevant log output

[rtimms]

Some of the terms in your rhs evaluate on edges in the finite volume scheme (gradients) and other on nodes (divergences of gradients). Note you can set a coordinate system for the spatial variable r_n = pybamm.SpatialVariable("r_n", domain=["negative particle"], coord_sys="spherical polar") so you dont need to explicitly work out the coordinate transforms. It would be simpler to write your equation as a divergence of a flux.