Unable to locate NewtonIntegrator within O.engines.

Unable to locate MWE.py :-)

[1]https://www.yade-dem.org/wiki/Howtoask

As the error says, you are using a >0 value for InsertionSortCollider.verletDist. Whereas the source code [*] requires (for an efficient approximate collision detection ?) NewtonIntegrator to belong to O.engines, which you decided against.

[*] see this if-block at https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/common/InsertionSortCollider.cpp#L393 --- sorry I do not warranty this version is exactly the one of your YADE 2018.02b version

Thanks rcaulk and jduriez,

Sorry for getting back late.

I met this problem by:
###
Firstly run script1 to generate sample and save it:
O.save('sample.yade.gz')
Then using script2 to reload the sample: O.load('sample.yade.gz'),
Then, carrying out triaxial test in script2.
In script2, when I define the O.engines as follow, it shows error:
###Following engines shows error###
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
 PyRunner(command='stopIfDamaged()',iterPeriod=1000),
 NewtonIntegrator(damping=0.4),
]

When I change the engines into the following one, it works well. The difference between the error one and good one is the position of NewtonIntegrator(damping=0.4).

##Following O.engines works well###
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
 NewtonIntegrator(damping=0.4),
    PyRunner(command='stopIfDamaged()',iterPeriod=1000),
]

I don't know why but it works.

Cheers,

Leonard

> I don't know why but it works.

for future reference, could you please post a complete script causing this strange problem?
Thanks
Jan

Hi Jan,
Thanks for your suggestion.
Here are the two scripts based on Yade 2018.02b:
#########
The following is the script1 which is for generating sample (takes around 30 seconds):
# -*- coding: utf-8 -*-
#*************************************************************************
from __future__ import division
from yade import pack, plot
import math
import numpy as np
import timeit

########## this script is modified from https://gitlab.com/yade-dev/trunk/blob/master/examples/triax-tutorial/script-session1.py#L142

utils.readParamsFromTable(lowerR=10.0,upperR=10.0)
from yade.params import table

num_spheres=2000# number of spheres
targetPorosity = 0.5 #the porosity we want for the packing
compFricDegree = 30 # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=-0.01 # loading rate (strain rate)
damp=0.4 # damping coefficient
stabilityThreshold=0.001 # we test unbalancedForce against this value in different loops (see below)
young=5e7 # contact stiffness

confinement=100e3

mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07)

## create materials for spheres and plates
MatWall=O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
MatSand = O.materials.append(CohFrictMat(isCohesive=False,young=young,poisson=0.5,frictionAngle=radians(30),
           density=2650.0,normalCohesion=1e6, shearCohesion=1e6,label='sand'))
## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()

sp.makeCloud(mn,mx,-1,0.3333,num_spheres,True, 0.95,seed=1)

O.bodies.append([sphere(center,rad,material='sand') for center,rad in sp])

Gl1_Sphere.quality=3

triax=TriaxialStressController(
 maxMultiplier=1.+6e5/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e4/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask = 7,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
 newton,
]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0

#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=triax.goal2=triax.goal3=-confinement

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1]
  if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001:
 break

print "### Isotropic completed ###"

import sys
while triax.porosity>targetPorosity:
 # we decrease friction value and apply it to all the bodies and contacts
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()

 O.run(500,1)

print "### state 2 Reach target porosity completed ###"

print 'top', -triax.stress(triax.wall_top_id)[1], 'right',-triax.stress(triax.wall_right_id)[0],'front',-triax.stress(triax.wall_front_id)[2]

triax.internalCompaction=False
triax.stressMask = 7
triax.goal1=triax.goal2=triax.goal3=-confinement
triax.max_vel=0.001

O.save('sample.yade.gz')
#####################################################
The following is script2 which is for triaxial test, it illustrates the problem.
##
# -*- coding: utf-8 -*-
#*************************************************************************
from yade import pack, plot
import math
import numpy as np

O.load('sample.yade.gz')
Gl1_Sphere.quality=3
confinement=100e3

triax=TriaxialStressController(
 ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
 ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
 thickness = 0,
 stressMask = 5,
 internalCompaction=False, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=0.4)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
    PyRunner(command='stopIfDamaged()',iterPeriod=1000),
 newton,
]

Gl1_Sphere.stripes=0
yade.qt.Controller()
yade.qt.View()
setContactFriction(radians(30))
triax.stressMask = 5
triax.goal2=-0.015
triax.goal1=-confinement
triax.goal3=-confinement

from yade import plot

### a function saving variables
def history():
 plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
   Volumetric_strian=triax.strain[0]+triax.strain[1]+triax.strain[2],
   s11=-triax.stress(triax.wall_right_id)[0],
   s22=-triax.stress(triax.wall_top_id)[1],
   s33=-triax.stress(triax.wall_front_id)[2],
   dev=-triax.stress(triax.wall_top_id)[1]-confinement,
   i=O.iter)

def stopIfDamaged():
 if -triax.strain[1]>0.15:
  O.pause()
  plot.saveDataTxt('data')

if 1:
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=1000,command='history()',label='recorder')]+[PyRunner(iterPeriod=100,command='stopIfDamaged()',label='checkDamage')]+O.engines[5:7]

plot.plots={'e22':('s11','s22','s33'),'e22 ':('dev'),'e22 ':('Volumetric_strian')}
plot.labels={'e22':'$\epsilon_{a}$','s11':'$\sigma_{11}$','s22':'$\sigma_{22}$','s33':'$\sigma_{33}$'}
plot.plot()
########
Run script2 after script1, I have:
In [1]: FATAL /build/yade-fDuCoe/yade-2018.02b/core/ThreadRunner.cpp:30 run: Exception occured:
InsertionSortCollider.verletDist>0, but unable to locate NewtonIntegrator within O.engines.

You may switch the position of newton with PyRunner in the O.engines (in script2) to see the difference. For me, the error is solved by this way.

Cheers,

Leonard

Thanks for the complete script.
Next time please post it in the OP, we could save all the thread :-)

> if 1:
> O.engines=O.engines[0:5]+...+O.engines[5:7]

just print O.engines before and after ;-)

cheers
Jan

One more comment. Between O.engines=[...] and this error-prone O.engines=O.engines[0:5]+...+O.engines[5:7] you have no relevant code (like running), you can put everything directly in the first O.engines=[...]
Jan

Thanks Jan for your comment.

Cheers,

Leonard