Picking the right BroadPhase

[tstucky]

Here's the list of possible BroadPhases we have to chose from

[BroadPhase<DEM::DiscreteElement, DEM::DiscreteElement>] : Between2SetOf
[BroadPhase<DEM::DiscreteElement, DEM::DiscreteElement>] : Hash
[BroadPhase<DEM::DiscreteElement, DEM::DiscreteElement>] : LCM
[BroadPhase<DEM::DiscreteElement, DEM::DiscreteElement>] : Raw
[BroadPhase<DEM::DiscreteElement, DEM::SupportShape>] : SupportShape
[BroadPhase<DEM::DiscreteElement, DEM::Tool>] : Tool
[BroadPhase<Physic::Body, Physic::Body>] : AABBTree
[BroadPhase<Physic::Body, Physic::Body>] : Hash
[BroadPhase<Physic::Body, Physic::Body>] : Lcm
[BroadPhase<Physic::Body, Physic::Body>] : Raw
[BroadPhase<Physic::Body, Physic::Ground>] : Raw
[BroadPhase<Physic::Particle, Physic::Ground>] : Raw
[BroadPhase<SPH::Particle, SPH::Particle>] : Lcm

With many of them missing documentation it is not at all clear to me which makes the most sense depending on a given simulation scenario.

I'll describe my scenario, and if somebody can please assist me in picking the proper BroadPhases for each collisional pair that would be a tremendous help. Since I am operating on a work laptop, I am aiming for speed over simulation accuracy.

My simulation needs to estimate the forces acting on a scoop as it digs through unconsolidated granular material that sits inside of a rectangular trench. The scoop geometry is described by a .off file and is imported into the sim as a DEM::Tool. The trench is built with several different Physic::Ground using the Rectangle shape. And the trench is filled with ~3000 DEM::DiscreteElements. See the attached image for reference.



Things that should collide:
  - DEM::DiscreteElement with DEM::Tool
  - DEM::DiscreteElement with Physic::Ground
  - DEM::DiscreteElement with DEM::DiscreteElement
Things that should not collide:
  - DEM::Tool with Physic::Ground - Since the DEM::Tool takes a linear scoop path that intersects Physic::Ground

A major point I am hung up on is that there appears to be no optimized BroadPhase for Physic::Ground. All that is available is Raw. The SupportShape BroadPhase looks to be far more optimized, but it does not appear I can build open shapes with SupportShape like my open trench.

[Cédric Hubert]

You are right, support shapes cannot be open. They are used to build discrete domains and to measure things on it, such as their length, radius, etc. during a simulation.

The Raw and SupportShape broad phases work similarly for DiscreteElement to SupportShape and DiscreteElement to Ground contact. The SupportShape broad phase simply handles inside and outside cases.

For these specific cases, I'm afraid that the broad phases have not really been improved. I will add it to the TODO list ;-)

For DiscreteElement to DiscreteElement contact, you have more choices:
 - Raw is just prohibitive (it checks possible contacts between all discrete elements !)
 - LCM "draws" a regular grid around the domain, which leads to very fast contact detection when the domain is compact
 - if not compact, you can use the AABBTree broad phase, which avoid building a huge LCM with plenty of voids in it
 - Hash is another fast broad phase detection method for non compact domains, but I don't know it well

As a side note, keep in mind that a DiscreteElement is a Body with a spherical shape. Thus, thanks to polymorphism, from the broad phase point of view, a DiscreteElement is a Body.

Hope this will help !

[tstucky]

It seems using the Raw BroadPhase really doesn't work for this scenario either. An anomalous force occurs that quickly pushes DE's to one side. Here is what I'm seeing after only 1000 iterations with a 1e-7s timesteps. Only force is from the Particle/Ground BroadPhase.

Iteration 0


Iteration 0 with X force shown


Iteration 0 with Z force shown


(Y force remained zero for all DE's)

Iteration 1000 with displacement norm shown


Here is what the trench geometry looks like without discrete elements. It is a rectangular concave, formed out of Physic::Ground Rectangles.

[Damien André]

hello, to fix this issue you can add the following PlugIn just after the ReadDomain plugIn.

<PlugIn Id="_Delete" What="DEM::SupportShape" ID="Boundary"/>

Because, the supportshape is also a ground.

[tstucky]

In the simulation results from my last post the support shape had been delete. I did notice I was forgetting to add mass to the discrete elements, but giving them mass hasn't change the behavior, only slowed it down a bit.

[tstucky]

Maybe you can answer for me if there is anything that looks amiss in my inp file? Attached is also the gdd (and lgdd in case your version is different enough from mine) file of the starting domain, but I cannot upload the scoop.off, so you will have to use a stand-in if you wish to run it yourself.

ScoopScrape.inp

CompactDomain.lgdd

CompactDomain.gdd

[Damien André]

Hello, sorry I forgot to make you a reply. You made a little mistake in your inp file. Instead of :

  <!-- build trench -z face -->
    <PlugIn Id="_NewFrame" Center="(0,0,-0.075)" Quat="(0,1,0)(90)" ID="F5"/>

Please use :

  <!-- build trench -z face -->
    <PlugIn Id="_NewFrame" Center="(0,0,-0.075)" Quat="(0,1,0)(-90)" ID="F5"/>

Hope this help, Damien.

[tstucky]

It is not clear to me why, but that has fixed it. Thank you!

Is there a proper and improper direction for Physic::Ground Rectangles to face?

[tstucky]

I still have difficulties with this sim. Now the particles are not flying to one direction, they are instead all blasting violently out of the open top of the rectangular container. I suspect this is due to the domain generation finishing with discrete elements still penetrating the walls of its their support shape, so when the domain is placed inside a series of Physic::Ground rectangles that match exactly the dimensions of its SupportShape, the collisions cause it to explode.

I have attempted several different fixes:
 1. Relaxing their energy with either dampening, zeroing out their velocities every 1,000 iterations, or a combination of both.
       * This likely did not lower the average penetration enough, because when it was turned off particles still achieved far too high kinetic energies. Perhaps it was not allowing the particles to achieve enough displacement to get sufficiently far from the Physic::Ground walls.
 2. Increasing the WallStiffDecreaseCoeff in Relax up to 0.99.
       * Not much change in wall penetration.
 3. Using the same gravity force, Stiffness, and Restitution values in the domain generation input as I use in the simulation input.
       * Doing this made the domain cooker script go on forever.
 4. A modification of #3, only the Relax plugin was disabled with the Process attribute, which allowed the simulation to complete.
       * This method produced the least wall penetration of any of the previous methods; however, when loading it into the simulation that replaces the SupportShape with a series of Physic::Ground objects, the elements will still quickly gain kinetic energy and explode out of the top.

Is there a proper way to do this I am just missing? I am new to discrete element method, and it is possible I am missing something that is obvious to others.

[Damien André]

Hello,
If you skip the relax step.... explosions are expected !
Because you built a domain where discrete elements are highly interpenetrated. I will make a patch today....

Let's keep updated.

[Damien André]

Hello, I have just implemented a patch.

Please update your granoo version with "svn update" and recompile it.
Now, add this option to your CookingBook.inp

 
<PlugIn Id="BuildBoundarySet" Update="False"/>

It should untrigger the unexpected rotation of the bounding shape and (I hope) solve your problem.
Cheers, Damien.

[tstucky]

You're right that turning off the Relax algorithm is a bad solution; however, as I express in item #3 of my last message, leaving Relax on and with my desired SupportShape stiffness set, the values shown by DisplayInfo (DE number, Cardinality, Compacity, and Average Penetration) would converge, but the simulation would not end, and therefore the BuildBoundarySet step was never reached and a valid domain (that I can add mass to) was never achieved.

The SupportShape stiffness used is only 1e5 N/m, which doesn't seem too high to me.

EDIT: Just to clarify, my difficulty seems to be with element-wall penetration. The way the CookingBook.inp file is initially configured, you wind up with a domain that has a lot of penetration of the SupportShape. I am wondering how I can lower that penetration, so that when I replace the support shape with a series of Physic::Grounds of higher stiffness, that starting wall penetration doesn't cause the whole domain to explode.

[Damien André]

Hello, it's strange, can you share with me a simplified version of your simulation that reproduces your bug ?

[tstucky]

Would you like a simplified version of the domain generation inp, the simulation inp, or both?

[tstucky]

Here is a simplified version of the input file that reproduces the convergent domain properties, but doesn't end the iteration for at least 50,000 iterations (when I manually terminated it). I tested this on revision 1911 of GranOO (I haven't installed your update yet).

domain_generation_stiff.inp

Thanks so much for your help!

[Damien André]

Hello, please replace the line :

    <PlugIn Id="_ManageCollision" Between="DiscreteElement/SupportShape"
      UniqueId="Wall" RecordContactIn="WallShell"
      BroadPhase="CustomSupportShape" SupportShapeID="Boundary" Mode="Inside"
      CallBack="Standard" Stiffness="1e5" RestitutionCoeff="0.1" />
    
    <PlugIn Id="_ManageCollision" Between="DiscreteElement/DiscreteElement"
      UniqueId="LCM" RecordInterpenetrationValues="Yes"
      BroadPhase="CustomLcm" 
      CallBack="Standard" Stiffness="1e5" RestitutionCoeff="0.1" />

by these ones :

    <PlugIn Id="_ManageCollision" Between="DiscreteElement/SupportShape" UniqueId="Wall"
  	  CallBack="Custom"
  	  BroadPhase="CustomSupportShape" RecordContactIn="WallShell"/>
    
    <PlugIn Id="_ManageCollision" Between="DiscreteElement/DiscreteElement" UniqueId="LCM"
    	    CallBack="Custom"
    	    BroadPhase="CustomLcm" RecordInterpenetrationValues="Yes"/>

Here we use the "Custom" callback because it is triggered to the "Relax" plugin.

Hope this help, Damien.

[tstucky]

I have already tried to produce domains using exactly the those _ManageCollision calls, and the domain it produces explodes when placed inside of a container because the DE-wall penetration is too high.

Here is another set of input files to illustrate this. Please let me know if you are able to reproduce the effect.

First: the domain generation, which follows exactly the method you showed.
domain_generation.inp

Second: the simulation, which starts to explode after about 10,000 iterations.
scoop_sim.inp

[Damien André]

Hello, please find below a working simulation :

<?xml version="1.0" encoding="UTF-8" standalone="no"?>

<GranOO Version="2.0" OutDir="ScoopScrape"
  Verbose="No" DisplayProgression="True"
  ThreadNumber="6" 
  TimeStep="1e-5"
  TotIteration="50000" >  
<!-- 
TotTime="1e-2" "TotTime="11.5"  TimeStep="1e-3"
  TotIteration="1"
-->

  <PreProcessing>
    <PlugIn Id="_ReadDomain"
      FileName="TrenchDomain/CompactDomain.gdd"
      Set="Elements" />

    <PlugIn Id="_SetContinuousDensity" Density="916.7" SupportShape="Boundary"/>
    <!-- Delete support shape. No longer needed now that we have mass. -->
    <PlugIn Id="_Delete" What="DEM::SupportShape" />

    <!-- build support surfaces -->
      <!-- build trench -y face -->
    <PlugIn Id="_NewFrame" Center="(0,-0.03,0)" Quat="(0,0,1)(90)" ID="F1"/>
    <PlugIn Id="_NewGround"
      Type="Rectangle" DimY="0.6" DimZ=".15" FrameID="F1" ID="G1" Set="Trench"/>

      <!-- build trench -x face -->
    <PlugIn Id="_NewFrame" Center="(-0.304,0,0)" Quat="(0,0,1)(0)" ID="F3"/>
    <PlugIn Id="_NewGround"
      Type="Rectangle" DimY=".05" DimZ=".15" FrameID="F3" ID="G3" Set="Trench"/>

      <!-- build trench +x face -->
    <PlugIn Id="_NewFrame" Center="(0.304,0,0)" Quat="(0,0,1)(180)" ID="F4"/>
    <PlugIn Id="_NewGround"
      Type="Rectangle" DimY=".05" DimZ=".15" FrameID="F4" ID="G4" Set="Trench"/>

      <!-- build trench -z face -->
    <PlugIn Id="_NewFrame" Center="(0,0,-0.08)" Quat="(0,1,0)(-90)" ID="F5"/>
    <PlugIn Id="_NewGround"
      Type="Rectangle" DimY="0.05" DimZ="0.6" FrameID="F5" ID="G5" Set="Trench"/>
    
      <!-- build trench +z face -->
    <PlugIn Id="_NewFrame" Center="(0,0,0.08)" Quat="(0,1,0)(90)" ID="F6"/>
    <PlugIn Id="_NewGround"
      Type="Rectangle" DimY="0.05" DimZ="0.6" FrameID="F6" ID="G6" Set="Trench"/>

      <!-- build trench flooring -x side -->
    <PlugIn Id="_NewFrame" Center="(-0.4,0.025,0)" Quat="(0,0,1)(90)" ID="F7"/>
    <PlugIn Id="_NewGround"
      Type="Rectangle" DimY="0.2" DimZ="0.55" FrameID="F7" ID="G7" Set="Trench"/>
    
      <!-- build trench flooring +x side -->
    <PlugIn Id="_NewFrame" Center="(0.4,0.025,0)" Quat="(0,0,1)(90)" ID="F8"/>
    <PlugIn Id="_NewGround"
      Type="Rectangle" DimY="0.2" DimZ="0.55" FrameID="F8" ID="G8" Set="Trench"/>

      <!-- build trench flooring -z side -->
    <PlugIn Id="_NewFrame" Center="(0,0.025,-0.175)" Quat="(0,0,1)(90)" ID="F9"/>
    <PlugIn Id="_NewGround"
      Type="Rectangle" DimY="0.6" DimZ="0.2" FrameID="F9" ID="G9" Set="Trench"/>

      <!-- build trench flooring +z side -->
    <PlugIn Id="_NewFrame" Center="(0,0.025,0.175)" Quat="(0,0,1)(90)" ID="F10"/>
    <PlugIn Id="_NewGround"
      Type="Rectangle" DimY="0.6" DimZ="0.2" FrameID="F10" ID="G10" Set="Trench"/>

  </PreProcessing>
  
  <Processing>
  	<PlugIn Id="_Check" />
  	<PlugIn Id="_ClearLoad" />

    <!-- between DE and ground -->
    <PlugIn Id="_ManageCollision"
      Between="Body/Ground" Set1="Elements" Set2="Trench"
      BroadPhase="Raw"
      NarrowPhase="WithShape"
      CallBack="Standard" Stiffness="1e5" RestitutionCoeff="0.1"
      RegularizationType="piecewise"
      DryFrictionSlope="1e5"
      StaticDryFrictionCoeff="0.5" DynamicDryFrictionCoeff="0.8"/>

    <!-- between DE and DE -->
    <PlugIn Id="_ManageCollision" ExcludeInteraction="False"
      Between="DiscreteElement/DiscreteElement" 
      BroadPhase="LCM" Set="Elements"
      CallBack="Standard" Stiffness="1e5" RestitutionCoeff="0.1"
      RegularizationType="piecewise"
      DryFrictionSlope="1e5"
      StaticDryFrictionCoeff="0.5" DynamicDryFrictionCoeff="0.8" />
    
    <!-- Europan surface gravity -->
    <PlugIn Id="_ApplyGravity" X="0." Y="-1.315" Z="0." />

    <PlugIn Id="_IntegrateAcceleration" Linear="Yes" Angular="Yes" />
    
    <!-- Save frame -->
    <PlugIn Id="_SaveDomain" IterLoop="1000" />
  </Processing>

  <PostProcessing>
  </PostProcessing>
</GranOO>

I give some margins to your bounding shapes because at the beginning of the simulation, the discrete elements highly interpenetrate the planes.

In my opinion, you can apply two strategies to overcome this difficulty :
 - increase the number of discrete elements (but it is bad for calculation time)
 - applying an intermediate step before your simulation : you let your domain relaxing with bounding shapes placed with margins and when your domain is stable you record it. After that, you start your real simulation with this domain instead of the one that comes directly from the cooking algorithm.

Hope this help, Damien.

[tstucky]

Thanks so much for all the help Damien!

Eventually increasing the number of DE's is desired, but in this early debugging phase, I am trying to keep the number of DE's and computation time low.

Since the second method could be automated I will try to implement that in the future.

Is there any reason we cannot use stiffness and the Relax plugin in the domain generation input file? If not, I would like to look into the algorithm myself and see if this can be done. Relaxing the domain with stiff walls would be a desirable simplification to the overall workflow.

[Damien André]

You are welcome :)

tstucky wrote

Is there any reason we cannot use stiffness and the Relax plugin in the domain generation input file? If not, I would like to look into the algorithm myself and see if this can be done. Relaxing the domain with stiff walls would be a desirable simplification to the overall workflow.

The relaxing algorithm of the cooker decreases gradually the wall stiffness. While the wall stiffness tend to zero, interpenetrations between discrete elements and walls are expected especially for domains with a small number of discrete elements. So this strategy does not fit with your problem.

Another approach could be decreasing the discrete element radii while keeping the wall stiffnesses to high value for preventing interpenetration.

For your problem, I recommend to use intermediate simulation with the following workflow :

Cooker (with bounding dimensions a little lower than the targetted one) ---> Relaxing sim (with the right bounding dimension) --> Final sim

If you prefer to implement a new relaxing strategy in the cooking algorithm, it should be nice if you can share with us your developments !
Kind regards, Damien.