Closed jfkonecn closed 3 years ago
coveralls
commented
4 years ago Coverage increased (+0.4%) to 89.171% when pulling ce8116f27efcb7b0f76586960bc997c981627cca on jfkonecn:feature/issue-42 into cba02552ea1a6947fcaa8407d4553e42e248e77c on CalebBell:master.
coveralls
commented
4 years ago Coverage increased (+0.3%) to 89.148% when pulling 926413d86a212acaf3ff99e26256a249eb039450 on jfkonecn:feature/issue-42 into cba02552ea1a6947fcaa8407d4553e42e248e77c on CalebBell:master.
CalebBell
commented
4 years ago Hi, I've given a lot of thought to your PR. You are clearly a talented coder - you figured out the python nuances of packaging, wrote good tests where you felt they were necessary, and wrote code the way you felt code is supposed to be written. There are things I am trying to move away from - like reading data from files instead of including the constants in the program, which is responsible for thermo's load time being measured in seconds. I am avoiding SciPy and Numpy lately in favor of custom functions, as their speed tends to favor arrays with lots of numbers - but in chemical engineering, our arrays are often smaller and work faster with pure Python.
The PTVEntry class is something sort of needed for this type of an addition - the current structure of the library is in favor of doing everything in the Chemical and Mixture classes, but I am moving in the development version into classes I am calling "Phase" (properties of a specific phase), "Flasher" (does the solving of equilibrium equations), and "EquilibriumState" (stores the results).
The biggest thing I want to talk to you about is speed. IAPWS's math is slow enough that most simulations will be bound by its calculation speed. Admittedly it has great accuracy - but it comes at a cost. A lot of that cost can be avoided by programming for speed. This is harder in Python than most languages, but still possible. I have done a bit of work on IAPWS95; here is what I think is the fastest way to implement IAPWS formulas. The IAPWS97 implementation has fewer coefficients and also has backward equations by the idea could be the same for it, and its time is still spent doing math.
I use CoolProp as a reference model for speed - I hope you know of it?
from CoolProp.CoolProp import PropsSI
%timeit PropsSI('DMASS', 'T', T, 'P', P_spec, 'water')
There here is a basic pure-python implementation of just one of the many functions, but the only one needed to solve a PVT point.
cis = [None, None, None, None, None, None, None, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 3, 3, 3, 3, 4, 6, 6, 6, 6, None, None, None]
dis = [1, 1, 1, 2, 2, 3, 4, 1, 1, 1, 2, 2, 3, 4, 4, 5, 7, 9, 10, 11, 13,
15, 1, 2, 2, 2, 3, 4, 4, 4, 5, 6, 6, 7, 9, 9, 9, 9, 9, 10, 10, 12,
3, 4, 4, 5, 14, 3, 6, 6, 6, 3, 3, 3]
tis = [-0.5, 0.875, 1., 0.5, 0.75, 0.375, 1., 4., 6., 12., 1., 5., 4., 2.,
13., 9., 3., 4., 11., 4., 13., 1., 7., 1., 9., 10., 10., 3., 7.,
10., 10., 6., 10., 10., 1., 2., 3., 4., 8., 6., 9., 8., 16., 22.,
23., 23., 10., 50., 44., 46., 50., 0., 1., 4.]
nis = [0.12533547935523E-1, 0.78957634722828E1, -0.87803203303561E1,
0.31802509345418, -0.26145533859358, -0.78199751687981E-2,
0.88089493102134E-2, -0.66856572307965, 0.20433810950965,
-0.66212605039687E-4, -0.19232721156002, -0.25709043003438,
0.16074868486251, -0.40092828925807E-1, 0.39343422603254E-6,
-0.75941377088144E-5, 0.56250979351888E-3, -0.15608652257135E-4,
0.11537996422951E-8, 0.36582165144204E-6, -0.13251180074668E-11,
-0.62639586912454E-9, -0.10793600908932, 0.17611491008752E-1,
0.22132295167546, -0.40247669763528, 0.58083399985759,
0.49969146990806E-2, -0.31358700712549E-1, -0.74315929710341,
0.47807329915480, 0.20527940895948E-1, -0.13636435110343,
0.14180634400617E-1, 0.83326504880713E-2, -0.29052336009585E-1,
0.38615085574206E-1, -0.20393486513704E-1, -0.16554050063734E-2,
0.19955571979541E-2, 0.15870308324157E-3, -0.16388568342530E-4,
0.43613615723811E-1, 0.34994005463765E-1, -0.76788197844621E-1,
0.22446277332006E-1, -0.62689710414685E-4, -0.55711118565645E-9,
-0.19905718354408, 0.31777497330738, -0.11841182425981,
-0.31306260323435E2, 0.31546140237781E2, -0.25213154341695E4,
-0.14874640856724, 0.31806110878444]
alphas = [20., 20., 20.]
betas = [150., 150., 250., 0.3, 0.3]
gammas = [1.21, 1.21, 1.25]
epsilons = [1., 1., 1.]
ais = [3.5, 3.5]
bis = [0.85, 0.95]
Bis = [0.2, 0.2]
Cis = [28., 32.]
Dis = [700., 800.]
Ais = [0.32, 0.32]
def d_psi_d_delta(i, tau, delta):
'''
>>> d_psi_d_delta(0, 647.096/647., 358./322.)
-4.411951793785948
'''
psi = exp(-Cis[i]*(delta-1)**2 - Dis[i]*(tau-1)**2)
ans = -2*Cis[i]*(delta-1)*psi
return ans
def d_Delta_bd_delta(i, tau, delta):
''''
>>> d_Delta_bd_delta(1, 647.096/647., 358./322.)
7.931159558108671e-06
'''
theta = (1-tau) + Ais[i]*((delta-1)**2)**(1./(2*betas[i+3]))
Delta = theta**2 + Bis[i]*((delta-1)**2)**ais[i]
_d_Delta_d_delta = (delta - 1)*(
Ais[i]*theta*2/betas[i+3]*((delta-1.)**2 )**(1./(2*betas[i+3])-1.)
+ 2*Bis[i]*ais[i]*((delta-1)**2)**(ais[i]-1))
ans = bis[i]*Delta**(bis[i]-1)*_d_Delta_d_delta
return ans
def dAddelta_res(tau, delta):
'''
>>> dAddelta_res(647.096/647., 358./322.)
-0.714012024371285
'''
phir = 0.0
for i in range(7):
phir += nis[i]*dis[i]*delta**(dis[i]-1.0)*tau**tis[i]
for i in range(7,51):
phir += nis[i]*exp(-delta**cis[i])*(delta**(dis[i]-1.0)*tau**tis[i])*(dis[i]-cis[i]*delta**cis[i])
for i in range(51, 54):
phir += nis[i]*delta**dis[i]*tau**tis[i]*exp(
-alphas[i-51]*(delta-epsilons[i-51])**2.0 - betas[i-51]*(tau-gammas[i-51])**2.0)*(
dis[i]/delta - 2.0*alphas[i-51]*(delta-epsilons[i-51]))
for i in range(2):
theta = (1.0-tau) + Ais[i]*((delta-1.0)**2.0)**(1.0/(2.0*betas[i+3]))
psi = exp(-Cis[i]*(delta-1.0)**2.0 - Dis[i]*(tau-1.0)**2.0)
Delta = theta**2.0 + Bis[i]*((delta-1.0)**2.0)**ais[i]
_d_psi_d_delta = d_psi_d_delta(i, tau, delta)
_d_Delta_bd_delta = d_Delta_bd_delta(i, tau, delta)
phir += nis[i+54]*(Delta**bis[i]*(psi + delta*_d_psi_d_delta)
+ _d_Delta_bd_delta*psi*delta)
return phir
That function is pretty slow because of all the math. But it turns out a lot of the powers, etc. are repeated, and indexing the coefficients take time - and we can do a little better. Do you know if SymPy, the mathematical library? It can assist here.
from sympy import *
tau, delta = symbols('tau, delta')
new_expr = dAddelta_res(tau, delta)
(cse(new_expr, optimizations='basic'))This yields the massive expression:
([(x0, tau**0.5),
(x1, tau**1.0),
(x2, delta**1.0),
(x3, delta**2.0),
(x4, delta**3.0),
(x5, tau**4.0),
(x6, delta - 1),
(x7, exp(-delta)),
(x8, x6*x7),
(x9, tau**6.0),
(x10, x1*x7),
(x11, x2*(delta - 2)),
(x12, x5*x7),
(x13, tau**2.0),
(x14, x4*x7*(delta - 4)),
(x15, tau**13.0),
(x16, delta**4.0),
(x17, tau**9.0),
(x18, delta**6.0),
(x19, tau**3.0),
(x20, delta**8.0),
(x21, x20*x5),
(x22, delta**9.0),
(x23, delta**2),
(x24, 2*x23),
(x25, exp(-x23)),
(x26, tau**7.0*x25),
(x27, x2*(x23 - 1)),
(x28, x25*x27),
(x29, tau**10.0),
(x30, x25*x29),
(x31, delta**3),
(x32, exp(-x31)),
(x33, x4*(x23 - 2)),
(x34, delta**5.0),
(x35, x34*(x23 - 3)),
(x36, delta**6),
(x37, exp(-x36)),
(x38, x34*(x36 - 1)),
(x39, x37*x38),
(x40, tau**50.0*x37),
(x41, x22*x25*(x23 - 5)),
(x42, tau**8.0*x25),
(x43, 3*x31),
(x44, x4*(x43 - 4)),
(x45, tau**23.0*x32),
(x46, x24 - 9),
(x47, x25*x46),
(x48, x20*x47),
(x49, delta**4),
(x50, delta - 1.0),
(x51, x50**2.0),
(x52, -20.0*x51),
(x53, x31*(-40.0*delta + 40.0 + 3/delta)),
(x54, x53*exp(x52 - 219.615*(0.826446280991736*tau - 1)**2.0)),
(x55, -tau),
(x56, 0.2*x51**3.5 + (0.32*x51**1.66666666666667 + x55 + 1.0)**2.0),
(x57, 32.0*x51),
(x58, (tau - 1.0)**2.0),
(x59, 800.0*x58),
(x60, x6**2),
(x61, (tau - 1)**2),
(x62, delta*x6),
(x63, 28.0*x51),
(x64, 700.0*x58),
(x65, x60**1.66666666666667),
(x66, 0.2*x60**3.5 + (x55 + 0.32*x65 + 1)**2),
(x67,
x62*(1.4*x60**2.5 + (-2.13333333333333*tau + 0.682666666666667*x65 + 2.13333333333333)*(x50**2)**0.666666666666667))],
[-3.6582165144204e-7*delta**10.0*x12*(delta - 11) + 3.277713668506e-5*delta**11.0*x42*(x23 - 6) + 1.3251180074668e-12*delta**12.0*x15*x7*(delta - 13) + 0.00012537942082937*delta**13.0*x29*(2*x49 - 7)*exp(-x49) + 6.2639586912454e-10*delta**14.0*x10*(delta - 15) - 0.0234599255063943*tau**0.375*x3 - 0.52291067718716*tau**0.75*x2 + 7.8957634722828*tau**0.875 + 0.25709043003438*tau**5.0*x11*x7 - 1.1537996422951e-9*tau**11.0*x22*x7*(delta - 10) + 6.6212605039687e-5*tau**12.0*x8 - 0.130840847171433*tau**16.0*x3*x32*(x31 - 1) - 0.034994005463765*tau**22.0*x32*x44 + 1.19434310126448*tau**44.0*x39 - 1.90664983984428*tau**46.0*x39 + 0.63605018690836*x0*x2 - 0.035222982017504*x1*x28 + 0.0352357972408536*x1*x4 - 0.0083326504880713*x1*x48 + 31.546140237781*x1*x54 - 8.7803203303561*x1 + 0.19232721156002*x10*x11 - 0.16074868486251*x12*x3*(delta - 3) + 0.040092828925807*x13*x14 + 0.029052336009585*x13*x48 - 3.9343422603254e-7*x14*x15 + 7.5941377088144e-6*x16*x17*x7*(delta - 5) - 0.4780732991548*x16*x30*(x24 - 5) - 0.022446277332006*x16*x45*(x43 - 5) - 0.44264590335092*x17*x28 - 0.00031740616648314*x17*x41 - 0.00056250979351888*x18*x19*x7*(delta - 7) - 0.014180634400617*x18*x30*(x24 - 7) - 0.0099938293981612*x19*x25*x33 - 0.038615085574206*x19*x48 + 0.0016554050063734*x20*x42*x46 + 0.020393486513704*x21*x47 + 1.5608652257135e-5*x21*x7*(delta - 9) - 0.041055881791896*x25*x35*x9 + 0.062717401425098*x26*x33 + 0.10793600908932*x26*(x24 - 1) + 0.80495339527056*x27*x30 - 0.58083399985759*x3*x30*(x24 - 3) + 1.67133355696935e-9*x3*x40*(2*x36 - 1) + 1.48631859420682*x30*x33 + 0.27272870220686*x30*x35 + 0.71047094555886*x38*x40 - 0.0039911143959082*x41*x9 + 0.076788197844621*x44*x45 - 2521.3154341695*x5*x53*exp(x52 - 390.625*(0.8*tau - 1)**2.0) + 0.66856572307965*x5*x8 - 31.306260323435*x54 - 0.14874640856724*x56**0.85*(-56.0*x62*exp(-28.0*x60 - 700.0*x61) + exp(-x63 - x64)) + 0.31806110878444*x56**0.95*(-64.0*x62*exp(-32.0*x60 - 800.0*x61) + exp(-x57 - x59)) - 0.126434447282154*x66**(-0.15)*x67*exp(-x63 - x64) + 0.302158053345218*x66**(-0.05)*x67*exp(-x57 - x59) - 0.20433810950965*x8*x9 + 0.012533547935523/x0])
With some work, this can be turned into code:
def dAddelta_res_fast(tau, delta):
tau05 = x0 = tau**0.5
tau_eighth = tau**0.125
tau_quarter = tau_eighth*tau_eighth
x1 = tau
x2 = delta
delta2 = x3 = x23 = delta*delta
delta3 = x4 = x31 = delta*x3 # delta 3
delta4 = x16 = x49 = x4*delta #delta**4.0 # delta 4
delta5 = x34 = delta*delta4
delta6 = x18 = x36 = delta*delta5
delta8 = x20 = delta6*delta2
delta9 = x22 = delta*delta8
x6 = delta - 1.0
x7 = exp(-delta)
x8 = x6*x7
tau2 = x13 = tau*tau
tau3 = x19 = tau*tau2
tau4 = x5 = tau*tau3
tau6 = x9 = tau3*tau3
tau7 = tau6*tau
tau9 = x17 = tau6*tau3
tau10 = x29 = tau*tau9
tau12 = tau10*tau2
tau13 = x15 = tau10*tau3
tau23 = tau10*tau13
x10 = x1*x7
x11 = x2*(delta - 2.0)
x12 = x5*x7
x14 = x4*x7*(delta - 4.0)
x21 = x20*x5
x24 = x23 + x23
x25 = exp(-x23)
x26 = tau*tau6*x25
x27 = x2*(x23 - 1.0)
x28 = x25*x27
x30 = x25*x29
x32 = exp(-x31)
x33 = x4*(x23 - 2.0)
x35 = x34*(x23 - 3.0)
x37 = exp(-x36)
x38 = x34*(x36 - 1.0)
x39 = x37*x38
x40 = tau23*tau23*tau4*x37
x41 = x22*x25*(x23 - 5.0)
x42 = tau4*tau4*x25
x43 = 3.0*x31
x44 = x4*(x43 - 4.0)
x45 = tau23*x32
x46 = x24 - 9.0
x47 = x25*x46
x48 = x20*x47
x50 = delta - 1.0
x51 = x50*x50
dn1_2_23 = x51**(2.0/3.0)
x52 = -20.0*x51
x53 = x31*(-40.0*delta + 40.0 + 3.0/delta)
c54 = (tau - 1.21)
x54 = x53*exp(x52 - 150.0*c54*c54)
x55 = -tau
y50 =(0.32*x51*dn1_2_23 + x55 + 1.0)
x56 = 0.2*x51**3.5 + y50*y50
x57 = 32.0*x51
x58 = x61 = (tau - 1.0)*(tau - 1.0)
x59 = 800.0*x58
x60 = x51 #= x6**2
x62 = delta*x6
x63 = 28.0*x51
x64 = 700.0*x58
x65 = x60**1.66666666666667
c100 = (x55 + 0.32*x65 + 1.0)
x66 = 0.2*x60**3.5 + c100*c100
c51 = exp(-x57 - x59)
c52 = exp(-x63 - x64)
x67 = x62*(1.4*x60**2.5 + (-2.13333333333333*tau + 0.682666666666667*x65 + 2.13333333333333)*dn1_2_23)
return (-3.6582165144204e-7*delta9*delta*x12*(delta - 11.0) + 3.277713668506e-5*delta9*delta2*x42*(x23 - 6)
+ 1.3251180074668e-12*delta9*delta3*x15*x7*(delta - 13.0)
+ 0.00012537942082937*delta9*delta4*x29*(x49 + x49 - 7.0)*exp(-x49)
+ 6.2639586912454e-10*delta9*delta5*x10*(delta - 15.0) - 0.0234599255063943*tau_quarter*tau_eighth*x3
- 0.52291067718716*tau_quarter*tau05*x2 + 7.8957634722828*tau_quarter*tau05*tau_eighth + 0.25709043003438*tau4*tau*x11*x7
- 1.1537996422951e-9*tau*tau10*x22*x7*(delta - 10.0) + 6.6212605039687e-5*tau12*x8
- 0.130840847171433*tau9*tau7*x3*x32*(x31 - 1.0) - 0.034994005463765*tau10*tau12*x32*x44
+ 1.19434310126448*tau23*tau12*tau9*x39 - 1.90664983984428*tau23*tau23*x39 + 0.63605018690836*x0*x2
- 0.035222982017504*x1*x28 + 0.0352357972408536*x1*x4 - 0.0083326504880713*x1*x48
+ 31.546140237781*x1*x54 - 8.7803203303561*x1 + 0.19232721156002*x10*x11
- 0.16074868486251*x12*x3*(delta - 3.0) + 0.040092828925807*x13*x14 + 0.029052336009585*x13*x48
- 3.9343422603254e-7*x14*x15 + 7.5941377088144e-6*x16*x17*x7*(delta - 5.0)
- 0.4780732991548*x16*x30*(x24 - 5.0) - 0.022446277332006*x16*x45*(x43 - 5.0)
- 0.44264590335092*x17*x28 - 0.00031740616648314*x17*x41 - 0.00056250979351888*x18*x19*x7*(delta - 7.0)
- 0.014180634400617*x18*x30*(x24 - 7.0) - 0.0099938293981612*x19*x25*x33
- 0.038615085574206*x19*x48 + 0.0016554050063734*x20*x42*x46 + 0.020393486513704*x21*x47
+ 1.5608652257135e-5*x21*x7*(delta - 9.0) - 0.041055881791896*x25*x35*x9 + 0.062717401425098*x26*x33
+ 0.10793600908932*x26*(x24 - 1.0) + 0.80495339527056*x27*x30 - 0.58083399985759*x3*x30*(x24 - 3.0)
+ 1.67133355696935e-9*x3*x40*(x36 + x36 - 1.0) + 1.48631859420682*x30*x33 + 0.27272870220686*x30*x35
+ 0.71047094555886*x38*x40 - 0.0039911143959082*x41*x9 + 0.076788197844621*x44*x45
- 2521.3154341695*x5*x53*exp(x52 - 250.0*(tau - 1.25)**2.0) + 0.66856572307965*x5*x8
- 31.306260323435*x54 - 0.14874640856724*x56**0.85*(-56.0*x62*exp(-28.0*x60 - 700.0*x61)
+ c52) + 0.31806110878444*x56**0.95*(-64.0*x62*exp(-32.0*x60 - 800.0*x61)
+ c51) - 0.126434447282154*x66**(-0.15)*x67*c52
+ 0.302158053345218*x66**(-0.05)*x67*c51 - 0.20433810950965*x8*x9 + 0.012533547935523/x0)
I just have a basic secant, derivative free solver from another of my libraries (also implemented for speed) to compare the performance of the two functions:
from math import exp
P_spec = 5e6
T = 600.0
tau = 647.096/T
rhoc_inv = 1.0/322.0
def to_solve(rho):
R = 461.51805
delta = rho*rhoc_inv
dAddelta_res_val = dAddelta_res_fast(tau, delta)
P_calc = (1.0 + dAddelta_res_val*delta)*rho*R*T
err = P_calc - P_spec
return err
def to_solve_slow(rho):
R = 461.51805
delta = rho*rhoc_inv
dAddelta_res_val = dAddelta_res(tau, delta)
P_calc = (1.0 + dAddelta_res_val*delta)*rho*R*T
err = P_calc - P_spec
return err
As you can see the faster version is 5x faster than the slow version. The faster version is also faster than CoolProp, although CoolProp does a whole lot more in the function call so it is not a fair comparison. If the faster code is ran in PyPy, the function call takes 10.5 us - which I believe is about as fast as any implementation in any language.
I hope this is helpful; it's just my thoughts.
jfkonecn
commented
4 years ago Have you tried running a profiler on your code? You can isolate the method calls that are giving you the most problems. https://docs.python.org/3.8/library/profile.html
Are you sure that's strictly file loading which is the problem? It seems like you, as well as my pull request, load one file at a time when you don't have to. You could put them on different threads. Python has async and await built into 3.5 and above to make threading easier. https://docs.python.org/3/library/asyncio-task.html https://docs.python.org/3/library/asyncio-stream.html I know this would require dropping 2.7 support, but python no longer supports it so I don't think it's a bad idea to stop supporting it as well especially if you're missing out really nice language features.
I haven't done a test with a profiler yet, but my guess would be that JSON would be faster to parse and less bug prone since python does all of the parsing for you. https://docs.python.org/3/library/json.html You could easily make a csv to JSON file converter so that you can keep the human readable csv and have the JSON for your code.
The biggest advantage you have to keeping massive amounts data in files is that it's easier to debug. It would be very easy to make a mistake when putting all of that data into code. It's good that you have unit tests to help you, but you could potential spend hours searching the code to find the break.
In my pull request I load the files right away, but I don't have to do this. I can wait to load the file until I absolutely need it.
Do you need to shave hundreds of microseconds off of execution time? If the users of this code base aren't going to notice a difference then I'm not sure adding potentially more difficultly in debugging is worth it. I don't think it's a big deal to keep the third party libraries, unless we see a performance boost that users will notice. It's one less area of code you have to worry about.
If you're think about using something like PTVEntry in the rest of your code then I think it needs to be a data transfer object, meaning the class only has get and set property methods. By doing this you aren't tied to a single implementation of object creation. It would be like how I have the interpolate_entry, __gibbs_method, vapor_method and __region3_density functions in my PR. PTVEntry can now be created anywhere in the library and all of the calling methods know they are entirely responsible for the values of each and every property.
resolves #42
This is a rough draft so feel free to be picky.
I created a new "PVTEntry" object to return all the steam information.
Do you happen to have an existing object which serves the same purpose? I didn't see one, but there's a good chance I missed it.