handling of concentrations not getting negative

[dion-w]

Hello everyone,
yesterday i checked out the new version of PyBaMM. In my setup i usually encounter thick eledtrodes with high amounts of elctrolyte depletion (see #1982). For this reason i introduced back in version 22.3 the electrolyte concentration as a heaviside step function which sets every concentration smaller 0 to 1e-10.
I wanted to introduce this in version 22.10 in base_electrolyte_diffusion.py in _get_standard_domain_concentration_variables but somehow the pybamm.EqualHeaviside changes the shape of the array.

A minimum working example is the following:

import pybamm

model = pybamm.lithium_ion.DFN()

sim = pybamm.Simulation(model)
sim.solve([0, 3600])
sim.plot()

While having the follwoling line introduced after line 64 in base_electrolyte_diffusion.py in _get_standard_domain_concentration_variables:

c_e_k = 1e-10 + (c_e_k - 1e-10) * pybamm.EqualHeaviside(0, c_e_k)

So why do you think EqualHeaviside changes the shape of the array? Or is this behaviour of "concentrations must not be negative" already handled somehow differently?

[valentinsulzer]

I think what's happening is that you're changing the variable that goes into the grad(c_e) term, which means it can't find the correct boundary conditions associated with that variable.

We now have some functions that make c_e positive before it goes into specific functions, look for pybamm.settings.tolerances in lithium_ion_parameters.py. This seems to help with allowing the solver to finish, though c_e itself still goes a little bit negative

[valentinsulzer]

Usually when you get shape errors that are off by 2, it means it didn't find the correct boundary conditions (since boundary conditions add a ghost node on each side)

[dion-w]

how can i add boundary conditions then? my basic attemps on envoking it as a method on the variable or function from pybamm. didn't work.