[ANN] PyThermo.jl, an interface to the Thermo library for thermophysical properties

Package Announcements
1

PyThermo may not be the fastest thermophysical property library in the land, but thanks to the tireless efforts of Caleb Bell (the author of the underlying Python library, Thermo), it’s quite thorough.

Data is available for a vast number of chemical species, with Unitful accessor functions for the most common properties. Additional properties can be queried via the underlying PyObject, and should be assumed to have units of J, kg, m, s, K unless otherwise specified in Thermo’s documentation.

Water:

julia> H2O = Species("water", T = 40u"°C")
Species(water, 313.1 K, 1.013e+05 Pa)

julia> density(H2O)
997.0326202881853 kg m^-3

julia> H2O.calculate(T = 1000, P = 3000) # K, Pa

julia> density(H2O)
0.006500222660286359 kg m^-3

Air:

julia> air = Mixture(["N2" => 0.78, "O2" => 0.21, "Ar" => 0.01], P = 1u"atm")
Mixture({N2: 0.78, O2: 0.21, Ar: 0.01}, 298.1 K, 1.013e+05 Pa)

julia> R_specific(air)
287.0055796266994 J kg^-1 K^-1

Vodka:

julia> vodka = Mixture(["water", "ethanol"], Vfls = [.6, .4]) # liquid volume frac
Mixture({water: 0.83, ethanol: 0.17}, 298.1 K, 1.013e+05 Pa)

julia> vodka.Pc # critical pressure
1.933817214110451e7

1,1,1,3,5,5,5-heptamethyltrisiloxane:

julia> siloxane = Species("1,1,1,3,5,5,5-heptamethyltrisiloxane", T=500u"°F")
Species(1,1,1,3,5,5,5-heptamethyltrisiloxane, 533.1 K, 1.013e+05 Pa)

julia> soundspeed(siloxane)
142.92248228212017 m s^-1

julia> siloxane.Cp
1853.573089850276
13 Likes
2

Nice. How does Thermo compare to CoolProp?

3

As a big hairball of a Python library, it’s a few times slower than CoolProp.jl (which is just a wrapper around the CoolProp C++ library), but includes more properties & many more species (see here for CoolProp’s, or here for a subset of Thermo’s).

using PyThermo, CoolProp

function ρ_TP(species::Species, T, P)
    species.calculate(T = T, P = P)
    density(species)
end

function ρ_TP(cp_handle, T, P)
    CoolProp.AbstractState_update(cp_handle,CoolProp.get_input_pair_index("PT_INPUTS"), P, T)
    CoolProp.AbstractState_keyed_output(cp_handle,CoolProp.get_param_index("D"))
end

julia> thermo_H2O = Species("water");

julia> CP_H2O = CoolProp.AbstractState_factory("HEOS", "Water");

julia> @time [ρ_TP(thermo_H2O, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];
  3.598842 seconds (448.88 k allocations: 22.063 MiB, 0.38% gc time, 1.86% compilation time)

julia> @time [ρ_TP(CP_H2O, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];
  0.538447 seconds (115.41 k allocations: 4.366 MiB, 6.34% compilation time)

Thermo also allows the CoolProp Python binding to be used for certain properties:

julia> thermo_H2O.Psat
17217.65709618894

julia> thermo_H2O.VaporPressure.set_user_methods("COOLPROP")

julia> thermo_H2O.Psat
17213.15165587802
1 Like
4

Thanks. I will have to take a deeper look. With all the Symbolics discussion going on, it is beginning to feel that a pure julia solution would be nice.

Paulo

2 Likes
5

Yes, a similarly-featured Julia package would be great. There are a number of existing thermo-oriented Julia packages, but they’re all kinda pointing in different directions, with very different takes on API.

1 Like
6

hi, im the main developer of ThermoModels.jl. basically im porting all the relevant functionality to OpenSAFT.jl (on Juliacon 2021, chat in Zulip). What i learned abut implementing EoS in julia is the following:

  • AD is a game changer, thermodynamic models are basically combinations of partial derivatives, being able to automatically differenciate though models allows to use wild methods (for example, using a virial aproximation as a first guess for gas volume).
  • its necesary to stablish a standard on how to specify a thermodynamic state. how do we say to a model, i want to calculate this with fixed enthalpy and temperature, or fixed temperature and on a saturation line. one of my attempts to do this is in ThermoState.jl using julia dispatch to select the correct specification. (i used this on MoistAir.jl and WaterIF97.jl ,as a proof of concept on how such system could be implemented). decoupling the state from the model is important, as allows to calculate multiple states in parallel. For example, in a rachford rice iteration, the liquid and volume root finders could be executed at the same time, or in global optimization methods for multiphase equilibria, we could just throw a bunch of processes on a cluster to calculate everithing
  • with general derivatives, the only things necessary to define a EoS is the base function (namely, gibbs or helmholtz energy) and a minimun, lower bound volume. this lower bound volume is used to calculate liquid volumes. this is specially important when defining SAFT equations, as we dont have the rich amount of correlations that critical models have, the idea is that the function and the lower bound volume should be enough for defining models, not caring about starting points for iterations, as those could in theory be derived by the information stored in the function and its behaviour.
  • data is necessary, thats where the Thermo python package hits the jackpot, it doesnt matter if you have the most advanced equations, if you don’t have data to use those equations.

On a related theme, have you seen GitHub - CalebBell/chemicals: chemicals: Chemical database of Chemical Engineering Design Library (ChEDL) , its basically the database of models + helper functions, a good approach could be to use Thermo or Chemicals to query the results and julia to calculate properties. for example, this issue numba_vectorized reports TypeError · Issue #28 · CalebBell/chemicals · GitHub is basically solved by julia broadcasting

3 Likes
7

CoolProp seems to use the same general approach when calculating derivatives, Pure and Pseudo-Pure fluid properties — CoolProp 6.4.1 documentation, so in theory, we can just plug in the helmholtz model defined on coolprop and make it work on julia.

8

Yeah, I don’t think it makes much sense to try to duplicate Caleb’s massive property/coefficient databases in Julia - we can stick to nailing the computationally-intensive stuff and implementing clean, Julian interfaces. ThermoState seems like it could be a solid foundation for thermo efforts in Julia - have you done any benchmarks versus CoolProp?

1 Like
10

ThermoState.does not compete with CoolProp, instead proposes a common property interface api available to exchange between models. A dificulty on comparing packages is that we need the same equation of state to compare.a PC SAFT requires that requires an internal fixed point calculation on each property call, is about 100 times slower than a Peng Robinson eos, for example.

What could be done with a ThermoState object is to pass it between different models, for example, for speed, maybe a cubic suffices when simulating a compressor, but you need a more complete EoS when simulating a distillation tower. you just pass the corresponding state and the eos uses the appropiate equations to obtain all the needed properties.
the apropiate benchmark package is OpenSAFT, as it has component data and eos models. both packages have an IAPWS95 implementation, we can compare those as a starting point. from CoolProp docs, it seems that IAPWS95 is the standard when you calculate with water.

Pkg.add("OpenSAFT#database") #latest branch
using OpenSAFT
water = IAPWS95()
OpenSAFT.volume(water,101325,300)

on CoolProp (using CoolProp.jl)

Pkg.clone("https://github.com/vimalaad/CoolProp.jl.git")
Pkg.build("CoolProp") # to download the latest binaries
using CoolProp
PropsSI("D", "T", 300, "P", 101325, "Water")
11

reading more about coolprop, they use a concept named extended corresponding states (in their transport properties), where they use a cubic (SRK, but could be others) to obtain a reduced temperature and volume with respect to a reference fluid (in their case R134a) and then evaluate those reduced parameters on the multiparameter reference EOS. this is also known as SPUNG methodology: Evaluation of SPUNG* and Other Equations of State for Use in Carbon Capture and Storage Modelling - ScienceDirect . its a very handy technique to obtain more accurate properties than resorting to a cubic alone, specially on the 2-phase zone.
this methodology was recently incorporated in the package, via the SPUNG meta-eos, that takes a reference model and a shape (cubic at the moment) model.

1 Like
12

i was going to publish a benchmark inmediately after that comment, but i found my IAPWS95 could be optimized. after the rewrite i accidentally changed one sign of one constant, and it was a pain to find (IAPWS95 comes with some tests for implementations, really a life saver). after that those are the results, using the benchmark on volume calculation:

Pythermo (i dont know what eos uses here)

julia> @time [ρ_TP(thermo_H2O, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];
#4.271689 seconds (2.22 M allocations: 130.980 MiB, 1.10% gc time, 33.34% compilation time

julia> @benchmark [ρ_TP(thermo_H2O, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];

BenchmarkTools.Trial: 
  memory estimate:  14.73 MiB
  allocs estimate:  320005
  --------------
  minimum time:     2.903 s (0.00% GC)   
  median time:      2.994 s (0.00% GC)   
  mean time:        2.994 s (0.00% GC)   
  maximum time:     3.085 s (0.00% GC)   
  --------------
  samples:          2
  evals/sample:     1

Coolprop (IAPWS95)

julia> @time [ρ_TP(CP_H2O, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];
  0.538447 seconds (115.41 k allocations: 4.366 MiB, 6.34% compilation time)

julia> @benchmark [ρ_TP(CP_H2O, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)]
BenchmarkTools.Trial: 
  memory estimate:  859.66 KiB
  allocs estimate:  50005
  --------------
  minimum time:     696.510 ms (0.00% GC)
  median time:      700.885 ms (0.00% GC)
  mean time:        701.333 ms (0.00% GC)
  maximum time:     708.771 ms (0.00% GC)
  --------------
  samples:          8
  evals/sample:     1

OpenSAFT (IAPWS95):

using OpenSAFT #development branch
function ρ_TP(species::TT, T, P) where TT<:OpenSAFT.EoSModel
  v = OpenSAFT.volume(species,P,T)
  return 1/v
end

OpenSAFT_EOS = IAPWS95()

julia> @time [ρ_TP(w, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];
  0.830912 seconds (451.73 k allocations: 25.016 MiB, 1.16% gc time, 7.03% compilation time)
julia> @benchmark [ρ_TP(w, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)]
BenchmarkTools.Trial:
  memory estimate:  20.48 MiB
  allocs estimate:  368673
  --------------
  minimum time:     765.410 ms (0.00% GC)
  median time:      777.859 ms (0.00% GC)
  mean time:        779.380 ms (0.21% GC)
  maximum time:     796.346 ms (0.00% GC)
  --------------
  samples:          7
  evals/sample:     1

same equation, same results

julia> OpenSAFT.sat_pure(w,330) 
(17213.151668037026, 1.8294254890028448e-5, 0.15862466665274236)

so, near CoolProp. i think that, on one component, they just use a fitted vapor pressure curve and look if the pressure-temperature pair is above or below said line. whereas here there are 2 aditional helmholtz evaluations to select the correct phase. a good thing choosing a functional api, is that you can parallelize (from my laptop, 4 threads):

julia>  @time Threads.@threads for T in LinRange(280, 330, 100)           for P in LinRange(8e4, 12e5, 100)
               ρ_TP(w, T, P)
           end
       end
  0.312199 seconds (433.66 k allocations: 24.927 MiB, 16.72% gc time, 27.31% compilation time)

if we only need evaluations on multiple points, we can just throw a cluster and evaluate points simultaneously. if you know your compounds, you can use Cubics:

julia> cubic = SRK(["water"])
SRK{BasicIdeal} with 1 component:
 "water"
Contains parameters: a, b, acentricfactor, Tc, Pc, Mw

julia> @benchmark [ρ_TP(cubic, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)]
BenchmarkTools.Trial: 
  memory estimate:  18.08 MiB
  allocs estimate:  300005
  --------------
  minimum time:     19.961 ms (0.00% GC)
  median time:      21.105 ms (0.00% GC)
  mean time:        23.406 ms (11.20% GC)
  maximum time:     33.276 ms (22.45% GC)
  --------------
  samples:          214
  evals/sample:     1

(obviusly there is a speed improvement here, as volume roots of cubics have an analytical solution)

3 Likes