Add in-plane current to 2D models

[rtimms]

At the moment you can't access the in-plane current in the current collectors

You can add the variable to get the gradient of the potential

import pybamm

options = {"cell geometry": "pouch", "thermal": "x-lumped", "dimensionality": 2, "current collector": "potential pair"}
model = pybamm.lithium_ion.SPMe(options=options)
phi_s_cn = model.variables["Negative current collector potential"]
grad_phi_s_cn = pybamm.grad(phi_s_cn)
model.variables.update(
    {
        "Gradient of negative current collector potential": grad_phi_s_cn,
    }
)
experiment = pybamm.Experiment([
    "Discharge with 5 C for 1 minute"
])
sim = pybamm.Simulation(model, experiment=experiment)
sim.solve()

but when you try to access it you get the error

NotImplementedError: Wrong number or type of arguments for overloaded function 'new_Function'.
  Possible prototypes are:
    Function()
    Function(str)
    Function(Function)
    Function(str,[SX],[SX],dict)
    Function(str,[MX],[MX],dict)
    Function(str,dict:SX,[str],[str],dict)
    Function(str,dict:MX,[str],[str],dict)
    Function(str,[SX],[SX],[str],[str],dict)
    Function(str,[MX],[MX],[str],[str],dict)
  You have: '(str,[MX],[NoneType])'

you can do a workaround, but it's ugly

# manually discretise gradient of potential 
grad_phi_s_cn_disc = sim._disc.process_symbol(pybamm.grad(phi_s_cn))
# extract t and y
t = sim.solution.t
y = sim.solution.y.full()
# evaluate at initial time (for example)
grad_phi_s_cn_eval = grad_phi_s_cn_disc.evaluate(t=t[0],y=y[:,0])
# these will both be column vectors (can reshape back to matrix if required)
grad_phi_s_cn_y = grad_phi_s_cn_eval[:, 0]
grad_phi_s_cn_z = grad_phi_s_cn_eval[:, 1]
# define current in y and z direction (aspect ratio delta appears due to scaling)
sigma_cn = sim._parameter_values.evaluate(model.param.sigma_cn)
delta = sim._parameter_values.evaluate(model.param.delta)
i_cn_y = -sigma_cn * delta * grad_phi_s_cn_y
i_cn_z = -sigma_cn * delta * grad_phi_s_cn_z
# rescale to dimensoinal current in [A/m2]
i_scale = sim._parameter_values.evaluate(model.param.i_typ)
i_cn_y_dim = i_cn_y * i_scale
i_cn_z_dim = i_cn_z * i_scale

and I don't think it's giving quite the right thing at the moment. Need to look into it.

[rtimms]

Note: we do have a test that the gradient is doing the right thing in 2D here