usnistgov / REFPROP-issues

A repository solely used for reporting issues with NIST REFPROP
26 stars 13 forks source link

Inconsistency between GUI and API for R454C #501

Open crlaugh opened 2 years ago

crlaugh commented 2 years ago

When calculating the specific enthalpies along the saturation lines for R454C.MIX, I again seem to observe some inconsistencies between the low-level API and the GUI. Based on the issue the other day, I am setting up the splines, as can be seen in the following code.

I am using PyCall in Julia to load ctREFPROP; my code is as follows:

# Initialization
using Revise
using PyCall

ctRP = pyimport("ctREFPROP.ctREFPROP")
RP = ctRP.REFPROPFunctionLibrary(ENV["RPPREFIX"])
RP.SETPATHdll(ENV["RPPREFIX"])

PVal = 3640

# Using low-level interface
resMix = RP.SETMIXdll("R454C.MIX", "HMX.BNC", "DEF")
zVal = resBub1[1]
msg = RP.SATSPLNdll(zVal)
resBub = RP.SATPdll(PVal, zVal, 1)
hBubVal = RP.THERM2dll(resBub[1], resBub[2], zVal)[3]
resDew = RP.SATPdll(PVal, zVal, 2)
hDewVal = RP.THERM2dll(resDew[1], resDew[3], zVal)[3]

# Sweeping through the two-phase region
xVec = -0.1:0.005:1.1
hVec = hBubVal .+ xVec.*(hDewVal - hBubVal)

# Legacy API
dVec = zeros(size(hVec))
for (ih, hVal) in enumerate(hVec)
    dVec[ih] = RP.PHFLSHdll(PVal, hVal, zVal)[2]
end 

This is visible in the following plot: image

I can see the same behavior when using the high-level API, though the outputs produced at the erroneous points are -9.99e6. This makes some sense to me, and it appears that the iterations don't converge at these points - the ierr values at these points are nonzero.

However, when I enter the same points into the GUI, I find much smoother behavior - even if I directly specify some of the erroneous points.
image

Is there anything else that I should be doing to get smoother behavior out of the REFPROP API calls?

I am using REFPROP 10.0.0.97 on Win10.

ianhbell commented 2 years ago

Where is that pressure on the phase diagram? If near the pmax of the isopleth, it wouldn't surprise me you are running into trouble, mixtures are hard. Do all the saturated calls work?

crlaugh commented 2 years ago

The critical pressure is 4.318 MPa, so it's not all that close to the critical point - I'm not sure how to calculate PMax of the isopleth. The saturated calls (with SATP) do work at this pressure.

crlaugh commented 2 years ago

I totally get that mixtures are hard, so it wouldn't surprise me if I just can't calculate the density at these particular state points. What I don't understand is why it seems like the GUI provides much smoother output. Does the GUI perturb the (P,h) values if the underlying flash calculations don't converge?

ianhbell commented 2 years ago

No, but the composition is changing as you move through the two-phase region, so the GUI might call SATSPLNdll at intermediate compositions

@EricLemmon , any idea?

ianhbell commented 2 years ago

Are the densities coming out the same between the GUI and the DLL?

crlaugh commented 2 years ago

They are not - good catch. For (PVal, hVal) = (3640 kPa, 29781 J/mol), RP.REFPROPdll("R454C.MIX", "PH", "D", DEFAULT, 0, 1, PVal, round(hVec[1]), [])[2][1] yields 8.489983168451976.

However the GUI outputs 8.0379 mol/L: image

I get the same densities from the high-level API and the legacy API.

ianhbell commented 2 years ago

Reference state then. Make sure you call SETREFdll to specify the desired reference state. That should have no impact on the iterations, but ???

crlaugh commented 2 years ago

How should I call SETREFdll to get the same references as the GUI?

ianhbell commented 2 years ago

The code in the XLA gives a clue:

        'Use the following line to calculate enthalpies and entropies with a reference state based on the currently defined mixture (and thus calculations will be the same as those produced from the graphical interface for mixtures).
        'See the file manual.txt in your REFPROP\Fortran directory for further information on the use of the SETREF subroutine.
           i = 2 'Set to 1 or 2 (see the SETREF routine in the file SETUP.FOR) or the information on the "Welcome" sheet.
           hFld = "": hIn = "Flags": hOut = "SETREF": Call REFPROPdll(hFld, hIn, hOut, 0&, 0&, i, 0#, 0#, z(1), Output(1), hUnits, iUnits, x(1), y(1), x3(1), q, ierr, herr, 10000&, 255&, 255&, 255&, 255&)
crlaugh commented 2 years ago

Thanks - I'll try this.

crlaugh commented 2 years ago

Where is this code from in the REFPROP distribution?

ianhbell commented 2 years ago

Which piece?

crlaugh commented 2 years ago

The XLA code.

ianhbell commented 2 years ago

https://github.com/usnistgov/REFPROP-wrappers/tree/master/wrappers/Excel

crlaugh commented 2 years ago

In looking at this code, it's not clear how to do this - it looks like there are a bunch of placeholders in this code, and I don't quite see how to map this to the API. @EricLemmon, can you provide some guidance?

ianhbell commented 2 years ago

Here is one example: https://nbviewer.org/github/usnistgov/REFPROP-wrappers/blob/master/wrappers/python/notebooks/Tutorial.ipynb#Setting-Up-Predefined-Mixture

ianhbell commented 2 years ago

The same thing can be done via REFPROPdll function, but if you are in Julia/Python, then I would just use the SETREFdll function directly.

crlaugh commented 2 years ago

I set up SETMIXdll using resMix = RP.SETMIXdll("R454C.MIX", "HMX.BNC", "DEF") and reran the calculation. As far as I can tell, there are significant differences in the densities between the GUI and the API. I have used both the high-level and legacy APIs to get the same output:

PVal = 3640
hVal = 31448

resMix = RP.SETMIXdll("R454C.MIX", "HMX.BNC", "DEF")
zVal = resMix[3]
msg = RP.SATSPLNdll(zVal)

RP.REFPROPdll("R454C.MIX", "PH", "D", DEFAULT, 0, 1, PVal, hVal, [])[2][1]
# 6.201493346592637

RP.PHFLSHdll(PVal, hVal, zVal)[2]
# 6.201493346592637

If I use the GUI at this same state point, I get the following:

image

Of course, I can believe the composition is different at this point in the two-phase region, but it's not clear to me how to calculate that for specification to PHFLSHdll or REFPROPdll.

Since that point is in the two-phase region, I also calculated the density in the liquid region, and also got discrepancies:

PVal = 3640
hVal = 29781

RP.REFPROPdll("R454C.MIX", "PH", "D", DEFAULT, 0, 1, PVal, hVal, [])[2][1]
# 8.489983168451976

RP.PHFLSHdll(PVal, hVal, zVal)[2]
# 8.489923355503572

whereas the GUI generates

image

?!

ianhbell commented 2 years ago

I don't see the call to SETREFdll with the mixture flags?

crlaugh commented 2 years ago

Thanks for your patience - I had used SETMIXdll, but not SETREFdll. Now that I used both of them, the densities between the GUI and both of the APIs match.

resMix = RP.SETMIXdll("R454C.MIX", "HMX.BNC", "DEF")
resREF = RP.SETREFdll("DEF", 2, resMix[3], 0, 0, 0, 0)
zVal = resMix[3]
msg = RP.SATSPLNdll(zVal)

Unfortunately, both REFPROPdll and PHFLSHdll still have nonconvergence issues that I don't seem to see in the GUI. Any other thoughts?

image

crlaugh commented 2 years ago

Actually, one of the points that I tried (P,h = 3640 kPa, 36452.21 J/mol) did not converge, but when I entered a second point (P,h = 3640 kPa, 36294.97 J/mol), the REFPROP GUI changed the input point to 36291 and converged. So presumably the GUI cannot calculate either of these state points, but there's some logic that makes it shift the input a bit so that it can return an answer.

crlaugh commented 2 years ago

I think that this answers my question - now the GUI and the API are providing consistent results, and just not converging at the same set of points because mixture calculations are hard.

EricLemmon commented 2 years ago

Which version of the GUI are you using?

EricLemmon commented 2 years ago

Sorry, I meant which version of the DLL?

crlaugh commented 2 years ago

10.0.0.97

EricLemmon commented 2 years ago

Ian, can you send him the latest one? I don't think it will fix the non-convergent issues, but let's see.

crlaugh commented 2 years ago

@ianhbell Do you have my personal email to try out the latest beta? Let me know if not.

ianhbell commented 2 years ago

I don't think I have a clear way of getting you the code for the latest one. @EricLemmon and I will have to iterate a bit on our side.