Commit 8a32c1e8 authored by MARTIN Antoine's avatar MARTIN Antoine

made k0 isotropic at each step+extraction of HH local strains

parent ac3dd0ff
import numpy as np
from numba import jit, float64, int64, prange
from numba import jit, float64, int64, prange, generated_jit
import matplotlib.pyplot as plt
import time
from scipy.optimize import fsolve
......@@ -9,6 +9,19 @@ from greenOperatorConductionNumba import *
#from constitutiveLaw import *
from microstructure import *
@jit(nopython=True)
def rec_mult(t,position,pshifts,A0,d,stride):
if stride==1:
t[int(position[0]),int(position[1])]=np.dot(t[int(position[0]),int(position[1])],A0)#il faudrait trouver une alternative pour les autres dimensions...
else:
for i in prange(2**d):
nextPosition = position+stride/2*pshifts[i]
rec_mult(t,nextPosition,pshifts,A0,d,stride/2)
@jit
def nullifyField(field,N):
#z=np.zeros(d)
......@@ -60,7 +73,7 @@ def getGreenMatrix(M,k0,gammaLocal,shifts):
i+=1
return matrix
def resetReference(N,k0,greenMatrices,phaseProps,phaseInvDeltaProp):
def resetReference(N,k0,greenMatrices,phaseProps,phaseInvDeltaProp,ind,prevk0):
M=2*np.ones(N.shape,dtype=np.int64)
if k0.__class__ == np.ndarray:
refProp = k0
......@@ -69,13 +82,21 @@ def resetReference(N,k0,greenMatrices,phaseProps,phaseInvDeltaProp):
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
if ind==0:
depth = 0
while M[0]<=N[0]:
gammaLocal = getLocalGreenNumba(M,k0)
greenMatrices[depth] = getGreenMatrix(M,k0,gammaLocal,shifts)
M*=2
depth+=1
else:
greenM = greenMatrices*prevk0[0,0]/k0[0,0]
depth = 0
while M[0]<=N[0]:
greenMatrices[depth]=greenM[depth]
M*=2
depth+=1
return(refProp)
@jit(nopython=True)#(float64[:,:](int64,float64[:,:],float64[:,:],float64[:,:]))
......@@ -87,7 +108,9 @@ def getEquivalentProperty(d, localProperties, refProp, greenMatrix):
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)
A=np.zeros((2**d,d,d),dtype=np.float64)
#MT=np.zeros((2**d,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):
......@@ -95,15 +118,18 @@ def getEquivalentProperty(d, localProperties, refProp, greenMatrix):
tau = np.linalg.solve(matrix,rhs)
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
return avgMT.dot(np.linalg.inv(avgT))#refProp+avgT.dot(np.linalg.inv(avgMT))
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))
@jit(nopython=True)#(float64[:,:](int64,int64[:],int64[:],int64,float64[:,:],float64[:],float64[:,:],float64[:,:,:]))
def recursion(depth, M, pshifts, position, stride, phasemap, mapStrides, phaseProps, refProp, greenMatrices):
def recursion(depth, M,N, pshifts, position, stride, phasemap, mapStrides, phaseProps, refProp, greenMatrices,t_A):
d=len(M)
localProperties = np.zeros(greenMatrices[0].shape)
#i=0
......@@ -112,27 +138,33 @@ def recursion(depth, M, pshifts, position, stride, phasemap, mapStrides, phasePr
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)
localProperties[d*i:d*(i+1),d*i:d*(i+1)] = recursion(depth+1,M*2,N,pshifts,nextPosition,stride/2, phasemap, mapStrides, phaseProps, refProp, greenMatrices,t_A)
else:
n = position+stride*pshifts[i]#np.array(shift) #tuple() constructor not supported by numba
nflat = np.sum(n*mapStrides)
nflat = int(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])
equivalentProperty,A = getEquivalentProperty(d, localProperties, refProp, greenMatrices[depth])
for i in prange(2**d):
nextPosition = position+stride*pshifts[i]
rec_mult(t_A,nextPosition,pshifts,A[i],d,stride)
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):
def recursionParallel(depth, M, N, pshifts, position, stride, phasemap, mapStrides, phaseProps, refProp, greenMatrices,t_A):
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
localProperties[d*i:d*(i+1),d*i:d*(i+1)] = recursion(depth+1,M*2,N,pshifts,nextPosition,stride/2, phasemap, mapStrides, phaseProps, refProp, greenMatrices,t_A)#do not use tuple as arguments in parallel mode
equivalentProperty,A = getEquivalentProperty(d, localProperties, refProp, greenMatrices[depth])
for i in prange(2**d):
nextPosition = position+stride*pshifts[i]
rec_mult(t_A,nextPosition,pshifts,A[i],d,stride)
return equivalentProperty,t_A
def hierachicHomogenize(N,phasemap,phaseProps):
d=len(N)
......@@ -140,37 +172,50 @@ def hierachicHomogenize(N,phasemap,phaseProps):
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)
pshifts = np.zeros((2**d,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)
ind=0
refProp = resetReference(N,k0,greenMatrices,phaseProps,phaseInvDeltaProp,ind,k0)
print("delta\treset(s)\trecurrence(s)")
delta = 1
while delta > 1e-3:
M=2*np.ones(N.shape,dtype=np.int64)
t_A=np.zeros(tuple(d*[N[0]] + 2*[d]))
def rec_fill(t):
if t[0].shape==(d,d):
for ii in range(N[0]):
t[ii]=np.eye(d)
else:
for ii in range(N[0]):
rec_fill(t[ii])
rec_fill(t_A)
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)
start = time.time()
invdk,t_A = recursionParallel(depth, M, N, pshifts, position, stride, phasemap, mapStrides, phaseInvDeltaProp, refProp, greenMatrices,t_A)
kest = refProp+np.linalg.inv(invdk)
end = time.time()
recursionTime = end - start
prevk0 = k0
k0 = (kest+kest.transpose())/2#np.diag(np.diag(kest))
k1 = (kest+kest.transpose())/2#np.diag(np.diag(kest))
k0=(k1[0,0]+k1[1,1])/2*np.eye(d)
ind+=1
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)
refProp = resetReference(N,k0,greenMatrices,phaseProps,phaseInvDeltaProp,ind,prevk0)
end = time.time()
resetTime = end-start
#print(k0)
print("%e\t%s\t%s" % (delta,resetTime,recursionTime))
return kest
print(kest)
#print(kest,k0)
#print("%e\t%s\t%s" % (delta,resetTime,recursionTime))
return kest,t_A
#print(kest)
khom.append(kest)
......@@ -182,10 +227,10 @@ def hierachicHomogenize(N,phasemap,phaseProps):
N0 = 128
d=3
N0 = 16
d=2
N=N0*np.ones(d,dtype=np.int64)
radius = 8
radius = 4
shape = {2:"disc",3:"sphere"}
V = pi*radius**2 if d==2 else 4*pi/3*radius**3
k1 = 1.
......@@ -211,19 +256,23 @@ for model in ["A","B"]:
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)
plt.imshow(microstructure.phasemap,origin='lower')
plt.show()
khh,t_A_hh = 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)))
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), 'ab')
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()
......
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