Integration with Symbolics.jl

[longemen3000]

at the moment, due to how we define functions and the fact that our models need to support ForwardDiff.jl, we have the ability, in theory, to generate fully symbolic versions of an EoS. for a proof of concept, this can be done (if we remove the NaN-safe protections):

using Clapeyron, Symbolics
model = UNIFAC(["water","ethanol"])
@variables v0 T0, x1, x2
x = [x1,x2]
Ge = Clapeyron.excess_gibbs_free_energy(model,v0,T0,x)
dGe = Symbolics.gradient(Ge,x)
act = Symbolics.simplify(exp.(dGe) ./(Clapeyron.R̄ * T0))

that gives the expression for the activity coefficient in terms of T0 and x (v0 is ignored).

 (0.12027235504272604exp(8.31446261815324T0*((-2.4561x1*(2.4561 / (2.4561x1 + 2.6616x2) + (-6.032427210000001x1) / ((2.4561x1 + 2.6616x2)^2) + (-2.19256047x2*exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-1.73916441x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.60543088x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((2.4561x1) / (2.4561x1 + 2.6616x2) + (0.7081x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 +
2.6616x2) + (0.8927x2*exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / (2.4561x1 + 2.6616x2)) - 2.4561log((2.4561x1) / (2.4561x1 + 2.6616x2) + (0.7081x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)
+ (0.8927x2*exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / (2.4561x1 + 2.6616x2)) - x2*((1.7689((-1.73916441x2) / ((2.4561x1 + 2.6616x2)^2) + (2.4561exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (-2.60543088x2) / ((2.4561x1 + 2.6616x2)^2) + (-6.032427210000001x1*exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.19256047x2*exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((0.7081x2) / (2.4561x1 + 2.6616x2) + (1.0608x2) / (2.4561x1 + 2.6616x2) + (0.8927x2*exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (2.4561x1*exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)) + (0.8927((-2.19256047x2) / ((2.4561x1 + 2.6616x2)^2) + (2.4561exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (-6.032427210000001x1*exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-1.73916441x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.60543088x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((0.8927x2) / (2.4561x1 + 2.6616x2) + (0.7081x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (2.4561x1*exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)))) + 8.31446261815324T0*((x1*(1.5106856162723488x1 + 1.9853131596298699x2)*(1.5106856162723488 / (1.5106856162723488x1 + 1.9853131596298699x2) + (-2.282171031212166(x1 + x2)) / ((1.5106856162723488x1 + 1.9853131596298699x2)^2))) / (1.5106856162723488(x1 + x2)) + (x2*(1.5106856162723488x1 + 1.9853131596298699x2)*(1.9853131596298699 / (1.5106856162723488x1 + 1.9853131596298699x2) + (-2.999184034049054(x1 + x2)) / ((1.5106856162723488x1 + 1.9853131596298699x2)^2))) / (1.9853131596298699(x1 + x2)) + (21.2870187x1*(2.4561x1 + 2.6616x2)*(4.25740374 / (1.7334(2.4561x1 + 2.6616x2)) + (-10.456609325814(1.7334x1 + 2.4951999999999996x2)) / (3.0046755600000004((2.4561x1 + 2.6616x2)^2)))) / (2.4561(1.7334x1 + 2.4951999999999996x2)) + (33.206121599999996x2*(2.4561x1 + 2.6616x2)*(4.6136174400000005 / (2.4951999999999996(2.4561x1 + 2.6616x2)) + (-16.311511052352(1.7334x1 + 2.4951999999999996x2)) / (6.2260230399999985((2.4561x1 + 2.6616x2)^2)))) / (2.6616(1.7334x1 + 2.4951999999999996x2)) + 12.2805log((2.4561(1.7334x1 + 2.4951999999999996x2)) / (1.7334(2.4561x1 + 2.6616x2))) + log((1.5106856162723488(x1 + x2)) / (1.5106856162723488x1 + 1.9853131596298699x2))))) / T0
 (0.12027235504272604exp(8.31446261815324T0*((-2.4561x1*((0.7081exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (0.8927exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / (2.4561x1 + 2.6616x2) + (-6.53715576x1) / ((2.4561x1 + 2.6616x2)^2) + (-2.3760103200000002x2*exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-1.8846789599999998x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.82342528x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((2.4561x1) / (2.4561x1 + 2.6616x2) + (0.7081x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((3.6156T0 - 1391.3 - 0.001144(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (0.8927x2*exp((801.9 + 0.007514(T0^2) - 3.824T0) / T0)) / (2.4561x1 + 2.6616x2)) - 0.8927(log((0.8927x2) / (2.4561x1 + 2.6616x2) + (0.7081x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (2.4561x1*exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)) - log(0.37571385632702137(0.8927 + 1.7689exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)))) - 1.7689(log((0.7081x2) / (2.4561x1 + 2.6616x2) + (1.0608x2) / (2.4561x1
+ 2.6616x2) + (0.8927x2*exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (2.4561x1*exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)) - log(0.37571385632702137(1.7689 + 0.8927exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)))) - x2*((1.7689(0.7081 / (2.4561x1 +
2.6616x2) + 1.0608 / (2.4561x1 + 2.6616x2) + (-1.8846789599999998x2) / ((2.4561x1 + 2.6616x2)^2) + (0.8927exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0))
/ (2.4561x1 + 2.6616x2) + (-2.82342528x2) / ((2.4561x1 + 2.6616x2)^2) + (-6.53715576x1*exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.3760103200000002x2*exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((0.7081x2) / (2.4561x1 + 2.6616x2) + (1.0608x2) / (2.4561x1 + 2.6616x2) + (0.8927x2*exp((4.746T0 - 1606.0 - 0.0009181(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (2.4561x1*exp((17.253 - 0.8389T0 - 0.0009021(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)) + (0.8927(0.8927 / (2.4561x1 + 2.6616x2) + (-2.3760103200000002x2) / ((2.4561x1 + 2.6616x2)^2) + (0.7081exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (-6.53715576x1*exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-1.8846789599999998x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2) + (-2.82342528x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / ((2.4561x1 + 2.6616x2)^2))) / ((0.8927x2) / (2.4561x1
+ 2.6616x2) + (0.7081x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0)) / (2.4561x1 + 2.6616x2) + (1.0608x2*exp((4.674T0 - 2777.0 - 0.001551(T0^2)) / T0))
/ (2.4561x1 + 2.6616x2) + (2.4561x1*exp((8.673T0 - 1460.0 - 0.01641(T0^2)) / T0)) / (2.4561x1 + 2.6616x2)))) + 8.31446261815324T0*((x1*(1.5106856162723488x1 + 1.9853131596298699x2)*(1.5106856162723488 / (1.5106856162723488x1 + 1.9853131596298699x2) + (-2.999184034049054(x1 + x2)) / ((1.5106856162723488x1 +
1.9853131596298699x2)^2))) / (1.5106856162723488(x1 + x2)) + (x2*(1.5106856162723488x1 + 1.9853131596298699x2)*(1.9853131596298699 / (1.5106856162723488x1 + 1.9853131596298699x2) + (-3.941468341799537(x1 + x2)) / ((1.5106856162723488x1 + 1.9853131596298699x2)^2))) / (1.9853131596298699(x1 + x2)) + (21.2870187x1*(2.4561x1 + 2.6616x2)*(6.12846072 / (1.7334(2.4561x1 + 2.6616x2)) + (-11.331505794384002(1.7334x1 + 2.4951999999999996x2)) / (3.0046755600000004((2.4561x1 + 2.6616x2)^2)))) / (2.4561(1.7334x1 + 2.4951999999999996x2)) + (33.206121599999996x2*(2.4561x1 + 2.6616x2)*(6.641224319999999 / (2.4951999999999996(2.4561x1 + 2.6616x2)) + (-17.676282650111997(1.7334x1 + 2.4951999999999996x2)) / (6.2260230399999985((2.4561x1 + 2.6616x2)^2)))) / (2.6616(1.7334x1 + 2.4951999999999996x2)) + 13.308log((2.6616(1.7334x1 + 2.4951999999999996x2)) / (2.4951999999999996(2.4561x1 + 2.6616x2))) + log((1.9853131596298699(x1 + x2)) / (1.5106856162723488x1 + 1.9853131596298699x2))))) / T0

In practice, we need additional support to make that a seamless experience. my vision is that:
Clapeyron.activity_coefficient(model,v,T,z) just works and returns the expression above. for that, we would need:

• Register nan-safe functions (and its derivatives) (tier 0: support primal functions)
• Register the differentials so they return symbolic derivatives instead of ForwardDiff ones (tier 1:support bulk properties)

i don't know if we can/should support higher tiers. any commentary about that it is appreciated

[ysyecust]

I encountered this error when trying your code.
{
"name": "LoadError",
"message": "LoadError: AssertionError: @variables expects a tuple of expressions or an expression of a tuple (@variables x y z(t) v[1:3] w[1:2,1:4] or @variables x y z(t) v[1:3] w[1:2,1:4] k=1.0)
in expression starting at c:\Users\YSY_e\Documents\TOP\BestProcessSimulation\equation\solver\jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_W0sZmlsZQ==.jl:3",
"stack": "LoadError: AssertionError: @variables expects a tuple of expressions or an expression of a tuple (@variables x y z(t) v[1:3] w[1:2,1:4] or @variables x y z(t) v[1:3] w[1:2,1:4] k=1.0)
in expression starting at c:\Users\YSY_e\Documents\TOP\BestProcessSimulation\equation\solver\jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_W0sZmlsZQ==.jl:3

Stacktrace:
[1] _parse_vars(macroname::Symbol, type::Type, x::Tuple{Symbol, Expr}, transform::Function)
@ Symbolics C:\Users\YSY_e\.julia\packages\Symbolics\2UpZj\src\variable.jl:134
[2] _parse_vars(macroname::Symbol, type::Type, x::Tuple{Symbol, Expr})
@ Symbolics C:\Users\YSY_e\.julia\packages\Symbolics\2UpZj\src\variable.jl:80
[3] var"@variables"(source::LineNumberNode, module::Module, xs::Vararg{Any})
@ Symbolics C:\Users\YSY_e\.julia\packages\Symbolics\2UpZj\src\variable.jl:377
[4] eval
@ .\boot.jl:429 [inlined]
[5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:2571
[6] #invokelatest#2
@ .\essentials.jl:1043 [inlined]
[7] invokelatest
@ .\essentials.jl:1040 [inlined]
[8] (::VSCodeServer.var"#217#218"{VSCodeServer.NotebookRunCellArguments, String})()
@ VSCodeServer c:\Users\YSY_e\.vscode\extensions\julialang.language-julia-1.105.2\scripts\packages\VSCodeServer\src\serve_notebook.jl:24
[9] withpath(f::VSCodeServer.var"#217#218"{VSCodeServer.NotebookRunCellArguments, String}, path::String)
@ VSCodeServer c:\Users\YSY_e\.vscode\extensions\julialang.language-julia-1.105.2\scripts\packages\VSCodeServer\src\repl.jl:276
[10] notebook_runcell_request(conn::VSCodeServer.JSONRPC.JSONRPCEndpoint{Base.PipeEndpoint, Base.PipeEndpoint}, params::VSCodeServer.NotebookRunCellArguments)
@ VSCodeServer c:\Users\YSY_e\.vscode\extensions\julialang.language-julia-1.105.2\scripts\packages\VSCodeServer\src\serve_notebook.jl:13
[11] dispatch_msg(x::VSCodeServer.JSONRPC.JSONRPCEndpoint{Base.PipeEndpoint, Base.PipeEndpoint}, dispatcher::VSCodeServer.JSONRPC.MsgDispatcher, msg::Dict{String, Any})
@ VSCodeServer.JSONRPC c:\Users\YSY_e\.vscode\extensions\julialang.language-julia-1.105.2\scripts\packages\JSONRPC\src\typed.jl:67
[12] serve_notebook(pipename::String, debugger_pipename::String, outputchannel_logger::Base.CoreLogging.SimpleLogger; error_handler::var"#5#10"{String})
@ VSCodeServer c:\Users\YSY_e\.vscode\extensions\julialang.language-julia-1.105.2\scripts\packages\VSCodeServer\src\serve_notebook.jl:147
[13] top-level scope
@ c:\Users\YSY_e\.vscode\extensions\julialang.language-julia-1.105.2\scripts
otebook
otebook.jl:35"
}

[longemen3000]

Hello,

That seems more like a failure in the @variables macro than in the Clapeyron code itself, can you share the notebook, or the cell that generated the error?

[ysyecust]

Hello,

That seems more like a failure in the @variables macro than in the Clapeyron code itself, can you share the notebook, or the cell that generated the error?

using Clapeyron, Symbolics
model = UNIFAC(["water","ethanol"])
@variables v0 T0, x1, x2
x = [x1,x2]
Ge = Clapeyron.excess_gibbs_free_energy(model,v0,T0,x)
dGe = Symbolics.gradient(Ge,x)
act = Symbolics.simplify(exp.(dGe) ./(Clapeyron.R̄ * T0))

Thank you for your replay, I also have a question about “How to obtain the derivative of the fugacity calculation ？”

namelist = ["ethane","ethylene","propylene","methane","propane"]
numcomp = length(namelist)
model = SRK(namelist;alpha=SoaveAlpha)
fugacity_coefficient(model,802384,233,[0.2,0.2,0.2,0.2,0.2,],phase=:l)

I want to get dFugdT(size=n) dFugdP(size=n) and dFugdN (size =n*n) n=5

[longemen3000]

The first error seems to be here @variables v0 T0, x1, x2 , it should be @variables v0 T0 x1 x2

For the derivatives of fugacity there is Clapeyron.∂lnϕ∂n∂P∂T(model,p,T,z), returns lnϕ, ∂lnϕ∂n, ∂lnϕ∂P, ∂lnϕ∂T, V. You can see the function here:

Clapeyron.jl/src/methods/fugacity_coefficient.jl

Line 89 in cd0a7d1

function ∂lnϕ∂n∂P∂T(model::EoSModel, p, T, z=SA[1.],cache = ∂lnϕ_cache(model,p,T,z,Val{true}()); phase=:unknown, vol0=nothing)