Dimension error when setting variable domains for custom model

[nril0001]

I am trying to model a coupled electrochemical reaction using the BaseModel class. The reaction contains two surface-confined species (denoted by *) and two species in solution:

Ox* + e- <-> Red*
Red* + S -> Ox* + P

Currently, I have an issue with getting the equations simulating the surface-confined species and the species in solution to interact with each other. These are the model variables:

c_Ox = pybamm.Variable("O(surf) [non-dim]")
c_s = pybamm.Variable("S(soln) [non-dim]", domain="solution")
c_p = pybamm.Variable("P(soln) [non-dim]", domain="solution")
i = pybamm.Variable("Current [non-dim]")

and the PDEs:

    # PDEs 
    model.rhs = {
        c_Ox: i_f + cat_con, #i_f is the echem contribution, cat_con is chemical contribution
        i: (i_f - cat_con - i)/(Cdl * Ru), # current divided by non-dim "s" obtained from Cdl and Ru term
        c_s: pybamm.div(pybamm.grad(c_s)), 
        c_p: pybamm.div(pybamm.grad(c_p))
    }

And here is what is setting up the boundary and initial conditions as well as the spatial variables and domain geometry:

    # Setting boundary and initial conditions
    model.boundary_conditions = {
        c_Ox: {
            "right": (pybamm.Scalar(0), "Dirichlet"), #0 - no Ox should be present at xmax
            "left": (i_f, "Neumann"),
        },
        c_s: {
            "right": (pybamm.Scalar(1), "Dirichlet"),
            "left": (cat_con, "Neumann"),                    
        },

        c_p: {
            "right": (pybamm.Scalar(0), "Dirichlet"),   #always be starting with no product
            "left": (-cat_con, "Neumann"),                 
        } 
    }

    model.initial_conditions = {
        c_Ox: pybamm.Scalar(1),
        i: Cdl * Eapp.diff(pybamm.t), 
        c_s: (CS_d/Ctot),
        c_p: (CP_d/Ctot),
    }

    # set spatial variables and domain geometry
    x = pybamm.SpatialVariable('x', domain="solution")
    model.geometry = pybamm.Geometry({
        "solution": {
                x: {
                    "min": pybamm.Scalar(0),
                    "max": pybamm.Scalar(10)  
                }
        }

        })

I believe my issue might be with the domains. However, when I tried adding a domain to the c_Ox state variables or even just putting them in the same domain as c_s and c_p I got a dimension error:

test_shape
raise pybamm.ShapeError("Cannot find shape (original error: {})".format(e))

ShapeError: Cannot find shape (original error: dimension mismatch)

[brosaplanella]

Hi Can you share a minimum working example we can run at our end? At the moment some components of the equations (like i_f or cat_con) are not defined. There is certainly an issue with domains & boundary conditions though, as c_Ox has no domain but has a boundary condition.

[nril0001]

I have never encountered that before and just to be sure there isn't an error with the code I uploaded I checked and confirmed that it runs how it has been on previous occasions

[brosaplanella]

Which PyBaMM version, Python version and operating system are you using?

[nril0001]

it should be the latest pybamm, python 3.10 and windows 11

[brosaplanella]

Ok, I managed to remove the error I was getting on my end but, weirdly, then it runs fine with no errors. I suspect that the two errors might be related actually. Can you try defining Eapp = pybamm.Scalar(0) (and commenting out the current definition of Eapp) and see if the error disappears? If so, then we can focus on fixing it.

[nril0001]

Hi brosaplanella I have been able to fix the issue on my end, but now I have an issue with attempting to implement capacitative current since it requires a derivative to calculate correctly (see new discussion). Thank you for attempting to help however

[nril0001]

Hi brosaplanella. I was wondering if you could explain how you got rid of the error you mentioned above

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

[brosaplanella]

I just defined Eapp = pybamm.Scalar(0) and removed the old definition to get rid of the error I was getting (different from yours) but then it turns out the code ran fine. As the two errors might be related, I would suggest you to do the same and see if your error disappears. If so, then we can focus on fixing the error.