search

usnistgov / REFPROP-wrappers

Wrappers around NIST REFPROP for languages such as Python, MATLAB, etc.
195 stars 127 forks source link

Calling the same routine again shows no convergence #595

Closed smoe38 closed 6 months ago

smoe38 commented 6 months ago

Hello,

I am successfully using my Python program to calculate processes with different medium compositions. Now I have encountered a strange behaviour and I cannot find a solution. See sample code below:

RP, z = set_fluid(medium_components, medium_composition)

enthalpy = some_function(RP, z, pressure, temperature)

iteration = iteration_function(RP, z, pressure, temperature)

Description

In the set_fluid function I define the medium with RP = REFPROPFunctionLibrary(environ['RPPREFIX']) and so on. In the second function I calculate some single values like enthalpy, entropy, ... In the third function I calculate a process iteratively.

So this works really well and I have calculated some single component fluids and also some mixtures with 5 components in the 2-phase region. Now I have tried another mixture of ["CO2", "H2O", "O2", "N2", "Ar", "He"] and I encounter the following behaviour:

But when I do this:

RP, z = set_fluid(medium_components, medium_composition)
y = RP.REFPROPdll("", "PT", "H", 21, 0, 0, p_in, T_in, z)
enthalpy = some_function(RP, z, pressure, temperature)

iteration = iteration_function(RP, z, pressure, temperature)

I get a value for enthalpy and it seems plausible. But I cannot get it to work so that this also works in the iterative function as well: 

RP, z = set_fluid(medium_components, medium_composition)
y = RP.REFPROPdll("", "PT", "H", 21, 0, 0, p_in, T_in, z)
enthalpy = some_function(RP, z, pressure, temperature)
y = RP.REFPROPdll("", "PT", "H", 21, 0, 0, p_in, T_in, z)
iteration = iteration_function(RP, z, pressure, temperature)

When I call y = RP.REFPROPdll("", "PT", "H", 21, 0, 0, p_in, T_in, z) the second time after enthalpy, I get an error again. So for the same input properties I get a value at first but on the second call I get an error.

I do not think this is due to the mixture rather than a 2-phase calculation. It seems that the dll is reset in a strange way. By the way, I use iMass = 0 and iFlag = 1 to call SATSPLN. But I also tried with iFlag = 0 and it does not change.

The tip for calling Refpropdll before the function I found in issue #461.

Is there any fix or workaround to "activate" the 2-phase-calculation in the whole scripts?

Thank you in advance.

Versions

REFPROP Version: v10.0
Operating System and Version: Windows 10 64 bit Access Method: Python

ianhbell commented 6 months ago

Can you provide a complete, runnable example where you experience this problem? Otherwise I don't see anything obvious in your code

smoe38 commented 6 months ago

Ok, so I could reproduce it a bit easier, I think. See example below:

def set_fluid(medium_components, medium_composition):
    # Import necessary functions
    from os import environ
    from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary

    # Calculate composition if given composition != 1.0
    z = []
    if sum(medium_composition) != 1.0:
        for i in range(len(medium_composition)):
            z.append(medium_composition[i] / sum(medium_composition))
    else:
        z = medium_composition

    # Refprop: Load function library
    RP = REFPROPFunctionLibrary(environ['RPPREFIX'])
    RP.SETPATHdll(environ['RPPREFIX'])
    # Set active fluid
    hFld = str(';'.join(medium_components))
    RP.SETFLUIDSdll(hFld)
    return RP, z

def calc_crit_point(RP, z):

    # Set unit base for Refprop
    MASS_BASE_SI = RP.GETENUMdll(0, "MASS BASE SI").iEnum

    # Calculate crit point and add to dict
    crit_point = {}
    for k in ["TC", "PC", "DC"]:
        x = RP.REFPROPdll("", "", k, MASS_BASE_SI, 0, 1, 0, 0, z)
        crit_point.update({k + ' [' + x.hUnits+']:': x.Output[0]})

    return crit_point

# Process data
p_in = 86 * 1e5 #Pa
T_in = 20 + 273.15 #K

# Medium components and composition
x_CO2 = 0.771
x_H2O = 0.1594
x_O2 = 0.0154
x_N2 = 0.0381
x_Ar = 0.016
x_He = 0.0001

medium_components = ["CO2", "H2O", "O2", "N2", "Ar", "He"]
medium_composition = [x_CO2, x_H2O, x_O2, x_N2, x_Ar, x_He]

# Function to set fluid
RP, z = set_fluid(medium_components, medium_composition)

# Part1
rho_in = RP.REFPROPdll("", "PT", "D", 21, 0, 0, p_in, T_in, z)
print(rho_in.Output[0], rho_in.herr)

# Part2
crit_point = calc_crit_point(RP, z)

# Part3
rho_in = RP.REFPROPdll("", "PT", "D", 21, 0, 0, p_in, T_in, z)
print(rho_in.Output[0], rho_in.herr)

So to reproduce my problem: With the code above I get two result prints. When only Part1 is commented I get -9999990.0 [TPFL2 error 226] 2-phase iteration did not converge.. When Part2 is commented I get 2 results prints.

When you try this out with another composition:

# Medium components and composition
x_CO2 = 0.913
x_H2O = 0.004
x_O2 = 0.018
x_N2 = 0.045
x_Ar = 0.019

medium_components = ["CO2", "H2O", "O2", "N2", "Ar"]
medium_composition = [x_CO2, x_H2O, x_O2, x_N2, x_Ar]

You will only get a results when nothing is in comments.

ianhbell commented 6 months ago

Ah this is a subtle one but I know what is going on. The call to critical point routines internally calculates the saturation spline, and that spline should be also on for the next call, even though you didn't request it. Here is code that demonstrates that is what is going on:

def set_fluid(medium_components, medium_composition):
    # Import necessary functions
    from os import environ
    from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary

    # Calculate composition if given composition != 1.0
    z = []
    if sum(medium_composition) != 1.0:
        for i in range(len(medium_composition)):
            z.append(medium_composition[i] / sum(medium_composition))
    else:
        z = medium_composition

    # Refprop: Load function library
    RP = REFPROPFunctionLibrary(environ['RPPREFIX'])
    RP.SETPATHdll(environ['RPPREFIX'])
    # Set active fluid
    hFld = str(';'.join(medium_components))
    RP.SETFLUIDSdll(hFld)
    return RP, z

def calc_crit_point(RP, z):

    # Set unit base for Refprop
    MASS_BASE_SI = RP.GETENUMdll(0, "MASS BASE SI").iEnum

    # Calculate crit point and add to dict
    crit_point = {}
    for k in ["TC", "PC", "DC"]:
        x = RP.REFPROPdll("", "", k, MASS_BASE_SI, 0, 1, 0, 0, z)
        crit_point.update({k + ' [' + x.hUnits+']:': x.Output[0]})

    return crit_point

# Process data
p_in = 86 * 1e5 #Pa
T_in = 20 + 273.15 #K

# Medium components and composition
x_CO2 = 0.771
x_H2O = 0.1594
x_O2 = 0.0154
x_N2 = 0.0381
x_Ar = 0.016
x_He = 0.0001

medium_components = ["CO2", "H2O", "O2", "N2", "Ar", "He"]
medium_composition = [x_CO2, x_H2O, x_O2, x_N2, x_Ar, x_He]

# Function to set fluid
RP, z = set_fluid(medium_components, medium_composition)

# Part1
rho_in = RP.REFPROPdll("", "PT", "D", 21, 0, 0, p_in, T_in, z)
print(rho_in.Output[0], rho_in.herr)

# Part2
crit_point = calc_crit_point(RP, z)

# This line wipes the saturation spline
RP.SATSPLNdll([0]*5)

# Now we programatically extract the nodes of the isopleth following https://nbviewer.org/github/usnistgov/REFPROP-wrappers/blob/master/wrappers/python/notebooks/Tutorial.ipynb
Nnodes = RP.REFPROPdll("","SPLINENODES","",0,0,0,0,0,z).iUCode
data = []
for i in range(1,Nnodes+1):
    rho_vap,T,p,rho_eq,z0,z1 = RP.REFPROPdll("","SPLINENODES","",0,0,i,0,0,z).Output[0:6]
    data.append(dict(rho_vap = rho_vap,T=T,p=p,z0=z0,z1=z1))
print(len(data), 'spline nodes are present currently')

# Part3
rho_in = RP.REFPROPdll("", "PT", "D", 21, 0, 0, p_in, T_in, z)
print(rho_in.Output[0], rho_in.herr)

If you comment out the call to SATSPLNdll then you get the old behavior, but if you don't you get the second call to always work.

The saturation spline usually helps, but not always, and this is an example where it makes things worse.

A related comment is that REFPROP has no notion of multiple critical points or solid phases, and you might be running into that in this application.

smoe38 commented 6 months ago

Thank you Ian for your quick responses. Unfortunately nothing changes when I comment the call to SATSPLNdll, I get the exact same output:

652.9213042437344 
319 spline nodes are present currently
652.9213042437344

Did I understand correctly, that I should "empty" the splines after I calculate the critical point? If I call later with iFlag = 1? In my program sometimes it works better with iFlag = 0, sometimes with iFlag = 1. I think this depends on the mixture and composition (mixing rules etc.)?

Do you have any tip for calculation in 2-phase-region?

ianhbell commented 6 months ago

That is correct - the call to SATSPLNdll is exactly the fix, and you uncommented the fix :)

In this case, yes. The saturation splines are a mixed blessing. Sometimes they help, sometimes not. I have found I need to experiment. No, I do not have any further recommendations about two-phase inputs. You might have to do your own iteration on your side. Better behaved are usually saturation calls, but I know you might not be at saturation, so that's not very helpful guidance.

smoe38 commented 6 months ago

If the call to SATSPLNdll is active or not does not change the output for me, it is exactly the same. There are still 319 nodes present.

ianhbell commented 6 months ago

What precise version of REFPROP do you have? Obtained via the RPVersion method?

ianhbell commented 6 months ago

Does the call to SATSPLNdll give you an error message on your side if you print the returned value?

smoe38 commented 6 months ago

I have Version 10.0 I get this output: SATSPLNdlloutput(ierr=312, herr='[SATSPLN error 312] sum(z)=0')

ianhbell commented 6 months ago

Sorry it should have been:

# This line wipes the saturation spline
print(RP.FLAGSdll("Splines off", 1))

see https://refprop-docs.readthedocs.io/en/latest/DLL/high_level.html#f/_/FLAGSdll

smoe38 commented 6 months ago

Thank you, this works.

ianhbell commented 6 months ago

Closed because I don't have any magic bullets or workarounds to offer