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

**URL:** <https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310>\
**Category:** Package Announcements\
**Created:** [March 31, 2021, 4:00pm UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310 "2021-03-31T16:00:13Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [March 31, 2021, 4:00pm UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310/1 "2021-03-31T16:00:13Z")

</div>

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](https://github.com/CalebBell/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](https://thermo.readthedocs.io/thermo.chemical.html).

Water:

```julia
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
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
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
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

```

> **[GitHub - stillyslalom/PyThermo.jl: Julia interface to Thermo...](https://github.com/stillyslalom/PyThermo.jl)**
>
> Julia interface to Thermo (github.com/CalebBell/thermo), a Python library for thermodynamics - GitHub - stillyslalom/PyThermo.jl: Julia interface to Thermo (github.com/CalebBell/thermo), a Python l...

---

<div class="post-metadata">

**Author:** ![Paulo\_Jabardo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paulo_jabardo/32/3196_2.png) [@Paulo\_Jabardo](https://discourse.julialang.org/u/Paulo_Jabardo)\
**Post date:** [March 31, 2021, 5:04pm UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310/2 "2021-03-31T17:04:02Z")

</div>

Nice. How does Thermo compare to CoolProp?

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [March 31, 2021, 6:11pm UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310/3 "2021-03-31T18:11:44Z")

</div>

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](http://coolprop.sourceforge.net/fluid_properties/PurePseudoPure.html#list-of-fluids) for CoolProp’s, or [here](https://github.com/CalebBell/thermo/blob/master/thermo/Critical%20Properties/Yaws%20Collection.tsv) for a subset of Thermo’s).

```julia
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
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
julia> thermo_H2O.Psat
17217.65709618894

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

julia> thermo_H2O.Psat
17213.15165587802

```

---

<div class="post-metadata">

**Author:** ![Paulo\_Jabardo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paulo_jabardo/32/3196_2.png) [@Paulo\_Jabardo](https://discourse.julialang.org/u/Paulo_Jabardo)\
**Post date:** [March 31, 2021, 8:20pm UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310/4 "2021-03-31T20:20:58Z")

</div>

> [@stillyslalom](#):
>
> many

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

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [March 31, 2021, 10:04pm UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310/5 "2021-03-31T22:04:56Z")

</div>

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.

- [Thermodynamics.jl](https://github.com/CliMA/Thermodynamics.jl): properties of water, ice, and moist air for climate modeling
- [Psychro.jl](https://github.com/termolivre/Psychro.jl): properties of moist air
- [Thermodynamics.jl](https://github.com/doomphoenix-qxz/Thermodynamics.jl): ideal-gas, Peng-Robinson, and Helmholtz equations of state in skeletal form. Not updated since ~Julia v0.7.
- [Airship.jl](https://github.com/ordovician/Airship.jl): atmospheric buoyancy & ideal gases
- [ModiaMedia.jl](https://github.com/ModiaSim/ModiaMedia.jl): perhaps the most feature-complete; based on [Modelica.Media](https://doc.modelica.org/Modelica%203.2.3/Resources/helpDymola/Modelica_Media.html#Modelica.Media)
- [ThermoModels.jl](https://github.com/longemen3000/ThermoModels.jl): saturation & phase equilibria with several popular equations of state
- [ThermofluidQuantities.jl](https://github.com/UCLAMAEThreads/ThermofluidQuantities.jl): Unitful extensions for thermophysical properties, plus models for ideal gases & several liquids

![](https://sea2.discourse-cdn.com/julialang/images/transparent.png)

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [March 31, 2021, 11:02pm UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310/6 "2021-03-31T23:02:47Z")

</div>

hi, im the main developer of ThermoModels.jl. basically im porting all the relevant functionality to [OpenSAFT.jl](https://github.com/ypaul21/OpenSAFT.jl) (on Juliacon 2021, chat in [Zulip](https://julialang.zulipchat.com/#narrow/stream/265161-OpenSAFT.2Ejl)). 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](https://github.com/longemen3000/ThermoState.jl) using julia dispatch to select the correct specification. (i used this on [MoistAir.jl](https://github.com/longemen3000/MoistAir.jl) and [WaterIF97.jl](https://github.com/longemen3000/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)](https://github.com/CalebBell/chemicals) , 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](https://github.com/CalebBell/chemicals/issues/28) is basically solved by julia broadcasting

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [March 31, 2021, 11:05pm UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310/7 "2021-03-31T23:05:43Z")

</div>

CoolProp seems to use the same general approach when calculating derivatives, [Pure and Pseudo-Pure fluid properties — CoolProp 6.4.1 documentation](http://www.coolprop.org/fluid_properties/PurePseudoPure.html#thermodynamic-properties-of-fluid), so in theory, we can just plug in the helmholtz model defined on coolprop and make it work on julia.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [April 1, 2021, 12:40am UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310/8 "2021-04-01T00:40:38Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [April 1, 2021, 1:06am UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310/10 "2021-04-01T01:06:53Z")

</div>

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.

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

```

on CoolProp (using CoolProp.jl)

```julia
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")

```

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [April 24, 2021, 12:21am UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310/11 "2021-04-24T00:21:08Z")

</div>

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](https://www.sciencedirect.com/science/article/pii/S1876610212010648) . 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.

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [April 24, 2021, 7:00am UTC](https://discourse.julialang.org/t/ann-pythermo-jl-an-interface-to-the-thermo-library-for-thermophysical-properties/58310/12 "2021-04-24T07:00:29Z")

</div>

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
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
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):

```julia
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
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
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
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
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)