Skip to content

Issue 4884 unify hysteresis #4893

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 29 commits into
base: develop
Choose a base branch
from
Open

Conversation

rtimms
Copy link
Contributor

@rtimms rtimms commented Mar 5, 2025

Description

Fixes #5038
Fixes #5081
Fixes #5021

Unifies the hysteresis models in PyBaMM and adds heat generation due to hysteresis. In particular:

  • Adds a new BaseHysteresisOpenCircuitPotential class that sets variables for the lithiation and delithiation OCP, as well as the average (simply the mean), and the hysteresis voltage (H = U_lith - U_delith).
  • Makes sure all models, even CurrentSigmoidOpenCircuitPotential, are expressed in terms of a hysteresis state variable -1<h<1
  • Allows the initial hysteresis state to be a function of position through the electrode
  • Allows the hysteresis decay rates of the Axen and Wycisk models to be a function of stoichiometry and temperature
  • Adds a heat source term in each active material phase Q_hys = i_vol * (U - U_eq) where i_vol is the volumetric interfacial current density, U is the OCP (i.e. includes hysteresis), and U_eq is the "equilibrium OCP"
  • Renames the notebook differential-capacity-hysteresis-state.ipynb tohysteresis-state-models.ipynb and extends it to include a comparison of all of the existing hysteresis models

Related #4884

What this doesn't do from #4884:

  • Reformulate the Wycisk model to use $\gamma = K \cdot \left(\frac{\frac{\partial q_\mathrm{vol}}{\partial U}}{\left.\frac{\partial q_\mathrm{vol}}{\partial U}\right|_\mathrm{ref}}\right)^{-n}$. This would be breaking as it requires rescaling parameters. Is this a breaking change worth making? Changes to parameters have, in the past, been problematic/easy for users to miss.
  • Express all the models simply as $\frac{\partial h}{\partial t} = \gamma\cdot\left(\frac{i_\mathrm{surf}a_\mathrm{vol}}{\phi_\mathrm{act}Fc_\mathrm{max}}\right)\cdot \left(1 - \mathrm{sgn}(i_\mathrm{surf})h\right)$ using the base class, and then specify \gamma in the subclasses. This would be an easy change to make, but I'm not sure it is much harder to just specify the whole rhs in the base class and it makes the code more readable (even if there is some repetition). The way it currently is, you see the whole ODE in the set_rhs method and don't need to chase back to the base class. Happy to take input on this.

Breaking

  • Need to explicitly give the equilibrium, delithiation, and lithiation OCPs when using a hysteresis model. E.g. you must provide all three of "Negative electrode OCP [V]", "Negative electrode delithiation OCP [V]", and "Negative electrode lithiation OCP [V]"

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

@rtimms rtimms requested a review from a team as a code owner March 5, 2025 09:33
Copy link

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

@rtimms rtimms marked this pull request as draft March 5, 2025 09:33
@rtimms
Copy link
Contributor Author

rtimms commented Mar 5, 2025

@ejfdickinson would appreciate your feedback on this

Copy link

codecov bot commented Mar 5, 2025

Codecov Report

Attention: Patch coverage is 98.76543% with 3 lines in your changes missing coverage. Please review.

Project coverage is 99.11%. Comparing base (c0528c2) to head (4d716d7).
Report is 1 commits behind head on develop.

Files with missing lines Patch % Lines
.../one_state_differential_capacity_hysteresis_ocp.py 96.22% 2 Missing ⚠️
...open_circuit_potential/one_state_hysteresis_ocp.py 97.22% 1 Missing ⚠️
Additional details and impacted files
@@             Coverage Diff             @@
##           develop    #4893      +/-   ##
===========================================
- Coverage    99.12%   99.11%   -0.02%     
===========================================
  Files          304      305       +1     
  Lines        23573    23606      +33     
===========================================
+ Hits         23366    23396      +30     
- Misses         207      210       +3     

☔ 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.

@ejfdickinson
Copy link

ejfdickinson commented Mar 11, 2025

@rtimms I'll look over it.

Relating to #4884 incomplete items currently excluded from the PR:

  • If the breaking change to introduce a reference value is unappetising (fair enough), I recommend that we consider this a "specification bug fix" where we redefine from

    $\left(\frac{\partial q_\mathrm{vol}}{\partial U}\right)^{-n}$

    to

    $\left(\frac{\frac{\partial q_\mathrm{vol}}{\partial U}}{C_\mathrm{ref,vol}}\right)^{-n}$

    where $C_\mathrm{ref,vol}$ is (implicitly) hardcoded = 1 F/m3. Then the equation has legitimate unit syntax but the implementation in code doesn't need altering.

  • Current approach (avoid inheriting an equation using $\gamma$ from the base class) is fine imo.

  • In my view we should prefer a solution with future-proofed flexibility for the true equilibrium OCP to be unequal to the mean of the two branches. Not doing this now is just storing another breaking change for the future, isn't it? From my opinion and in terms of my use cases, I'd be ok to suck up a breaking change that impacts the lumped heat source. Of course, it'd be good to check this against a general policy.

Copy link

@ejfdickinson ejfdickinson left a comment

Choose a reason for hiding this comment

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

In addition to specific comments on the docs, one additional comment:

  • Should we use this as an opportunity to move from arbitrary author-name definitions to more descriptive inputs options, e.g.

"Axen" -> "one-state"
"Wycisk" -> "one-state differential capacity"

and similarly for module names? Although they may be the original reports, I don't think that "Axen" or "Wycisk" are really that intuitively descriptive for what are actually quite general hysteresis descriptions.

@rtimms
Copy link
Contributor Author

rtimms commented Jun 23, 2025

To do:

Fix #5038

and

-Should we use this as an opportunity to move from arbitrary author-name definitions to more descriptive inputs options, e.g.
"Axen" -> "one-state"
"Wycisk" -> "one-state differential capacity"

@rtimms rtimms marked this pull request as ready for review June 27, 2025 16:40
Copy link
Contributor

@aabills aabills left a comment

Choose a reason for hiding this comment

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

Overall looks fine. Forgive my ignorance, but Wycisk/Axen are specific cases of the "new models" right? In that case, should they be kept as options to make this a non-breaking change? In any case, I don't think they should be both kept as subclasses and deleted as options -- I'd be surprised if there are any users that use submodels directly rather than going through the options, so this effectively just creates non-accessible code.

"$$ U = \\text{sigmoid}\\left(-\\frac{KI}{Q_{\\text{cell}}}\\right) U_{\\text{delith}} + \\text{sigmoid}\\left(\\frac{KI}{Q_{\\text{cell}}}\\right) U_{\\text{lith}}\n",
"$$\n",
"\n",
"Where $K$ is a fitting parameter (in the current implementation this is hard-coded to 100, as in the original publication), $I$ is the cell current, and $Q_{\\text{cell}}$ is the cell capacity. To simplify the comparison with the other models, we can rewrite this expression in terms of the hysteresis state variable $h$ given by\n",
Copy link
Contributor

Choose a reason for hiding this comment

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

didn't look at this all that closely but I think it's important here to make the things which are functions of time functions of time in the documentation

if (isinstance(ocp_option, str) and ocp_option == "Wycisk") or (
isinstance(ocp_option, tuple) and "Wycisk" in ocp_option
):
raise pybamm.OptionError(
Copy link
Contributor

Choose a reason for hiding this comment

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

is there a good reason to throw an error rather than a deprecation warning here? This could be less breaking

if (isinstance(ocp_option, str) and ocp_option == "Axen") or (
isinstance(ocp_option, tuple) and "Axen" in ocp_option
):
raise pybamm.OptionError(
Copy link
Contributor

Choose a reason for hiding this comment

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

here too


variables.update(
{
f"{Domain} electrode {phase_name}equilibrium open-circuit potential [V]": ocp_surf,
Copy link
Contributor

Choose a reason for hiding this comment

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

I think it's going to be pretty confusing for users to have this variable...

self.rhs[h] = dhdt


class WyciskOpenCircuitPotential(
Copy link
Contributor

Choose a reason for hiding this comment

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

these can't really be accessed without a lot of work right? What's the point in keeping it?

Copy link
Contributor

Choose a reason for hiding this comment

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

What if this implemented the old model, since the new model is slightly different? And then remove the option errors. That would make this a non-breaking change.

@@ -146,11 +150,27 @@ def _get_standard_coupled_variables(self, variables):
# Reversible electrochemical heating
dUdT_p = variables[f"Positive electrode {phase}entropic change [V.K-1]"]
Q_rev_p += a_j_p * T_p * dUdT_p
# Hysteresis electrochemical heating
if ocp_option == "single":
Copy link
Contributor

Choose a reason for hiding this comment

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

I think this should cite something, it's not necessarily obvious to most engineers that it should be included.

Copy link
Member

@valentinsulzer valentinsulzer left a comment

Choose a reason for hiding this comment

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

Looks good, a couple of minor comments for the notebook.
I also prefer the previous names (Axen/Wycisk) over the new names which are a real mouthful, but not a strong opinion

Copy link
Member

Choose a reason for hiding this comment

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

unfinished sentence "We add some"

Copy link
Member

Choose a reason for hiding this comment

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

Could also do with a sentence at the end describing the results

@kratman
Copy link
Contributor

kratman commented Jul 2, 2025

I am looking into the docs failure now, probably something pretty minor

@kratman
Copy link
Contributor

kratman commented Jul 2, 2025

I am looking into the docs failure now, probably something pretty minor

I pushed a temporary fix for the failing citation so that review could continue, I will push a in a real fix once I get it working again

@kratman
Copy link
Contributor

kratman commented Jul 2, 2025

Docs are fixed now. I also updated the description to make sure all tickets get closed

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
6 participants