now the program flow is working

This commit is contained in:
Martin Diehl 2014-07-08 16:09:02 +00:00
parent a1162e6fb4
commit fdb4f2ffea
1 changed files with 47 additions and 30 deletions

View File

@ -6,6 +6,8 @@ from scipy.optimize import curve_fit
from scipy.linalg import svd from scipy.linalg import svd
import threading,time,os,subprocess,shlex import threading,time,os,subprocess,shlex
import damask import damask
geomName='20grains16x16x16'
popt1=np.ones(6,'d')
def execute(cmd,dir='./'): def execute(cmd,dir='./'):
@ -30,11 +32,11 @@ def Hill48(x, F,G,H,L,M,N):
def vonMises(x, S_y): def vonMises(x, S_y):
sv=np.zeros(0,'d') sv=np.zeros(0,'d')
for i in xrange(np.shape(x)[1]): for i in xrange(np.shape(x)[0]):
U, l, Vh = svd(np.array(x[:,i]).reshape(3,3)) U, l, Vh = svd(np.array(x[i,:]).reshape(3,3))
sv = np.append(sv,l) sv = np.append(sv,l)
sv = sv.reshape(3,np.shape(x)[1]) sv = sv.reshape(np.shape(x)[0],3)
ooo= (sv[2,:]-sv[1,:])**2+(sv[1,:]-sv[0,:])**2+(sv[0,:]-sv[2,:])**2-2*S_y**2 ooo= (sv[:,2]-sv[:,1])**2+(sv[:,1]-sv[:,0])**2+(sv[:,0]-sv[:,2])**2-2*S_y**2
return ooo return ooo
#--------------------------------------------------------------------------------------------------- #---------------------------------------------------------------------------------------------------
@ -85,13 +87,15 @@ class Criterion(object):
def fit(self,stress): def fit(self,stress):
try: try:
popt, pcov = curve_fit(Hill48, stress, np.zeros(np.shape(stress)[1])) spopt1, pcov = curve_fit(Hill48, stress, np.zeros(np.shape(stress)[1]),p0=popt1)
popt = popt1
print 'Hill 48', popt print 'Hill 48', popt
print Hill48(stress,popt[0],popt[1],popt[2],popt[3],popt[4],popt[5])
except Exception as detail: except Exception as detail:
print detail print detail
pass pass
popt, pcov = curve_fit(vonMises, stress, np.zeros(np.shape(stress)[1])) #popt, pcov = curve_fit(vonMises, stress, np.zeros(np.shape(stress)[0]))
print 'von Mises', popt #print 'von Mises', popt
#--------------------------------------------------------------------------------------------------- #---------------------------------------------------------------------------------------------------
@ -118,39 +122,53 @@ def doSim(delay,thread):
s.acquire() s.acquire()
me=getLoadcase() me=getLoadcase()
print('generating loadcase for sim %i from %s'%(me,thread)) if not os.path.isfile('%s.load'%me):
f=open('%s.load'%me,'w') print('generating loadcase for sim %s from %s'%(me,thread))
f.write(myLoad.getNext(me)) f=open('%s.load'%me,'w')
f.close() f.write(myLoad.getNext(me))
f.close()
print('starting simulation %i from %s'%(me,thread)) s.release()
s.release() else: s.release()
execute('DAMASK_spectral -l %i -g 20grains16x16x16'%me) s.acquire()
if not os.path.isfile('%s_%s.spectralOut'%(geomName,me)):
print('starting simulation %i from %s'%(me,thread))
s.release()
execute('DAMASK_spectral -l %i -g %i'%(me,geomName))
else: s.release()
s.acquire() s.acquire()
print('startin post processing for sim %i from %s'%(me,thread)) if not os.path.isfile('./postProc/%s_%s.txt'%(geomName,me)):
print('starting post processing for sim %i from %s'%(me,thread))
s.release()
execute('postResults --cr f,p %i_%i.spectralOut'%(geomName,me))
execute('addCauchy ./postProc/%i_%i.txt'%(geomName,me))
execute('addStrainTensors -l -v ./postProc/%i_%i.txt'%(geomName,me))
execute('addMises -s Cauchy -e ln(V) ./postProc/%i_%i.txt'%(geomName,me))
else: s.release()
s.acquire()
print('reading values for sim %i from %s'%(me,thread))
s.release() s.release()
execute('postResults --cr f,p 20grains16x16x16_%i.spectralOut'%me)
execute('addCauchy ./postProc/20grains16x16x16_%i.txt'%me)
execute('addStrainTensors -l -v ./postProc/20grains16x16x16_%i.txt'%me)
execute('addMises -s Cauchy -e ln(V) ./postProc/20grains16x16x16_%i.txt'%me)
refFile = open('./postProc/20grains16x16x16_%i.txt'%me) refFile = open('./postProc/20grains16x16x16_%i.txt'%me)
refFile = open('./postProc/20grains16x16x16_1.txt')
table = damask.ASCIItable(refFile) table = damask.ASCIItable(refFile)
table.head_read() table.head_read()
for l in ['Mises(ln(V))','1_Cauchy']: for l in ['Mises(ln(V))','1_Cauchy']:
if l not in table.labels: print '%s not found'%l if l not in table.labels: print '%s not found'%l
while table.data_read(): while table.data_read():
if float(table.data[table.labels.index('Mises(ln(V))')]) > 0.002: if float(table.data[table.labels.index('Mises(ln(V))')]) > 0.002:
yieldStress = np.array(table.data[table.labels.index('1_Cauchy'):table.labels.index('9_Cauchy')+1],'d').reshape(3,3) yieldStress = np.array(table.data[table.labels.index('1_Cauchy'):table.labels.index('9_Cauchy')+1],'d')
s.acquire() s.acquire()
print('startin fitting for sim %i from %s'%(me,thread)) print('starting fitting for sim %i from %s'%(me,thread))
global stressAll global stressAll
stressAll=np.append(stressAll,yieldStress.reshape(9)) stressAll=np.append(yieldStress,stressAll)
stressAll=stressAll.reshape(9,len(stressAll)//9) for i in range(np.shape(stressAll.reshape(len(stressAll)//9,9).transpose())[1]):
myFit.fit(stressAll) print i, stressAll.reshape(len(stressAll)//9,9).transpose()[:,i]
myFit.fit(stressAll.reshape(len(stressAll)//9,9).transpose())
s.release() s.release()
def getLoadcase(): def getLoadcase():
@ -169,13 +187,13 @@ def converged():
# main # main
minN_simulations=20 minN_simulations=20
maxN_simulations=20 maxN_simulations=14
N_simulations=0 N_simulations=0
s=threading.Semaphore(1) s=threading.Semaphore(1)
scale = 0.02 scale = 0.02
incs = 10 incs = 10
duration = 10 duration = 10
stressAll=np.zeros(0,'d') stressAll=np.zeros(0,'d').reshape(0,0)
myLoad = Loadcase() myLoad = Loadcase()
myFit = Criterion('Hill48') myFit = Criterion('Hill48')
@ -189,5 +207,4 @@ for i in range(N_threads):
for i in range(N_threads): for i in range(N_threads):
t[i].join() t[i].join()
print "Exiting Main Thread" print "Exiting Main Thread"