Changed fundamental SEI variable from thickness to concentration

[DrSOKane]

Description

Previously, lithium was not conserved if both SEI growth and loss of active material were enabled. To fix this, the fundamental variable to be solved for is now SEI concentration, not thickness. This has the added bonus of simplifying the SEI on cracks model.

Fixes #3006

Type of change

Please add a line in the relevant section of CHANGELOG.md to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

• New feature (non-breaking change which adds functionality)
• Optimization (back-end change that speeds up the code)
• Bug fix (non-breaking change which fixes an issue)

Key checklist:

• No style issues: $ pre-commit run (see CONTRIBUTING.md for how to set this up to run automatically when committing locally, in just two lines of code)
• All tests pass: $ python run-tests.py --all
• The documentation builds: $ python run-tests.py --doctest

You can run unit and doctests together at once, using $ python run-tests.py --quick.

Further checks:

• Code is commented, particularly in hard-to-understand areas
• Tests added that prove fix is effective or that feature works

[codecov]

Codecov Report

Attention: 2 lines in your changes are missing coverage. Please review.

Comparison is base (063a383) 99.59% compared to head (4a80658) 99.58%.

Files	Patch %	Lines
...mm/models/submodels/interface/sei/sei_thickness.py	97.64%	2 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #3171      +/-   ##
===========================================
- Coverage    99.59%   99.58%   -0.01%     
===========================================
  Files          257      258       +1     
  Lines        20802    20822      +20     
===========================================
+ Hits         20718    20736      +18     
- Misses          84       86       +2     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[valentinsulzer]

Thanks, this is great! Exactly what I was envisioning.
Could you post some "before and after" plots of SEI thickness and concentration vs time just to confirm that the results are unchanged (except where the lithium conservation bugs are fixed)?
Are there any tests you can think of adding to verify that lithium is now conserved (these tests should fail on develop but pass on this branch)?

[valentinsulzer]

Thanks, test looks good, need to fix the failing example though

[DrSOKane]

Hi Tino, the test is NOT good as it passed even on develop. I changed it so it now uses the loss of lithium due to loss of active material variable, but after much experimenting I found this variable was calculated using the incorrect assumption that the average of a product is the same as the product of averages. I managed to fix this, but the test still fails!

Here are three plots, all of which were calculated using the following code:

model = pybamm.lithium_ion.DFN({
    "SEI": "reaction limited",
    "SEI porosity change": "true",
    "lithium plating": "partially reversible",
    "particle mechanics": ("swelling and cracking", "swelling only"),
    "SEI on cracks": "true",
    "loss of active material": "stress and reaction-driven",
})
param = pybamm.ParameterValues("OKane2022")
var_pts = {"x_n": 5, "x_s": 5, "x_p": 5, "r_n": 30, "r_p": 20}
solver = pybamm.CasadiSolver(mode="safe")
sim = pybamm.Simulation(
    model, parameter_values=param, var_pts=var_pts, solver=solver
)
solution = sim.solve([0, 3500])
t = solution["Time [s]"].entries
particles = solution["Total lithium in particles [mol]"].entries
neg = solution["Loss of lithium due to loss of active material in negative electrode [mol]"].entries
pos = solution["Loss of lithium due to loss of active material in positive electrode [mol]"].entries
side = solution["Total lithium lost to side reactions [mol]"].entries
plt.figure()
plt.plot(t,particles+side+neg+pos)
plt.xlabel("Time [s]")
plt.ylabel("Lithium in particles + lithium lost [mol]")
plt.show()

The third plot is a massive improvement on the first two but still doesn't meet the requirement of 12 decimal places.

Regarding the failing example, I have no idea how to fix it! It's saying Loss of lithium to SEI [mol] has a domain [current collector]. I ran it in debug mode, set a checkpoint in _get_standard_total_concentration_variables and every time the checkpoint was triggered, the variable in question had a domain of [], as expected for a yz_averaged variable.

Give yourself a gold star for reading this essay!

[RuiheLi]

A bit more information on the glithces Simon found regarding Loss of lithium to SEI [mol] has a domain [current collector].
I find that: (1) this only occurs in SPM model, for SPMe or DFN it is fine; (2) this may relate to how PyBaMM handle solution.
If we print the type of Loss of lithium to SEI [mol], which is Q_sei in def _get_standard_total_concentration_variables under pybamm/models/submodels/interface/sei/base_sei.py, then it shows that it is 'pybamm.expression_tree.scalar.Scalar'> or 'pybamm.expression_tree.binary_operators.Multiplication'>

however, when I do type(sol["Loss of lithium to SEI [mol]"]), it shows that it is pybamm.solvers.processed_variable.ProcessedVariable

[rtimms]

@DrSOKane, thanks for this! have you had any luck in debugging the notebook? would be good to get this in for the next release

[brosaplanella]

Might be worth to merge develop into this branch first, as we have improved how notebooks are tested, which might help debug.

[agriyakhetarpal]

Note: that also means they need to be written in a slightly different way, i.e. in a literate programming sense rather than as a Python script which they were earlier converted to at the time of testing

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[DrSOKane]

Update: It's worse than I thought. This bug also exists for the lithium plating model! If you run either SEI or lithium-plating with "x-average side reactions" set to true, the relevant loss of lithium variable has the domain "current collector" when it shouldn't.