Crash when using CPM model

Hello Colin,

as the error says, sigmaT[0] is NaN.

I usually get these errors when I have for some reason time step larger
than the critical one and the simulation goes to instable regime and soon
some value goes to NaN. Could this be your case? As you mentioned
"scaling", did you also correctly modify the time step?

If this is not the case, we will find another reason :-)

cheers
Jan

2014-08-05 14:37 GMT+02:00 Colin Power <<email address hidden>
>:

> New question #252656 on Yade:
> https://answers.launchpad.net/yade/+question/252656
>
> Hi, I have a simulation that uses randomDensePack() and the CPM model to
> pack about 5000 spheres. On my initial testing, I appended about 250
> spheres using a script and was able to run the simulation to completion.
> When I scaled up the simulation I got this error after about 370000
> iterations:
>
> FATAL
> /tmp/buildd/yadedaily-1.10.0-72-d9ab58c~precise/pkg/dem/ConcretePM.cpp:399
> go: Verification `!isnan(sigmaT[0])' failed!
> FATAL
> /tmp/buildd/yadedaily-1.10.0-72-d9ab58c~precise/pkg/dem/ConcretePM.cpp:399
> go: in interaction #25639+#36271
> terminate called without an active exception
> Aborted (core dumped)
>
> I took a quick glance at the file linked below but was not able to figure
> out my problem.
>
> https://github.com/yade/trunk/blob/1997c194c0aa759cae101a3dd0a559fcf049b29f/pkg/dem/ConcretePM.cpp
>
> Thank You,
> Colin Power
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

I can try running my simulation with a smaller time step and see what happens.

The successful simulation was ran with time step 5.67819191123e-05, while the failed one used 2.6296979111e-06. The sphere radii in the first simulation are very closely distributed around 0.18. The failed simulation uses radius of 0.05 and a fuzz of 0.6666. I used 0.5*PWaveTimeStep() to achieve these time steps.

0.5*PWaveTimeStep() should be safe enough.. Do you have also something else
in your simulation (e.g. some facets or walls), or only spheres?
Jan

2014-08-05 15:51 GMT+02:00 Colin Power <<email address hidden>
>:

> Question #252656 on Yade changed:
> https://answers.launchpad.net/yade/+question/252656
>
> Colin Power posted a new comment:
> I can try running my simulation with a smaller time step and see what
> happens.
>
> The successful simulation was ran with time step 5.67819191123e-05,
> while the failed one used 2.6296979111e-06. The sphere radii in the
> first simulation are very closely distributed around 0.18. The failed
> simulation uses radius of 0.05 and a fuzz of 0.6666. I used
> 0.5*PWaveTimeStep() to achieve these time steps.
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Yes, there are multiple facets in the simulation and there are also are multiple material types. I am also using a buoyancy model that is outlined in this example: https://github.com/yade/trunk/blob/f50c914fd20b550e06512a774022653467b2e041/examples/clumps/apply-buoyancy-clumps.py

I have pasted my script below. The ccore_yade module just contains helper functions for file management.

import sys, time, random, os, gts, math
from yade import ymport
from yade import pack
from yade import qt
import ccore_yade

# Simulation Parameters
nbIter=6000000
filename = 'seabed'
volumeName = 'keel'
iterPeriod = 100

folder = ccore_yade.newSimulationFolder(filename)
print "Simulation is being stored in: " + folder
filepath = folder+'/'+filename

#used for randomDensePack, experimenting with this.
memoizeDb='keel_packings.sqlite'

#define material properties:
shearModulus = 3.6e9
poissonRatio = 0.28
youngModulus = 2*shearModulus*(1+poissonRatio)
angle = radians(36) #friction angle in radians
angle_r = radians(36) #friction angle in radians
rho_p = 917 #density of particles
rho_f = 1029 #density of the fluid, rho_f > rho_p = floating particles
rho_r = 1602 #density of seafloor rock
rho_s = 8050 #density of steel
v_iw = 0 #velocity of increasing water-level
gravity = -9.80665 #m/s
damping = 0.3 #ratio
radius = 0.05 #m
seabed_velocity = 0.05 #m/s

#define water boundaries:
boundaryMin = Vector3(-6, -6, -0.14)
boundaryMax = Vector3(10, 10, 10)
waterLevel = boundaryMin[2]
t0 = O.time

#define colors:
sphereColor = (0.0,0.0,0.05) #dirty yellow
waterColor = (.2,.2,.7) #blue

#define materials:
id_Mat=O.materials.append(CpmMat(young=youngModulus, poisson=poissonRatio, density=rho_p, frictionAngle=angle, sigmaT=5.5e3, relDuctility=30, epsCrackOnset=1e10, isoPrestress=0, damLaw=0) )
Mat=O.materials[id_Mat]

id_Mat=O.materials.append(FrictMat(young=50e6, poisson=0.3, density=rho_r, frictionAngle=angle_r))
rock=O.materials[id_Mat]

#steel materials have been seperated for filtering in Paraview
steel_Mat1=O.materials.append(FrictMat(young=50e6, poisson=0.3, density=rho_s, frictionAngle=angle))
steel1=O.materials[steel_Mat1]

steel_Mat2=O.materials.append(FrictMat(young=50e6, poisson=0.3, density=rho_s, frictionAngle=angle))
steel2=O.materials[steel_Mat2]

steel_Mat3=O.materials.append(FrictMat(young=50e6, poisson=0.3, density=rho_s, frictionAngle=angle))
steel3=O.materials[steel_Mat3]

# Import mesh files to create facets
seabed = O.bodies.append(ymport.gmsh(meshfile=filename+'.mesh', shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity, color=(0.75,0.35,0), material=rock, wire=False))
leftWall = O.bodies.append(ymport.gmsh(meshfile='leftWall.mesh', shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity, color=(0.8,0.8,0.85), material=steel1, wire=False))
rightWall = O.bodies.append(ymport.gmsh(meshfile='rightWall.mesh', shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity, color=(0.8,0.8,0.85), material=steel2, wire=False))
topWall = O.bodies.append(ymport.gmsh(meshfile='topWall.mesh', shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity, color=(0.8,0.8,0.85), material=steel3, wire=False))

# Create water level indicator
idsWaterFacets = []
idsWaterFacets.append(O.bodies.append(facet( \
   [ boundaryMin, [boundaryMax[0], boundaryMin[1], boundaryMin[2]], [boundaryMax[0], boundaryMax[1], boundaryMin[2]] ], \
   fixed=True, material=FrictMat(young=0), color=waterColor, wire=False))) #no interactions will appear
idsWaterFacets.append(O.bodies.append(facet( \
   [ [boundaryMax[0], boundaryMax[1], boundaryMin[2]], [boundaryMin[0], boundaryMax[1], boundaryMin[2]], boundaryMin ], \
   fixed=True, material=FrictMat(young=0), color=waterColor, wire=False))) #no interactions will appear

surface=gts.read(open(volumeName+'.gts'))
volume=pack.inGtsSurface(surface)

O.bodies.append(pack.randomDensePack(predicate=volume, radius=radius, material=Mat, cropLayers=0, rRelFuzz=0.6666, spheresInCell=50000, memoizeDb=None, useOBB=False, memoDbg=False, color=None, returnSpherePack=False))

#Define timestep
O.dt=0.5*PWaveTimeStep()

ccore_yade.recordTimeStep(filepath, iterPeriod)

O.engines=[
 ForceResetter(),

 InsertionSortCollider([ Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
  [Ip2_Mat_CpmMat_CpmPhys(cohesiveThresholdIter=10000), Ip2_FrictMat_CpmMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGeom_CpmPhys_Cpm(omegaThreshold=1, yieldEllipseShift=0, yieldSurfType=0)]
 ),

 NewtonIntegrator(damping=damping, gravity=[0, 0, gravity], label='integrator'),
    TranslationEngine(dead=False, ids=seabed, translationAxis=[-1,0,0], velocity=seabed_velocity, label='seabed_movement'),

    #Moveable walls for future simulations
    #TranslationEngine(dead=False, ids=leftWall, translationAxis=[0,1,0], velocity=0.1, label='leftWall_movement'),
    #TranslationEngine(dead=False, ids=rightWall, translationAxis=[0,-1,0], velocity=0.1, label='rightWall_movement'),
    #TranslationEngine(dead=False, ids=topWall, translationAxis=[0,0,-1], velocity=0.1, label='topWall_movement'),

    # Save data for Paraview
    VTKRecorder(fileName=filepath+'_vtk', recorders=['all', 'cpm'], iterPeriod=iterPeriod),

    PyRunner(command='O.pause()', iterPeriod=nbIter),
    PyRunner(iterPeriod=25, command='applyBuoyancy()', label='buoLabel'),
    PyRunner(iterPeriod=10000, command='print O.iter')
]

print "Start simulation: " + filename
qt.View() # launch Yade with the 3D view
O.pause() # launch Yade with the simulation paused

# Based on https://github.com/yade/trunk/blob/f50c914fd20b550e06512a774022653467b2e041/examples/clumps/apply-buoyancy-clumps.py
def applyBuoyancy():
 global waterLevel
 waterLevel = (O.time-t0)*v_iw + boundaryMin[2]
 for b in O.bodies:
  zMin = 1e9
  zMax = -1e9
  if b.isStandalone and isinstance(b.shape,Sphere):
   rad = b.shape.radius
   zMin = b.state.pos[2] - rad
   dh = min((waterLevel - zMin), 2*rad) #to get sure, that dh is not bigger than 2*radius
  elif b.isClump: #determine rad, zMin and zMax for clumps:
   for ii in b.shape.members.keys():
    pos = O.bodies[ii].state.pos
    zMin = min(zMin,pos[2]-O.bodies[ii].shape.radius)
    zMax = max(zMax,pos[2]+O.bodies[ii].shape.radius)
   #get equivalent radius from clump mass:
   rad = ( 3*b.state.mass/(4*pi*O.bodies[b.shape.members.keys()[0]].mat.density) )**(1./3.)
   #get dh relative to equivalent sphere, but acting when waterLevel is between clumps z-dimensions zMin and zMax:
   dh = min((waterLevel - zMin)*2*rad/(zMax - zMin), 2*rad)
  else:
   continue
  if dh > 0:
   F_buo = -1*(pi/3)*dh*dh*(3*rad - dh)*rho_f*integrator.gravity # = -V*rho*g
   O.forces.addF(b.id, F_buo, permanent=True)

If some spheres have multiple contacts with facets, then maybe PWave timestep is irrelevant. Just an idea.
If you want to be sure it is not a timestep problem you can use automatic timestep with GlobalStiffnessTimeStepper.
Bruno

Hi Colin,
sorry for late reply. I tried your script, but to run it, I also
need ccore_yade module. Could you please send it in email or place it
somewhere to the internet and send a link?
thanks
Jan

2014-08-05 16:51 GMT+02:00 Colin Power <<email address hidden>
>:

> Question #252656 on Yade changed:
> https://answers.launchpad.net/yade/+question/252656
>
> Colin Power posted a new comment:
> Yes, there are multiple facets in the simulation and there are also are
> multiple material types. I am also using a buoyancy model that is
> outlined in this example:
>
> https://github.com/yade/trunk/blob/f50c914fd20b550e06512a774022653467b2e041/examples/clumps
> /apply-buoyancy-clumps.py
>
> I have pasted my script below. The ccore_yade module just contains
> helper functions for file management.
>
> import sys, time, random, os, gts, math
> from yade import ymport
> from yade import pack
> from yade import qt
> import ccore_yade
>
> # Simulation Parameters
> nbIter=6000000
> filename = 'seabed'
> volumeName = 'keel'
> iterPeriod = 100
>
> folder = ccore_yade.newSimulationFolder(filename)
> print "Simulation is being stored in: " + folder
> filepath = folder+'/'+filename
>
> #used for randomDensePack, experimenting with this.
> memoizeDb='keel_packings.sqlite'
>
> #define material properties:
> shearModulus = 3.6e9
> poissonRatio = 0.28
> youngModulus = 2*shearModulus*(1+poissonRatio)
> angle = radians(36) #friction angle in
> radians
> angle_r = radians(36) #friction angle in
> radians
> rho_p = 917 #density of particles
> rho_f = 1029 #density of the
> fluid, rho_f > rho_p = floating particles
> rho_r = 1602 #density of seafloor rock
> rho_s = 8050 #density of steel
> v_iw = 0 #velocity of
> increasing water-level
> gravity = -9.80665 #m/s
> damping = 0.3 #ratio
> radius = 0.05 #m
> seabed_velocity = 0.05 #m/s
>
> #define water boundaries:
> boundaryMin = Vector3(-6, -6, -0.14)
> boundaryMax = Vector3(10, 10, 10)
> waterLevel = boundaryMin[2]
> t0 = O.time
>
> #define colors:
> sphereColor = (0.0,0.0,0.05) #dirty yellow
> waterColor = (.2,.2,.7) #blue
>
> #define materials:
> id_Mat=O.materials.append(CpmMat(young=youngModulus, poisson=poissonRatio,
> density=rho_p, frictionAngle=angle, sigmaT=5.5e3, relDuctility=30,
> epsCrackOnset=1e10, isoPrestress=0, damLaw=0) )
> Mat=O.materials[id_Mat]
>
> id_Mat=O.materials.append(FrictMat(young=50e6, poisson=0.3, density=rho_r,
> frictionAngle=angle_r))
> rock=O.materials[id_Mat]
>
> #steel materials have been seperated for filtering in Paraview
> steel_Mat1=O.materials.append(FrictMat(young=50e6, poisson=0.3,
> density=rho_s, frictionAngle=angle))
> steel1=O.materials[steel_Mat1]
>
> steel_Mat2=O.materials.append(FrictMat(young=50e6, poisson=0.3,
> density=rho_s, frictionAngle=angle))
> steel2=O.materials[steel_Mat2]
>
> steel_Mat3=O.materials.append(FrictMat(young=50e6, poisson=0.3,
> density=rho_s, frictionAngle=angle))
> steel3=O.materials[steel_Mat3]
>
> # Import mesh files to create facets
> seabed = O.bodies.append(ymport.gmsh(meshfile=filename+'.mesh',
> shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity,
> color=(0.75,0.35,0), material=rock, wire=False))
> leftWall = O.bodies.append(ymport.gmsh(meshfile='leftWall.mesh',
> shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity,
> color=(0.8,0.8,0.85), material=steel1, wire=False))
> rightWall = O.bodies.append(ymport.gmsh(meshfile='rightWall.mesh',
> shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity,
> color=(0.8,0.8,0.85), material=steel2, wire=False))
> topWall = O.bodies.append(ymport.gmsh(meshfile='topWall.mesh',
> shift=Vector3(0, 0, 0), scale=1.00, orientation=Quaternion.Identity,
> color=(0.8,0.8,0.85), material=steel3, wire=False))
>
> # Create water level indicator
> idsWaterFacets = []
> idsWaterFacets.append(O.bodies.append(facet( \
> [ boundaryMin, [boundaryMax[0], boundaryMin[1],
> boundaryMin[2]], [boundaryMax[0], boundaryMax[1], boundaryMin[2]] ], \
> fixed=True, material=FrictMat(young=0),
> color=waterColor, wire=False))) #no interactions will appear
> idsWaterFacets.append(O.bodies.append(facet( \
> [ [boundaryMax[0], boundaryMax[1],
> boundaryMin[2]], [boundaryMin[0], boundaryMax[1], boundaryMin[2]],
> boundaryMin ], \
> fixed=True, material=FrictMat(young=0),
> color=waterColor, wire=False))) #no interactions will appear
>
> surface=gts.read(open(volumeName+'.gts'))
> volume=pack.inGtsSurface(surface)
>
> O.bodies.append(pack.randomDensePack(predicate=volume, radius=radius,
> material=Mat, cropLayers=0, rRelFuzz=0.6666, spheresInCell=50000,
> memoizeDb=None, useOBB=False, memoDbg=False, color=None,
> returnSpherePack=False))
>
> #Define timestep
> O.dt=0.5*PWaveTimeStep()
>
> ccore_yade.recordTimeStep(filepath, iterPeriod)
>
> O.engines=[
> ForceResetter(),
>
> InsertionSortCollider([ Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
> [Ip2_Mat_CpmMat_CpmPhys(cohesiveThresholdIter=10000),
> Ip2_FrictMat_CpmMat_FrictPhys()],
> [Law2_ScGeom_FrictPhys_CundallStrack(),
> Law2_ScGeom_CpmPhys_Cpm(omegaThreshold=1, yieldEllipseShift=0,
> yieldSurfType=0)]
> ),
>
> NewtonIntegrator(damping=damping, gravity=[0, 0, gravity],
> label='integrator'),
> TranslationEngine(dead=False, ids=seabed, translationAxis=[-1,0,0],
> velocity=seabed_velocity, label='seabed_movement'),
>
> #Moveable walls for future simulations
> #TranslationEngine(dead=False, ids=leftWall, translationAxis=[0,1,0],
> velocity=0.1, label='leftWall_movement'),
> #TranslationEngine(dead=False, ids=rightWall,
> translationAxis=[0,-1,0], velocity=0.1, label='rightWall_movement'),
> #TranslationEngine(dead=False, ids=topWall, translationAxis=[0,0,-1],
> velocity=0.1, label='topWall_movement'),
>
> # Save data for Paraview
> VTKRecorder(fileName=filepath+'_vtk', recorders=['all', 'cpm'],
> iterPeriod=iterPeriod),
>
> PyRunner(command='O.pause()', iterPeriod=nbIter),
> PyRunner(iterPeriod=25, command='applyBuoyancy()', label='buoLabel'),
> PyRunner(iterPeriod=10000, command='print O.iter')
> ]
>
> print "Start simulation: " + filename
> qt.View() # launch Yade with the 3D view
> O.pause() # launch Yade with the simulation paused
>
> # Based on
> https://github.com/yade/trunk/blob/f50c914fd20b550e06512a774022653467b2e041/examples/clumps/apply-buoyancy-clumps.py
> def applyBuoyancy():
> global waterLevel
> waterLevel = (O.time-t0)*v_iw + boundaryMin[2]
> for b in O.bodies:
> zMin = 1e9
> zMax = -1e9
> if b.isStandalone and isinstance(b.shape,Sphere):
> rad = b.shape.radius
> zMin = b.state.pos[2] - rad
> dh = min((waterLevel - zMin), 2*rad) #to get
> sure, that dh is not bigger than 2*radius
> elif b.isClump: #determine rad,
> zMin and zMax for clumps:
> for ii in b.shape.members.keys():
> pos = O.bodies[ii].state.pos
> zMin =
> min(zMin,pos[2]-O.bodies[ii].shape.radius)
> zMax =
> max(zMax,pos[2]+O.bodies[ii].shape.radius)
> #get equivalent radius from clump mass:
> rad = (
> 3*b.state.mass/(4*pi*O.bodies[b.shape.members.keys()[0]].mat.density)
> )**(1./3.)
> #get dh relative to equivalent sphere, but acting
> when waterLevel is between clumps z-dimensions zMin and zMax:
> dh = min((waterLevel - zMin)*2*rad/(zMax - zMin),
> 2*rad)
> else:
> continue
> if dh > 0:
> F_buo = -1*(pi/3)*dh*dh*(3*rad -
> dh)*rho_f*integrator.gravity
> # = -V*rho*g
> O.forces.addF(b.id, F_buo, permanent=True)
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi, I have an altered scrip that caused the same crash. This time it does not include facets, it took me about 2 hours for the crash to occur at around iteration 190000.

The script uses a .gts file for the packing so I have copied the two files to pastebin.

http://pastebin.com/gSrUxhzQ
http://pastebin.com/8Wquq33c

This time the error message was:

FATAL /tmp/buildd/yadedaily-1.11.0-4-90640cc~trusty/pkg/dem/ConcretePM.cpp:399 go: Verification `!isnan(sigmaT[0])' failed!
FATAL /tmp/buildd/yadedaily-1.11.0-4-90640cc~trusty/pkg/dem/ConcretePM.cpp:399 go: in interaction #1636+#3172
terminate called without an active exception
Aborted (core dumped)

Thank You,
Colin Power