lanl / singularity-eos

Performance portable equations of state and mixed cell closures
https://lanl.github.io/singularity-eos/
BSD 3-Clause "New" or "Revised" License
27 stars 8 forks source link

Fix bulk modulus calculation in spiner EOS #370

Closed pdmullen closed 6 months ago

pdmullen commented 6 months ago

Fixes a (potential) bug in Spiner EOS bulk modulus calculation.

PR Summary

Despite the bulk modulus being a required table in a SpinerDependsRhoT or SpinerDependsRhoSie EOS file, it is overwritten by functions fixBulkModulus_ and calcBMod_, respectively. Inside these functions, the bulk modulus formulae assume dpdrho to be $\partial P / \partial \rho \; \vert{T}$. However, I believe (please check me) that spiner expects that the tables provided correspond to $\partial P / \partial \rho \; \vert{e}$. This PR changes the bulk modulus formula to reflect this. There is this trick that eos_spiner plays wherein if $\partial P / \partial e \; \vert{\rho} < 0$, it modifies the bulk modulus calculation. I tried to preserve this behavior by defining $\partial P / \partial \rho \; \vert{T}$ in terms of thermodynamic derivatives (please double check me...).

PR Checklist

If preparing for a new release, in addition please check the following:

Yurlungur commented 6 months ago

Thanks for catching this @pdmullen ! I'm going to trigger tests on Darwin.

chadmeyer commented 6 months ago

Looking at this now. I think this might be OK, but I want to make sure (and perhaps also ensure that this is what we want).

Let's agree that $c^2 = \left.\frac{\partial P}{\partial \rho}\right|_{e}+\frac{P}{\rho^2}\left.\frac{\partial P}{\partial e}\right|_{\rho}$. Is that what is being done in both of these cases? I ask because one of these cases is DependsRhoT and the other one is DependsRhoSie, so I'm not sure what was tabulated (is it the same between them?).

Is that what it is? I'm thinking so, but I might need to dig into Sesame2Spiner to be sure.

Yurlungur commented 6 months ago

Is that what it is? I'm thinking so, but I might need to dig into Sesame2Spiner to be sure.

I had the same confusion. I thought Sesame2Spiner was tabulating $dP/d\rho |_T$ but it turns out it's tabulating $dP/d\rho |_e$.

jhp-lanl commented 6 months ago

Is that what it is? I'm thinking so, but I might need to dig into Sesame2Spiner to be sure.

I had the same confusion. I thought Sesame2Spiner was tabulating dP/dρ|T but it turns out it's tabulating dP/dρ|e.

I think the another issue is that Sesame2Spiner gets density-energy derivatives for density-temperature tables via thermodynamic identities and then uses another thermodynamic identity on them to get the bulk modulus. I would propose (for the future) to just ask for the bulk modulus directly via the EOS_BSt_DT table.

This would be advantageous for two reasons: 1) it avoids roundoff due to unnecessary multiplications and additions and 2) it can leverages the effort that the EOSPAC folks put into making sure that EOSPAC returns the most accurate bulk modulus possible.

pdmullen commented 6 months ago

Is that what it is? I'm thinking so, but I might need to dig into Sesame2Spiner to be sure.

I had the same confusion. I thought Sesame2Spiner was tabulating dP/dρ|T but it turns out it's tabulating dP/dρ|e.

I think the another issue is that Sesame2Spiner gets density-energy derivatives for density-temperature tables via thermodynamic identities and then uses another thermodynamic identity on them to get the bulk modulus. I would propose (for the future) to just ask for the bulk modulus directly via the EOS_BSt_DT table.

This would be advantageous for two reasons: 1) it avoids roundoff due to unnecessary multiplications and additions and 2) it can leverages the effort that the EOSPAC folks put into making sure that EOSPAC returns the most accurate bulk modulus possible.

I'd add that it feels a bit weird for SpinerEOS to demand a bulk modulus field in tables when it actually just recomputes them from thermodynamic derivatives. Maybe a future MR could optionally enroll "bulk modulus correction/fixing". That way if a user wants to use bulk moduli as directly computed from EOSPAC, they could set correct_bulk_mod = false (or similar).

jhp-lanl commented 6 months ago

I think this looks ready for merging

Yurlungur commented 6 months ago

@jhp-lanl @pdmullen agreed regarding using the bulk modulus from eospac. I'll make an issue.