[Bug]: Wycisk model differential capacity variable error

[Daniel-Nicolae23]

PyBaMM Version

25.4.2

Python Version

3.12.7

Describe the bug

One bug and one discussion question really.

The bug: I'm unable to access the differential capacity variable after simulating a model with the Wycisk option:

sol["Negative electrode secondary differential capacity [A.s.V-1]"]
# or sol["Negative electrode differential capacity [A.s.V-1]"] for single phase anode behaves the same

The numpy boolean negative, the `-` operator, is not supported, use the `~` operator or the logical_not function instead.

I tried to do some debugging by looking at the WyciskOpenCircuitPotential class. Printing dQdU from where it's defined at:

dU = self.phase_param.U(sto_surf, T_bulk).diff(sto_surf)
dQdU = Q_mag / dU

prints something that's probably right at the model initialisation:

0.0002777777777777778 * yz-average(x-average(Secondary: Negative electrode active material volume fraction)) * Negative electrode thickness [m] * Electrode width [m] * Electrode height [m] * Secondary: Maximum concentration in negative electrode [mol.m-3] * Faraday constant [C.mol-1] 
/ 
(Secondary: Negative electrode OCP [V] + (yz-average(Ambient temperature [K]) - Reference temperature [K]) * Secondary: Negative electrode OCP entropic change [V.K-1] + 1e-06 * (-(1e-10 <= minimum(boundary value(X-averaged negative secondary particle concentration [mol.m-3]) / Secondary: Maximum concentration in negative electrode [mol.m-3], 0.9999999999)) * (boundary value(X-averaged negative secondary particle concentration [mol.m-3]) / Secondary: Maximum concentration in negative electrode [mol.m-3] <= 0.9999999999) / ((maximum(minimum(boundary value(X-averaged negative secondary particle concentration [mol.m-3]) / Secondary: Maximum concentration in negative electrode [mol.m-3], 0.9999999999), 1e-10)) ** 2.0) + -(1e-10 <= minimum(boundary value(X-averaged negative secondary particle concentration [mol.m-3]) / Secondary: Maximum concentration in negative electrode [mol.m-3], 0.9999999999)) * (boundary value(X-averaged negative secondary particle concentration [mol.m-3]) / Secondary: Maximum concentration in negative electrode [mol.m-3] <= 0.9999999999) / ((-1.0 + maximum(minimum(boundary value(X-averaged negative secondary particle concentration [mol.m-3]) / Secondary: Maximum concentration in negative electrode [mol.m-3], 0.9999999999), 1e-10)) ** 2.0)))

Or at least the numerator is for sure correct, and I'll trust the denominator too 😅.

Later, in the set_rhs(), I see we're only extracting a component from this variable:

dQdU = dQdU.orphans[0]

and using that inside the ODE parameter. Printing this new dQdU here just gives the numerator from above, i.e. the phase capacity:
0.0002777777777777778 * yz-average(x-average(Secondary: Negative electrode active material volume fraction)) * Negative electrode thickness [m] * Electrode width [m] * Electrode height [m] * Secondary: Maximum concentration in negative electrode [mol.m-3] * Faraday constant [C.mol-1]

So the question: are we not removing the stoichiometry dependency of the $k(z)$ rate by using this? What is the purpose of the orphans[0] line in set_rhs()?

I tried commenting it, hoping to put back the full dQdU in the ODE. If I do that, I get exactly the same error message as above, but this time it's when running sim.solve().

Steps to Reproduce

import pybamm

model = pybamm.lithium_ion.SPMe({
    "working electrode": "both",
    "open-circuit potential": (("single", "Wycisk"), "single"),
    "particle phases": ("2", "1"),
})

parameters_DCHS = pybamm.ParameterValues("Chen2020_composite")
lithiation_ocp = parameters_DCHS["Secondary: Negative electrode lithiation OCP [V]"]
delithiation_ocp = parameters_DCHS["Secondary: Negative electrode delithiation OCP [V]"]

def ocp_avg(sto):
    return (lithiation_ocp(sto) + delithiation_ocp(sto)) / 2

parameters_DCHS.update({
    "Secondary: Initial hysteresis state in negative electrode": 1,
    "Secondary: Negative particle hysteresis decay rate": 1,
    "Secondary: Negative particle hysteresis switching factor": 2,
    "Secondary: Negative electrode OCP [V]": ocp_avg,
}, check_already_exists=False)


experiment = pybamm.Experiment([
    "Rest for 1 second",
    "Discharge at 1C for 3 minutes"
])
sim = pybamm.Simulation(model, parameter_values=parameters_DCHS, experiment=experiment, solver=pybamm.IDAKLUSolver())
sol = sim.solve(calc_esoh=False)

sol["Negative electrode secondary differential capacity [A.s.V-1]"].entries

Relevant log output