usnistgov / REFPROP-issues

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

REFPROP accuracy #187

Closed JHRhoral closed 4 years ago

JHRhoral commented 5 years ago

Hello,

I am developer. I checked the results of calculation reduced Helmholtz free energy and their derivatives according to ISO 20765-2:2015 (GERG-2008). I was checked results for all gases from page 57 in ISO 20765-2:2015 and from page 44 in AGA Report No. 8 Part 2:2017. Result I was checked with REFROP library and work written by Eric W. Lemmon Applied Chemicals and Materials Division National Institute of Standards and Technology (NIST) Boulder, Colorado, USA. version 2.0 of routines for the calculation of thermodynamic properties from the AGA 8 Part 2 GERG-2008 equation of state. April, 2017 - hereinafter NIST-GERG. Before calculation using REFPROP library it was necessary to add corresponding coefficients and constants from ISO 20765-2:2015 to *.FLD files of cause. Problem For example when calculating reduced free Helmholtz energy for gas Gas 3 arise significant differences in calculation with library REFPROP and calculation with work NIST-GERG. My results come out as well according to NIST-GERG. Because in REFPROP installation are source files in FORTRAN, I was debug it and I found that in line 128 in file MIX_HMX.FOR is used limit to the total number of adders to 55. Becouse in calculation sum have more that 56 adders last sums is not realised and the result from REFPROP library is other. Results for gasses Gulf coast, Amarillo, Ekofisk High N2, High CO2-N2, Gas 1 Gas 2 from AGA 8 Part 2 GERG-2008 is right.

REFPROP version 10 Windowes 8.1 64-bit CZE Access Method: FORTRAN original source files from installation and Intel Parallel Studion XE

ianhbell commented 5 years ago

I take it you did not call the GERG08 function: https://refprop-docs.readthedocs.io/en/latest/DLL/legacy.html#f/_/GERG08dll prior to calling your code? If you had done so, the results should have agreed with GERG-08 essentially perfectly.

REFPROP implements more accurate pure fluid models for the constituent fluids by default, but you can force the less-accurate GERG-2008 models by the appropriate functions (as indicated above).

There is also an option in the REFPROP GUI to enable the use of the GERG 2008 model.

JHRhoral commented 5 years ago

Hello, there is my code to tests and validation of my calculations. The REFPROP library is used as static library.

    program REFPROP_JHR_ST_Console
    implicit none

!DEC$ ATTRIBUTES DLLIMPORT :: OutputDebugString
!DEC$ ATTRIBUTES STDCALL, ALIAS : '_OutputDebugStringA@4' :: OUTPUTDEBUGSTRING

INTERFACE
    INTEGER*4 FUNCTION OUTPUTDEBUGSTRING (lpOutputString)
        !DEC$ ATTRIBUTES REFERENCE :: lpOutputString
        CHARACTER*(*)   lpOutputString
    END FUNCTION OUTPUTDEBUGSTRING
END INTERFACE

    ! Variables
    !character*10000 :: hFiles = 'METHANE.FLD;NITROGEN.FLD;CO2.FLD;ETHANE.FLD;PROPANE.FLD;BUTANE.FLD;ISOBUTAN.FLD;PENTANE.FLD;IPENTANE.FLD;HEXANE.FLD;HEPTANE.FLD;OCTANE.FLD'
    character*10000 :: hFiles = 'METHANE.FLD;NITROGEN.FLD;CO2.FLD;ETHANE.FLD;PROPANE.FLD;BUTANE.FLD;ISOBUTAN.FLD;PENTANE.FLD;IPENTANE.FLD;HEXANE.FLD;HEPTANE.FLD;OCTANE.FLD;OXYGEN.FLD;CO.FLD'
    integer :: ncomp = 14

    character*255 ::hFmix = 'HMX.BNC'
    character :: hrf*3 = 'DEF'
    integer :: ierr = 0
    character*255 :: herr
    integer:: itau = 0, idel = 0
    !double precision :: T = 180.0d+000, D = 19.677996111770913D+000, P = 10000.0d+000
    double precision :: T = 220.0d+000, D = 16.840970176229725D+000, P = 11000.0d+000

    double precision :: Tout = 0.0d+000, Pout = 0.0d+000
    double precision :: phi00 = 0.0d+000

    !character*10000 :: hMixNme = 'n:\REFPROP\Data\GAS 1 ISO.MIX'
    character*10000 :: hMixNme = 'n:\\HelmholtzFreeEnergy\\Knihovny\\GAS 3.MIX'    
    double precision :: z(20)
    double precision :: phi(20)

    character*255 :: hFlagGerg = 'GERG 2008'
    character*255 :: hFlagDebug = 'Debug'
    character*255 :: hFlagDir = 'Dir search'
    integer :: jFlag = 1
    integer:: kFlag = 0

    double precision ::  a = 0.0
    double precision ::  g = 0.0

    character :: hpth*255 = 'c:\Program Files (x86)\REFPROP\FLUIDS'

    integer :: iFlag = 1

    character :: htype*3 = 'EOS', hmix*3 = 'HMX', hcomp*60 = 'FEQ', hFiles2*1000, hFmix2*1000, hrf2*1000

    INTEGER*4 Debug, i1, i2, i3

    double precision :: Dl = 0.0d+000, Dv = 0.0d+000
    double precision :: x(20), y(20)
    double precision :: q= 0.0d+000,e= 0.0d+000,h= 0.0d+000,s= 0.0d+000,Cv= 0.0d+000,Cp= 0.0d+000,w= 0.0d+000
    double precision:: resDeriv = 0.0d+000, tau = 1.6069943273625296d+000, delta = 2.6261618807721607d+000;

    double precision PHI0, PHIHMX

    ! Body of REFPROP_JHR_ST_Console
    print *, 'Hello World'
    call FLAGS (hFlagGerg,jFlag,kFlag,ierr,herr)
    call FLAGS (hFlagDebug,jFlag,kFlag,ierr,herr)
    call FLAGS (hFlagDir,jFlag,kFlag,ierr,herr)    
    call SETPATH (hpth)          
    call GERG08(ncomp,iFlag,ierr,herr)    
    call SETMOD (ncomp,htype,hmix,hcomp,ierr,herr) 

      i1=INDEX(hFiles,CHAR(0))-1
      i2=INDEX(hFmix, CHAR(0))-1
      i3=INDEX(hrf,   CHAR(0))-1
      if (i1.lt.0) i1=LEN(hFiles)
      if (i2.lt.0) i2=LEN(hFmix)
      if (i3.lt.0) i3=LEN(hrf)
      hFiles2=hFiles(:i1)
      hFmix2 =hFmix (:i2)
      hrf2   =hrf   (:i3)

    call SETUP0 (ncomp,hFiles2,hFmix2,hrf2,ierr,herr)    
    call SETMIXTURE (hMixNme,z,ierr)    

    Debug = OUTPUTDEBUGSTRING('%%%')

    resDeriv = PHIHMX(itau,idel,tau,delta,z)

    !call AG (T,D,z,a,g)
    phi00 = PHI0(itau,idel,T,D,z)

    !call FUGCOF (T,D,z,phi,ierr,herr)

    !call TPFLSH (T,P,z,D,Dl,Dv,x,y,q,e,h,s,Cv,Cp,w,ierr,herr)

    !call ABFLSH ('TP',T,P,z,iFlag,Tout,Pout,D,Dl,Dv,x,y,q,e,h,s,Cv,Cp,w,ierr,herr)

    end program REFPROP_JHR_ST_Console
JHRhoral commented 5 years ago

Hello, I tried increase nmsav parameter to 58 and result of calculation from chenged REFPROP library and from routines preparedby NIST according to AGA 8 Part 2 GERG-2008 equation of state April, 2017 (NIST-GERG) is the same now.

parameter (nmsav=58) !Number of mixture models allowed (KW0-KW9, KWA-KWZ, etc.)

ianhbell commented 5 years ago

Interesting.... If you don't make that change, what happens?

Did you compile REFPROP as a static library using my REFPROP-cmake package: https://github.com/usnistgov/REFPROP-cmake ?

JHRhoral commented 5 years ago

What happens? The absolute value of the difference between the Helmholtz energy calculation by the REFPROP library and the work (NIST-GERG) is about 1.0e-5. Increasing nmsav to 58 is less than 1.0e-10. The absolute value of the difference between the calculation of the first pressure derivative with respect to density by the REFPROP library and the work (NIST-GERG) is about 1.0e-2. Increasing nmsav to 58 is less than 1.0e-8.

The static library used for the tests was created in Visual Studio 2019 via Intel Parallel Studio XE. I created a dynamic library using the same tools, I ran the same tests by program in C++ with it, and the test results are the same as using the static library.

With regards Ing. Jan Horník

ianhbell commented 5 years ago

Seems like this might be a bug. I've passed this on to @EricLemmon to investigate

EricLemmon commented 5 years ago

I am very impressed with the depth you searched into the code to find this bug! It was not obvious and the change in calculated values was small, something that most people would have ignored. It took me quite a bit of time to verify that the error was real, find a good solution to fix the problem, and obtain identical numbers between the Refprop code and the AGA8 code. I appreciate the work you put into this and for reporting it.

The bug is fixed now. I would suggest you change the 58 to 100 because I noticed that even 58 could be exceeded. I am currently working on a number of other bug fixes and do not have a version that I am comfortable with passing on, but please feel free to email me directly from time to time to see if the new code is ready: Eric.Lemmon@nist.gov