Bug #923929 “OMP force container fails to reset forces on clump ...” : Bugs : Yade

Revision history for this message

Bruno Chareyre (bruno-chareyre) wrote on 2012-01-30:

#1

Workarounded in bzr3009.
Would still be nice to really fix.

Changed in yade:
importance:	Undecided → Wishlist
status:	New → Confirmed

Revision history for this message

Anton Gladky (gladky-anton) wrote on 2012-04-13:

#2

Hi Bruno,

I was trying to reproduce the bug last night, but could not.
Please, confirm and close the bug, or send a test-script.

Thanks.

Anton

Changed in yade:
status:	Confirmed → Incomplete

Revision history for this message

Chiara Modenese (chiara-modenese) wrote on 2012-04-13: Re: [Yade-dev] [Bug 923929] Re: OMP force container fails to reset forces on clump members

#3

Hi,
I have been using clumps for a while but never found a problem with
that. I am now using a release which is a bit older than the r3009 where
this bug seems to be fixed. I never used gravity engine too. Could you
please tell me how to reproduce the bug with a simple script?

Thank you.
Chiara

On 30/01/2012 21:42, Chareyre wrote:
> Workarounded in bzr3009.
> Would still be nice to really fix.
>
> ** Changed in: yade
> Importance: Undecided => Wishlist
>
> ** Changed in: yade
> Status: New => Confirmed
>

Revision history for this message

Bruno Chareyre (bruno-chareyre) wrote on 2012-04-13:

#4

It is fixed with the code below, which introduces a loop in Newton just to apply 0-force on each clump.
This is how I made the bug disapear, but I would not really call this a fix...
We can keep this in wishlist.

+#ifdef YADE_OPENMP
+void NewtonIntegrator::ensureSync()
+{
+ if (syncEnsured) return;
+ YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
+ if(b->isClump()) continue;
+ scene->forces.addForce(b->getId(),Vector3r(0,0,0));
+ } YADE_PARALLEL_FOREACH_BODY_END();
+ syncEnsured=true;
+}
+#endif
+
void NewtonIntegrator::action()
{
+ #ifdef YADE_OPENMP
+ //prevent https://bugs.launchpad.net/yade/+bug/923929
+ ensureSync();
+ #endif

Revision history for this message

Bruno Chareyre (bruno-chareyre) wrote on 2012-04-13:

#5

For testing the bug, first comment the ensureSync() line. Second, ask Klaus's script. I couldn't find it.

Revision history for this message

Bruno Chareyre (bruno-chareyre) wrote on 2012-11-28:

#6

Download full text (3.5 KiB)

Ok, the fix was in fact not sufficient, as shown by François Kneib with the script below.
An improved fix is in rev. fc20180e55ae1bb774fdce

To see the bug, run the script below after commenting
https://github.com/yade/trunk/blob/fc20180e55ae1bb774fdce750300660af7bad76a/pkg/dem/NewtonIntegrator.cpp#L101

#!/usr/local/bin/yade-trunk -x
# -*- coding: utf-8 -*-
# encoding: utf-8
from yade import pack,ymport,export,geom,bodiesHandling, plot
import math
from yade import qt
import numpy as np

O.periodic=False
damp=0.1
compFricDegree=30.
finalFricDegree=30.
length=1.2
height=1.8
width=1.2
thickness=0.00001
radius=1.3*0.07
EY=1e6 #Young modulus
stressGoal=30
stressTolerance=0.01
vol=length*height*width

#O.cell.hSize=Matrix3( length, 0, 0,
# 0, 3.*height, 0,
# 0, 0, width)

O.materials.append(CohFrictMat(isCohesive=True,density=1e3,young=EY,poisson=0.3,momentRotationLaw=1,frictionAngle=radians(finalFricDegree),normalCohesion=10e10,shearCohesion=10e10,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=1e3,young=EY,poisson=0.5,momentRotationLaw=1,frictionAngle=radians(compFricDegree),normalCohesion=10e7,shearCohesion=10e7,label='sphereMat'))

upBox = utils.box( center=(0,2*height+thickness/2.0,0),
     extents=(length*100.0,thickness/2.0,width*100.0),
     fixed=1,
     wire=False,
     material='boxMat')

lowBox = utils.box( center=(0,height-thickness/2.0,0),
     extents=(length*100.0,thickness/2.0,width*100.0),
     fixed=1,
     wire=False,
     material='boxMat')

O.bodies.append([upBox,lowBox])
sp=pack.SpherePack()
c1=pack.SpherePack([((0,0,0),radius),((radius,0,0),radius)])
sp.makeClumpCloud((0.,height+thickness/2.,0.),(length,2.*height-thickness/2.,width),[c1],periodic=False)
bds=sp.toSimulation(material='sphereMat')

O.usesTimeStepper=False
O.dt=1e-4

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()],allowBiggerThanPeriod=True),
InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label='ip2')],[Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]),
# GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.7,defaultDt=utils.PWaveTimeStep()),
NewtonIntegrator(damping=damp,gravity=[0.,-10.,0.],exactAsphericalRot=0),
]

def plotStress():
stress = getStress(vol)
plot.addData(
  iter=O.iter,
  stress00=(stress[0,0]),
  stress01=(stress[0,1]),
  stress02=(stress[0,2]),
  stress10=(stress[1,0]),
  stress11=abs(stress[1,1]),
  stress12=(stress[1,2]),
  stress20=(stress[2,0]),
  stress21=(stress[2,1]),
  stress22=abs(stress[2,2]),
  f_total0=-O.forces.f(0)[0]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]),
  f_total1=-O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]),
  f_total2=-O.forces.f(0)[2]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]),
  gamma=O.bodies[0].state.pos[0]/O.bodies[0].state.pos[1],
  volume=vol,
  #UBF=utils.unbalancedForce(),
  y=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])

#O.engines = O.engines+[PyRunner(command='plotStress()',iterPeriod=100)]
plot.plots={'iter':('stress11','stress22',None,'volume')}
#plot.plot()

phase=0
def main...

Ok, the fix was in fact not sufficient, as shown by François Kneib with the script below.
An improved fix is in rev. fc20180e55ae1bb774fdce

To see the bug, run the script below after commenting
https://github.com/yade/trunk/blob/fc20180e55ae1bb774fdce750300660af7bad76a/pkg/dem/NewtonIntegrator.cpp#L101

#!/usr/local/bin/yade-trunk -x
# -*- coding: utf-8 -*-
# encoding: utf-8
from yade import pack,ymport,export,geom,bodiesHandling, plot
import math
from yade import qt
import numpy as np

O.periodic=False
damp=0.1
compFricDegree=30.
finalFricDegree=30.
length=1.2
height=1.8
width=1.2
thickness=0.00001
radius=1.3*0.07
EY=1e6 #Young modulus
stressGoal=30
stressTolerance=0.01
vol=length*height*width

#O.cell.hSize=Matrix3(	length,	0,			0,
#						0,		3.*height,	0,
#						0,		0,			width)

O.materials.append(CohFrictMat(isCohesive=True,density=1e3,young=EY,poisson=0.3,momentRotationLaw=1,frictionAngle=radians(finalFricDegree),normalCohesion=10e10,shearCohesion=10e10,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=1e3,young=EY,poisson=0.5,momentRotationLaw=1,frictionAngle=radians(compFricDegree),normalCohesion=10e7,shearCohesion=10e7,label='sphereMat'))

upBox = utils.box( 	center=(0,2*height+thickness/2.0,0),
					extents=(length*100.0,thickness/2.0,width*100.0),
					fixed=1,
					wire=False,
					material='boxMat')
					
lowBox = utils.box(	center=(0,height-thickness/2.0,0),
					extents=(length*100.0,thickness/2.0,width*100.0),
					fixed=1,
					wire=False,
					material='boxMat')

O.bodies.append([upBox,lowBox])
sp=pack.SpherePack()
c1=pack.SpherePack([((0,0,0),radius),((radius,0,0),radius)])
sp.makeClumpCloud((0.,height+thickness/2.,0.),(length,2.*height-thickness/2.,width),[c1],periodic=False)
bds=sp.toSimulation(material='sphereMat')

O.usesTimeStepper=False
O.dt=1e-4

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()],allowBiggerThanPeriod=True),
	InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label='ip2')],[Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]),
#	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.7,defaultDt=utils.PWaveTimeStep()),
	NewtonIntegrator(damping=damp,gravity=[0.,-10.,0.],exactAsphericalRot=0),
]

def plotStress():
	stress = getStress(vol)
	plot.addData(
		iter=O.iter,
		stress00=(stress[0,0]),
		stress01=(stress[0,1]),
		stress02=(stress[0,2]),
		stress10=(stress[1,0]),
		stress11=abs(stress[1,1]),
		stress12=(stress[1,2]),
		stress20=(stress[2,0]),
		stress21=(stress[2,1]),
		stress22=abs(stress[2,2]),
		f_total0=-O.forces.f(0)[0]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]),
		f_total1=-O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]),
		f_total2=-O.forces.f(0)[2]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]),
		gamma=O.bodies[0].state.pos[0]/O.bodies[0].state.pos[1],
		volume=vol,
		#UBF=utils.unbalancedForce(),
		y=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])

#O.engines = O.engines+[PyRunner(command='plotStress()',iterPeriod=100)]
plot.plots={'iter':('stress11','stress22',None,'volume')}
#plot.plot()

phase=0
def main():
	global phase
	if phase==0:
		if(abs(getStress(vol)[1,1])<stressGoal-stressTolerance):growParticles(1.01,1,1)
		elif(abs(getStress(vol)[1,1])>stressGoal+stressTolerance):growParticles(0.99,1,1)
		elif(unbalancedForce()<1e-3):phase+=1
	if phase==1:
		phase+=1
		O.pause()
#O.engines = O.engines+[PyRunner(command='main()',iterPeriod=100)]

#O.run(10000,1)
O.saveTmp()

Revision history for this message

Reza Housseini (reza-housseini-3) wrote on 2013-12-17:

#7

Was there some progress with syncing the forces? My forces are still accumulating to NaN with the changes Bruno suggested.

Revision history for this message

Bruno Chareyre (bruno-chareyre) wrote on 2013-12-19:

#8

No progress since #6.
It should not be a problem for users. Why are you asking?

Revision history for this message

Reza Housseini (reza-housseini-3) wrote on 2013-12-23:

#9

I wrote a ClumpsFactory extension and for some values there seems to be a data race. After each iteration the forces are adapted from wrong clumps and are multiplied by a factor 9 (?). When I use

scene->forces.addForce(id, f);

instead of

scene->forces.getForceUnsynced(id) += f;

I don't encounter this problem. Therefore I concluded there is a flaw in the synchronization mechanism.

I also don't see how you can guarantee that there's no data race. STL containers aren't thread safe if you are reading and writing at the same time.

Revision history for this message

Bruno Chareyre (bruno-chareyre) wrote on 2013-12-23:

#10

addForce() is thread-safe because each thread is computing its own partial sum of forces. The total is sync'ed only once a step by NewtonIntegrator.
getForceUnsynced(id)+= is not thread safe, OTOH, all threads would access the same value. It should not be used except in a few exceptions.

Revision history for this message

Reza Housseini (reza-housseini-3) wrote on 2013-12-31:

#11

I'm still wondering about some parts in this code.

1. Why will adding a vector of zeros [1] helps resetting the forces of clumps?
2. Why is the collected force of the clumps members added to the ForceContainer [2] ...
3. and why is this force then again retaken from the ForceContainer? [3]

There is no need to update the forces in the ForceContainer, because there will be reset to zero at the next step. The force is only used to calculate the acceleration in [4].

Can someone please clarify these points?

[1] https://github.com/yade/trunk/blob/master/pkg/dem/NewtonIntegrator.cpp#L93
[2] https://github.com/yade/trunk/blob/master/pkg/dem/NewtonIntegrator.cpp#L150
[3] https://github.com/yade/trunk/blob/master/pkg/dem/NewtonIntegrator.cpp#L157
[4] https://github.com/yade/trunk/blob/master/pkg/dem/NewtonIntegrator.cpp#L178

Revision history for this message

Bruno Chareyre (bruno-chareyre) wrote on 2014-01-02:

#12

> 1. Why will adding a vector of zeros [1] helps resetting the forces of clumps?
Because the size of the force container is adjusted dynamically. Adding
a zero force to body "id" instructs the container that "id" is a valid
index, then it will be resetted, otherwise it is not known by the
container. This is weird, of course.

> 2. Why is the collected force of the clumps members added to the ForceContainer [2] ...
Why not? The container is supposed to contain total forces on bodies.
Thus, we need one force per clump. The line you are linking is
incrementing by the forces from each member of the clump.

> 3. and why is this force then again retaken from the ForceContainer? [3]
Because "total += f" (L150) does not mean total==f.
If we need the total we have to get it.
There can be forces applied directly on the clump (user defined forces,
for instance). That is why "total" and "f" may be different.

> There is no need to update the forces in the ForceContainer, because
> there will be reset to zero at the next step. The force is only used to
> calculate the acceleration in [4].
I don't understand this point. It doesn't seem correct to state that the
force is "only used in [4]".
If I write this line is a script, I will be using the force as well, but
this is not something you can anticipate looking at the c++ code:
f=O.forces.f[clumpId]
So we better have them defined correctly.

I hope it helps.

Bruno

Revision history for this message

Bruno Chareyre (bruno-chareyre) wrote on 2014-09-15:

#13

The bug disappeared, fixed in [1] most probably.
Reverting bzr3009, which is no longer needed [2].

[1] https://github.com/yade/trunk/commit/4ea76ad6e47ac5074a389ad61712a0840e8560a5
[2] https://github.com/yade/trunk/commit/4a52ebdd25d57ee348aec6375a65ff388ffd97be

Changed in yade:
status:	Incomplete → Fix Released

Yade

OMP force container fails to reset forces on clump members

Bug Description

Other bug subscribers

Remote bug watches