setPermF broken?

Asked by Daniel Kiracofe on 2020-05-13

I have a user-defined force that I need to apply to a clump. The force changes with time, but rather slowly relative to the time step. I noticed that the PyRunner to calculate and apply this force was taking 20 - 30% of the entire simulation time, so I thought to recalculate the force only every 100 iterations or so. I was previously using this to apply the force

O.forces.addF(ClumpId, ( 0, externalForce , 0) )

I switched it to this:

O.forces.setPermF(ClumpId, ( 0, externalForce , 0) )

The result is that the clump in question does not move. It's as if zero force was applied to it. I checked O.forces.permF(ClumpId), and it gives me the expected force, but the body does not move. I also tried addF with permanent=True, also adding the force directly to a single member of the clump instead of the clump directly. None of these work. Any suggestions? Or anybody can point me a to working example of setPermF?

I am using recent daily build (20200511-3819~5bf8512~buster1).

Here's the key part of the source:

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()], verletDist=verlet_dist ),
 InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
   [Ip2_FrictMat_FrictMat_FrictPhys()],
   [Law2_ScGeom_FrictPhys_CundallStrack()] ),
 timestepper,
 PyRunner(command='swing()',iterPeriod=100),
 newton,
   PyRunner(command='addPlotData()',iterPeriod=plotDataPeriod),

]

def swing():
        global theta

        theta = 2 * math.pi * ( f_sweep_start * O.time + f_sweep_rate * O.time**2)
        extF = F_amp * math.sin(theta)

        kF = -k * ( wall_state.pos[1]-x0 )
        cF = -c * ( wall_state.vel[1] )

# O.forces.addF( ClumpId, ( 0, extF + kF + cF , 0) ) # works
        O.forces.setPermF(ClumpId, ( 0, extF + kF +cF , 0) ) # does nothing?

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bruno Chareyre
Solved:
2020-05-16
Last query:
2020-05-16
Last reply:
2020-05-15
Jan Stránský (honzik) said : #1

Hello,

> Here's the key part of the source:
> Any suggestions?

please provide more information:
- a complete code reproducing your problem (trying to keep it minimal) [1]
(The code provided is not really "key" - engines are mostly irrelevant, computations inside swing are not relevant either).
- how you start yade (e.g. -j option)

> Or anybody can point me a to working example of setPermF?

a MWE [1]:
###
O.bodies.append([sphere((0,0,z),.1) for z in range(100)])
O.forces.setPermF(0,(1,0,0))
O.run(10,True)
print(O.bodies[0].state.pos)
###

works nicely with "yade", gives zero motion with "yade -j 3".
So it seems there is a bug in setPermF for parallel run, compare [2] and [3].

@Bruno? [4]

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/core/ForceContainerSerial.cpp#L45
[3] https://gitlab.com/yade-dev/trunk/-/blob/master/core/ForceContainerParallel.cpp#L104
[4] https://gitlab.com/yade-dev/trunk/-/commit/da054fee5fe98ff2ebe46d45dfdf839d5b2bb9eb

Daniel Kiracofe (kiracodl) said : #2

Jan, I appreciate the quick reply. Here are some more details.

I was originally running with "yade -j 4". I tried it with just "yade", but that did not change anything.

The original code was quite complicated. I tried to create a stripped down version of my code to produce the minimum possible bug report. I stripped it down to code #1 below. In that code, I create 4 walls which are clumped together. if I apply the setPermF to the clump, it does not work. If I apply the setPermF to an individual wall that is a member of the clump, it does work.

However, in my original code, neither one of these worked. So I went back to my original code again, and started stripping it down again to try to get the minimum possible bug report for a single body that is not part of a clump. I got code #2 below. In this code, I create a single wall, and then either 0, 1 or a bunch of spheres. If I create a bunch of spheres, then setPermF (on the wall) does not work. If I create 0 or 1 sphere, then it does work. Not quite sure what is going on here.

So, unless I'm doing something wrong or misunderstanding, it looks like there may be three separate bugs with setPermF
1) the parallel vs single-threaded issue raised above
2) the clump vs single body issue demonstrated in code #1 below
3) the still unknown issue demonstrated in code #2 below

Please let me know if you need any more details.

Daniel

#########################################
#code #1
#in this version, setPermF works if I apply to one of the members of the clump
#but does not work if apply it directly to the clump

###############################################3
box_x = 0.00025
box_y = 0.0013
box_z = box_y
f_sweep_start = 45.5 #
f_sweep_stop = 52 #Hertz
f_sweep_rate = 1 #Hertz/second
out_file_name = "TestOutput.txt"
F_amp = 4e-6
grav_y=0
damp = 0.000
ball_gen_method=1

#################################################

import qt, plot
qt.View() #open the controlling and visualization interfaces

finalFricDegree = atan(0.25) #contact friction during the deviatoric loading

young2 = 200e9 #the walls' Young's modulus, same as the spheres
rho = 8230 #

mn = Vector3(0, box_y, 0)
mx = Vector3(box_x,2*box_y, box_z) #corners of the initial packing
thick = 0.010 # the thickness of the walls

nat_freq_hz = 50.0
nat_period = 1 / nat_freq_hz
nat_freq_rad = nat_freq_hz * 2 * math.pi

#GENERATING BODIES #####

O.dt = nat_period / 100
plotDataPeriod = 10

#special 4 walls case for 2D because of all of the blocked DOF
expand=1.5

walls_density = 8000

wallIds = []

O.materials.append(FrictMat(young=young2,poisson=0.7,frictionAngle=0, density=walls_density , label='walls'))

# yade.utils.box( center extents= half sizes!
topWall = yade.utils.box( (box_x/2,2*box_y+thick/2,box_z/2), (box_x/2*expand,thick/2,box_z/2*expand), material='walls')
topWallId = O.bodies.append(topWall)
wallIds.append(topWallId)

botWall = yade.utils.box( (box_x/2,box_y-thick/2,box_z/2), (box_x/2*expand,thick/2, box_z/2*expand), material='walls')
botWallId = O.bodies.append(botWall)
wallIds.append(botWallId)

rightWall = yade.utils.box( (box_x+thick/2,1.5*box_y,box_z/2), (thick/2,box_y/2*expand, box_z/2*expand), material='walls')
rightWallId = O.bodies.append(rightWall)
wallIds.append(rightWallId)

leftWall = yade.utils.box( (-thick/2,1.5*box_y,box_z/2), (thick/2,box_y/2*expand, box_z/2*expand), material='walls')
leftWallId = O.bodies.append(leftWall)
wallIds.append(leftWallId)

wallClump=O.bodies.clump(wallIds)

for i in wallIds:
        O.bodies[i].state.blockedDOFs = 'xzXYZ' #fixme is this needed with clumping?

#DEFINING ENGINES ####

newton=NewtonIntegrator(damping=damp, gravity=(0,grav_y,0))

O.engines=[
 ForceResetter(),
 PyRunner(command='swing()',iterPeriod=100),
 newton,
   PyRunner(command='addPlotData()',iterPeriod=plotDataPeriod),
]

def swing():
        force = 5
        #this works
        #O.forces.addF( wallClump, ( 0, force , 0) )

        #this does not work
        O.forces.setPermF( wallClump , ( 0, force , 0) )

        #this does work
        #O.forces.setPermF( topWallId , ( 0, force , 0) )

def addPlotData():
        position=O.bodies[topWallId].state.pos[1]

        plot.addData(freq=O.time , position=position)

plot.plots={'freq':('position')}
plot.plot()

###################################################################################
###################################################################################
#code #2
from yade import pack

#this variable selection between three cases
#ball_gen_method=0, generates one wall and zero spheres, then setPermF on the wall works
#ball_gen_method=1, generates one wall and one spheres, then setPermF on the wall works
#ball_gen_method=2, generates one wall and many spheres, then setPermF on the wall does NOT work
ball_gen_method=1

###############################################3
box_x = 0.00025
box_y = 0.0013
f_sweep_start = 45.5 #Hertz
f_sweep_stop = 52 #Hertz
f_sweep_rate = 1 #Hertz/second
periodic_bc = 2
particle_dia = 2.5e-5
grav_y=0
damp = 0.000

box_z = particle_dia*1.01

#################################################
##### DEFINING VARIABLES AND MATERIALS ######
#################################################

import qt, plot
qt.View() #open the controlling and visualization interfaces

finalFricDegree = atan(0.25) #contact friction during the deviatoric loading

young1 = 200e9 #the spheres' Young's modulus, Pa. AFRL published papers used Inco718. young's modulus about 29e6 psi at room temp, or 200 GPa
young2 = 200e9 #the walls' Young's modulus, same as the spheres
rho = 8230 #generic inco 718, kg/m^3

mn = Vector3(0, box_y, 0)
mx = Vector3(box_x,2*box_y, box_z) #corners of the initial packing
thick = 2*particle_dia # the thickness of the walls

dia = [particle_dia, particle_dia * 1.001 ] #these are the points on a distribution, but for now we basically have them all the same size
phi = [0.0, 1.0]

O.materials.append(FrictMat(young=young1,poisson=0.7,frictionAngle=finalFricDegree,density=rho,label='spheres'))

#GENERATING BODIES #####

ballIds = ()
if (ball_gen_method==2):
        #for this setting, setPermF does not work
        sp=pack.SpherePack()
        sp.makeCloud(minCorner=mn,maxCorner=mx,periodic=False,psdSizes=dia,psdCumm=phi,distributeMass=True,seed=1)

        ballIds = O.bodies.append([sphere(center,rad,material='spheres',color=(0,0,1)) for center,rad in sp])
elif (ball_gen_method == 1):
        #only only 1 single ball.
        #with this setting setPermF does work
        center = ( (mn[0]+mx[0])/2, (mn[1]+mx[1])/2, (mn[2]+mx[2])/2)
        rad = 0.01
        ballIds = O.bodies.append(sphere(center,rad,material='spheres',color=(0,0,1)))
elif (ball_gen_method == 0):
        #turn off all spheres. this is only for debugging
        # for this setting, setPermF does work
        print('this space intentionally left blank')

timestepper = GlobalStiffnessTimeStepper(active=False,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8)
O.dt = 2e-9

plotDataPeriod = 10
expand=1.5

walls_density = 8000

O.materials.append(FrictMat(young=young2,poisson=0.7,frictionAngle=0, density=walls_density , label='walls'))

# yade.utils.box( center extents= half sizes!
topWall = yade.utils.box( (box_x/2,2*box_y+thick/2,box_z/2), (box_x/2*expand,thick/2,box_z/2*expand), material='walls')
topWallId = O.bodies.append(topWall)

x0 = O.bodies[topWallId].state.pos[1]

newton=NewtonIntegrator(damping=damp, gravity=(0,grav_y,0))

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()] ),
 timestepper,
 PyRunner(command='swing()',iterPeriod=100),
 newton,
   PyRunner(command='addPlotData()',iterPeriod=plotDataPeriod),
]

def swing():
        F = 100

#this works, always
# O.forces.addF(topWallId , ( 0, F , 0))
#works when there is only 0 or 1 sphere, but not when there are a bunch
        O.forces.setPermF(topWallId , ( 0, F , 0) )

def addPlotData():
        position=O.bodies[topWallId].state.pos[1]-x0

        plot.addData(freq=O.time , position=position)

        if (O.time > 8):
                O.pause()

plot.plots={'freq':('position')}
plot.plot()

Robert Caulk (rcaulk) said : #3

> a MWE
>works nicely with "yade", gives zero motion with "yade -j 3".
>So it seems there is a bug in setPermF for parallel run, compare [2] and [3].

Hey Jan,

I tried to reproduce your bug and I was unable to. In both cases the body moves (running yadedaily latest), but in single core the body moves much more than in parallel. Something is certainly wrong with one or the other. I modified your MWE a bit to track some stats:

########################################
O.bodies.append([sphere((0,0,z),.1) for z in range(10)])
O.forces.setPermF(0,(1,0,0))

def printStats():
 print('pos',O.bodies[0].state.pos)
 print('permF',O.forces.permF(0))
 print('force',O.forces.f(0))
 print('vel',O.bodies[0].state.vel)

O.engines=O.engines+[PyRunner(dead=0,iterPeriod=1000,command='printStats()',label='stats')]

newton.damping=0

O.run(200000,True)
########################################

Observations:

* The single core run accelerates the body much more quickly than the multi core run, perhaps the permanent force is being added over and over for single core runs, while remaining constant in multi core runs. I believe this was the issue that Bruno's patch [1] intended to fix for multi core runs. I am unsure why he did not apply the patch to the single core force container.
* O.forces.f(0) reports 0 (unexpected) in multicore runs and 1 (expected value) in single core runs.

[1] https://gitlab.com/yade-dev/trunk/-/commit/da054fee5fe98ff2ebe46d45dfdf839d5b2bb9eb

Jérôme Duriez (jduriez) said : #4

An issue maybe ? ;-)

https://gitlab.com/yade-dev/trunk/-/issues/156 with my observations (more or less yours, guys)

Chareyre (bruno-chareyre-9) said : #5

I suspect fixing one problem recently could have cause that new one. Thanks
for reporting, please provide the info you can. I'll be back.
Bruno

Le jeu. 14 mai 2020 00:55, Jan Stránský <
<email address hidden>> a écrit :

> Question #690691 on Yade changed:
> https://answers.launchpad.net/yade/+question/690691
>
> Status: Open => Answered
>
> Jan Stránský proposed the following answer:
> Hello,
>
> > Here's the key part of the source:
> > Any suggestions?
>
> please provide more information:
> - a complete code reproducing your problem (trying to keep it minimal) [1]
> (The code provided is not really "key" - engines are mostly irrelevant,
> computations inside swing are not relevant either).
> - how you start yade (e.g. -j option)
>
> > Or anybody can point me a to working example of setPermF?
>
> a MWE [1]:
> ###
> O.bodies.append([sphere((0,0,z),.1) for z in range(100)])
> O.forces.setPermF(0,(1,0,0))
> O.run(10,True)
> print(O.bodies[0].state.pos)
> ###
>
> works nicely with "yade", gives zero motion with "yade -j 3".
> So it seems there is a bug in setPermF for parallel run, compare [2] and
> [3].
>
> @Bruno? [4]
>
> cheers
> Jan
>
> [1] https://www.yade-dem.org/wiki/Howtoask
> [2]
> https://gitlab.com/yade-dev/trunk/-/blob/master/core/ForceContainerSerial.cpp#L45
> [3]
> https://gitlab.com/yade-dev/trunk/-/blob/master/core/ForceContainerParallel.cpp#L104
> [4]
> https://gitlab.com/yade-dev/trunk/-/commit/da054fee5fe98ff2ebe46d45dfdf839d5b2bb9eb
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

A dumb mis-placement of that line 196 [1] in a parallel loop. The result was actually non-deterministic.
That was nasty bug... thanks.
The fix is in branch fixPermForcesMultithread [2], hopefully merged soon.
Cheers
Bruno

[1] https://gitlab.com/yade-dev/trunk/-/merge_requests/481/diffs#f91335edf66ecbb442af1573f3370b076ad3f416_196_196
[2] https://gitlab.com/yade-dev/trunk/-/tree/fixPermForcesMultithread

Daniel Kiracofe (kiracodl) said : #7

Bruno, I have compiled branch fixPermForcesMultithread and verified that it does fix the bug I was experiencing.

Thank you all for the help!

Daniel

Daniel Kiracofe (kiracodl) said : #8

Thanks Bruno Chareyre, that solved my question.