Skip to content

Improve numerical stability of exchange current density and OCPs#5371

Merged
MarcBerliner merged 30 commits intomainfrom
regularise-functions
Feb 5, 2026
Merged

Improve numerical stability of exchange current density and OCPs#5371
MarcBerliner merged 30 commits intomainfrom
regularise-functions

Conversation

@MarcBerliner
Copy link
Member

@MarcBerliner MarcBerliner commented Feb 2, 2026

Description

Exchange current densities include sqrt(c_e), sqrt(c_s), and sqrt(c_max - c_s) terms that blow up when concentrations -> 0/c_max. Inverse BV kinetics do arcsinh(j/j0), which also blows up when j0 is small or causes divide by zero errors for j0=0. These cause solver failures in edge cases

This PR adds regularised versions of two expressions:

  • RegPower(x, a): approximates |x|^a * sign(x) but stays linear near zero (based on modelica's regPow)
  • Arcsinh2(a, b): computes arcsinh(a/b) while avoiding blowup/division by zero errors as b -> 0

Major changes:

  • The exchange current density now automatically replaces all Sqrt/Power -> RegPower to make them more numerically stable at small concentrations
  • The barrier function for the OCP is changed from 1e-6/sto to a more complicated logarithmic expression, whose parameters were tuned to provide a barrier of 1 mV at sto = 0.001 and 1000 mV at sto = 0. This is a significant from the previous version because we no longer have an asymptote as sto -> 0. Some calculations that relied on this assumption, like the ESOH solver, will have slightly different answers in this extreme limit. Most usecases will not change

Type of change

Please add a line in the relevant section of CHANGELOG.md to document the change (include PR #)

Important checks:

Please confirm the following before marking the PR as ready for review:

  • No style issues: nox -s pre-commit
  • All tests pass: nox -s tests
  • The documentation builds: nox -s doctests
  • Code is commented for hard-to-understand areas
  • Tests added that prove fix is effective or that feature works

@MarcBerliner MarcBerliner changed the title Regularise exchange current density and inverse kinetics Improve numerical stability of exchange current density and OCPs Feb 2, 2026
@codecov
Copy link

codecov bot commented Feb 3, 2026

Codecov Report

❌ Patch coverage is 91.76030% with 22 lines in your changes missing coverage. Please review.
✅ Project coverage is 98.32%. Comparing base (c18ea21) to head (66a0155).
⚠️ Report is 1 commits behind head on main.

Files with missing lines Patch % Lines
src/pybamm/expression_tree/functions.py 86.04% 18 Missing ⚠️
src/pybamm/expression_tree/binary_operators.py 90.90% 2 Missing ⚠️
...rc/pybamm/expression_tree/operations/regularise.py 98.30% 1 Missing ⚠️
...nterface/kinetics/inverse_kinetics/base_inverse.py 0.00% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main    #5371      +/-   ##
==========================================
+ Coverage   98.03%   98.32%   +0.28%     
==========================================
  Files         327      328       +1     
  Lines       28690    28941     +251     
==========================================
+ Hits        28126    28455     +329     
+ Misses        564      486      -78     

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

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@MarcBerliner MarcBerliner marked this pull request as ready for review February 3, 2026 17:08
@MarcBerliner MarcBerliner requested a review from a team as a code owner February 3, 2026 17:08
Copy link
Member

@BradyPlanden BradyPlanden left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @MarcBerliner - this is looking good. Just a few questions / comments

BradyPlanden
BradyPlanden previously approved these changes Feb 5, 2026
Copy link
Member

@BradyPlanden BradyPlanden left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks great, thanks @MarcBerliner! Just one non-blocking comment

Co-authored-by: Brady Planden <55357039+BradyPlanden@users.noreply.github.com>
@MarcBerliner MarcBerliner merged commit 86531fb into main Feb 5, 2026
24 of 26 checks passed
@MarcBerliner MarcBerliner deleted the regularise-functions branch February 5, 2026 14:57
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants