Commit 31874bd1 authored by MARTIN Antoine's avatar MARTIN Antoine

add optimization on k0 + curves k_insulating/k_hom

parent 376a9274
......@@ -23,57 +23,93 @@ def rec_fill(t,N):
rec_fill(t[ii],N)
def A_CA(tau_field,N,k0,ki):
def Gamma_Tau(tau_field,N):
d=len(N)
k0=np.eye(d)
filter_level=2
Eps_f=np.zeros(tuple(N)+(d,d))
for j in range(d):
tau = np.zeros(tuple(N)+(d,))
for i in np.ndindex(tuple(N)):
tau[i]=tau_field[i][:,j]
it=0
while it<=0:
field,field_fourier,fft,ifft,tupleK,frequencies,filters,tupleFilters = initializeGreen(N,filter_level=filter_level)
nullifyField(field,tuple(N))
Eps_field = np.zeros(tuple(N)+(d,))
for i in np.ndindex(tuple(N)):
Eps_field[i+(j,)]=1.
for i in np.ndindex(tuple(N)+(d,)):
field[i] = tau[i]
field=operate_field(field,field_fourier,fft,ifft,tupleK,N,frequencies,filter_level,filters,tupleFilters,k0)
for i in np.ndindex(tuple(N)):
Eps_field[i]+=field[i]
for i in np.ndindex(tuple(N)):
tau[i]=np.dot(ki[i]-k0,Eps_field[i])
it+=1
print(Eps_field[0,0,0],it)
# plt.figure()
# plt.imshow(Eps_field[:,:,0],origin='lower')
# plt.show()
# plt.close()
for i in np.ndindex(tuple(N)+(d,)):
Eps_f[i+(j,)]=Eps_field[i]
return Eps_f
def Sig(Eps_field,A_field,k0,N,ki,khh):
d=len(N)
Sigma_field=np.zeros(tuple(N)+(d,d))
field,field_fourier,fft,ifft,tupleK,frequencies,filters,tupleFilters = initializeGreen(N,filter_level=filter_level)
nullifyField(field,tuple(N))
Eps_field = np.zeros(tuple(N)+(d,))
for i in np.ndindex(tuple(N)+(d,)):
field[i] = tau_field[i]
field=operate_field(field,field_fourier,fft,ifft,tupleK,N,frequencies,filter_level,filters,tupleFilters,k0)
for i in np.ndindex(tuple(N)):
Sigma_field[i]=np.dot(np.dot(k0,Eps_field[i]-A_field[i])+np.dot(ki[i],A_field[i]),np.linalg.inv(khh))
Eps_field[i]+=field[i]
# plt.figure()
# plt.imshow(Sigma_field[:,:,0,0],origin='lower')
# plt.imshow(Eps_field[:,:,0],origin='lower')
# plt.show()
# plt.close()
return Sigma_field
return -Eps_field
# def Sig(Eps_field,A_field,k0,N,ki,khh):
# d=len(N)
# Sigma_field=np.zeros(tuple(N)+(d,d))
# for i in np.ndindex(tuple(N)):
# Sigma_field[i]=np.dot(np.dot(k0,Eps_field[i]-A_field[i])+np.dot(ki[i],A_field[i]),np.linalg.inv(khh))
# # plt.figure()
# # plt.imshow(Sigma_field[:,:,0,0],origin='lower')
# # plt.show()
# # plt.close()
# return Sigma_field
def k_energetic(tau_field,A_field,N,k0,ki,khh):
def alpha_CA(eps_HH,ki,N,j):
d=len(N)
Eps_field=A_CA(tau_field,N,k0,ki)
Sigma_field=Sig(Eps_field,A_field,k0,N,ki,khh)
kCA=np.zeros((d,d))
kSAinv=np.zeros((d,d))
sig_HH=np.zeros(tuple(N)+(d,))
Id=np.eye(d)
for i in np.ndindex(tuple(N)):
sig_HH[i]=np.dot(ki[i],eps_HH[i])
Gamma_eps_HH=Gamma_Tau(eps_HH,N)
eps_HH_inf=np.zeros(tuple(N)+(d,))
for i in np.ndindex(tuple(N)):
eps_HH_inf[i]=Id[:,j]+Gamma_eps_HH[i]
eps_HH_star=-Gamma_Tau(sig_HH,N)
a1=0.
a2=0.
for i in np.ndindex(tuple(N)):
a1+=np.dot(eps_HH_star[i],np.dot(ki[i],eps_HH_inf[i]))/N[0]**d
a2+=np.dot(eps_HH_star[i],np.dot(ki[i],eps_HH_star[i]))/N[0]**d
return -a2/a1,eps_HH_star,eps_HH_inf
def alpha_SA(eps_HH,eps_inf,eps_star,ki,N):
d=len(N)
sig_HH=np.zeros(tuple(N)+(d,))
for i in np.ndindex(tuple(N)):
kCA+=np.dot((Eps_field[i]).transpose(),np.dot(ki[i],Eps_field[i]))/N[0]**d
kSAinv+=np.dot((Sigma_field[i]).transpose(),np.dot(np.linalg.inv(ki[i]),Sigma_field[i]))/N[0]**d
return kCA,np.linalg.inv(kSAinv)
sig_HH[i]=np.dot(ki[i],eps_HH[i])
delta_eps=np.zeros(tuple(N)+(d,))
sig_star=np.zeros(tuple(N)+(d,))
for i in np.ndindex(tuple(N)):
delta_eps[i]=eps_inf[i]-eps_HH[i]
sig_star[i]=eps_star[i]+sig_HH[i]
a1=0.
a2=0.
for i in np.ndindex(tuple(N)):
a1+=np.dot(delta_eps[i],np.dot(np.linalg.inv(ki[i]),delta_eps[i]))/N[0]**d
a2+=np.dot(delta_eps[i],np.dot(np.linalg.inv(ki[i]),sig_star[i]))/N[0]**d
return -a2/a1,delta_eps,sig_star
def k_energetic(A_field,N,ki,khh):
d=len(N)
kCA=np.zeros((d,d))
kinv=np.zeros((d,d))
kSA=np.zeros((d,d))
for j in range(d):
eps_HH=np.zeros(tuple(N)+(d,))
for i in np.ndindex(tuple(N)):
eps_HH[i]=A_field[i][:,j]
al_CA,eps_HH_star,eps_HH_inf=alpha_CA(eps_HH,ki,N,j)
print(al_CA)
eps_CA=np.zeros(tuple(N)+(d,))
sig_SA=np.zeros(tuple(N)+(d,))
for i in np.ndindex(tuple(N)):
eps_CA[i]=eps_HH_inf[i]+1/al_CA*eps_HH_star[i]
al_SA,delta_eps,sig_star=alpha_SA(eps_HH,eps_HH_inf,eps_HH_star,ki,N)
print(al_SA)
for i in np.ndindex(tuple(N)):
sig_SA[i]=al_SA*delta_eps[i]+sig_star[i]
for i in np.ndindex(tuple(N)):
kCA[j,j]+=np.dot(eps_CA[i],np.dot(ki[i],eps_CA[i]))/N[0]**d
kinv[j,j]+=np.dot(sig_SA[i],np.dot(np.linalg.inv(ki[i]),sig_SA[i]))/N[0]**d
kSA=np.dot(khh,np.dot(np.linalg.inv(kinv),khh.transpose()))
return kCA,kSA
......@@ -26,6 +26,7 @@ def Ref(tau_field,N,k0,ki):
d=len(N)
filter_level=2
Eps_f=np.zeros(tuple(N)+(d,d))
prevEpsf=np.zeros(tuple(N)+(d,d))
for j in range(d):
tau = np.zeros(tuple(N)+(d,))
for i in np.ndindex(tuple(N)):
......@@ -34,7 +35,7 @@ def Ref(tau_field,N,k0,ki):
err=1.
preverr=1.
krefj=k0[:,j]
while err>1e-3 or preverr>1e-3:
while err>1e-4 or preverr>1e-4:
print(it)
field,field_fourier,fft,ifft,tupleK,frequencies,filters,tupleFilters = initializeGreen(N,filter_level=filter_level)
nullifyField(field,tuple(N))
......@@ -56,15 +57,16 @@ def Ref(tau_field,N,k0,ki):
for i in np.ndindex(tuple(N)):
krefj+=np.dot(ki[i],Eps_field[i])/N[0]**d
err=np.abs((krefj[j]-prevkrefjj)/prevkrefjj)
print(err)
print(krefj)
for i in np.ndindex(tuple(N)+(d,)):
Eps_f[i+(j,)]=Eps_field[i]
prevEpsf[i+(j,)]=prevEps[i]
for i in np.ndindex(tuple(N)):
Sig_field[i]=np.dot(ki[i],Eps_field[i])
#plt.figure()
#plt.imshow(Sig_field[:,:,0],origin='lower')
#plt.show()
#plt.close()
plt.figure()
plt.imshow(Sig_field[:,:,0],origin='lower')
plt.show()
plt.close()
kref=np.zeros((d,d))
for i in np.ndindex(tuple(N)):
kref+=np.dot(ki[i],Eps_f[i])/N[0]**d
......
......@@ -230,14 +230,14 @@ N0 = 128
d=2
N=N0*np.ones(d,dtype=np.int64)
radius = 8
a=14
b=3
#shape = {2:"disc",3:"sphere"}
shape = {2:"ellipse",3:"ellipsoïd"}
a=30
b=1
shape = {2:"disc",3:"sphere"}
#shape = {2:"ellipse",3:"ellipsoïd"}
V = pi*radius**2 if d==2 else 4*pi/3*radius**3
k1 = 1.
k2 = 0.001
for model in ["A"]:
for model in ["A","B"]:
if model == "A":
phaseProps = np.array([k1,k2],dtype=np.float64)
elif model == "B":
......@@ -247,16 +247,16 @@ for model in ["A"]:
kCA_list=[]
kSA_list=[]
kref_list=[]
for target_phi in range(2,99,5):
#ncenter = int(-np.log(1-target_phi/100.)*N0**d/V)
ncenter = int(-np.log(1-target_phi/100.)*N0**d/np.pi/a/b)
for target_phi in range(5,99,5):
ncenter = int(-np.log(1-target_phi/100.)*N0**d/V)
#ncenter = int(-np.log(1-target_phi/100.)*N0**d/np.pi/a/b)
centers = np.random.rand(ncenter,d)*N0
#radii = np.ones(ncenter)*radius
radii = np.ones(ncenter)*radius
ai=np.ones(ncenter)*a
bi=np.ones(ncenter)*b
thetai = np.random.rand(ncenter,)*np.pi
#microstructure = booleanSpheres(centers,radii,N)
microstructure = booleanEllipses(centers,ai,bi,thetai,N)
microstructure = booleanSpheres(centers,radii,N)
#microstructure = booleanEllipses(centers,ai,bi,thetai,N)
#microstructure = booleanSpheres(np.array([[0,0],[N0/3,N0/3],[3*N0/4,N0/2]]),[N0/3,N0/3,N0/3],N)
mapStrides = np.array([N0**(d-i-1) for i in range(d)],dtype=np.int64)
phasemap = microstructure.phasemap.flatten() #TODO : flatten to allow scalar index (see numba issue with tuple index)
......@@ -264,31 +264,30 @@ for model in ["A"]:
frac.append(phi)
print("target phi %f, achieved phi %f, num. centers %d" % (target_phi/100.,phi,ncenter))
#microstructure = square(N0/2)
# plt.figure()
# plt.imshow(microstructure.phasemap,origin='lower')
# plt.show()
k0,khh,t_A_hh = hierachicHomogenize(N,phasemap,phaseProps)
k0=((1-phi)*phaseProps[0]+phi*phaseProps[1])*np.eye(len(N))
ki=np.zeros(tuple(N)+(d,d))
for i in np.ndindex(tuple(N)):
if microstructure.phasemap[i]==1:
ki[i]=phaseProps[1]*np.eye(d)
else:
ki[i]=phaseProps[0]*np.eye(d)
tau_hh=np.zeros(tuple(N)+(d,d))
for i in np.ndindex(tuple(N)):
tau_hh[i]=np.dot(ki[i]-k0,t_A_hh[i])
# tau_hh=np.zeros(tuple(N)+(d,d))
# for i in np.ndindex(tuple(N)):
# tau_hh[i]=np.dot(ki[i]-k0,t_A_hh[i])
# t11=np.sum(tau_hh[:,:,0,0])/N0**d
# t12=np.sum(tau_hh[:,:,0,1])/N0**d
# t21=np.sum(tau_hh[:,:,1,0])/N0**d
# t22=np.sum(tau_hh[:,:,1,1])/N0**d
# print(k0+np.array([[t11,t12],[t21,t22]]))
kCA,kSA=k_energetic(tau_hh,t_A_hh,N,k0,ki,khh)
kCA,kSA=k_energetic(t_A_hh,N,ki,khh)
k0=np.mean(phaseProps)*np.eye(len(N))
kref=Ref(tau_hh,N,k0,ki)
print(khh,kCA,kSA,kref)
tau=np.zeros(tuple(N)+(d,d))
for i in np.ndindex(tuple(N)):
tau[i]=np.dot(ki[i]-k0,np.eye(d))
kref=Ref(tau,N,k0,ki)
print(khh,kCA,kSA)#,kref)
khom.append(khh)
kCA_list.append(kCA)
kSA_list.append(kSA)
......@@ -299,46 +298,46 @@ for model in ["A"]:
# plt.imshow(tau_hh[:,:,1,0],origin='lower')
# plt.show()
# plt.close()
savefile=open("boolean_%s_contrast%dmodel%s.txt" % (shape[d],int(k1/k2),model), 'ab')
savefile=open("boolean_%s_contrast%dmodel%s2.txt" % (shape[d],int(k1/k2),model), 'ab')
np.savetxt(savefile,np.array([[phi,np.trace(khh)/d]]))
savefile.close()
savefile=open("boolean_%s_contrast%dmodel%sSA.txt" % (shape[d],int(k1/k2),model), 'ab')
savefile=open("boolean_%s_contrast%dmodel%sSA2.txt" % (shape[d],int(k1/k2),model), 'ab')
np.savetxt(savefile,np.array([[phi,np.trace(kSA)/d]]))
savefile.close()
savefile=open("boolean_%s_contrast%dmodel%sCA.txt" % (shape[d],int(k1/k2),model), 'ab')
savefile=open("boolean_%s_contrast%dmodel%sCA2.txt" % (shape[d],int(k1/k2),model), 'ab')
np.savetxt(savefile,np.array([[phi,np.trace(kCA)/d]]))
savefile.close()
savefile=open("boolean_%s_contrast%dmodel%sref.txt" % (shape[d],int(k1/k2),model), 'ab')
np.savetxt(savefile,np.array([[phi,np.trace(kref)/d]]))
savefile.close()
#savefile=open("boolean_%s_contrast%dmodel%sref2.txt" % (shape[d],int(k1/k2),model), 'ab')
#np.savetxt(savefile,np.array([[phi,np.trace(kref)/d]]))
#savefile.close()
Kmoritanaka = lambda x,k1,k2,d : ((1-x)*k1+x*k2*d*k1/((d-1)*k1+k2))/((1-x)+x*d*k1/((d-1)*k1+k2))
eqKselfconsistent = lambda k,x,k1,k2,d : 1e6 if k <0 else ((1-x)*(k1-k)*k*d/((d-1)*k+k1)+x*(k2-k)*k*d/((d-1)*k+k2))
Kselfconsistent= lambda x,k1,k2,d : fsolve(eqKselfconsistent,(k1+k2)/2,args=(x,k1,k2,d))
x=np.linspace(0,1,101)
plt.plot(x,Kmoritanaka(x,k1,k2,d),label="Hashin-Shtrikman")
plt.plot(x,Kmoritanaka(1-x,k2,k1,d),label="Hashin-Shtrikman")
plt.plot(x,[Kselfconsistent(f2,k1,k2,d) for f2 in x],label="Self-Consistent")
# x=np.linspace(0,1,101)
# plt.plot(x,Kmoritanaka(x,k1,k2,d),label="Hashin-Shtrikman")
# plt.plot(x,Kmoritanaka(1-x,k2,k1,d),label="Hashin-Shtrikman")
# plt.plot(x,[Kselfconsistent(f2,k1,k2,d) for f2 in x],label="Self-Consistent")
resA = np.loadtxt("boolean_%s_contrast%dmodel%s.txt" % (shape[d],int(k1/k2),"A"))
plt.plot(resA[:,0],resA[:,1],'k.',label="hierarchical num. self-consistent A")
resASA = np.loadtxt("boolean_%s_contrast%dmodel%sSA.txt" % (shape[d],int(k1/k2),"A"))
plt.plot(resASA[:,0],resASA[:,1],'k.',color='red',label="hierarchical num. self-consistent A SA")
resACA = np.loadtxt("boolean_%s_contrast%dmodel%sCA.txt" % (shape[d],int(k1/k2),"A"))
plt.plot(resACA[:,0],resACA[:,1],'k.',color='green',label="hierarchical num. self-consistent A CA")
resrefA = np.loadtxt("boolean_%s_contrast%dmodel%sref.txt" % (shape[d],int(k1/k2),"A"))
plt.plot(resrefA[:,0],resrefA[:,1],'k+',color='red',label="A ref")
# resA = np.loadtxt("boolean_%s_contrast%dmodel%s2.txt" % (shape[d],int(k1/k2),"A"))
# plt.plot(resA[:,0],resA[:,1],'k.',label="hierarchical num. self-consistent A")
# resASA = np.loadtxt("boolean_%s_contrast%dmodel%sSA2.txt" % (shape[d],int(k1/k2),"A"))
# plt.plot(resASA[:,0],resASA[:,1],'k.',color='red',label="hierarchical num. self-consistent A SA")
# resACA = np.loadtxt("boolean_%s_contrast%dmodel%sCA2.txt" % (shape[d],int(k1/k2),"A"))
# plt.plot(resACA[:,0],resACA[:,1],'k.',color='green',label="hierarchical num. self-consistent A CA")
# resrefA = np.loadtxt("boolean_%s_contrast%dmodel%sref2.txt" % (shape[d],int(k1/k2),"A"))
# plt.plot(resrefA[:,0],resrefA[:,1],'k+',color='red',label="A ref")
plt.legend(loc="upper right")
plt.title(r"boolean %s contrast = %d" % (shape[d],int(k1/k2)))
plt.xlabel("volume fraction of insulating phase")
plt.ylabel(r"$K_{hom} / K_{conductive \, phase}$")
plt.savefig("boolean_%s_contrast%d_model%s.pdf" % (shape[d],int(k1/k2),"A"))
plt.show()
# plt.legend(loc="upper right")
# plt.title(r"boolean %s contrast = %d" % (shape[d],int(k1/k2)))
# plt.xlabel("volume fraction of insulating phase")
# plt.ylabel(r"$K_{hom} / K_{conductive \, phase}$")
# plt.savefig("boolean_%s_contrast%d_model%s2.pdf" % (shape[d],int(k1/k2),"A"))
# plt.show()
plt.close()
# plt.close()
# x=np.linspace(0,1,101)
......@@ -346,19 +345,68 @@ plt.close()
# plt.plot(x,Kmoritanaka(1-x,k2,k1,d),label="Hashin-Shtrikman")
# plt.plot(x,[Kselfconsistent(f2,k1,k2,d) for f2 in x],label="Self-Consistent")
# resB = np.loadtxt("boolean_%s_contrast%dmodel%s.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resB[:,0],resB[:,1],'k.',label="hierachical num. self-consistent B")
# resBSA = np.loadtxt("boolean_%s_contrast%dmodel%sSA.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resBSA[:,0],resBSA[:,1],'k.',color='red',label="hierachical num. self-consistent B SA")
# resBCA = np.loadtxt("boolean_%s_contrast%dmodel%sCA.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resBCA[:,0],resBCA[:,1],'k.',color='green',label="hierachical num. self-consistent B CA")
# resrefB = np.loadtxt("boolean_%s_contrast%dmodel%sref.txt" % (shape[d],int(k1/k2),"B"))
# resB = np.loadtxt("boolean_%s_contrast%dmodel%s2.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resB[:,0],resB[:,1],'k.',label="hierarchical num. self-consistent B")
# resBSA = np.loadtxt("boolean_%s_contrast%dmodel%sSA2.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resBSA[:,0],resBSA[:,1],'k.',color='red',label="hierarchical num. self-consistent B SA")
# resBCA = np.loadtxt("boolean_%s_contrast%dmodel%sCA2.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resBCA[:,0],resBCA[:,1],'k.',color='green',label="hierarchical num. self-consistent B CA")
# resrefB = np.loadtxt("boolean_%s_contrast%dmodel%sref2.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resrefB[:,0],resrefB[:,1],'k+',color='red',label="B ref")
# plt.legend(loc="upper right")
# plt.title(r"boolean %s contrast = %d" % (shape[d],int(k1/k2)))
# plt.xlabel("volume fraction of insulating phase")
# plt.ylabel(r"$K_{hom} / K_{conductive \, phase}$")
# plt.savefig("boolean_%s_contrast%d_model%s.pdf" % (shape[d],int(k1/k2),"B"))
# plt.savefig("boolean_%s_contrast%d_model%s2.pdf" % (shape[d],int(k1/k2),"B"))
# plt.show()
# plt.close()
x=np.linspace(0,1,101)
plt.plot(x,k2/Kmoritanaka(x,k1,k2,d),label="Hashin-Shtrikman")
plt.plot(x,k2/Kmoritanaka(1-x,k2,k1,d),label="Hashin-Shtrikman")
plt.plot(x,[k2/Kselfconsistent(f2,k1,k2,d) for f2 in x],label="Self-Consistent")
resA = np.loadtxt("boolean_%s_contrast%dmodel%s2.txt" % (shape[d],int(k1/k2),"A"))
plt.plot(resA[:,0],k2/resA[:,1],'k.',label="hierarchical num. self-consistent A")
resASA = np.loadtxt("boolean_%s_contrast%dmodel%sSA2.txt" % (shape[d],int(k1/k2),"A"))
plt.plot(resASA[:,0],k2/resASA[:,1],'k.',color='red',label="hierarchical num. self-consistent A SA")
resACA = np.loadtxt("boolean_%s_contrast%dmodel%sCA2.txt" % (shape[d],int(k1/k2),"A"))
plt.plot(resACA[:,0],k2/resACA[:,1],'k.',color='green',label="hierarchical num. self-consistent A CA")
# resrefA = np.loadtxt("boolean_%s_contrast%dmodel%sref2.txt" % (shape[d],int(k1/k2),"A"))
# plt.plot(resrefA[:,0],k2/resrefA[:,1],'k+',color='red',label="A ref")
plt.legend(loc="upper right")
plt.title(r"boolean %s contrast = %d" % (shape[d],int(k1/k2)))
plt.xlabel("volume fraction of insulating phase")
plt.ylabel(r"$K_{insulating \, phase} / K_{hom}$")
plt.savefig("boolean_%s_contrast%d_model%s_k_insulating_phase2.pdf" % (shape[d],int(k1/k2),"A"))
plt.show()
plt.close()
plt.plot(x,k2/Kmoritanaka(x,k1,k2,d),label="Hashin-Shtrikman")
plt.plot(x,k2/Kmoritanaka(1-x,k2,k1,d),label="Hashin-Shtrikman")
plt.plot(x,[k2/Kselfconsistent(f2,k1,k2,d) for f2 in x],label="Self-Consistent")
resB = np.loadtxt("boolean_%s_contrast%dmodel%s2.txt" % (shape[d],int(k1/k2),"B"))
plt.plot(1-resB[:,0],k2/resB[:,1],'k.',label="hierachical num. self-consistent B")
resBSA = np.loadtxt("boolean_%s_contrast%dmodel%sSA2.txt" % (shape[d],int(k1/k2),"B"))
plt.plot(1-resBSA[:,0],k2/resBSA[:,1],'k.',color='red',label="hierachical num. self-consistent B SA")
resBCA = np.loadtxt("boolean_%s_contrast%dmodel%sCA2.txt" % (shape[d],int(k1/k2),"B"))
plt.plot(1-resBCA[:,0],k2/resBCA[:,1],'k.',color='green',label="hierachical num. self-consistent B CA")
# resrefB = np.loadtxt("boolean_%s_contrast%dmodel%sref2.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resrefB[:,0],k2/resrefB[:,1],'k+',color='red',label="B ref")
plt.legend(loc="upper right")
plt.title(r"boolean %s contrast = %d" % (shape[d],int(k1/k2)))
plt.xlabel("volume fraction of insulating phase")
plt.ylabel(r"$K_{insulating \, phase} / K_{hom}$")
plt.savefig("boolean_%s_contrast%d_model%s_k_insulating_phase2.pdf" % (shape[d],int(k1/k2),"B"))
plt.show()
plt.close()
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment