[Bug]: <, > operators being evaluated as booleans in functions

[HOLL95]

PyBaMM Version

24.1

Python Version

3.12

Describe the bug

Attempting to define a function in rhs using < and > operators results in the following error

   return left - right
           ~~~~~^~~~~~~
TypeError: numpy boolean subtract, the `-` operator, is not supported, use the bitwise_xor, the `^` operator, or the logical_xor function instead.

removing these operators removes the error; presumably they are causing the expression to be evaluated as a Boolean

Steps to Reproduce

import pybamm
R=8.314459848
T=209
F= 96485.3328959
E_0 = R * T / F
model = pybamm.BaseModel()
theta = pybamm.Variable("O(surf) [non-dim]")
i = pybamm.Variable("Current [non-dim]")
E_start=0.4/E_0
E_reverse=-0.4/E_0
alpha=0.5
Ru=100
k0=100
E0=0.1/E_0
Cdl=0.5
v=1
deltaE=0.15/E_0
t_reverse=(E_reverse-E_start)
omega=10
D=7.2e-6
param = pybamm.ParameterValues({
                "Diffusion Constant [cm2 s-1]": D,
                "Faraday Constant [C mol-1]": F,
                "Gas constant [J K-1 mol-1]": R,
                "Electrode Area [cm2]": 0.07,
                "Temperature [K]": T,
                "Voltage frequency [rad s-1]": omega,
                "Voltage start [V]": E_start,
                "Voltage reverse [V]": E_reverse,
                "Voltage amplitude [V]":deltaE,
                "Scan Rate [V s-1]": v,
                "Electrode Coverage [mol cm2]": 6.5e-10,
            })


Edc_forward = -pybamm.t
Edc_backwards = pybamm.t - 2*t_reverse
Eapp = E_start + \
            (pybamm.t <= t_reverse) * Edc_forward + \
            (pybamm.t > t_reverse) * Edc_backwards + \
            deltaE * pybamm.sin(omega * pybamm.t)
Eeff = Eapp - i * Ru

i_f = k0 * ((1 - theta) * pybamm.exp((1-alpha) * (Eeff - E0))
            - theta * pybamm.exp(-alpha * (Eeff - E0))
            )


model.rhs = {
    theta: i_f ,
    i: 1/(Cdl * Ru) * (i_f + Cdl * Eapp.diff(pybamm.t) - i),
}


model.algebraic = {
}


model.boundary_conditions = {
    theta: {
        "right": (pybamm.Scalar(1), "Dirichlet"),
        "left": (i_f/D, "Neumann"),
    }
}

model.initial_conditions = {
    theta: pybamm.Scalar(1),
    i: Cdl * Eapp.diff(pybamm.t),
}


x = pybamm.SpatialVariable('x', domain=["solution"])
geometry = {
    "solution": {
            x: {
                "min": pybamm.Scalar(0),
                "max": pybamm.Scalar(20)
            }}
}
var_pts = {
    x: 100
}

submesh_types = {
    "solution": pybamm.MeshGenerator(
        pybamm.Exponential1DSubMesh, {'side': 'left'}
    )
}

mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
spatial_methods = {
    "solution": pybamm.FiniteVolume()
}


model.variables = {
    "Current [non-dim]": i,
    "O(surf) [non-dim]": theta,
    "Applied Voltage [non-dim]": Eapp,
}



param.process_model(model)
geometry = model.geometry
param.process_geometry(geometry)



disc = pybamm.Discretisation(mesh, spatial_methods)
disc.process_model(model)

Relevant log output

No response