Commit ac3dd0ff authored by Jeremy BLEYER's avatar Jeremy BLEYER

Add diffusion files from Francois

parent c301e824
import time
import numpy as np
from math import pi
from numba import jit, prange, njit
import pyfftw
@jit
def discrete_frequency(k,n,h):
if 2*k>n:
k-=n
return 2*pi*k/(n*h)#2*np.tan(pi*k/n)/h
@jit
def initialize_frequencies(N,K,H,filter_level):
d=len(N)
frequencies = np.zeros((d,max(K),filter_level))
for i in range(d):
for ki in range(K[i]):
for p in range(filter_level):
kip = ki+N[i]*p
frequencies[i][ki][p] = discrete_frequency(kip,N[i]*filter_level,H[i])
return frequencies
@jit
def initialize_filters(N,K,n):
d=len(N)
filters = np.zeros((d,max(K),n))
for i in range(d):
for ki in range(K[i]):
for p in range(n):
z=ki/float(N[i])+p
if z==0:
filters[i][ki][p] = 1
else:
filters[i][ki][p] = (np.sin(pi*z)/(n*np.sin(pi*z/n)))**2
return filters
def initializeGreen(N,filter_level=2):
d = len(N)
numComp=d
H=1./N
K=N.copy()
K[-1] = N[-1]//2+1
Npadded = N.copy()
Npadded[-1] = K[-1]*2
full_array = pyfftw.n_byte_align_empty(np.append(Npadded,numComp), pyfftw.simd_alignment,'float64')
field = full_array[...,:N[-1],:]
field_fourier = full_array.ravel().view('complex128').reshape(np.append(K,numComp))
fft = pyfftw.FFTW(field, field_fourier, axes=range(d))
ifft = pyfftw.FFTW(field_fourier, field, axes=range(d),direction='FFTW_BACKWARD')
frequencies = initialize_frequencies(N,K,H,filter_level)
filters = initialize_filters(N,K,filter_level)
tupleK=tuple(K[1:])
tupleFilters = np.zeros((filter_level**d,d),dtype=np.int64)#tuple(d*[filter_level])
i=0
for shift in np.ndindex(tuple(d*[filter_level])):
tupleFilters[i] = np.array(shift)
i+=1
return field,field_fourier,fft,ifft,tupleK,frequencies,filters,tupleFilters
@jit
def get_matrix_fourier_filtered(k,N,frequencies,filter_level,filters,tupleFilters,k0):
"""
k : multi-index in frequency grid
"""
d=len(k)
Gamma = np.zeros((d,d),dtype=np.complex128)
q = np.empty(d)
#for p in np.ndindex(tupleFilters):
# kp = k+N*np.array(p)
for j in range(len(tupleFilters)):
p = tupleFilters[j]
kp = k+N*p
f = 1
for i in range(d):
q[i] = frequencies[i][k[i]][p[i]]
f*=filters[i][k[i]][p[i]]
g=np.outer(q,q)
qk0q = q.dot((k0).dot(q))
if not(qk0q==0):# or nullify_frequency(kp/filter_level,N)):
Gamma+=f*g/(qk0q)
return Gamma
@jit
def operate_fourier(tau,k,N,frequencies,filter_level,filters,tupleFilters,k0):
"""
tau : d-dimensionnal array of polarization
k : multi-index in frequency grid
"""
return -get_matrix_fourier_filtered(k,N,frequencies,filter_level,filters,tupleFilters,k0).dot(tau)
@jit(parallel=False)
def operate_fourier_field(x,y,tupleK,N,frequencies,filter_level,filters,tupleFilters,k0):
"""
x : input, field of comp-dimensionnal array
y : output, field of comp-dimensionnal array (operates in place if x = y)
tupleK : a tuple of fourier grid dimensions
"""
for kx in prange(N[0]):
for kyz in np.ndindex(tupleK):
k = (kx,)+kyz
y[k] = operate_fourier(x[k],np.array(k),N,frequencies,filter_level,filters,tupleFilters,k0)
#do not use tuple as arguments in parallel mode -> remove tupleFIlters
#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):
#start = time.time()
fft()
#end = time.time()
#t1 = end-start
#start = time.time()
operate_fourier_field(xFourier,yFourier,tupleK,N,frequencies,filter_level,filters,tupleFilters,k0)
#end = time.time()
#t2 = end-start
#start = time.time()
ifft()
#end = time.time()
#t3 = end-start
#print("%s\t%s\t%s" % (t1,t2,t3))
"""
TODO: optimize operate_fourier_field,
300 times slower than FFT
"""
import numpy as np
from numba import jit, float64, int64, prange
import matplotlib.pyplot as plt
import time
from scipy.optimize import fsolve
from greenOperatorConductionNumba import *
#from discretizationGrid import *
#from constitutiveLaw import *
from microstructure import *
@jit
def nullifyField(field,N):
#z=np.zeros(d)
for n in np.ndindex(N):
field[n] = np.zeros(len(N))
def getLocalGreenNumba(M,k0):
d=len(M)
filter_level=2
field,field_fourier,fft,ifft,tupleK,frequencies,filters,tupleFilters = initializeGreen(M,filter_level=filter_level)
gammaLocal = np.zeros(tuple(d*[3] + 2*[d]))
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)
for m in np.ndindex(tuple(d*[3])):
s = tuple([m[j]-1 for j in range(d)])
gammaLocal[s][i]+=field[s]
return gammaLocal
"""
def getLocalGreen(M,k0):
d = len(M)
grid = discretizationGrid(M)
gamma = greenOperatorConduction(grid,k0,discretization="truncated",filter_level=2)
gammaLocal = np.zeros(tuple(d*[3] + 2*[gamma.numComp]))
for i in range(gamma.numComp):
for n in np.ndindex(tuple(M)):
gamma.field[n] = np.zeros(d)
gamma.field[tuple(d*[0]+[i])] = 1
gamma.operate_field()
for m in np.ndindex(tuple(d*[3])):
s = tuple([m[j]-1 for j in range(d)])
gammaLocal[s][i]+=gamma.field[s]
return gammaLocal
"""
def getGreenMatrix(M,k0,gammaLocal,shifts):
d=len(M)
numDOF = d*2**d
matrix = np.zeros((numDOF,numDOF))
i = 0
for I in np.ndindex(shifts):
j = 0
for J in np.ndindex(shifts):
shift=tuple(np.array(I)-np.array(J))
matrix[d*i:d*(i+1),d*j:d*(j+1)] -= gammaLocal[shift]
j+=1
i+=1
return matrix
def resetReference(N,k0,greenMatrices,phaseProps,phaseInvDeltaProp):
M=2*np.ones(N.shape,dtype=np.int64)
if k0.__class__ == np.ndarray:
refProp = k0
else:
refProp = k0*np.eye(len(N))
for i in range(len(phaseProps)):
kloc=phaseProps[i]
phaseInvDeltaProp[i] = np.linalg.inv(np.eye(d)*kloc-refProp)
depth = 0
shifts = tuple(d*[2])
while M[0]<=N[0]:
gammaLocal = getLocalGreenNumba(M,k0)
greenMatrices[depth] = getGreenMatrix(M,k0,gammaLocal,shifts)
M*=2
depth+=1
return(refProp)
@jit(nopython=True)#(float64[:,:](int64,float64[:,:],float64[:,:],float64[:,:]))
def getEquivalentProperty(d, localProperties, refProp, greenMatrix):
#Tk = array of tensors obtained by inversion of HS matrix
#Ck-1 = Cref + avg(Tk):inv(avg(Mk:Tk))
matrix = localProperties+greenMatrix
#print(matrix)
numDOF = len(matrix)#d*2**d
avgT = np.zeros((d,d),dtype=np.float64)
avgMT = np.zeros((d,d),dtype=np.float64)
rhs = np.zeros(numDOF,dtype=np.float64)
for i in range(d):
macroGradient = np.array([1. if i==j else 0. for j in range(d)],dtype=np.float64)
for cell in range(2**d):
rhs[d*cell:d*(cell+1)] = macroGradient
tau = np.linalg.solve(matrix,rhs)
eps = localProperties.dot(tau)
for cell in range(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
return avgMT.dot(np.linalg.inv(avgT))#refProp+avgT.dot(np.linalg.inv(avgMT))
@jit(nopython=True)#(float64[:,:](int64,int64[:],int64[:],int64,float64[:,:],float64[:],float64[:,:],float64[:,:,:]))
def recursion(depth, M, pshifts, position, stride, phasemap, mapStrides, phaseProps, refProp, greenMatrices):
d=len(M)
localProperties = np.zeros(greenMatrices[0].shape)
#i=0
#for shift in np.ndindex(shifts):#tuple(d*[2])):
for i in prange(2**d):
nextPosition = position+stride*pshifts[i]#np.array(shift)
if stride > 1:
#localProperties[d*i:d*(i+1),d*i:d*(i+1)] = np.linalg.inv(recursion(depth+1,M*2,shifts,nextPosition,stride/2, phasemap, mapStrides, phaseProps, refProp, greenMatrices)-refProp)
localProperties[d*i:d*(i+1),d*i:d*(i+1)] = recursion(depth+1,M*2,pshifts,nextPosition,stride/2, phasemap, mapStrides, phaseProps, refProp, greenMatrices)
else:
n = position+stride*pshifts[i]#np.array(shift) #tuple() constructor not supported by numba
nflat = np.sum(n*mapStrides)
#kloc=phaseProps[phasemap[nflat]]
#localProperties[d*i:d*(i+1),d*i:d*(i+1)] = np.linalg.inv(np.eye(d)*kloc-refProp)
localProperties[d*i:d*(i+1),d*i:d*(i+1)] = phaseProps[phasemap[nflat]]
# i+=1
equivalentProperty = getEquivalentProperty(d, localProperties, refProp, greenMatrices[depth])
return equivalentProperty
@jit(nopython=True,parallel=False)#very minor speed up, communication problem ?
def recursionParallel(depth, M, pshifts, position, stride, phasemap, mapStrides, phaseProps, refProp, greenMatrices):
d=len(M)
numDOF = d*2**d
localProperties = np.zeros((numDOF,numDOF))#np.zeros(greenMatrices[0].shape)
for i in prange(2**d):
nextPosition = position+stride*pshifts[i]
localProperties[d*i:d*(i+1),d*i:d*(i+1)] = recursion(depth+1,M*2,pshifts,nextPosition,stride/2, phasemap, mapStrides, phaseProps, refProp, greenMatrices)#do not use tuple as arguments in parallel mode
equivalentProperty = getEquivalentProperty(d, localProperties, refProp, greenMatrices[depth])
return equivalentProperty
def hierachicHomogenize(N,phasemap,phaseProps):
d=len(N)
phaseInvDeltaProp = np.zeros((len(phaseProps),d,d))
numDOF = d*2**d
greenMatrices=np.zeros((int(np.log2(N[0])),numDOF,numDOF),dtype=np.float64)
shifts = tuple(d*[2])
pshifts = np.zeros((d**2,d),dtype=np.int64)
i=0
for shift in np.ndindex(shifts):
pshifts[i] += np.array(shift)
i+=1
k0=np.mean(phaseProps)*np.eye(len(N))
refProp = resetReference(N,k0,greenMatrices,phaseProps,phaseInvDeltaProp)
print("delta\treset(s)\trecurrence(s)")
delta = 1
while delta > 1e-3:
M=2*np.ones(N.shape,dtype=np.int64)
depth = 0
position=np.zeros(N.shape,dtype=np.int64)
stride = N0/2
start = time.time()
invdk = recursionParallel(depth, M, pshifts, position, stride, phasemap, mapStrides, phaseInvDeltaProp, refProp, greenMatrices)
kest = refProp+np.linalg.inv(invdk)
end = time.time()
recursionTime = end - start
prevk0 = k0
k0 = (kest+kest.transpose())/2#np.diag(np.diag(kest))
delta = np.linalg.norm(k0-prevk0)/np.linalg.norm(k0)
start = time.time()
#refProp = resetReference(N,k0,greenMatrices)
refProp = resetReference(N,k0,greenMatrices,phaseProps,phaseInvDeltaProp)
end = time.time()
resetTime = end-start
#print(k0)
print("%e\t%s\t%s" % (delta,resetTime,recursionTime))
return kest
print(kest)
khom.append(kest)
phi=np.sum(phasemap)*1./N0**d
frac.append(phi)
print("voigt %f" % ((1-phi)*k1+phi*k2))
print("reuss %f" % (1./((1-phi)/k1+phi/k2)))
N0 = 128
d=3
N=N0*np.ones(d,dtype=np.int64)
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"]:
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):
ncenter = int(-np.log(1-target_phi/100.)*N0**d/V)
centers = np.random.rand(ncenter,d)*N0
radii = np.ones(ncenter)*radius
microstructure = booleanSpheres(centers,radii,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)
phi=np.sum(phasemap)*1./N0**d
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, cmap="Greys")
plt.show(block=False)
"""
khh = hierachicHomogenize(N,phasemap,phaseProps)
print(khh)
khom.append(khh)
print("voigt %f" % ((1-phi)*k1+phi*k2))
print("reuss %f" % (1./((1-phi)/k1+phi/k2)))
Khom.append(np.array(khom))
Frac.append(np.array(frac))
savefile=open("boolean_%s_contrast%dmodel%s.txt" % (shape[d],int(k1/k2),model), 'ab')
np.savetxt(savefile,np.array([[f,np.trace(kh)/d] for (f,kh) in zip(frac,khom)]))
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")
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.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.show()
from numba import jitclass,int16,int8,jit
import numpy as np
#@jitclass
#class microstructure():
# def __init__(self):
# pass
#
# def get_phase(self):
# pass
@jitclass([('a',int16)])
class square(object):
def __init__(self,a):
self.a=a
def get_phase(self,n):
for ni in n:
if ni>=self.a:
return 0
return 1
class booleanSpheres(object):
def __init__(self,centers,radii,N):
self.centers=centers
self.radii=radii
self.N = N
self.d = len(N)
self.phasemap = np.zeros(tuple(N),dtype=np.int8)
self.initializePhaseMap()
def initializePhaseMap(self):
for i,center in enumerate(self.centers):
r = self.radii[i]
c = np.rint(center).astype(int)
for s in np.ndindex(tuple(self.d*[2*int(r)+1])):
ds = np.array(s)-int(r)
n = c-ds
if (n-center+0.5).dot(n-center+0.5)<= r*r:
self.phasemap[tuple(n % self.N)]=1
def get_phase(self,n):
return self.phasemap[n]
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