Having previously discussed Excess thermodynamic properties, excess energy models, the NRTL model specifically, and VLE, we now proceed to study LLE.
for array structures and operations
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
from scipy.constants import R
print(R)
8.3144598
def stdv(v): #will standardize an array, transforming all values to between 0 and 1
w=(v-np.min(v))/(np.max(v)-np.min(v))
return w
from numba import jit
@jit
def Gamma(T,c_x,q_alpha, q_A):
q_tau = q_A/T
q_G = np.exp(-(q_alpha*q_tau))
l_D = ((1/((q_G.T) @ c_x)).T)
q_E = (q_tau*q_G) * l_D
gamma = np.exp(((q_E+(q_E.T))-(((q_G * l_D) * (c_x.T)) @ (q_E.T))) @ c_x)
return gamma
Parameters from Renon and Prausnitz, 1969.
# Ethyl acetate (1) + water (2) + ethanol (3)
alpha12 = 0.4
alpha23 = 0.3
alpha13 = 0.3
# 6 binary Aij parameters
Dg12 = 1335 * 4.184 #J/K
Dg21 = 2510 * 4.184 #J/K
Dg23 = 976 * 4.184 #J/K
Dg32 = 88 * 4.184 #J/K
Dg13 = 301 * 4.184 #J/K
Dg31 = 322 * 4.184 #J/K
we will assemble the parameters in a matrix structure compatible with the current imlementation of the model
#assemble matrix with regressed parameters Dg_i,j, according to the model all diagonal terms are zero
Dg = np.array([[0, Dg12, Dg13],
[Dg21, 0, Dg23],
[Dg31, Dg32, 0]])
#assemble symmetric matrix alpha
alpha = np.array([[0, alpha12, alpha13],
[alpha12, 0, alpha23],
[alpha13, alpha23, 0]])
A= Dg/R
T = 273+70 #K
Now we start discussing the Liq-Liq equilibria flash calculation
it works as follows:
The degrees of Freedom for a flash calculation are temperature - T, pressure - P and global composition - z
|
Equilibriuma criteria
$$\mu_{i}^{\rm {L1}}=\mu_{i}^{\rm {L2}}, \forall i$$Devised algorithm after analytical simplification of repeated contributions:
$$x_{i}^{\rm {L1}} \gamma_{i}^{\rm {L1}}= x_{i}^{\rm {L2}} \gamma_{i}^{\rm {L2}}, \forall i$$based on Rachford and Rice (1952) approach to LV equilibrium
Ref: Rachford HH, Rice JD, Petroleum Trans AIME 1952;195:327-328
def ELLflash_explicit(Z,beta0,K0,MODEL,qolog):
#copy intial guesses to local variables
beta=beta0.copy()
K=K0.copy()
converged = 0 # flag indicating not converged yet
if qolog:
itlogi=0
itlogv=np.ndarray(2,dtype=object)
itlogv[0]=np.zeros(100*100)
itlogv[1]=np.zeros([100*100,3])
#initialize local arrays
a = np.zeros(3)
XI = np.zeros(3)
XII = np.zeros(3)
looped2=0
looped3=0
#managing loop for fixed point point algorithm
looping2 = 0
while ( (looping2 == 0) # means do at least once
or ((ResK)>0.0001 # convergence criteria
and looping2 < 100) ): # max number of iterations
prevK = K.copy()
#print(np.may_share_memory(K,prevK))
#Solve RachfordRice equation with a newton raphson loop
looping3 = 0
while ( (looping3 == 0) # means do at least once
or ((ResBeta)>0.0001 # convergence criteria
and looping3 < 100) ): #max number of iterations
prevbeta = beta.copy()
for i in range(3):
a[i] = ( (K[i]-1) / (1+beta*(K[i]-1)) )
R=0 #residue funtion (F)
for j in range(3):
R+=Z[j]*a[j]
J=0 #jacobian (Fprime)
for j in range(3):
J-=Z[j]*a[j]*a[j]
beta -= R/J #NR step
if beta[0,0]<0:
beta[0,0]=0
if 1<beta[0,0]:
beta[0,0]=1
ResBeta = abs(beta[0,0]-prevbeta[0,0])
looping3 += 1
if qolog:
#print(looping3,beta0,prevbeta[0,0],beta[0,0],ResBeta)
itlogv[0][itlogi] = beta[0,0]
itlogv[1][itlogi,0:3] = K[0:3,0]
itlogi += 1
looped3+= looping3
for i in range(3):
XI[i] = Z[i] * (1./(1+beta * (K[i]-1))) #Rachford-Rice composition updating scheme
gamaI = MODEL(XI/np.sum(XI)) #enforce normalization in the MODEL call
for i in range(3):
XII[i] = K[i]*XI[i]
gamaII = MODEL(XII/np.sum(XII)) #enforce normalization in the MODEL call
for i in range(3):
K[i] = (gamaI[i])/(gamaII[i]) #definition of K
ResK=0
for i in range(3):
ResK += np.sqrt((prevK[i]-K[i])**2)
looping2 += 1
looped2 += looping2
#print('Calculation has ended')
if (looping2 < 100 and looping3 < 100 and (0.<=beta[0,0]<=1.) ):
converged = 1
if qolog:
return XI, XII, beta, converged, looped2, looped3, itlogv
else:
return XI, XII, beta, converged, looping2, looped3
#test
#exact ad-hoc NRTL solution
#Ztest = np.array([[ 0.10526316, 0.84210526, 0.05263158]]).T
#XItest = np.array([[ 0.02616646, 0.93009914, 0.0437344 ]]).T
#XIItest = np.array([[ 0.51520944, 0.38604626, 0.0987443 ]]).T
#BETAtest = np.array([[ 0.16173772]])
#appproximate
Ztest = np.array([[ 0.1, 0.85, 0.05]]).T
XItest = np.array([[ 0.025, 0.93, 0.045 ]]).T
XIItest = np.array([[ 0.5, 0.39, 0.11 ]]).T
BETAtest = np.array([[ 0.16]])
Ktest=1/(XItest/XIItest)
beta0=BETAtest
k0=Ktest
Z=Ztest
MODEL = lambda x: Gamma(T,x,alpha,A)
ans=ELLflash_explicit(Z,beta0,k0,MODEL,True)
print('at T', T, 'Z', Z.T)
print('with ig of k0', k0.T,'beta0', beta0)
print('solution is')
print('xi',ans[0].T)
print('xii',ans[1].T)
print('beta',ans[2])
print('converged?',ans[3], '-- 1 means yes')
at T 343 Z [[ 0.1 0.85 0.05]] with ig of k0 [[ 20. 0.41935484 2.44444444]] beta0 [[ 0.16]] solution is xi [ 0.02558642 0.93241587 0.04199771] xii [ 0.52329144 0.38118861 0.09551994] beta [[ 0.14951341]] converged? 1 -- 1 means yes
plt.plot(stdv(ans[6][0][:ans[5]]))
plt.plot(stdv(ans[6][1][:ans[5],0]))
plt.plot(stdv(ans[6][1][:ans[5],1]))
plt.plot(stdv(ans[6][1][:ans[5],2]))
plt.xlim(0,ans[5])
plt.show()
%timeit ans=ELLflash_explicit(Z,beta0,k0,MODEL,False)
100 loops, best of 3: 8.45 ms per loop
def ELLflash(Z,beta0,K0,MODEL): #ELLflash_linalg
beta=beta0.copy()
K=K0.copy()
converged = 0
looped2 = 0
looped3 = 0
looping2 = 0
while ( (looping2 == 0) # means do at least once
or ((ResK)>0.0001 # convergence criteria
and looping2 < 100) ): # max number of iterations
prevK = K.copy()
looping3 = 0
while ( (looping3 == 0)
or (ResBeta>0.0001 # convergence criteria
and looping3 < 100) ): #max number of iterations
prevbeta = beta.copy()
a = (1/(1+beta*(K-1))) * (K-1)
beta = beta+(1/((Z.T*a.T) @ a)) @ (Z.T @ a)
if beta[0,0]<0:
beta[0,0]=0
if 1<beta[0,0]:
beta[0,0]=1
ResBeta = beta[0,0]-prevbeta[0,0]
looping3 += 1
looped3 += looping3
XI = (1./(1+beta * (K-1))) * Z
gamaI = MODEL(XI/np.sum(XI))
XII = K*XI
gamaII = MODEL(XII/np.sum(XII))
K = (gamaI)/(gamaII)
ResK = np.linalg.norm(K-prevK)
looping2 += 1
looped2 += looping2
#print('Calculation has ended')
if (looping2 < 100 and looping3 < 100 and ( 0. <= beta[0,0] <= 1.) ):
converged = 1
return XI, XII, beta, converged, looped2, looped3
# no linalg, close to solution
%timeit ELLflash_explicit(Z,beta0,k0,MODEL,False)
ans=ELLflash_explicit(Z,beta0,k0,MODEL,False)
print(ans[3],ans[2])
# no linalg, a little far from solution
%timeit ELLflash_explicit(Z,beta0*beta0,k0/Z,MODEL,False)
ans=ELLflash_explicit(Z,2*beta0,k0/Z,MODEL,False)
print(ans[3],ans[2])
# linalg, close to solution
%timeit ELLflash(Z,beta0,k0,MODEL)
ans=ELLflash(Z,beta0,k0,MODEL)
print(ans[3],ans[2])
# linalg, a little far from solution
%timeit ELLflash(Z,2*beta0,k0/Z,MODEL)
ans=ELLflash(Z,2*beta0,k0/Z,MODEL)
print(ans[3],ans[2])
#approximately half the time here
# close to solution, 1.6ms drops to 0.9
# a little far from solution, 6.4ms drops to 3,4
100 loops, best of 3: 8.37 ms per loop 1 [[ 0.14951341]] 100 loops, best of 3: 13.3 ms per loop 1 [[ 0.14951384]] The slowest run took 213.01 times longer than the fastest. This could mean that an intermediate result is being cached. 1 loop, best of 3: 3.97 ms per loop 1 [[ 0.14951426]] 100 loops, best of 3: 5.71 ms per loop 1 [[ 0.14951415]]
presented by Abreu, C. R. A., 2010 - COBEQ as "a mdified version of Ohanomah e Thompson (1984)"
#Estimativas Iniciais
def iguess(Z,MODEL):
#print(Z)
gama0 = MODEL(Z/np.sum(Z))
Xlin = 1/(gama0.T @ Z) * (gama0*Z)
gamalin = MODEL(Xlin/np.sum(Xlin))
Klin = (1./gamalin) * gama0
KlinMin = np.min(Klin)
KlinMax = np.max(Klin)
NT = (KlinMax - KlinMin)
nI0 = (KlinMax - Klin) * Z
nII0 = (Klin-KlinMin) * Z
beta0 = np.array([[ 1/(KlinMax-KlinMin) * np.sum(nII0) ]])
###
XI0 = 1/(np.sum(nI0)) * nI0
gamaI0=MODEL(XI0/np.sum(XI0))
XII0 = 1/(np.sum(nII0)) * nII0
gamaII0=MODEL(XII0/np.sum(XII0))
K0 = (gamaI0)/(gamaII0)
#K0 = (XII0+1e-5)/(XI0+1e-5)
#print(XI0, XII0, gamaI0, gamaII0, K0)
return beta0, K0
#test
#exact ad-hoc NRTL solution
#Ztest = np.array([[ 0.10526316, 0.84210526, 0.05263158]]).T
#XItest = np.array([[ 0.02616646, 0.93009914, 0.0437344 ]]).T
#XIItest = np.array([[ 0.51520944, 0.38604626, 0.0987443 ]]).T
#BETAtest = np.array([[ 0.16173772]])
#appproximate
Ztest = np.array([[ 0.1, 0.85, 0.05]]).T
#XItest = np.array([[ 0.025, 0.93, 0.045 ]]).T
#XIItest = np.array([[ 0.5, 0.39, 0.11 ]]).T
igtest = iguess(Ztest,MODEL)
BETAtest = igtest[0]
Ktest=igtest[1]
beta0=BETAtest
k0=Ktest
Z=Ztest
MODEL = lambda x: Gamma(T,x,alpha,A)
ans=ELLflash_explicit(Z,beta0,k0,MODEL,True)
print('at T', T, 'Z', Z.T)
print('with ig of k0', k0.T,'beta0', beta0)
print('solution is')
print('xi',ans[0].T)
print('xii',ans[1].T)
print('beta',ans[2])
print('converged?',ans[3], '-- 1 means yes')
at T 343 Z [[ 0.1 0.85 0.05]] with ig of k0 [[ 53.51704507 0.08210211 1.89504205]] beta0 [[ 0.10903288]] solution is xi [ 0.02558676 0.93241488 0.04199836] xii [ 0.52333357 0.38114549 0.09552094] beta [[ 0.14950019]] converged? 1 -- 1 means yes
plt.plot(stdv(ans[6][0][:ans[5]]))
plt.plot(stdv(ans[6][1][:ans[5],0]))
plt.plot(stdv(ans[6][1][:ans[5],1]))
plt.plot(stdv(ans[6][1][:ans[5],2]))
plt.xlim(0,ans[5])
plt.show()
nspace = 20
Z2 = np.linspace(0,1,nspace)
Z3 = np.linspace(0,1,nspace)
Z2[0]=1e-9
Z3[0]=1e-9
Z2[-1]=1-1e-9
Z3[-1]=1-1e-9
rawZs = np.ndarray((nspace,nspace), dtype=object)
flags = np.zeros((nspace,nspace))
for i in range(nspace):
for j in range(nspace):
rawZs[i,j]=np.array([1-Z2[i]-Z3[j],Z2[i],Z3[j]])
if ( Z2[i] + Z3[j] <= 1 ):
flags[i,j] = 1
Zs = rawZs[np.where(flags==1)]
npts = Zs.shape[0]
print(npts)
Resultados_beta = np.zeros([npts])
Resultados_conv = np.zeros([npts])
Resultados_Z = np.zeros([npts,3])
Resultados_XI = np.zeros([npts,3])
Resultados_XII = np.zeros([npts,3])
Resultados_K0 = np.ndarray(npts, dtype = object)
Resultados_K = np.ndarray(npts, dtype = object)
Resultados_beta0 = np.zeros([npts])
sResultados_beta = np.zeros([npts])
sResultados_conv = np.zeros([npts])
sResultados_Z = np.zeros([npts,3])
sResultados_XI = np.zeros([npts,3])
sResultados_XII = np.zeros([npts,3])
sResultados_K0 = np.ndarray(npts, dtype = object)
sResultados_K = np.ndarray(npts, dtype = object)
sResultados_beta0 = np.zeros([npts])
MODEL = lambda x: Gamma(T,x,alpha,A)
l=0
for k in range(npts):
beta0, K0 = iguess(np.array([Zs[k]]).T,MODEL)
#print(beta0,K0)
#print(Zs[k],beta0,K0)
ans = ELLflash(np.array([Zs[k]]).T,beta0,K0,MODEL)
if (0<ans[2] and 1>ans[2] and ans[3]!=0):
#print(ans[0][0][0],ans[0][1][0],ans[0][2][0])
#print(ans[0][0:3,0])
Resultados_XI[l,:] = ans[0][0:3,0]#np.array([ans[0][0][0],ans[0][1][0],ans[0][2][0]])
Resultados_XII[l,:] = ans[1][0:3,0] #np.array([ans[1][0][0],ans[1][1][0],ans[1][2][0]])
Resultados_beta[l] = ans[2]
Resultados_conv[l] = ans[3]
Resultados_Z[l,:] = Zs[k]
Resultados_beta0[l] = beta0[0,0]
Resultados_K0[l] = K0
Resultados_K[l] = Resultados_XII[l,:]/Resultados_XI[l,:]
l+=1
#print('PASS', Zs[k], ans[0][0:3,0], ans[1][0:3,0], ans[2])
else:
#print('FAIL', k, Zs[k])
pass
print('done',l)
ndone=l-1
210 done 33
fResultados_XI = Resultados_XI[:][np.where(Resultados_conv!=0)]
fResultados_XII = Resultados_XII[:][np.where(Resultados_conv!=0)]
fResultados_beta = Resultados_beta[np.where(Resultados_conv!=0)]
fResultados_Z = Resultados_Z[:][np.where(Resultados_conv!=0)]
# barycentric coords: (a,b,c)
a=fResultados_XI[:,0]
b=fResultados_XI[:,1]
c=fResultados_XI[:,2]
d=fResultados_XII[:,0]
e=fResultados_XII[:,1]
f=fResultados_XII[:,2]
g=fResultados_Z[:,0]
h=fResultados_Z[:,1]
ii=fResultados_Z[:,2]
# plot mole fractions results in axis of ethanol vc water
plt.axis([0,1,0,1])
plt.scatter(a,c)
plt.scatter(d,f)
ZS=np.zeros([npts,3])
for i in range(npts):
ZS[i,:]=Zs[i]
plt.scatter(ZS[:,0],ZS[:,2],marker='o',s=80, facecolors='none', edgecolors='r')
plt.scatter(g,ii,marker='*')
#UNCOMMENT THIS LINE TO PLOT THE FAILED POINT
#failedpoint= [0.42105263 , 0.52631579 , 0.05263158]
#plt.scatter(failedpoint[0],failedpoint[2],marker='d')
plt.show()
# x crit = 0.2 aprox
#refiltrar a e b LEFT and RIGHT of the crit point
crit=.2
al=a[np.where(a<=crit)]
ar=a[np.where(a>crit)]
dl=d[np.where(d<=crit)]
dr=d[np.where(d>crit)]
bl=b[np.where(a<=crit)]
br=b[np.where(a>crit)]
el=e[np.where(d<=crit)]
er=e[np.where(d>crit)]
cl=c[np.where(a<=crit)]
cr=c[np.where(a>crit)]
fl=f[np.where(d<=crit)]
fr=f[np.where(d>crit)]
Ly = np.concatenate((al,dl),axis=0)
Lx = np.concatenate((cl,fl),axis=0)
Ry = np.concatenate((ar,dr),axis=0)
Rx = np.concatenate((cr,fr),axis=0)
sLy=np.sort(Ly)
sLx=np.sort(Lx)
sRy=np.sort(Ry)
sRx=-np.sort(-Rx)
print(a.shape,b.shape,d.shape,e.shape)
print(al.shape,bl.shape,dl.shape,el.shape)
print(ar.shape,br.shape,dr.shape,er.shape)
print(Ly.shape,Ry.shape,Lx.shape,Rx.shape)
(33,) (33,) (33,) (33,) (22,) (22,) (11,) (11,) (11,) (11,) (22,) (22,) (33,) (33,) (33,) (33,)
plt.axis([0,1,0,1])
plt.scatter(Ly,Lx)
plt.scatter(Ry,Rx)
plt.plot(sLy,sLx)
plt.plot(sRy,sRx)
for i in range(0,ndone):
plt.plot([a[i],d[i]],[c[i],f[i]])
plt.scatter(g,ii, marker='*')
plt.xlabel('mol. frac. water')
plt.ylabel('mol. frac. ethanol')
plt.title('LLE at '+str(T-273)+'degC')
plt.show()
#what happended to some points:
failedpoint= [0.42105263 , 0.52631579 , 0.05263158]
ig = iguess(np.array([failedpoint]).T,MODEL)
#print('ig')
#print(ig)
ans=ELLflash_explicit(np.array([failedpoint]).T,ig[0],(ig[1]),MODEL,True)
#print(ans)
#lets study those points in a different notebook!
print(ans[2])
print(ans[3],'0 means did not cnverge')
plt.plot(stdv(ans[6][0][:ans[5]]))
plt.plot(stdv(ans[6][1][:ans[5],0]))
plt.plot(stdv(ans[6][1][:ans[5],1]))
plt.plot(stdv(ans[6][1][:ans[5],2]))
plt.xlim(0,ans[5])
plt.show()
[[ 1.]] 0 0 means did not cnverge
#what happended to some points:
failedpoint= [0.42105263 , 0.52631579 , 0.05263158]
ig = iguess(np.array([failedpoint]).T,MODEL)
#print('ig')
#print(ig)
ans=ELLflash_explicit(np.array([failedpoint]).T,ig[0],(ig[1])/2,MODEL,True)
#print('ans')
#print(ans)
#print(ans[1])
print(ig[0],'vs',ans[2][0,0])
print(ig[1].T,'vs',(ans[1][:]/ans[0][:]))
print(ans[2])
print(ans[3],'1 means cnverge')
#lets study those points in a different notebook!
plt.plot(stdv(ans[6][0][:ans[5]]))
plt.plot(stdv(ans[6][1][:ans[5],0]))
plt.plot(stdv(ans[6][1][:ans[5],1]))
plt.plot(stdv(ans[6][1][:ans[5],2]))
plt.xlim(0,ans[5])
plt.xlabel(r'iterations')
plt.ylabel(r'variables')
plt.title(r'iterations history')
plt.savefig('9.png')
plt.show()
[[ 0.79986453]] vs 0.310193764201 [[ 0.41875185 4.62153782 1.02404211]] vs [ 0.03473608 2.8461674 0.41370058] [[ 0.31019376]] 1 1 means cnverge
testZ2 = 0.52631579
npts=100
Zs = np.ndarray(npts, dtype=object)
i=0
for testZ1 in np.linspace(0,testZ2,npts):
testZ3 = 1-testZ1-testZ2
Zs[i] = np.asarray([testZ1,testZ2,testZ3])
i+=1
Resultados_beta = np.zeros([npts])
Resultados_conv = np.zeros([npts])
Resultados_Z = np.zeros([npts,3])
Resultados_XI = np.zeros([npts,3])
Resultados_XII = np.zeros([npts,3])
Resultados_K0 = np.ndarray(npts, dtype = object)
Resultados_K = np.ndarray(npts, dtype = object)
Resultados_beta0 = np.zeros([npts])
sResultados_beta = np.zeros([npts])
sResultados_conv = np.zeros([npts])
sResultados_Z = np.zeros([npts,3])
sResultados_XI = np.zeros([npts,3])
sResultados_XII = np.zeros([npts,3])
sResultados_K0 = np.ndarray(npts, dtype = object)
sResultados_K = np.ndarray(npts, dtype = object)
sResultados_beta0 = np.zeros([npts])
MODEL = lambda x: Gamma(T,x,alpha,A)
l=0
ll=0
for k in range(npts):
beta0, K0 = iguess(np.array([Zs[k]]).T,MODEL)
#print(beta0,K0)
#print(Zs[k],beta0,K0)
ans = ELLflash(np.array([Zs[k]]).T,beta0,K0,MODEL)
if (0<ans[2] and 1>ans[2] and ans[3]!=0):
#print(ans[0][0][0],ans[0][1][0],ans[0][2][0])
#print(ans[0][0:3,0])
Resultados_XI[l,:] = ans[0][0:3,0]#np.array([ans[0][0][0],ans[0][1][0],ans[0][2][0]])
Resultados_XII[l,:] = ans[1][0:3,0] #np.array([ans[1][0][0],ans[1][1][0],ans[1][2][0]])
Resultados_beta[l] = ans[2]
Resultados_conv[l] = ans[3]
Resultados_Z[l,:] = Zs[k]
Resultados_beta0[l] = beta0[0,0]
Resultados_K0[l] = K0
Resultados_K[l] = Resultados_XII[l,:]/Resultados_XI[l,:]
l+=1
#print('PASS', Zs[k], ans[0][0:3,0], ans[1][0:3,0], ans[2])
else:
#print('FAIL', k, Zs[k])
pass
sResultados_conv[ll] = ans[3]
sResultados_Z[ll,:] = Zs[k]
sResultados_beta0[ll] = beta0[0,0]
sResultados_K0[ll] = K0
if ans[3]==1:
sResultados_beta[ll] = ans[2]
sResultados_K[ll] = ans[1][0:3,0]/ans[0][0:3,0]
else:
sResultados_beta[ll] = ans[2]*0
sResultados_K[ll] = ans[1][0:3,0]/ans[0][0:3,0]*0
ll+=1
print('done',l)
ndone=l-1
/home/segtovichisv/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:59: RuntimeWarning: invalid value encountered in true_divide
done 40
iterable = (sResultados_conv[l] for l in range(66,ll))
axis = np.fromiter(iterable, np.float)
print(axis)
plt.plot((axis+1),label='conv')
iterable = (sResultados_beta[l] for l in range(0,ll))
axis = np.fromiter(iterable, np.float)
print(axis)
FIRST_CONV = 66
iterable = (sResultados_beta[l] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.plot(axis,label='solution')
iterable = (sResultados_beta0[l] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.plot(axis,label='ig')
'''[[ 0.79986453]] vs 0.310193764201
[[ 0.41875185 4.62153782 1.02404211]] vs [ 0.03473608 2.8461674 0.41370058]'''
plt.scatter([13],[0.79986453],marker='*',label='bad_ig')
plt.scatter([13],[0.310193764201 ],label='failed_sol')
plt.scatter([13],[1-0.310193764201 ],marker='d',label='comp_failed_sol')
iterable = (sResultados_beta[l] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.plot(1-axis,label='comp_solution')
iterable = (sResultados_beta0[l] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.plot(1-axis,label='comp_ig')
plt.legend()
plt.show()
[ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [ 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 1. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0.9473989 0.92386607 0.90278913 0.88372061 0.86632151 0.85031202 0.83547414 0.82162759 0.80857401 0.79616533 0.78472324 0.77401123 0.76361938 0.75360881 0.74397174 0.73462545 0.72556555 0.71673969 0.70811277 0.69971917 0. 0.31691301 0.32501899 0.33303922 0.34100881 0.34891757 0.35678734 0.3645978 0.37235028 0.38001389 0.38755401 0.39492909 0.40207962 0.40894618 0.41546994 0.42160072 0.42729509 0.43252767 0.43729026 0.44158459 0.44542847]
iterable = (sResultados_conv[l] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
print(axis)
plt.semilogy(1000*(axis)+100,label='conv')
iterable = (sResultados_K[l][0] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.semilogy(axis,label='solution')
iterable = (sResultados_K0[l][0] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.semilogy(axis,label='ig')
'''[[ 0.79986453]] vs 0.310193764201
[[ 0.41875185 4.62153782 1.02404211]] vs [ 0.03473608 2.8461674 0.41370058]'''
plt.scatter([13],[0.41875185],marker='*',label='bad_ig')
plt.scatter([13],[0.03473608],label='failed_sol')
plt.scatter([13],[1/0.03473608 ],marker='d',label='inv_failed_sol')
iterable = (sResultados_K[l][0] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.plot(1/axis,label='inv_solution')
iterable = (sResultados_K0[l][0] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.plot(1/axis,label='inv_ig')
plt.legend()
plt.show()
[ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
/home/segtovichisv/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:25: RuntimeWarning: divide by zero encountered in true_divide
iterable = (sResultados_conv[l] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
print(axis)
plt.semilogy(1000*(axis)+100,label='conv')
iterable = (sResultados_K[l][1] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.semilogy(axis,label='solution')
iterable = (sResultados_K0[l][1] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.semilogy(axis,label='ig')
'''[[ 0.79986453]] vs 0.310193764201
[[ 0.41875185 4.62153782 1.02404211]] vs [ 0.03473608 2.8461674 0.41370058]'''
plt.scatter([13],[4.62153782],marker='*')
plt.scatter([13],[2.8461674 ],label='failed_sol')
plt.scatter([13],[1/2.8461674 ],marker='d',label='inv_failed_sol')
iterable = (sResultados_K[l][1] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.plot(1/axis,label='inv_solution')
iterable = (sResultados_K0[l][1] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
plt.plot(1/axis,label='inv_ig')
plt.legend()
plt.show()
[ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
/home/segtovichisv/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:24: RuntimeWarning: divide by zero encountered in true_divide
xiz = np.arange(FIRST_CONV,ll)
iterable = (sResultados_conv[l] for l in range(FIRST_CONV,ll))
caxis = np.fromiter(iterable, np.float)
#print(axis)
zaxis = np.ones(ll-FIRST_CONV)
maxis = np.ma.masked_where(caxis == 0, np.ones(ll-FIRST_CONV))
naxis = np.ma.masked_where(caxis == 1, np.ones(ll-FIRST_CONV))
mxiz = np.ma.masked_where(caxis == 0, xiz)
plt.semilogy()
plt.scatter(xiz,50*(zaxis),marker='o', s=80, facecolors='none', edgecolors='r')
plt.scatter(mxiz,50*(maxis),marker='x',label='converges')
plt.scatter(mxiz,50*(naxis),marker='o', s=80, facecolors='none', edgecolors='r',label='diverges')
iterable = (sResultados_K[l][2] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
maxis = np.ma.masked_where(caxis == 0, axis)
plt.semilogy(mxiz,maxis,label='solution')
iterable = (sResultados_K0[l][2] for l in range(FIRST_CONV,ll))
axis = np.fromiter(iterable, np.float)
maxis = np.ma.masked_where(caxis == 0, axis)
plt.semilogy(mxiz,maxis,label='initial guess')
'''[[ 0.79986453]] vs 0.310193764201
[[ 0.41875185 4.62153782 1.02404211]] vs [ 0.03473608 2.8461674 0.41370058]'''
plt.scatter([13+66],[1.02404211],marker='*')
plt.scatter([13+66],[0.41370058 ],label='failed_sol')
plt.scatter([13+66],[1/0.41370058 ],marker='d',label='inverse missed solution')
iterable = (sResultados_K[l][2] for l in range(FIRST_CONV,ll))
maxis = np.ma.masked_where(caxis == 0, axis)
axis = np.fromiter(iterable, np.float)
plt.semilogy(mxiz,1/maxis,label='inverse solution')
iterable = (sResultados_K0[l][2] for l in range(FIRST_CONV,ll))
maxis = np.ma.masked_where(caxis == 0, axis)
axis = np.fromiter(iterable, np.float)
plt.semilogy(mxiz,1/maxis,label='inverse initial guess')
#for i in range(34):
# print(mxiz[i],maxis[i])
plt.xlim(FIRST_CONV,ll-1)
plt.ylim(1e-1,1e2)
#plt.legend(loc=0)
plt.title(r'initial guess $\times$ solution')
plt.ylabel(r'$K_3$')
plt.xlabel(r'input index')
plt.savefig('fig7.png')
plt.show()
/home/segtovichisv/anaconda3/lib/python3.5/site-packages/matplotlib/scale.py:93: RuntimeWarning: invalid value encountered in less_equal mask = a <= 0.0
print( sResultados_K0[l][:] )
[[ 9.845702 ] [ 0.14959554] [ 1.21599902]]
### given K3 = 1, beta = .5
#grid K1 and K2 between .01 and 100
nspace = 5
K1s = np.logspace(-2,2,nspace)
K2s = np.logspace(-2,2,nspace)
Z = [0.42105263 , 0.52631579 , 0.05263158]
K3=1
beta0=np.array([[.5]])
###
npts = nspace*nspace
histK1 = np.ndarray([nspace,nspace], dtype=object)
histK2 = np.ndarray([nspace,nspace], dtype=object)
for k in range(nspace):
for l in range(nspace): # for each point in the grid,
# calculate the initial guess
K0 = np.asarray([[K1s[k], K2s[l], K3]]).T
#then calculate the flash
ans = ELLflash_explicit(np.array([Z]).T,beta0,K0,MODEL,True)
histK1[k,l] = np.fromiter( (ans[6][1][i,0] for i in range(100) ), dtype=np.float)
histK2[k,l] = np.fromiter( (ans[6][1][i,1] for i in range(100) ), dtype=np.float)
/home/segtovichisv/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:52: RuntimeWarning: invalid value encountered in true_divide
print(histK1[0,0])
print(histK2[0,0])
print(nspace)
plt.plot(histK1[0,0])
[ 0.01 0.01 0.86792087 1.09588658 1.09336134 1.08841374 1.08369072 1.07926302 1.0751117 1.07121577 1.06755599 1.06411489 1.0608766 1.05782663 1.05495177 1.05223993 1.04968005 1.04726197 1.04497636 1.04281465 1.04076891 1.03883186 1.03699674 1.03525731 1.0336078 1.03204283 1.03055742 1.02914696 1.02780711 1.02653387 1.02532349 1.02417247 1.02307753 1.02203563 1.02104389 1.02009964 1.01920035 1.01834366 1.01752734 1.01674932 1.01600763 1.01530042 1.01462595 1.01398257 1.01336874 1.01278299 1.01222394 1.01169029 1.01118081 1.01069432 1.01022972 1.00978597 1.00936208 1.0089571 1.00857015 1.00820037 1.00784698 1.0075092 1.00718632 1.00687764 1.00658252 1.00630034 1.00603051 1.00577246 1.00552566 1.0052896 1.00506381 1.00484781 1.00464118 1.00444349 1.00425434 1.00407336 1.00390018 1.00373446 1.00357586 1.00342408 1.00327882 1.00313979 1.00300672 1.00287934 1.00275741 1.00264069 1.00252895 1.00242199 1.00231958 1.00222154 1.00212767 1.00203779 1.00195174 1.00186935 1.00179046 1.00171491 1.00164257 1.0015733 1.00150697 1.00144345 1.00138262 1.00132436 1.00126857 1.00121514] [ 0.01 0.01 0.72458545 0.91732793 0.92380401 0.92800565 0.93187024 0.93548513 0.93887286 0.94205114 0.94503589 0.94784151 0.95048113 0.95296665 0.95530897 0.95751803 0.95960292 0.96157198 0.96343286 0.96519262 0.96685774 0.96843419 0.9699275 0.97134279 0.97268478 0.97395785 0.9751661 0.97631329 0.97740295 0.97843836 0.97942259 0.98035848 0.98124871 0.98209578 0.98290201 0.98366961 0.98440061 0.98509696 0.98576045 0.98639279 0.98699558 0.98757033 0.98811845 0.98864129 0.9891401 0.98961607 0.99007033 0.99050394 0.99091791 0.99131318 0.99169066 0.99205119 0.99239558 0.9927246 0.99303896 0.99333937 0.99362646 0.99390086 0.99416316 0.99441391 0.99465364 0.99488286 0.99510205 0.99531167 0.99551214 0.99570388 0.99588729 0.99606274 0.99623058 0.99639115 0.99654479 0.99669179 0.99683245 0.99696705 0.99709587 0.99721914 0.99733712 0.99745005 0.99755813 0.99766159 0.99776062 0.99785541 0.99794616 0.99803304 0.99811621 0.99819584 0.99827207 0.99834506 0.99841495 0.99848187 0.99854594 0.99860729 0.99866604 0.9987223 0.99877617 0.99882776 0.99887716 0.99892448 0.99896979 0.99901318] 5
[<matplotlib.lines.Line2D at 0x7f9f44aa2278>]
for k in range(nspace):
for l in range(nspace): # for each point in the grid,
plt.loglog(histK1[k,l],histK2[k,l])
plt.title(r'convergence path')
plt.xlabel(r'$K_1$')
plt.ylabel(r'$K_2$')
plt.savefig('fig8.png')
plt.show()
## para cada k1 ca2, calc RES, plot contourf