Request: extend cubic EOS to support van der Waals

[adowling2]

@eslickj How difficult is it to extend the cubic EOS routines here to support van der Waals (u = w = 0)? We are happy to open a PR in IDAES to test this once the cubic root finder is updated.

@agarciadiego @bbefort

[eslickj]

@adowling2, probably not to hard, but I got to thinking the way we have it setup is probably a little stupid. We should probably just leave that part of the calculation in Pyomo for maximum flexibility. I kind of regretted the setup right away. Recently I added a plain generic cubic root function. You give it the coefficients of z^3 + bz^2 + cz + d = 0 as function arguments and it gives you the root. If the root you want doesn't exist it gives you the one that does. That's nice for the property calculations, but non-smooth. I think the current cubic setup in IDAES does the same though. We can think about providing a sort of psuedo-root the is smoothish. I can do that most likely.

Check here: https://github.com/IDAES/idaes-ext/blob/main/src/cubic/cubic_roots.c, the functions you want are cubic_root_h(b, c, d), cubic_root_m(b, c, d), and cubic_root_l(b, c, d). If I can get those adopted across the board, maybe we could get rid of the more complicated nonsense. @andrewlee94 any opinion?

[adowling2]

@eslickj Yes, we looked through the cubic_roots.c file and got confused about how to specify u and w. The general, EOS agnostic approach you suggested sounds best.

[andrewlee94]

The advantage of the current cubic EoS external function is that it only needs the A and B parameters from the EoS as inputs, and handles everything else internally. This keeps things fairly neat, at least on the Pyomo side as it only involves two arguments (rather than needing to create expressions for the additional derived coefficients. The downside of course is that the external function needs to know how to calculate the b, c and d coefficients itself.

[adowling2]

A middle ground would be to extend the current approach to allow the user to provide u and w. That would keep it EOS-centric but flexible.

[eslickj]

@adowling2, that's certainly easily doable. I guess to keep from breaking anything, I can add more functions. As much as I would like to cut this down, I know the clutter has to go somewhere. I'll add the functions ASAP to keep things moving, and we can think about cleaning up later.

[eslickj]

I just added a new EOS option, that set u = w = 0. So 0=PR, 1=SRK, and 2=VDW. I should have a new release out tonight. That's was the least invasive change. We can plan to move to replacing the EOS argument with u and v later. That will require some more extensive work and coordination, although the changes are pretty straight forward. We could just implement new functions too, but I'd like to try to simplify a bit instead of making it more complicated.

[adowling2]

Thank you! I like this plan. We can get this to work with the current (now slightly extended) interface. At some point in the future, we can clean things up.

@agarciadiego

[andrewlee94]

The easiest and cleanest approach is to probably deprecate the old functions, but and convert them to a wrapper that coverts the current "EoS type" argument to the u and w values and then calls the general function.

Then we can move over gradually without breaking anything.

[adowling2]

@andrewlee94 I like your proposal. This would allow someone to add a new EOS (change u and w) without changing the c code, correct? They would just need to modify Python code in IDAES (if I understand your proposal correctly.)

[andrewlee94]

Yes, this would let us do any cubic EoS by just providing the relevant u and w values.

[andrewlee94]

To take this a step further, that general function could itself be a thin wrapper over a general cubic root finder. All it would need to do is take the A, B, u and w inputs and calculate the relevant b, c and d coefficients (`a=1) of the actual cubic equation (to save the Pyomo user from needing to create those expressions).

[eslickj]

@andrewlee94, do I member right that we aren't using the extended smooth root stuff? That's one thing that will cut this down a bit. With your last comment, I'm not clear on what you are saying. Are you saying the C function takes A, B, u, and w or the C function takes b, c, d and in Pyomo we calculated those from A, B, u, and w (in some general wrapper)?

[adowling2]

Is the cubic EOS root-finding independent of the extended smoothed root formulation?

[andrewlee94]

1. I don't think we are using the smooth stuff - the SmoothVLE formulation does all of that on the Pyomo side.
2. I am proposing 2 c-functions. A general cubic root finder which takes 4 arguments (a, b, c and d) and finds the roots, and a wrapper over this specifically for cubic EoS's which takes A, B, u and w and then calculated the actual cubic coefficients from that using the standard forms.

[eslickj]

Ok, that's basically how it works, just with the EOS arg instead of u and v. It's a little less straight forward than move part of it to Pyomo since we need to handle derivatives. The advantage of putting the A and B calculation in Pyomo is that we don't need to go from derivatives in terms of b, c, d to derivatives in terms of A and B in the C function. It's not a big deal, but it does make the C code harder to follow.

[andrewlee94]

@eslickj I had not thought about the partial derivatives - that is a good argument for keeping as much on the Pyomo side as possible (i.e. only have the basic cubic root finder in C). We should be able to write the wrapper functions in Python anyway, as it is just making Pyomo expressions at that point.

[eslickj]

@andrewlee94, we already have functions that take b, c, d so if we go that route migration should be pretty easy. We can just eventually delete the functions we don't need. It's not a big deal, but a few years in the future some new person taking over will thank us when they don't have to look at that mess.

[eslickj]

New binaries are here https://github.com/IDAES/idaes-ext/releases/tag/test-release, except Windows, which will be there soon. I don't really have a test for the van der Waals thing, so let me know if you have any trouble. Once we get a chance to test everything, I'll create a new official release.

[ksbeattie]

@adowling2, this is in the latest idaes-ext release, but just needs testing to confirm it is correct. Can you identify someone to confirm/test it?

[adowling2]

@agarciadiego @bbefort Do we have a good test case for the new van der Waals EOS support?

[adowling2]

@eslickj Please ping us if/when the new Windows binaries are posted.

[eslickj]

@adowling2. I think it's already in the latest binaries.

[eslickj]

I think this is done.