polyhedrals after collide to facet, throwing away strongly

Hi,
this is a typical symptom of too large time step (especially if decreasing stiffness helps), try to decrease it.
cheers
Jan

Thank you a lot for all your helps Jan Stransky

>> this is a typical symptom of too large time step
Do you think if I decrease time step, maybe this problem is solved?
for more particles, I set time step:10e-8.

 >> note: Most of polyhedrals have sharp angle, does this have effect on this problem?

Do you have any idea about this??

Actually I think contact model in polyhedrals has problem!As the two polyhedrals collide very slowly, their velocity and kinetic energy increase abnormally!!

>> this is a typical symptom of too large time step (especially if decreasing stiffness helps)
Unfortunately decreasing stiffness or time step did help at all. :(

Your help in this would be appreciated

> note: Most of polyhedrals have sharp angle, does this have effect on this problem?

no

> decreasing young module until 10e3 can help the situation
> Unfortunately decreasing stiffness or time step did help at all. :(

please be consistent. So does it help or not?

> for more particles, I set time step:10e-8.

Why? the value of time step should not depend on number of particles..

> Actually I think contact model in polyhedrals has problem!As the two polyhedrals collide very slowly, their velocity and kinetic energy increase abnormally!!

can you please provide a Minimal Working Example (e.g. two polhedrons) reproducing this? thanks

cheers
Jan

> please be consistent. So does it help or not?

("decreasing young module until 10e3 can help the situation"): this was a mistake and decreasing stiffness or time step did help at all.
If stiffness reaches less than 10e5, then particles are passed through facet.

> Why? the value of time step should not depend on number of particles..
Appropriate time step for simulation with all particles is 10e-8 to avoid exploding particles, in this test (for saving time) 10e-5 was selected for time step because particles number are small and they do not explode, whereas I checked, changing time step does not change situation (does not solve this problem).

Hello,

I have tried the script provided in OP, resulting in:
NameError: name 'ListVer2' is not defined
could you please provide a test script which we can try? thanks

cheers
Jan

In advance thank you for your time

> can you please provide a Minimal Working Example (e.g. two polhedrons) reproducing this?

Here is small data and input for my code, please copy in file:

*LWPOLYLINE
4.54097400033104,6.49184245861868
4.65891804001884,6.58769165827740
4.51425920185592,6.93490024048091
4.46361293495481,6.82361276486310
4.54097400033104,6.49184245861868
*LWPOLYLINE
3.38239216165079,4.01327469219457
3.43869334720801,4.27025377921128
3.31240868369778,4.27786275211975
3.18944272805255,4.21917403263651
3.17292105915418,4.00731315957292
3.38239216165079,4.01327469219457
*LWPOLYLINE
3.07472367219926,4.22146522059046
3.00951705252640,4.00710000000000
3.15872105915418,4.00710000000000
3.16994885128804,4.22031962661348
3.07472367219926,4.22146522059046
*LWPOLYLINE
3.19496730804472,4.54686225863838
3.10768290255748,4.56486506811422
3.08101342018444,4.41477399856679
3.07472367219926,4.23195436841659
3.17133336059188,4.25583338934358
3.19496730804472,4.54686225863838
*LWPOLYLINE
3.71267645811818,4.17335865185637
3.80626391854923,4.37990041469040
3.85709007230427,4.37990041469040
3.81834238444732,4.15093746143616
3.71267645811818,4.17335865185637
*LWPOLYLINE
3.86833315936705,7.18115411676252
3.75805619368548,6.85052980253022
3.96022603149273,6.83823476145616
3.98624257576362,7.12695113182557
3.86833315936705,7.18115411676252
*LWPOLYLINE
3.75705101349167,7.68673972685043
3.75705101349167,8.04757108283743
3.89622020287160,8.25921672832907
4.02762326650518,8.25921672832907
4.02762326650518,8.09425875012654
3.93578182400268,7.82502476121894
3.75705101349167,7.68673972685043
*LWPOLYLINE
4.29367658088641,7.54996433412856
4.25127773179602,7.75135995952029
4.17124793309971,7.75135995952029
4.11338686757211,7.59716107043463
4.21414905751508,7.54996433412856
4.29367658088641,7.54996433412856
*LWPOLYLINE
3.60265201447845,7.35103416471193
3.61010182622721,7.35103416471193
3.71828876180858,7.48918606800436
3.74275101349167,7.67253715703928
3.74275101349167,7.82185097460873
3.67394962688404,7.75228435593977
3.60265201447845,7.35103416471193
*LWPOLYLINE
3.75665696623051,7.66810154888463
3.94362927435906,7.81276307636047
4.02303237072992,7.79446112714408
3.89479310068892,7.42404447036142
3.73297934687869,7.49095982637576
3.75665696623051,7.66810154888463

Thanks for the working script. I have tried it. With O.dt=10e-5, I see clearly what was described.
Again, it is a typical symptom of too large time step (making the simulation unstable) and I don't believe it is not overcome with decreasing the time step to a **sufficient** value

> Unfortunately decreasing ... time step did help at all. :(

what values have you tried?

> Actually I think contact model in polyhedrals has problem!As the two polyhedrals collide very slowly, their velocity and kinetic energy increase abnormally!!

Can you please provide a script (two bodies as is described) reproducing this behavior?
(even the number of bodies is decreased, still the simulation is far from minimal to show and study such phenomenons in detail)

cheers
Jan

Thanks Jan for your time

> what values have you tried?
O.dt=10e-8

> Can you please provide a script (two bodies as is described) reproducing this behavior?

from yade import export,polyhedra_utils
import os
from yade import plot
import math
from yade import utils
import pylab
import matplotlib; matplotlib.rc('axes',grid=True)
from matplotlib import pyplot
import numpy as np
from numpy import *
from yade import export as expt

Dolomite = PolyhedraMat()
Dolomite.density = 2870 #kg/m^3
Dolomite.young = 24.36e9 #Pa
Dolomite.poisson = 0.2
Dolomite.frictionAngle = radians(55.12) #rad

Bound = PolyhedraMat()
Bound.density = 2870 #kg/m^3
Bound.young = 60e10 #Pa
Bound.poisson = 0.2
Bound.frictionAngle = radians(55.12) #rad

p1=((4.54097400033104,6.49184245861868,0.25),(4.65891804001884,6.58769165827740,0.25),(4.51425920185592,6.93490024048091,0.25),(4.46361293495481,6.82361276486310,0.25),(4.54097400033104,6.49184245861868,0.25),(4.54097400033104,6.49184245861868,-0.25),(4.65891804001884,6.58769165827740,-0.25),(4.51425920185592,6.93490024048091,-0.25),(4.46361293495481,6.82361276486310,-0.25),(4.54097400033104,6.49184245861868,-0.25))
p2=((3.60265201447845,7.35103416471193,0.25),(3.61010182622721,7.35103416471193,0.25),(3.71828876180858,7.48918606800436,0.25),(3.74275101349167,7.67253715703928,0.25),(3.74275101349167,7.82185097460873,0.25),(3.67394962688404,7.75228435593977,0.25),(3.60265201447845,7.35103416471193,0.25),(3.60265201447845,7.35103416471193,-0.25),(3.61010182622721,7.35103416471193,-0.25),(3.71828876180858,7.48918606800436,-0.25),(3.74275101349167,7.67253715703928,-0.25),(3.74275101349167,7.82185097460873,-0.25),(3.67394962688404,7.75228435593977,-0.25),(3.60265201447845,7.35103416471193,-0.25))

O.bodies.append(polyhedra_utils.polyhedra(Dolomite,v=p1,fixed=False, color=(0.1,0.5,0.2)))
O.bodies.append(polyhedra_utils.polyhedra(Dolomite,v=p2,fixed=False, color=(0.1,0.5,0.2)))

v1=((-2.4993,-0.01,2),(-2,-0.01,2),(-2.4993,-0.01,-2))
v2=((-2,-0.01,2),(-2,-0.01,-2),(-2.4993,-0.01,-2))
v3=((-1.99,-0.01,2),(-1.99,-0.01,-2),(-1,-0.01,2))
v4=((-1,-0.01,2),(-1,-0.01,-2),(-1.99,-0.01,-2))
v5=((-0.99,-0.01,2),(-0.99,-0.01,-2),(1,-0.01,2))
v6=((1,-0.01,2),(1,-0.01,-2),(-0.99,-0.01,-2))
v7=((1.01,-0.01,2),(1.01,-0.01,-2),(2,-0.01,2))
v8=((2,-0.01,2),(2,-0.01,-2),(1.01,-0.01,-2))
v9=((2.01,-0.01,2),(2.01,-0.01,-2),(3,-0.01,2))
v10=((3,-0.01,2),(3,-0.01,-2),(2.01,-0.01,-2))
v11=((3.01,-0.01,2),(3.01,-0.01,-2),(4.58,-0.01,2))
v12=((4.58,-0.01,2),(4.58,-0.01,-2),(3.01,-0.01,-2))
v13=((4.6162,0.0059,2),(5.25,0.93,2),(4.6162,0.0059,-2))
v14=((5.25,0.93,2),(5.25,0.93,-2),(4.6162,0.0059,-2))
v15=((5.25,0.94,2),(5.25,0.94,-2),(6.42,2.57,2))
v16=((6.42,2.57,2),(6.42,2.57,-2),(5.25,0.94,-2))
v17=((6.42,2.58,2),(6.42,2.58,-2),(6.7906,3.01,2))
v18=((6.7906,3.01,2),(6.7906,3.01,-2),(6.42,2.58,-2))
v19=((6.7906,3.0162,2),(6.7939,4,2),(6.7906,3.0162,-2))
v20=((6.7906,3.0162,-2),(6.7939,4,-2),(6.7939,4,2))
v21=((6.7953,4.0079,2),(7.02,5.6,2),(6.7953,4.0079,-2))
v22=((7.02,5.6,2),(7.02,5.6,-2),(6.7953,4.0079,-2))
v23=((7.02,5.61,2),(7.02,5.61,-2),(7.3,7.01,2))
v24=((7.3,7.01,2),(7.3,7.01,-2),(7.02,5.61,-2))
v25=((7.3,7.02,2),(7.3,7.02,-2),(7.54,8.51,2))
v26=((7.54,8.51,2),(7.54,8.51,-2),(7.3,7.02,-2))
v27=((7.54,8.52,2),(7.54,8.52,-2),(7.8,10,2))
v28=((7.8,10,2),(7.8,10,-2),(7.54,8.52,-2))
v29=((7.8,10.01,2),(7.8,10.01,-2),(8.14,12,2))
v30=((8.14,12,2),(8.14,12,-2),(7.8,10.01,-2))
v31=((8.14,12.01,2),(8.14,12.01,-2),(8.5,14,2))
v32=((8.5,14,2),(8.5,14,-2),(8.14,12,-2))
v33=((8.5,14.01,2),(8.5,14.01,-2),(8.75,15.5,2))
v34=((8.75,15.5,2),(8.75,15.5,-2),(8.5,14.01,-2))
v35=((8.75,15.51,2),(8.75,15.51,-2),(9.0615,16.9981,2))
v36=((9.0615,16.9981,2),(9.0615,16.9981,-2),(8.75,15.51,-2))
v37=((9.0621,17.0187,2),(9.09,19,2),(9.0621,17.087,-2))
v38=((9.09,19,2),(9.09,19,-2),(9.062,17.0187,-2))
v39=((9.09,19.01,2),(9.09,19.01,-2),(9.2006,20.9992,2))
v40=((9.2006,20.9992,2),(9.2006,20.9992,-2),(9.09,19.01,-2))
v41=((2.2009,21.0924,2),(2.2009,21.0924,-2),(2.22,19.01,2))
v42=((2.22,19.01,2),(2.22,19.01,-2),(2.2009,21.0924,-2))
v43=((2.22,19,2),(2.22,19,-2),(2.1538,17.043,2))
v44=((2.1538,17.043,2),(2.1538,17.043,-2),(2.22,19,-2))
v45=((2.0524,17.0104,2),(2.0524,17.0104,-2),(1.76,15.01,2))
v46=((1.76,15.01,2),(1.76,15.01,-2),(2.0524,17.0104,-2))
v47=((1.76,15,2),(1.76,15,-2),(1.42,13.01,2))
v48=((1.42,13.01,2),(1.42,13.01,-2),(1.76,15,-2))
v49=((1.42,13,2),(1.42,13,-2),(1.16,11.49,2))
v50=((1.16,11.49,2),(1.16,11.49,-2),(1.42,13,-2))
v51=((1.16,11.48,2),(1.16,11.48,-2),(0.8699,10.0876,2))
v52=((0.8699,10.0876,2),(0.8699,10.0876,-2),(1.16,11.48,-2))
v53=((1.1059,9.9916,2),(1.1059,9.9916,-2),(2.5,9.9927,2))
v54=((2.5,9.9927,2),(2.5,9.9927,-2),(1.1059,9.9916,-2))
v55=((2.51,9.9927,2),(2.51,9.9927,-2),(4.0026,9.9927,2))
v56=((2.51,9.9927,-2),(4.0026,9.9927,-2),(4.0026,9.9927,2))
v57=((4,9.94,2),(4,9.94,-2),(3.87,9,2))
v58=((3.87,9,2),(3.87,9,-2),(4,9.94,-2))
v59=((3.87,8.99,2),(3.87,8.99,-2),(3.69,8,2))
v60=((3.69,8,2),(3.69,8,-2),(3.87,8.99,-2))
v61=((3.69,7.99,2),(3.69,7.99,-2),(3.51,7,2))
v62=((3.51,7,2),(3.51,7,-2),(3.69,7.99,-2))
v63=((3.51,6.99,2),(3.51,6.99,-2),(3.34,5.99,2))
v64=((3.34,5.99,2),(3.34,5.99,-2),(3.51,6.99,-2))
v65=((3.34,5.98,2),(3.34,5.98,-2),(3.16,4.98,2))
v66=((3.16,4.98,2),(3.16,4.98,-2),(3.34,5.98,-2))
v67=((3.16,4.97,2),(3.16,4.97,-2),(2.95,4.0001,2))
v68=((2.95,4.0001,2),(2.95,4.0001,-2),(3.16,4.97,-2))
v69=((0,4,2),(0,4,-2),(2.9,3.9,2))
v70=((0,4,-2),(2.9,3.9,2),(2.9,3.9,-2))

O.bodies.append(utils.facet(v1,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v2,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v3,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v4,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v5,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v6,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v7,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v8,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v9,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v10,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v11,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v12,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v13,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v14,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v15,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v16,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v17,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v18,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v19,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v20,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v21,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v22,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v23,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v24,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v25,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v26,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v27,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v28,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v29,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v30,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v31,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v32,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v33,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v34,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v35,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v36,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v37,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v38,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v39,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v40,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v41,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v42,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v43,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v44,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v45,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v46,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v47,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v48,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v49,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v50,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v51,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v52,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v53,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v54,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v55,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v56,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v57,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v58,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v59,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v60,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v61,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v62,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v63,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v64,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v65,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v66,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v67,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v68,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v69,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v70,color=(1,0,1), material=Bound))

for b in O.bodies:
    if b.material is Dolomite:
        b.state.blockedDOFs='zXY'

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Facet_Aabb()],verletDist=0.05),
   InteractionLoop(
      [Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),Ig2_Facet_Polyhedra_PolyhedraGeom()],
      [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()]
   ),
   NewtonIntegrator(gravity=(0,-9.81,0),damping=0.2),
]

utils.calm()

O.dt=10e-5
O.saveTmp()
O.run()

>> what values have you tried?
>O.dt=10e-8

what about 1e-10? 1e-12? 1e-15? 1e-20?

> If stiffness reaches less than 10e5, then particles are passed through facet.

what are the striffness values of both polyhedra and facets?
You can decrease stiffness of polyhedrons and increase stiffness of facets artificially (e.g. 1e300, the harmonic average in contact law makes it ok)

> Can you please provide a script (two bodies as is described) reproducing this behavior?

thanks, but it is not really a good script:
- there is a lot of not needed facets
- the two polyhedrons may not interact at all. E.g. with
      Dolomite.young = 5e7
      Bound.young = 1e300
      O.dt=10e-5
   the simulation behaves quite OK (because the polyhedrons does not interact)

for some other values, the interaction polyhedrons interact and "bounces excessively".

> it is a typical symptom of too large time step (making the simulation unstable) and I don't believe it is not overcome with decreasing the time step to a **sufficient** value

this still holds. It seems O.dt=1e-8 is still too large.

cheers
Jan

I really appreciate your patience Jan.

> Can you please provide a script (two bodies as is described) reproducing this behavior?
> but it is not really a good script

I corrected the codes in question and two polyhedrons have interact (With using your suggestion for young modules (Dolomite.young = 5e7 Bound.young = 1e300) and O.dt=10e-6). But when there are more polygons, problem appears.

> You can decrease stiffness of polyhedrons and increase stiffness of facets artificially (e.g. 1e300, the harmonic average in contact law makes it ok)

You know, this works for two polyhedrons! When I checked with more polyhedrons, the problem was still there.

> what about 1e-10? 1e-12? 1e-15? 1e-20?
In fact, I have no idea because if problem is really solved with these time steps, simulation will be impossible for me because my original code have 920 particles and it takes a long time.

Yes, DEM is expensive.. there are several option you can try:
- improve current polyhedrons implementation
- try different shapes, like PFacet, PotentialParticle or PotentialBlock
- try totally different method (some event-driven particle simulation
- try different software

In any of the cases, I have no idea how much you can gain..

cheers
Jan

> - try different shapes, like PFacet, PotentialParticle or PotentialBlock

Thanks a lot Jan , You right.

I want to try clump, hope to be successful.

Again Thanks for all your helps.

I would like to emphasize that I really have no idea how much (if any) performance you can gain.
On the other hand, to get some basic feeling, it should not be too difficult and time consuming to try e.g. clumped PFacets (although I have no personal experience with them).
cheers
Jan

> I would like to emphasize that I really have no idea how much (if any) performance you can gain.
> On the other hand, to get some basic feeling, it should not be too difficult and time consuming to try e.g. clumped PFacets

Yes Jan, I understand what you mean. I'll try it.

Thanks for your advice.

Best regards,
Mahdeyeh

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Dear Jan Stransky
I have a question about:

> - try totally different method (some event-driven particle simulation

is there event-driven particle simulation in the YADE right now?

no..
cheers
Jan

Thanks Jan Stransky.

I appreciate for introducing if you know any software with the ability of event-driven particle simulation.

Best regards,
Mahdeyeh

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Hi,
sorry for late answer, but no, I don't know any..
Jan