Issue with pybamm.DefiniteIntegralVector causing ShapeError

[FaissalOSU]

Hello everyone,

I am working on building a submodel in PyBaMM. Everything in the submodel works as expected until I reach the following line:

pressure = pybamm.DefiniteIntegralVector(pressure_gradient)

This line causes the following error:

pybamm.expression_tree.exceptions.ShapeError: Cannot find shape
(original error: operands could not be broadcast together with shapes (19,1) (21,1))

I have double-checked, and all shapes and domains appear consistent.

Below is the relevant portion of the code for context:

Parameters
----------
param : parameter class
    The parameters to use for this submodel
domain : str
    The domain of the model either 'Negative' or 'Positive'
options: dict
    A dictionary of options to be passed to the model.
    See :class:`pybamm.BaseBatteryModel`
phase : str, optional
    Phase of the particle (default is "primary")

"""

def __init__(self, param, domain, options, phase):
    # domain_map = {
    #     "Positive": "positive",
    #     "Negative": "negative",
    #     "positive electrode": "positive",
    #     "negative electrode": "negative"
    # }
    # full_domain = domain_map.get(domain, domain)
    super().__init__(param, domain, options, phase)
    # self.domain = full_domain

def get_fundamental_variables(self):
    domain_param = self.domain_param
    domain, Domain = self.domain_Domain
    # phase_name = self.phase_name
    # phase_param = self.phase_param
    # Define electrode-level pressure variable (on electrode domain)
    pressure = pybamm.Variable(
        f"{Domain} particle surface pressure [Pa]",
        domain=self.domain,
    )

    stress_t_surf = pybamm.Variable(
        f"{Domain} particle surface tangential stress [Pa]",
        domain=self.domain,
    )

    stress_r_surf = pybamm.Variable(
        f"{Domain} particle surface radial stress [Pa]",
        domain=self.domain,
    )

    variables = {
        f"{Domain} particle surface pressure [Pa]": pressure,
        f"{Domain} particle surface tangential stress [Pa]": stress_t_surf,
        f"{Domain} particle surface radial stress [Pa]": stress_r_surf,
    }

    return variables

def get_coupled_variables(self, variables):
    domain_param = self.domain_param
    domain, Domain = self.domain_Domain
    # phase_name = self.phase_name
    # phase_param = self.phase_param
    c_s_rav = variables[
        f"R-averaged {domain} particle concentration [mol.m-3]"
    ]

    c_s_surf = variables[
        f"{Domain} particle surface concentration [mol.m-3]"
    ]
    T_xav = variables["X-averaged cell temperature [K]"]
    T = pybamm.PrimaryBroadcast(
        variables[f"{Domain} electrode temperature [K]"],
        [f"{domain} particle"],
    )
    eps_s = variables[
        f"{Domain} electrode active material volume fraction"
    ]

    # use a tangential approximation for omega
    c_0 = 50
    R0 = 8.31

    c_s_max = 30000  # mol/m^3, typical for graphite
    x_LiS = c_s_surf / c_s_max

    # Hardcoded material and geometric parameters
    R = 3e-6          # Particle radius in meters
    E = 15e9          # Young's modulus in Pascals
    nu = 0.3          # Poisson's ratio
    omega = 0.08      # partial molar volume change

    # Displacement and gradient
    # R0 / 3 * (c_s_rav - c_0)
    u = R * ((1 + (c_s_rav - c_0) * omega)**(1/3) - 1)
    w = u / R

    # Christensen stress equations
    stress_r_surf = (E / ((1 + nu) * (1 - 2 * nu))) * ((1 - nu) * w -
                                                       0.5 * w**2 + 2 * nu * (u / R - 0.5 * (u / R)**2))
    stress_t_surf = (E / ((1 + nu) * (1 - 2 * nu))) * \
        (nu * w - 0.5 * w**2 + (u / R - 0.5 * (u / R)**2))

    # pressure = 5*c_s_rav

    if self.domain == "negative":
        xi = pybamm.standard_spatial_vars.x_n
    elif self.domain == "positive":
        xi = pybamm.standard_spatial_vars.x_p

    # xi_broad = pybamm.SecondaryBroadcast(
    #     xi, self.domain)
    # print(stress_t_surf.shape)
    # print(stress_r_surf.shape)
    # # print(xi.shape)

    # xi_broad = pybamm.Broadcast(xi, self.domain)

    dsigma_r_dx = pybamm.grad(stress_r_surf)

    # dsigma_r_dx = pybamm.FullBroadcast(
    #     dsigma_r_dx, [f"{self.domain} electrode"])
    pressure_gradient = dsigma_r_dx + \
        (2 / xi) * (stress_r_surf - stress_t_surf)
    print(f"Domain streess{stress_r_surf.domains}")
    print(f"Domain grad{dsigma_r_dx.domains}")

    # print(f"The shape is: {dsigma_r_dx.evaluate_for_shape()}")
    # print(f"The shape is: {pressure_gradient.evaluate_for_shape()}")
    # print(f"The shape is: {stress_r_surf.evaluate_for_shape()}")
    # print(f"The shape is: {xi.evaluate_for_shape()}")

    # # Integrate to get pressure (assuming p=0 at surface ξ=1)

    # pressure_gradientB = pybamm.FullBroadcast(
    #     pressure_gradient, self.domain)
    # xiB = pybamm.SpatialVariable("x", domain=[self.domain])
    print()
    print(f"Domain pressure_gradient{pressure_gradient.domains}")
    print(f"Domain xi{xi.domains}")
    print("pressure_gradient shape:", pressure_gradient.shape_for_testing)
    print("stress_r_surf shape:", stress_r_surf.shape_for_testing)
    print("xi shape:", xi.shape_for_testing)

    # # print(self.domain)
    # pressure_gradient_broad = pybamm.PrimaryBroadcast(
    #     pressure_gradient, [self.domain])

    # # Integrate over electrode domain
    # xi_e = pybamm.SpatialVariable("x", domain=[self.domain])
    pressure = pybamm.DefiniteIntegralVector(pressure_gradient)

    variables.update(
        {
            f"{Domain} particle surface radial stress [Pa]": stress_r_surf,
            f"{Domain} particle surface tangential stress [Pa]": stress_t_surf,
            f"{Domain} particle surface pressure [Pa]": pressure,
        }
    )

    return variables

Has anyone encountered a similar issue with DefiniteIntegralVector? Any suggestions on how to resolve this shape mismatch while keeping the calculation fully symbolic would be greatly appreciated.

Thank you in advance!