RungeKuttaCashKarp54Integrator with gridConnections is causing errors

Asked by Rohit John on 2021-01-28

Hello,

I want to use the RungeKuttaCashKarp54Integrator engine in a simulation with cylinderConnections. I added viscous damping to the source code of Law2_ScGeom6D_CohFrictPhys_CohesionMoment, this caused the NewtonIntegrator to be unstable (The objects are exploding). I concluded this may be because of the integration scheme and opted to use RungeKuttaCashKarp54Integrator. However when I try to run the simulation with RungeKuttaCashKarp54Integrator, I get "RuntimeError: No IGeomDispatcher in engines or inside InteractionLoop".

The minimum working code is given below. Please help me to fix the issue

Kind regards,
Rohit John

##############################################################################
# encoding: utf-8
"An example showing bending beams."

from yade.gridpfacet import *

#### Parameter ####
L=10. # length of the beam
n=12 # number of nodes used to generate the beam
r=L/50. # radius of the beam element

#### Engines ####
integrator=RungeKuttaCashKarp54Integrator([
 ForceResetter(),
 GeneralIntegratorInsertionSortCollider([Bo1_GridConnection_Aabb()]),
 InteractionLoop(
  [Ig2_GridNode_GridNode_GridNodeGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False)],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 ),
 GravityEngine(gravity=Vector3(0,0,-9.81)),
])

#Tolerances can be set for the optimum accuracy
integrator.rel_err=1e-6;
integrator.abs_err=1e-6;

O.engines=[integrator,
]

#### Create materials and set different properties ####
O.materials.append(CohFrictMat(young=1e6,poisson=0.3,density=1e1,frictionAngle=radians(10),normalCohesion=1e40,shearCohesion=1e40,momentRotationLaw=True,label='mat1'))

#### Vertices ####
vertices1=[]

for i in range(0,n):
  vertices1.append( [i*L/n,0,0] )

#### Create cylinder connections ####
nodesIds=[]
cylIds=[]
cylinderConnection(vertices1,r,nodesIds,cylIds,color=[1,0,0],highlight=False,intMaterial='mat1')

#### Set boundary conditions ####
for i in range(0,1):
   O.bodies[nodesIds[i*n]].dynamic=False

#### Set a time step ####
O.dt=1e-05

#### Allows to reload the simulation ####
O.saveTmp()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Janek Kozicki
Solved:
2021-01-28
Last query:
2021-01-28
Last reply:
2021-01-28
Best Janek Kozicki (cosurgi) said : #1

See: https://gitlab.com/yade-dev/trunk/-/issues/171 it is not fixed. You can try the workaround mentioned in that issue, use flag dead=True.

> I added viscous damping to the source code of Law2_ScGeom6D_CohFrictPhys_CohesionMoment, this caused the NewtonIntegrator to be unstable

Quick note:
We can't really try since we don't have your specific damping model. You could push to a gitlab branch at some point if you want someone to be able to try the script with relevant code version.

Maybe you did that already, but before trying another integrator (or higher number precision ;) ), I would go through this checkList:
- what happens with super-small O.dt? (defining what is "super-small" needs to think about mass-viscosity-stiffness parameters)
- are the directions of forces, moment, velocities and spin consistent, with just 2-3 nodes after 2 iterations?

Cheers
Bruno

Rohit John (rohitkjohn) said : #3

Thanks Janek Kozicki, that solved my question.

Rohit John (rohitkjohn) said : #4

Dear Bruno,

> what happens with super-small O.dt? (defining what is "super-small" needs to think about mass-viscosity-stiffness parameters)
- It depends on the damping fraction and the material stiffness. I wrote the code such that the damping coefficient is a fraction of the stiffness. (Same fraction for normal, shear, twist and bending stiffness in Law2_ScGeom6D_CohFrictPhys_CohesionMoment). At small time steps it works, but if I want higher damping I must reduce the time step. For one simulation it had to be 1e-10. Even days were not enough to get some useful result.

> are the directions of forces, moment, velocities and spin consistent, with just 2-3 nodes after 2 iterations?
I checked the values and when the damping and dt were badly set, the forces and velocities grow exponentially, but when I reduce the dt or damping they "calm down".

From my understanding the integrator in YADE is similar to forward Euler which is unstable when there is damping. So I though changing the integrator might help. Janek Kozicki solution works but the RungeKuttaCashKarp54Integrator is extremely slow.

If possible could you check could the code found here: https://github.com/rohitkjohnedu/yade_laws .

Kind regards,
Rohit K. John

You are past-the-end of my checklist. :)
RungeKuttaCashKarp54Integrator makes sense then.

Last one maybe: in a working case (small damping), did you find consistency _at the structural scale_ in terms of damping the oscillations? It could be that there is a mismatch between what you input at the interactions scale and what you expect in terms of damping at the structural scale?

> is similar to forward Euler

Yes and no. It is a 2nd order time centered scheme from Newton's point of view. Problem is velocities are not known on step, they are know at mid-steps. So this wonderful scheme for mass-spring systems becomes a poor 1st order Euler when there are viscous terms. If you take mid-steps as your frame of reference and you think to what happens to velocities wrt. visous forces it is plain explicit Euler, you are right.

There are ways to evaluate on-step velocites (something similar to "velocity Verlet" [1]) which we don't use. It should improve your case, it needs some changes in the integrator and it needs two calculations of forces per time iteration.

> the RungeKuttaCashKarp54Integrator is extremely slow.

Because there are even more than two calculations of forces per time iteration, I guess?

Do you mean that despite larger timesteps it is still slower?

Bruno

[1] https://en.wikipedia.org/wiki/Verlet_integration

> It could be that there is a mismatch between what you input at the interactions scale and what you expect in terms of damping at the structural scale

I would say the bending viscosity of the interactions should be a multiple of ~N/L, do we agree on that?

> https://github.com/rohitkjohnedu/yade_laws .

Problem here is we can't see you changes.
From the top of my head, you need to do something like:
0. get sources with "git clone"
1. edit sources
2. "git commit -a -m 'that's my change'"
3. git push to gitlab/github whatever

This way we can see what was done specifically in step 1

Bruno

Rohit John (rohitkjohn) said : #7

Dear Bruno,

> Problem here is we can't see you changes.
- I see. I tried pushing the code into the YADE gitlab and I got this error
remote: You are not allowed to push code to this project.

So I pushed the source into my account. You can view the changes here.
https://github.com/rohitkjohnedu/my_yade/commit/908aef935ee1236cd5188dd3e20544cfeba00117

Kind regards,
Rohit K. John

Rohit John (rohitkjohn) said : #8

Dear Bruno,

>Last one maybe: in a working case (small damping), did you find consistency _at the structural scale_ in terms of damping the oscillations? It could be that there is a mismatch between what you input at the interactions scale and what you expect in terms of damping at the structural scale?

- I cannot say this with certainity because I did not find any literature with data that I can use to verify me results. But in the structural case I did find that the motion is damped. For example, in one simulation I had attached at weight to the tip of a cantilever (a forceEngine at the particle at the end) and let it reach equilibrium after deformation. Then the weight was released (forceEngine dead = True), the oscillations were similar to what we find in a damped oscillator.
The mismatch that you mentions is a bit unclear. Could you please elaborate or maybe give an example?

> I would say the bending viscosity of the interactions should be a multiple of ~N/L, do we agree on that?

- Could you please define N and L. Am I wrong to assume N is the number of gridNodes and L is the length of the cantilever? My first thought is that damping should be a material property so it should be constant. However, it could also depend on the shape of the cantilever, like the area moment of inertia perhaps. It is difficult to find good literature on damping in beams. Or I might looking in the wrong place or the wrong thing.

Kind regards,
Rohit K. John

Hi,

> the oscillations were similar to what we find in a damped oscillator

Cool. :)

> Am I wrong to assume N is the number of gridNodes and L is the length of the cantilever?

You are right.

Sorry, I wrote that comment right after reading another thread on a linked topic, where N and L were clearly defined. Overall, the thing is: interaction properties - if they are supposed to reflect a beam-like thing - necessarily depend on beam geometry and discretisation. It can't be conveyed by material properties alone. The "mismatch" I was mentioning was about that: material properties are not "per-unit length of a beam" properties and so someone can be off by orders of magnitudes if based on wrong expectations.

I couldn't check the github push in depth, yet, but thanks for clarity.
Let you report here if you have more insights on the issue.

Bruno

Rohit John (rohitkjohn) said : #10

Thanks for your advice Bruno, I shall put an update in this thread when I find something useful.

Kind regards,
Rohit K. John