Generic property package speed of sound

[eslickj]

We need speed of sound out of the generic property package. Not sure who will end up doing it, but this may save someone some time thinking about partial derivatives and https://www.osti.gov/servlets/purl/809947.

[andrewlee94]

Another option that has been brought up on a couple of equations would be to create an AD tool that can automatically generate these derivatives for us. As all thermophysical properties are interrelated, it is at least theoretically possible to write a model for a single thermophysical property (generally Gibbs energy, Helmholtz energy or specific volume), and then all other properties from there by differentiation.

If we could develop a tool that does that (semi-)automatically, that would be a huge advantage over other tools.

[eslickj]

@andrewlee94, Pyomo does have the Sympy module. We may be able to do something with that. I guess if we had some sort of AD tool, it would have to write Pyomo expressions, so I'm not sure if there is a good off-the-self solution. If we work out the derivatives manually you can do them implicitly so you don't have to solve anything, which is tedious but not hard. With AD I'm having trouble connecting the dots, but you're right it does seem possible and a great idea if we can make it work, especially it a full set of any derivative you may want just pop out like magic.

@jsiirola or @carldlaird do you guys have any thoughts on automatic calculation of thermodynamic partial derivatives? Like for speed of sound we need the partial derivative of P with respect to rho at constant T and partial P with respect to T at constant rho. In the case I care about both can be derived directly from a cubic equation of state, but that's not always the case, so the path to derivative may be more convoluted.

[dallan-keylogic]

@eslickj Based on my previous experience with Sympy, it wouldn't be my first choice for actual calculation of derivatives. A real AD tool evaluates derivatives much faster and scale much better. On the level of a thermo package, though, it might not be horrible. I've previously had some experience using CasADi for AD. It looks like the Pyomo DAE simulator has some methods for converting between Pyomo symbolic expressions and CasADi symbolic expressions.

@andrewlee94 @eslickj I am much more skeptical that taking derivatives is going to be a good way of deriving thermodynamic quantities, though. Maybe for a full equation of state things would work out. However, differentiating a signal makes it noisier, and just because one function is a good approximation for another doesn't mean its derivative would be a good approximation of the other's derivative. Things could get pretty hairy, I imagine, for, e.g., polynomial correlations giving the heat capacity.

[andrewlee94]

@dallan-keylogic Using AD/SD tools would be up to the model developer, and thus would not be a global solution but rather one for those cases where you know you have a rigorous thermo model behind everything.

The other thing to note is that because of the way Pyomo works, this evaluation would either only occur once per property during model construction, or during the model write stage - it is not something that would have to be called at each solver iteration (i.e. we are building expressions, which only need to be done once, not actual function evaluations).

[blnicho]

Pyomo does have a couple tools for taking derivatives of Pyomo expressions:

code: https://github.com/Pyomo/pyomo/blob/main/pyomo/core/expr/calculus/derivatives.py

examples: https://github.com/Pyomo/pyomo/blob/main/pyomo/core/tests/unit/test_symbolic.py or https://github.com/Pyomo/pyomo/blob/main/pyomo/core/tests/unit/test_derivs.py

These tools might not be the most efficient but it should be straightforward to test if the overhead of using them is worth not taking the derivatives by hand.

[eslickj]

This is happening in PR #460?

[andrewlee94]

I think PR #460 adds the basic hooks and the equation for the cubic case, but I think the other property packages would need to have equations added as well.

[andrewlee94]

I think the basic implementation of this is in place. We can address the need for equations in other EoS modules on a case-by-case basis.

[dallan-keylogic]

Just a note here in case it ever comes up again. I wanted to validate my implementation of partial molar quantities for ceos.py against derivatives calculated numerically. I found an equation for a partial molar quantity in terms of derivatives with respect to mole fractions in my thermo textbook.

def partial_molar(props,expr,comp_k):
        out = expr
        for comp_j in  props.component_list:
            if comp_j == comp_k:
                out += differentiate(expr,
                         wrt=props.mole_frac_phase_comp["Vap",comp_j])
            out -= (props.mole_frac_comp[comp_j]
                    * differentiate(expr,
                        wrt=props.mole_frac_phase_comp["Vap",comp_j]))
        return value(out)

However, in order for this method to work, the functions need to be explicit functions of all the mole fractions. If we use log mole fraction constraints of the form exp(log_mole_frac_comp[j]) == mole_frac_comp[j], the differentiation will not give a sensible answer. That rules out Gibbs energy and entropy. For enthalpy, I found a more subtle issue--- in ceos.py, the compressability factor Z is an implicit function of mole fraction, so the answers were wrong.

Future reader be warned!