Commit c303da6d authored by MARTIN Antoine's avatar MARTIN Antoine

add Champs_admissibles.py and Reference.py

parent 8a32c1e8
import numpy as np
from numba import jit, float64, int64, prange, generated_jit
import matplotlib.pyplot as plt
import time
from scipy.optimize import fsolve
from greenOperatorConductionNumba import *
from microstructure import *
@jit
def nullifyField(field,N):
for n in np.ndindex(N):
field[n] = np.zeros(len(N))
@jit
def rec_fill(t,N):
d=len(N)
if t[0].shape==(d,d):
for ii in prange(N[0]):
t[ii]=np.eye(d)
else:
for ii in prange(N[0]):
rec_fill(t[ii],N)
def A_CA(tau_field,N,k0,ki):
d=len(N)
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))
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):
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))
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)
import numpy as np
from numba import jit, float64, int64, prange, generated_jit
import matplotlib.pyplot as plt
import time
from scipy.optimize import fsolve
from greenOperatorConductionNumba import *
from microstructure import *
@jit
def nullifyField(field,N):
for n in np.ndindex(N):
field[n] = np.zeros(len(N))
@jit
def rec_fill(t,N):
d=len(N)
if t[0].shape==(d,d):
for ii in prange(N[0]):
t[ii]=np.eye(d)
else:
for ii in prange(N[0]):
rec_fill(t[ii],N)
def Ref(tau_field,N,k0,ki):
d=len(N)
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
err=1.
preverr=1.
krefj=k0[:,j]
while err>1e-3 or preverr>1e-3:
print(it)
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,))
Sig_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
preverr=err
prevkrefjj=krefj[j]
krefj=np.zeros((d,))
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)
for i in np.ndindex(tuple(N)+(d,)):
Eps_f[i+(j,)]=Eps_field[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()
kref=np.zeros((d,d))
for i in np.ndindex(tuple(N)):
kref+=np.dot(ki[i],Eps_f[i])/N[0]**d
return kref
......@@ -106,9 +106,9 @@ def operate_fourier_field(x,y,tupleK,N,frequencies,filter_level,filters,tupleFil
#if not enough, remove tupleK and use scalar index to flattened x,y
def operate_field(xFourier,yFourier,fft,ifft,tupleK,N,frequencies,filter_level,filters,tupleFilters,k0):
def operate_field(x,yFourier,fft,ifft,tupleK,N,frequencies,filter_level,filters,tupleFilters,k0):
#start = time.time()
fft()
xFourier=fft(x)
#end = time.time()
#t1 = end-start
#start = time.time()
......@@ -116,7 +116,7 @@ def operate_field(xFourier,yFourier,fft,ifft,tupleK,N,frequencies,filter_level,f
#end = time.time()
#t2 = end-start
#start = time.time()
ifft()
return ifft(yFourier)
#end = time.time()
#t3 = end-start
#print("%s\t%s\t%s" % (t1,t2,t3))
......
......@@ -5,10 +5,9 @@ import time
from scipy.optimize import fsolve
from greenOperatorConductionNumba import *
#from discretizationGrid import *
#from constitutiveLaw import *
from microstructure import *
from Champs_admissibles import *
from Reference import *
......@@ -36,7 +35,7 @@ def getLocalGreenNumba(M,k0):
for i in range(d):
nullifyField(field,tuple(M))
field[tuple(d*[0]+[i])] = 1
operate_field(field_fourier,field_fourier,fft,ifft,tupleK,M,frequencies,filter_level,filters,tupleFilters,k0)
field=operate_field(field,field_fourier,fft,ifft,tupleK,M,frequencies,filter_level,filters,tupleFilters,k0)
for m in np.ndindex(tuple(d*[3])):
s = tuple([m[j]-1 for j in range(d)])
gammaLocal[s][i]+=field[s]
......@@ -119,10 +118,10 @@ def getEquivalentProperty(d, localProperties, refProp, greenMatrix):
eps = localProperties.dot(tau)
for cell in range(2**d):
A[cell,:,i]=eps[d*cell:d*(cell+1)]
avgT[i]+=tau[d*cell:d*(cell+1)]
avgMT[i]+=eps[d*cell:d*(cell+1)]
avgT[i]/=2**d
avgMT[i]/=2**d
avgT[:,i]+=tau[d*cell:d*(cell+1)]
avgMT[:,i]+=eps[d*cell:d*(cell+1)]
avgT[:,i]/=2**d
avgMT[:,i]/=2**d
for cell in range(2**d):
A[cell]=np.dot(A[cell],np.linalg.inv(avgMT))
return avgMT.dot(np.linalg.inv(avgT)), A #refProp+avgT.dot(np.linalg.inv(avgMT))
......@@ -202,8 +201,8 @@ def hierachicHomogenize(N,phasemap,phaseProps):
end = time.time()
recursionTime = end - start
prevk0 = k0
k1 = (kest+kest.transpose())/2#np.diag(np.diag(kest))
k0=(k1[0,0]+k1[1,1])/2*np.eye(d)
k_ = (kest+kest.transpose())/2#np.diag(np.diag(kest))
k0=(k_[0,0]+k_[1,1])/2*np.eye(d)
ind+=1
delta = np.linalg.norm(k0-prevk0)/np.linalg.norm(k0)
start = time.time()
......@@ -214,7 +213,7 @@ def hierachicHomogenize(N,phasemap,phaseProps):
#print(kest,k0)
#print("%e\t%s\t%s" % (delta,resetTime,recursionTime))
return kest,t_A
return prevk0,kest,t_A
#print(kest)
khom.append(kest)
......@@ -227,24 +226,25 @@ def hierachicHomogenize(N,phasemap,phaseProps):
N0 = 16
N0 = 128
d=2
N=N0*np.ones(d,dtype=np.int64)
radius = 4
radius = 8
shape = {2:"disc",3:"sphere"}
V = pi*radius**2 if d==2 else 4*pi/3*radius**3
k1 = 1.
k2 = 0.001
Khom = []
Frac = []
for model in ["A","B"]:
for model in ["B"]:
if model == "A":
phaseProps = np.array([k1,k2],dtype=np.float64)
elif model == "B":
phaseProps = np.array([k2,k1],dtype=np.float64)
khom=[]
frac=[]
for target_phi in range(75,99,4):
kCA_list=[]
kSA_list=[]
kref_list=[]
for target_phi in [98,96,94,92,90,88,86,84,78,77,76,66,65]:
ncenter = int(-np.log(1-target_phi/100.)*N0**d/V)
centers = np.random.rand(ncenter,d)*N0
radii = np.ones(ncenter)*radius
......@@ -257,24 +257,52 @@ for model in ["A","B"]:
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()
# plt.figure()
# plt.imshow(microstructure.phasemap,origin='lower')
# plt.show()
khh,t_A_hh = hierachicHomogenize(N,phasemap,phaseProps)
print(khh)
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])
# 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)
k0=np.mean(phaseProps)*np.eye(len(N))
kref=Ref(tau_hh,N,k0,ki)
print(khh,kCA,kSA,kref)
khom.append(khh)
kCA_list.append(kCA)
kSA_list.append(kSA)
kref_list.append(kref)
print("voigt %f" % ((1-phi)*k1+phi*k2))
print("reuss %f" % (1./((1-phi)/k1+phi/k2)))
plt.figure()
plt.imshow(t_A_hh[:,:,0,0],origin='lower')
plt.show()
plt.close()
Khom.append(np.array(khom))
Frac.append(np.array(frac))
savefile=open("boolean_%s_contrast%dmodel%s.txt" % (shape[d],int(k1/k2),model), 'w')
np.savetxt(savefile,np.array([[f,np.trace(kh)/d] for (f,kh) in zip(frac,khom)]))
savefile.close()
# plt.figure()
# 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')
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')
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')
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()
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))
......@@ -284,14 +312,45 @@ 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="hierachical num. self-consistent A")
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")
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")
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.pdf" % (shape[d],int(k1/k2)))
plt.savefig("boolean_%s_contrast%d_model%s.pdf" % (shape[d],int(k1/k2),"A"))
plt.show()
plt.close()
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")
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"))
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.show()
plt.close()
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