yade-users team mailing list archive

[Mumić <question688652@xxxxxxxxxxxxxxxxxxxxx>]

Question #688652 on Yade changed:
https://answers.launchpad.net/yade/+question/688652

    Status: Solved => Open

Mumić is still having a problem:
Hello everyone,

I have managed to plot vary the values as described on my last post but
now I face another problem: "yade" and "yade-batch" yield different
results.

When I compute Poisson modulus with batch, the values obtained are akin
to those found on Šmilauer's thesis, but if I run the simulations
manually the values will differ.

My script can be found below:

###############################################
#!/usr/bin/python
# -*- coding: utf-8 -*-
from yade import plot,pack,timing
import time, sys, os, copy, numpy as np

savedir = './elastic_validation/'

# default parameters or from table
readParamsFromTable(
	young = 66e9,
	poisson = 0.1,

	sigmaT=1e99,
	#frictionAngle=atan(0.8),
	epsCrackOnset=1e99,
	relDuctility=1e99,
	intRadius=1.5,
	dtSafety=.1,
	damping=0.8,
	strainRateTension=.005,
	strainRateCompression=.5,
	setSpeeds=True,
	# 1=tension, 2=compression (ANDed; 3=both)
	doModes=1,

	specimenLength=.039850,
	specimenRadius=.019875,
	sphereRadius=0.75e-3,

	# isotropic confinement (should be negative)
	isoPrestress=0,
	noTableOk=True,
	
	n=100
)

from yade.params.table import *
if 'description' in list(O.tags.keys()): O.tags['id']=O.tags['id']+O.tags['description']

print(young, poisson)
# make geom; the dimensions are hard-coded here; could be in param table if desired
# z-oriented hyperboloid, length 20cm, diameter 10cm, skirt 8cm
# using spheres 7mm of diameter
concreteId=O.materials.append(
		CpmMat(
			young=young,
			#frictionAngle=frictionAngle,
			poisson=poisson,
			density=6420,
			sigmaT=sigmaT,
			relDuctility=relDuctility,
			epsCrackOnset=epsCrackOnset,
			#isoPrestress=isoPrestress,
			neverDamage=True,
			)
		)

#sps=SpherePack()

#sp=pack.randomDensePack(pack.inAlignedBox((0,0,0),(specimenLength,specimenLength,specimenLength)),
sp=pack.randomDensePack(pack.inCylinder((0,0,-.5*specimenLength),(0,0,.5*specimenLength),specimenRadius),
						spheresInCell=5000,
						radius=sphereRadius,
						memoizeDb='/tmp/calibCylinder3.sqlite',
						returnSpherePack=True)

sp.toSimulation(material=concreteId)

def getElementsCurrentPosition():
	return np.array(list(map(lambda x: [x.state.pos[0], x.state.pos[1], x.state.pos[2], x.id], O.bodies)))

ibpos = getElementsCurrentPosition()

def fetchElements(pos=(0,0,0), n=1):
	pos = np.array(pos)
	distances = np.sum(ibpos[:,:3]-pos, axis=1)**2
	return ibpos[distances.argsort()][:n,3].tolist()

def computePosition(n=1):
	elements = np.array(fetchElements(n=n)).astype(int)
	current_position = getElementsCurrentPosition()
	return current_position[elements]
	
allpos = computePosition(n)

bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.dt=dtSafety*PWaveTimeStep()
print('Timestep',O.dt)

mm,mx=[pt[axis] for pt in aabbExtrema()]
imm, imx = mm, mx
coord_25,coord_50,coord_75=mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm)
area_25,area_50,area_75=approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),],verletDist=.05*sphereRadius),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc')],
		[Ip2_CpmMat_CpmMat_CpmPhys()],
		[Law2_ScGeom_CpmPhys_Cpm()],
	),
	NewtonIntegrator(damping=damping,label='damper'),
	CpmStateUpdater(realPeriod=.5),
	UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
	PyRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
	PyRunner(realPeriod=1,command='stopIfDamaged()',label='damageChecker'),
]

# plot stresses in ¼, ½ and ¾ if desired as well; too crowded in the graph that includes confinement, though
plot.plots={'eps':('sigma')} #,'sigma.50')},'t':('eps')} #'sigma.25','sigma.50','sigma.75')}

O.saveTmp('initial');

O.timingEnabled=False

global mode
mode='tension' if doModes & 1 else 'compression'

def initTest():
	global mode
	print("init")
	if O.iter>0:
		O.wait();
		O.loadTmp('initial')
		print("Reversing plot data"); plot.reverseData()
	else: plot.plot(subPlots=False)
	strainer.strainRate=abs(strainRateTension) if mode=='tension' else -abs(strainRateCompression)
	try:
		from yade import qt
		renderer=qt.Renderer()
		renderer.dispScale=(1000,1000,1000) if mode=='tension' else (100,100,100)
	except ImportError: pass
	print("init done, will now run.")
	O.step(); # to create initial contacts
	# now reset the interaction radius and go ahead
	ss2sc.interactionDetectionFactor=1.
	is2aabb.aabbEnlargeFactor=1.

        O.run()

def stopIfDamaged():
	global mode, allpos, alldisp
	if O.iter<2 or 'sigma' not in plot.data: return # do nothing at the very beginning
	sigma,eps=plot.data['sigma'],plot.data['eps']
	allpos = np.append(allpos, computePosition(n), 1)
	print(allpos.shape)
	extremum=max(sigma) if (strainer.strainRate>0) else min(sigma)
	minMaxRatio=0.8 if mode=='tension' else 0.8
	if extremum==0: return
	import sys;	sys.stdout.flush()
	if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer.strain)>(5e-3 if isoPrestress==0 else 5e-2) or O.iter>5000:
		if mode=='tension' and doModes & 2: # only if compression is enabled
			numpy.savetxt(savedir+str(young)+' '+str(poisson)+' 1.txt', numpy.array([sigma, eps]).T, header='sigma\teps', delimiter='\t')
			numpy.savetxt(savedir+str(young)+' '+str(poisson)+' 3.txt', allpos, delimiter='\t')
			with open('./uniax.py', 'r') as f1:
				text = f1.read()
				with open(savedir+str(young)+' '+str(poisson)+'.py', 'w') as f2:
					f2.write(text)
			mode='compression'
			O.save('/tmp/uniax-tension.yade.gz')
			print("Saved /tmp/uniax-tension.yade.gz (for use with interaction-histogram.py and uniax-post.py)")
			print("Damaged, switching to compression... "); O.pause()
			# important! initTest must be launched in a separate thread;
			# otherwise O.load would wait for the iteration to finish,
			# but it would wait for initTest to return and deadlock would result
			import _thread; _thread.start_new_thread(initTest,())
			return
		else:
			numpy.savetxt(savedir+str(young)+' '+str(poisson)+' 2.txt', numpy.array([sigma, eps]).T, header='sigma\teps', delimiter='\t')
			numpy.savetxt(savedir+str(young)+' '+str(poisson)+' 4.txt', allpos, delimiter='\t')
			with open('./uniax.py', 'r') as f1:
				text = f1.read()
				with open(savedir+str(young)+' '+str(poisson)+'.py', 'w') as f2:
					f2.write(text)
			print("Damaged, stopping.")
			ft,fc=max(sigma),min(sigma)
			if doModes==3:
				print('Strengths fc=%g, ft=%g, |fc/ft|=%g'%(fc,ft,abs(fc/ft)))
			if doModes==2:
				print('Compressive strength fc=%g'%(abs(fc)))
			if doModes==1:
				print('Tensile strength ft=%g'%(abs(ft)))
			title=O.tags['description'] if 'description' in list(O.tags.keys()) else O.tags['params']
			print('Bye.')
			O.pause()
		
def addPlotData():
	yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,
		'sigma.25':forcesOnCoordPlane(coord_25,axis)[axis]/area_25+isoPrestress,
		'sigma.50':forcesOnCoordPlane(coord_50,axis)[axis]/area_50+isoPrestress,
		'sigma.75':forcesOnCoordPlane(coord_75,axis)[axis]/area_75+isoPrestress,
		})
initTest()
waitIfBatch()
###############################################

if you run this code on 'yade', it'll yield the poisson value of
0.16531931 but if I run it on 'yade-batch' the obtained value is
0.23966927.

Is there a cause for this difference, since the parameters are the same?

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.