Phase equilibria including derivatives

import feos
import numpy as np

Vapor pressure

n = 10_000_000
fit_params = ["m", "sigma", "epsilon_k"]

# order: m, sigma, epsilon_k, mu
parameters = np.array([[1.5, 3.4, 230.0, 2.3]] * n)
temperature = np.expand_dims(np.linspace(250.0, 400.0, n), 1)
eos = feos.EquationOfStateAD.PcSaftNonAssoc
%time feos.vapor_pressure_derivatives(eos, fit_params, parameters, temperature)

CPU times: user 1min 39s, sys: 654 ms, total: 1min 40s
Wall time: 2.3 s

(array([  99933.70173872,   99933.76369505,   99933.82565141, ...,
        4841855.06247647, 4841856.28251754, 4841857.50255883],
       shape=(10000000,)),
 array([[-1.04684676e+05,  1.37290974e+05, -2.94409278e+03],
        [-1.04684730e+05,  1.37291045e+05, -2.94409438e+03],
        [-1.04684784e+05,  1.37291116e+05, -2.94409599e+03],
        ...,
        [-3.06545385e+06,  2.31346293e+06, -8.79512373e+04],
        [-3.06545458e+06,  2.31346323e+06, -8.79512573e+04],
        [-3.06545531e+06,  2.31346353e+06, -8.79512774e+04]],
       shape=(10000000, 3)),
 array([ True,  True,  True, ...,  True,  True,  True], shape=(10000000,)))

Liquid density

n = 10_000_000
fit_params = ["m", "sigma", "epsilon_k"]

# order: m, sigma, epsilon_k, mu
parameters = np.array([[1.5, 3.4, 230.0, 2.3]] * n)
temperature = np.linspace(250.0, 400.0, n)
pressure = np.array([1e5] * n)
input = np.stack((temperature, pressure), axis=1)
eos = feos.EquationOfStateAD.PcSaftNonAssoc
%time feos.liquid_density_derivatives(eos, fit_params, parameters, input)

CPU times: user 1min 24s, sys: 716 ms, total: 1min 25s
Wall time: 2.1 s

(array([22.7191149 , 22.71911436, 22.71911382, ...,  0.03027266,
         0.03027266,  0.03027266], shape=(10000000,)),
 array([[-1.68941171e+01, -2.24206688e+01,  2.73083514e-02],
        [-1.68941167e+01, -2.24206685e+01,  2.73083541e-02],
        [-1.68941164e+01, -2.24206681e+01,  2.73083568e-02],
        ...,
        [ 9.49488124e-05,  9.60641306e-06,  7.80107221e-07],
        [ 9.49487960e-05,  9.60641033e-06,  7.80107128e-07],
        [ 9.49487795e-05,  9.60640761e-06,  7.80107035e-07]],
       shape=(10000000, 3)),
 array([ True,  True,  True, ...,  True,  True,  True], shape=(10000000,)))

Bubble points

n = 1_000_000
fit_params = ["k_ij"]
parameters = np.array([[
    # Substance 1: m, sigma, epsilon_k, mu
    1.5, 3.4, 230.0, 2.3,
    # Substance 2: m, sigma, epsilon_k, mu
    2.3, 3.5, 245.0, 1.4,
    # k_ij
    0.01
]] * n)
temperature = np.linspace(200.0, 406.0, n)
molefracs = np.array([0.5] * n)
pressure = np.array([1e5] * n)
input = np.stack((temperature, molefracs, pressure), axis=1)
eos = feos.EquationOfStateAD.PcSaftNonAssoc
%time feos.bubble_point_pressure_derivatives(eos, fit_params, parameters, input)

CPU times: user 11min 21s, sys: 122 ms, total: 11min 21s
Wall time: 11.3 s

(array([   5142.13808145,    5142.20830389,    5142.27852715, ...,
        3828278.20909898, 3828290.44546341, 3828302.68185363],
       shape=(1000000,)),
 array([[   40721.23744125],
        [   40721.73456205],
        [   40722.23168782],
        ...,
        [10996502.56325178],
        [10996528.1146768 ],
        [10996553.66610128]], shape=(1000000, 1)),
 array([ True,  True,  True, ...,  True,  True,  True], shape=(1000000,)))

Bubble points with cross-association

n = 100_000
fit_params = ["k_ij"]
parameters = np.array([[
    # Substance 1: m, sigma, epsilon_k, mu, kappa_ab, epsilon_k_ab, na, nb
    1.5, 3.4, 230.0, 2.3, 0.01, 1200.0, 1.0, 2.0,
    # Substance 2: m, sigma, epsilon_k, mu, kappa_ab, epsilon_k_ab, na, nb
    2.3, 3.5, 245.0, 1.4, 0.005, 500.0, 1.0, 1.0,
    # k_ij
    0.01,
]] * n)
temperature = np.linspace(200.0, 388.0, n)
molefracs = np.array([0.5] * n)
pressure = np.array([1e5] * n)
input = np.stack((temperature, molefracs, pressure), axis=1)
eos = feos.EquationOfStateAD.PcSaftFull
%time feos.bubble_point_pressure_derivatives(eos, fit_params, parameters, input)

CPU times: user 3min 11s, sys: 15.9 ms, total: 3min 11s
Wall time: 3.17 s

(array([7.34250975e+02, 7.34379943e+02, 7.34508930e+02, ...,
        2.34916950e+06, 2.34925463e+06, 2.34933975e+06], shape=(100000,)),
 array([[4.64318575e+03],
        [4.64395838e+03],
        [4.64473112e+03],
        ...,
        [7.29257342e+06],
        [7.29279667e+06],
        [7.29301993e+06]], shape=(100000, 1)),
 array([ True,  True,  True, ...,  True,  True,  True], shape=(100000,)))