pybamm.grad and (SpatialVariable or pybamm.laplacian) discretise differently?

[ZeliangSu]

PyBaMM Version

0.2.3

Python Version

3.6

Describe the bug

Hi, community,
Would like to implement an equation as the simple snippet below, however ran into a dimension problem, saying the shape of rhs and ics are not matching. Is it an expected behavior? Or something I got wrong?
Can somebody kindly explain and provide a correct one? Thanks a lot.

Steps to Reproduce

import pybamm
import numpy as np
model = pybamm.BaseModel()
c_s = pybamm.Variable('c_s', domain='positive particle')
r = pybamm.SpatialVariable('r', domain=['positive particle'], coord_sys="spherical polar")
model.variables = {'c_s': c_s}

model.rhs = {
    c_s: pybamm.laplacian(c_s) + 1/r * pybamm.grad(c_s) #  <=========== Here. grad() and (r or laplacian) discretises differently?
    # c_s: pybamm.laplacian(c_s) + 2/r * c_s.diff(r) #  then I provide my work around here using a diff(), but not sure if it's correct
}

geometry = { "positive particle": {r: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}}
submesh_types = {"positive particle": pybamm.Uniform1DSubMesh}
var_pts = {r: 30}  # 30 bins
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
spatial_methods = {"positive particle": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)

model.initial_conditions = {c_s: 1}
model.boundary_conditions = {c_s: {'left': (pybamm.Scalar(0), 'Neumann'), 'right': (pybamm.Scalar(2), 'Dirichlet')}}

disc.process_model(model)
solver = pybamm.CasadiSolver()
t_eval = np.linspace(0, 1, 600)
solution = solver.solve(model, t_eval)

Relevant log output

AssertionError: rhs and initial_conditions must have the same shape after discretisation but rhs.shape = (31, 1) and initial_conditions.shape = (30, 1) for variable 'c_s'.

[dion-w]

Maybe not from the PyBaMM side but from a more general viewpoint. The gradient calculates the derivative, so it calculates the changes between volume-elements. You can view it as the function being the values in the discretization elements, wherease the gradient calculates the weighted differences between the volume elements. The difference between the volume elements is naturally one dimension less than the actual values (example, three discretization elements where every elements has a temperature value, but there are only two values for the temperature difference/derivative/gradient namely the temperature difference between the first and the second, and the difference between the second and the third element).

Discretization-wise you could say: the gradient evaluates on the internal edges and has therefore one value less than the actual function and since the laplacian is the second derivative it also may have a different dimension. So i would suggest that you're checking your differential equations first. There maybe a functionality in PyBaMM though that helps you dealing with this problem on the software side but i don't know.

[brosaplanella]

Following on @dion-w's comment, note that your equation is not consistent as you are adding together pybamm.laplacian(c_s) (a scalar) with 1/r * pybamm.grad(c_s) (a vector). If you are trying to write the diffusion equation in cylindrical or spherical coordinates, note that you just need the Laplacian as PyBaMM will take care of the geometry.

BTW, I don't think it is related to your problem, but you are using a version of PyBaMM (0.2.3) which is almost two years old. I suggest updating it to a more recent version. In addition, note that from v21.9, Python 3.6 is not supported so you might need to update Python too.

This issue is better suited as a discussion, so I am transferring it there.

[ZeliangSu]

Thank you both,
@dion-w
'You can view it as the function being the values in the discretization elements...' I'm still trying to understand the Broadcast module and Discretisation one, and how they work together. To my surprise, you can see in my error that the bins changed from (30, 1) to (31, 1). ‘the gradient evaluates on the internal edges and has therefore one value...’ I don't know if his change is caused by a evaluate_on_edge.
@brosaplanella
'If you are trying to write the diffusion equation in cylindrical or spherical coordinates...' yes, that's exactly what I'm trying to do, trying to restreint it in 1D by ignoring the angular term. If I understood well, you are suggesting coord_sys="spherical polar" then c_s: pybamm.laplacian(c_s) without the rest, pybamm will do the transformation automatically? PyBaMM (0.2.3) updated now, thank you.

[valentinsulzer]

Here is the information on discretisation and broadcast