Aborted (core dumped) occured when I use the Hertz-Minlin contact law with non-zero krot and eta

Asked by Xin Li

python3.8: /builds/yade-dev/trunk/deb/yadedaily/pkg/common/InsertionSortCollider.cpp:778: bool yade::InsertionSortCollider::spatialOverlapPeri(yade::Body::id_t, yade::Body::id_t, yade::Scene*, yade::Vector3i&) const: Assertion `maxima[3 * id1 + axis] - minima[3 * id1 + axis] < .99 * dim' failed.
python3.8: /builds/yade-dev/trunk/deb/yadedaily/pkg/common/InsertionSortCollider.cpp:778: bool yade::InsertionSortCollider::spatialOverlapPeri(yade::Body::id_t, yade::Body::id_t, yade::Scene*, yade::Vector3i&) const: Assertion `maxima[3 * id1 + axis] - minima[3 * id1 + axis] < .99 * dim' failed.
Aborted (core dumped)
I am running a periodic triaxial test using hertz-mindlin contact law. However, the program abort suddenly as soon as I add the rolling stiffness coefficient and rolling plastic limit (as shown above). If these values are set to 0, the simulation can continue.
My code is as follows:
# encoding: utf-8

readParamsFromTable(
num=2500,
)

import glob, os
import numpy as np
from yade.params import table
from yade import pack, plot

isBatch=runningInBatch() #check if run in batch mode
if isBatch:
    print('Running: '+O.tags['description']) #mark the current simulation

num = table.num # number of soil particles
dScaling = 1e3 # density scaling to increase crticial time step
sigmac= -100000 # confining pressure
rate = 0.1 # strain rate
damp = 0.5 # damping coefficient
stabilityRatio = 1.e-3 # threshold for quasi-static condition
stressTolRatio = 1.e-3 # tolerance for stress goal

# corners to define specimen size
m_1=0.004*pow(num/2500,1/3)
m_2=0.008*pow(num/2500,1/3)
mn,mx=Vector3.Zero,Vector3(m_1,m_1,m_2)

# Soil sphere parameters
E=7e9 # micro Young's modulus
v=0.3 # micro Poisson's ratio
kr=0.5 # rolling/bending stiffness coefficient
eta=0.5 # rolling/bending plastic limit
mu = 26 # contact friction angle
rho = 2650*dScaling # soil density

# create materials
spMat = O.materials.append(
   FrictMat(young=E,poisson=v,frictionAngle=radians(mu),density=rho))

# create a cloud of ramdomly positioned spheres
O.periodic = True
O.cell.hSize = Matrix3(mx[0],0,0, 0,mx[1],0, 0,0,mx[2])
sp=pack.SpherePack()
sizes=[.00015,.000225,.0003,.0003625,.00045,.0006,0.00118]
cumm=[0,.215,.3807312826,.601,.8865351132,.9965177017,1]
sp.makeCloud(mn,mx,psdSizes=sizes,psdCumm=cumm,distributeMass=True,seed=1,num=num)

# load spheres to simulation
spIds=sp.toSimulation(material=spMat)

#define engines
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb()]),
   InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_MindlinPhys(
            krot=kr,
            ktwist=kr,
            eta=eta
        )],
   [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True)]
   ),
   PeriTriaxController(label='triax',
      goal=(sigmac,sigmac,sigmac),
      # whether they are strains or stresses
      stressMask=7,
      # type of servo-control
      dynCell=True,
      maxStrainRate=(10*rate,10*rate,10*rate),
      # wait until the unbalanced force goes below this value
      maxUnbalanced=stabilityRatio,
      relStressTol=stressTolRatio,
      doneHook='compactionFinished()'
   ),
   NewtonIntegrator(damping=damp,label='newton'),
   PyRunner(iterPeriod=500,command='addPlotData()',label='data')
]

O.dt=.3*PWaveTimeStep()

def compactionFinished():
    if unbalancedForce()<stabilityRatio:
        print(porosity())
        triax.doneHook='startShearing()'

def startShearing():
    # set the current cell configuration to be the reference one
    O.cell.trsf = Matrix3.Identity
    O.cell.velGrad = Matrix3.Zero
    setRefSe3()
    # set loading type: constant pressure in x,y, 30% compression in z
    triax.stressMask=3
    triax.globUpdate=1
    # allow faster deformation along x,y to better maintain stresses
    triax.maxStrainRate=(10*rate,10*rate,rate)
    triax.goal=(sigmac,sigmac,-.3)
    triax.maxUnbalanced=stabilityRatio
    triax.relStressTol=stressTolRatio
    triax.doneHook='Finished()'

def Finished():
    print(str(num)+'finished')
    O.pause()

def addPlotData():
    if triax.stressMask != 7:
        plot.addData(i=O.iter,t=O.realtime,s11=-triax.stress[0],\
                    s22=-triax.stress[1],s33=-triax.stress[2],e11=triax.strain[0],\
                    e22=triax.strain[1],e33=triax.strain[2],n=porosity(),\
                    V=O.cell.volume,M=triax.mass,Unb=unbalancedForece())
        if isBatch:
            dataName = 'triax2_'+ O.tags['description'] + '.txt'
        else:
            dataName = 'TriaxREV2_'+str(num)+'.txt'
        plot.saveDataTxt(dataName,vars=('s11','s22','s33','e11','e22','e33','n','V','M'))

O.run()
waitIfBatch()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Xin Li
Solved:
Last query:
Last reply:
Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#1

Hi,

please note that krot is not a problem. The simulation fails if ktwist is nonzero (but you prescribe 'kr' to both 'krot' and 'ktwist').
It works without setting 'ktwist' (e.g. comment out #ktwist=kr,)

Cheers,
Karol

PS Remember that version and system information is helpful too, and try to include it in your next questions.

Revision history for this message
Xin Li (portal3) said :
#2

Thank you very much Karol!

After canceling the setting of torsional stiffness of hertz-mindlin contact law, this simulation started to run. Besides, I am still trying to figure out why the value of ktwist can cause the error above. Do I need to change the Ubuntu system? My Ubuntu system is Ubuntu 21.10 and Yade 2021.01a using python 3.9.7. And, there is the same problem when the code is running in another computer (Ubuntu 20.04 LTS and yadedaily using python 3.8.10.

Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#3

Sorry, but I don't know how to help with 'ktwist'. I don't have any experience with twisting yet. Maybe you should start another topic.
My comment about system information was general. If you include it in your question, it may help someone to solve your problem.

Cheers,
Karol

Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#4

Hi,
It could be that a high torsional stiffness needs a smaller timestep, so that you get a classical numerical instability if you keep it the same.
GlobalStiffnessTimeStepper would take this into account, but PWaveTimeStep() doesn't.
Cheers
Bruno

Revision history for this message
Anton Gladky (gladky-anton) said :
#5

It would also be good if you compile yade with debug symbols, so we can get more information about the crash.

Could you also please create an issue on gitlab [1] with this information about this problem?

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

Thanks

Anton

Revision history for this message
Xin Li (portal3) said :
#6

Hi Bruno,

Thank you for your answer! I have changed the time step setting to GlobalStiffnessTimeStepper(), this program can continue to run. But, the time step suddenly drops to 0 after a few seconds. If I set the ktwist of Ip2_FrictMat_FrictMat_MindlinPhys() to 0, there is no problem. I am trying some other examples with Ip2_FrictMat_FrictMat_MindlinPhys(ktwist= non-zero value) to figure out if my system or yade version has bugs.

Revision history for this message
Xin Li (portal3) said :
#7

Hi Anton,

I have create an issue with the same information about this problem on gitlab.

Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#8

> the time step suddenly drops to 0

Zero exactly?

Revision history for this message
Xin Li (portal3) said :
#9

Hi Bruno,

> Zero exactly?

Yeah, actually from the 3D view of the simulation, no strain can be observed before timestep droped to 0 even though no errors are reported. But, I found that this was caused by small stable timestep like you said torsional stiffness needs a smaller timestep to maintian numerical stability. After reducing the timestepSafetyCoefficient, the simulation can run normally.
Really thank you for your answer! !

Revision history for this message
Klaus Thoeni (klaus.thoeni) said :
#10

I believe this was related to an issue in the GlobalStiffnessTimeStepper which was fixed here [1]. Please can you confirm that everything runs as expected with yadedaily so we can close the issue?

Thanks,
Klaus

[1] https://gitlab.com/yade-dev/trunk/-/merge_requests/618

Revision history for this message
Xin Li (portal3) said :
#11

Hi Klaus,

I reran the simulation with hertz-mindlin contact law (with rolling stiffness) and GlobalStiffnessTimeStepper(safetycoefficient=0.8), everything was ok.
Thank you!

Regards,
Xin