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

[mdmurbach]

Looks like we're missing the final term in equation 22 from the Paasch paper. This adds in B/2 to both the element code, docstring, and tests.

[coveralls]

Pull Request Test Coverage Report for Build 4226976993

• 1 of 1 (100.0%) changed or added relevant line in 1 file are covered.
• No unchanged relevant lines lost coverage.
• Overall coverage remained the same at 100.0%

---

Totals
Change from base Build 4226528233:	0.0%
Covered Lines:	848
Relevant Lines:	848

---

💛 - Coveralls

[BGerwe]

Thanks for fixing these. I'll have time to review this and #255 in the next day or two.

[BGerwe]

The changes LGTM, but digging into this unearthed other changes I think we should make.

[BGerwe]

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

[etrevis]

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}

[zhuhaotian519]

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

[BGerwe]

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)