Skip to content
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.

Sign up for GitHub

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

Jump to bottom

adds B/2 term to Paasch element, fixes #251 #256

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

adds B/2 term to Paasch element, fixes #251 #256

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5ed337b
adds B/2 term to Paasch element, fixes #251
mdmurbach Feb 20, 2023
5c36be8
update tests for new Paasch equation
mdmurbach Feb 20, 2023
6a72c28
break into multiple lines
mdmurbach Feb 20, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion impedance/models/circuits/elements.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -345,6 +345,7 @@ def T(p, f):
.. math::

Z = A\\frac{\\coth{\\beta}}{\\beta} + B\\frac{1}{\\beta\\sinh{\\beta}}
+ \\frac{B}{2}
Copy link
Collaborator

Choose a reason for hiding this comment

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

Oh boy, this one is a doozy. So, I haven't really read the paper to understand the full context of this model, but am I correct in assuming that A, B, a, and b are defined in order to reduce the number of fitting parameters (from 5 down to 4, since as written we'd have $\rho_1$, $\rho_2$, $k$, $K$, and $d$)?

I understand the point of that, but it does come at the cost of additional complexity while reading the code since the audience has to work through the conversion of A to $\rho_1$ and $\rho_2$, etc. Rewriting the code in terms of parameters as written in the paper would be (I think) a breaking change since it introduces an additional fitting parameter for the element, and changes the meaning of each parameter. I'm not super inclined towards this unless we first added a deprecation notice for a couple of minor version releases.

Regardless, I think there's some documentation we should add to this element. The first thing is that Eq. 22 in the paper is for $AZ(\omega)$, where A is the geometric area of the electrode. This means $Z(\omega)$ has units of $\Omega cm^{-2}$, but I believe all other elements are written with units $\Omega$. So I think we need to call that out. Maybe add something like (CAUTION: This element, in agreement with Eq. 22 Paasch et al., has units of Ohms cm ^-2. All other elements in this package )

It's also a bit confusing that the paper uses "A" as an area, but the code uses A as a lumped parameter. Maybe we can refactor A and B to M and N? I'm not attached to M and N, but ideally we use something that isn't already meaningful in electrochemistry, and isn't already used in the paper. (Kind of tough when the paper uses nearly the whole Latin and Greek alphabets!)

Another thing is that we have $\beta$ defined slightly differently than the paper, where $\beta_{code} = d \beta_{paper}$. That fact is a little bit hidden by how $\beta$ is defined in terms of lumped parameters a and b. Maybe we can add something like NOTE: \\beta defined here differed from Eq. 23 of Paasch et al. by a factor of d

Copy link
Contributor

@etrevis etrevis Mar 21, 2023 •
edited
Loading

Choose a reason for hiding this comment

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

I regularly use this element and I can confirm the current bundle of parameters makes interpretation and initial estimation/guess cumbersome.

If I had to refactor the code I would add electrode area 'A' and thickness 'd' as optional parameters that default to a value of 1. With both having a value of 1, the electronic and ionic resistivities $\rho_1$ , $\rho_2$ (unit $\Omega \cdot m$) would then assume the meaning of electronic and ionic resistance of the sample generating the dataset. Area and thickness are normalization values and should not be included in the fitting/optimization directly imho.

Something like this would have 4 open parameters ( $\rho_1$ , $\rho_2$, $k$, $K$ ) , coherence with the paper and a 'straightforward' physical interpretation:

$Z = \frac{d}{A} \left[ \frac{\rho_1^2 + \rho_2^2}{\rho_1 + \rho_2} \cdot \frac{\coth{\beta}}{\beta} + \frac{2\rho_1\rho_2}{\rho_1 + \rho_2} \cdot \frac{1}{\beta\sinh{\beta}}+ \frac{\rho_1\rho_2}{\rho_1 + \rho_2} \right]$ with $\beta = \left(\frac{k+i\omega}{K}\right)^\frac{1}{2}$

Suggested change
+ \\frac{B}{2}
+ Z = \\frac{d}{A} \left[ \\frac{\rho_1^2 + \rho_2^2}{\rho_1 \rho_2} \cdot \\frac{\\coth{\\beta}}{\\beta} + \\frac{2\rho_1\rho_2}{\rho_1 + \rho_2} \cdot \\frac{1}{\\beta\\sinh{\\beta}}+ \\frac{\rho_1\rho_2}{\rho_1 + \rho_2} \right]$ with $\beta = \left(\\frac{k+i\omega}{\omega_1}\right)^\\frac{1}{2}

Copy link

Choose a reason for hiding this comment

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

@etrevis Excellent work, could perhaps make it a custom circuit element?


where

Expand Down Expand Up @@ -377,7 +378,7 @@ def T(p, f):
else:
sinh.append(1e10)
Comment on lines 378 to 379
Copy link
Collaborator

Choose a reason for hiding this comment

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

I can't comment on the whole block but L374-379 make me a bit uneasy. L376 checks whether values of beta are larger than 100, but I believe it only checks the real part, even though beta will always be complex. Granted, the function only blows up when the real part is large, but I think we should make the check explicit.

Next, if beta.real is >100, we sub in 1e10 to avoid overflow errors. However, I think that fallback value is too low, and it's missing the imaginary component. For example, if we look at the overflow test on L89-92 of test_circuit_elements.py, we're using params=[1, 2, 50, 100] (I have no idea if those are reasonable or even physically possible).

If we use a frequency of 35 Hz, then beta = 74.15 +74.14j and np.sinh(beta) = (2.48e+31-7.58e+31j). If we instead use a frequency of 80 Hz then beta = 112.10+112.09j and np.sinh(beta) = 1e10. It seems weird that the output should be so much lower. I don't have a specific replacement number in mind, but I worry that 1e10 gets into the range of impedances people may actually be trying to fit?

Idk... pragmatically speaking 1e10 is 10GOhm so maybe it's not something to really worry about?

I also have problems with overflow in tanh that don't happen in the CI testing pipeline. I'm guessing that's because I'm on Windows, but the pipeline uses Unix? I've fixed this locally by replacing L374-379 with

    sinh, tanh = [], []
    for x in beta:
        # beta will be dtype np.complex128 by default, so real and imag components are each float64

        if x.real < 100:
            sinh.append(np.sinh(x))
            tanh.append(np.tanh(x))
        else:
            sinh.append(1e10)
            tanh.append(1 + 0j)


Z = A / (beta * np.tanh(beta)) + B / (beta * np.array(sinh))
Z = A / (beta * np.tanh(beta)) + B / (beta * np.array(sinh)) + B / 2
return Z


Expand Down
6 changes: 3 additions & 3 deletions impedance/tests/test_circuit_elements.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,9 +59,9 @@ def test_each_element():
(0.35594139 - 0.16491599j),
],
"T": [
(1.00041 - 0.00837309j),
(0.0156037 - 0.114062j),
(0.00141056 - 0.00141039j),
(1.10041251-0.00837309j),
(0.11560368-0.11406217j),
(0.10141056-0.00141039j),
],
"K": [
(0.099999842086579 - 0.000125663507704j),
Expand Down