error in porosity after erasing bodies

Bug #1774926 reported by azim on 2018-06-04
6
This bug affects 1 person
Affects Status Importance Assigned to Milestone
Yade
Undecided
Unassigned

Bug Description

Dear all
I had a problem in [1] about updating len(O.bodies) after erasing some particles.
after erasing I summoned utils.porosity and triax.porosity to check te porosity.
the former was OK but the latter was negative and it is wrong (porosity<0).
that's beacause triax.porosity (triax.boxVolume-triax.particlesVolume/triax.boxVolume) doesn't consider the erasing process and uses the particlesVolume equalse to the beggining of the simulation. sp.dim() , sp.center() and sp.aabb() have the same behaviour.
I'm using ubuntu 14.0d LTS and yadedaily 2018.02b-12211bed66~trusty.
[1] https://answers.launchpad.net/yade/+question/669845

Thanks in advance.
Azim

Script:

# -*- coding: utf-8 -*-
from yade import pack
num_spheres=1000
compFricDegree = 30
targetPorosity = 0.43
finalFricDegree = 30
rate=-0.02
damp=0.2
stabilityThreshold=0.01
young=5e6
mn,mx=Vector3(0,0,0),Vector3(1,1,1)
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()

clumps=False
if clumps:
 volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
 mean_rad = pow(0.09*volume/num_spheres,0.3333)
 c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
 sp.makeClumpCloud(mn,mx,[c1],periodic=False)
 sp.toSimulation(material='spheres')
 O.bodies.updateClumpProperties()
else:
 sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1)
 O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
triax=TriaxialStressController(
 maxMultiplier=1.0+2e4/young,
 finalMaxMultiplier=1.0+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]
triax.goal1=triax.goal2=triax.goal3=-10000

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
    break

triax.internalCompaction=False
for b in O.bodies:
 if b.id>5:
  if b.state.pos[0]<.5:
   O.bodies.erase(b.id)

for i in range(1,200):
     print "utils.porosity= ", porosity()
     print "TriaxialStressController.porosity= ", triax.porosity
     O.run(1,1)

Bruno Chareyre (bruno-chareyre) wrote :

Hi,
If solid volume is changed you should use TriaxialStressController.updatePorosity [1]
It would not be good to recalculate solid volume everytime by default.
Bruno

[1] https://yade-dem.org/doc/yade.wrapper.html?highlight=updateporosity#yade.wrapper.TriaxialStressController.updatePorosity

Changed in yade:
status: New → Invalid
azim (mirzavand) wrote :

Hi Bruno
thank you for your answer.
I looked at [1], I don't get it.
>>you should use TriaxialStressController.updatePorosity [1]
when and how do we use it ? it returns False!!
I mean, do we need to use it in a special position in the code?
Thanks
Azim

azim (mirzavand) wrote :

Hi Bruno,
Thank you.
I got it:

TriaxialStressController.updatePorosity=True
print "porosity", TriaxialStressController.porosity

thanks again for helping me.
Azim

To post a comment you must log in.
This report contains Public information  Edit
Everyone can see this information.

Other bug subscribers