Minor improvements

[yoelcortes]

This pull includes two minor updates:

1. Remove print statement for failed PR routine.
2. Use multiplication instead of division. This prevents division by zero error for a very marginal case with unreasonable inputs (which ultimately leads to a different exception).

Although the changes are small, they will help get some cases passing in BioSTEAM.
Thank you!

[CalebBell]

Hello Yoel,
I'm very happy to remove the print statement. I didn't think it would be hit, it remained from when I was refining the cubic solver. It is a very intricate piece of code and something I might clean up in the future as it is rather messy.

I don't understand what the issue with the numerical check is. Can you provide a test EOS object that triggers this? I also don't the mathematial validity of your proposed formula update. It doesn't check for the same thing. The existing check checks that

a_alpha/(b*(b + delta) + epsilon) is very small relative to P ( ~< 1e.16 P); your proposed check is not the same.

Sincerely,
Caleb

[yoelcortes]

Hi Caleb,

Thanks so much for reviewing! I dug deeper into why I was getting the errors.

1. Regarding a_alpha/(b*(b + delta) + epsilon), you are right, my proposed check is not the same. I cannot believe my algebra was so off!! I was unable to reproduce the divide-by-zero error I was getting a couple of months ago, so let's just forget about this until I can get a test case.

2. Regarding the print statement, I found that I can reproduce it if P=numpy.nan, which is not that useful to test for at lower levels. I was solving for the bubble point pressure and extremely large liquid fugacity coefficients resulted in this "nan" value.

On a related issue, I believe the liquid fugacity coefficient estimated in the test case below might be too large (up to exp(4345) for some chemicals). I do not fully understand the topology of liquid fugacities using an equation of state and how such a large number is possible. I was wondering if I could consult with you on this matter? I am happy to open this as a separate issue on github.

Test case:

from thermo import PR78MIX
hydrocarbon_chemicals = ['H2', 'Methane', 'Ethane', 'Propane', 'n-Butane']
eos_mix = PR78MIX(
    Tcs=[33.145, 190.564, 305.322, 369.89, 425.125], 
    Pcs=[1296400.0, 4599200.0, 4872200.0, 4251200.0, 3796000.0], 
    omegas=[-0.219, 0.01142, 0.0995, 0.1521, 0.201], 
    kijs=[[0.0, -0.0044, -0.0781, -0.1311, -0.397], 
          [-0.0044, 0.0, -0.0059, 0.0119, 0.0185], 
          [-0.0781, -0.0059, 0.0, 0.0011, 0.0089], 
          [-0.1311, 0.0119, 0.0011, 0.0, 0.0033], 
          [-0.397, 0.0185, 0.0089, 0.0033, 0.0]], 
    zs=[0.015151515151515152, 0.07575757575757576,
        0.30303030303030304, 0.30303030303030304, 
        0.30303030303030304], 
    T=150, 
    P=75102339158.62407
)
print(f'{eos_mix.phi_l=}')
eos_mix.fugacities()
print(f'{eos_mix.lnphis_l=}')

Outputs:

eos_mix.phi_l=1e+308
eos_mix.lnphis_l=[995.5175372103799, 1608.9970827712168, 2430.9889188004086, 3375.668731908805, 4345.471474543158]

Thank you!

[yoelcortes]

Hi Caleb,

I found a test case which gives the division by zero error. It is a trivial case and I think it should be handled at a outer level on my side. It's when the EOS has no components:

eos_mix = PR78MIX(
    Tcs=[], 
    Pcs=[], 
    omegas=[],
    zs=[], 
    T=298.15, 
    P=101325
)

No need to change anything, just following up.

Thanks,

[CalebBell]

Hello Yoel,
Cubic EOSs are not universally applicable. I've put a lot of effort into making the behavior at extreme conditions mathematically correct even if the results are not likely to be correct. In this case the 75,000 MPa you specify pushes the cubic into a region it predicts a liquid density 3226 times more than an ideal gas. This is mathematically accurate to the equation - you can check that the relative error of the volume is the same as machine precision with eos_mix.volume_error().
You can also see the difference in results comes as the covolume b which is constant is approached by the volume that gets eos_mix.V_l/eos_mix.b increasingly near 1, here 1.0003. The log fugacity coefficients are correct - eos_mix.lnphis_l = [995.5175372103799, 1608.9970827712168, 2430.9889188004086, 3375.668731908805, 4345.471474543158] and correctly calculated according to the peng-robinson equation.
However, as you notice those terms cannot be converted into fugacity or log fugacity coefficients - they are simply larger than what a float can fit.

I set P_MAX_FIXED = 1e9 Pa in the flasher architecture to prevent flash calculations in this fantasy region. Other types of equations of state are applicable to them but not cubic equations of state. Thermodynamics math really pushes floating point math to its limits!

Sincerely,
Caleb

[yoelcortes]

@CalebBell, thanks for the valuable insight! Super useful. Please feel free to either merge the PR to remove the print statement or take care of later on your side.

Thank you!