search

ClapeyronThermo / Clapeyron.jl

Clapeyron provides a framework for the development and use of fluid-thermodynamic models, including SAFT, cubic, activity, multi-parameter, and COSMO-SAC.
MIT License
194 stars 49 forks source link

Integration with Symbolics.jl #186

Open longemen3000 opened 1 year ago

longemen3000 commented 1 year ago

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:

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

ysyecust commented 1 month ago

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 commented 1 month ago

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 commented 1 month ago

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 commented 1 month ago

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:

https://github.com/ClapeyronThermo/Clapeyron.jl/blob/cd0a7d1d12f94221619941f371fb0d032ea409c4/src/methods/fugacity_coefficient.jl#L89