problem of erasing particle

Hi Guo,

You are removing particles while the simulation is still in process. You need to make the script wait, units the simulation pauses:

O.run(20, wait=True)

Cheers,
Karol

Dear Karol,

Thanks for your reply. How can I achieve it in the engine?

This code is not right.
def deletParticle():
    for bid in ParID:
        x=O.bodies[bid].state.pos[0]
        if (0.2<x<0.3):
            O.bodies.erase(bid)
            ParID.remove(bid)
            print(f'Particle {bid} is deleted.')

sp=pack.SpherePack()
sp.makeCloud([0,0,0],[0.5,0.5,0.5],num=10000)
ParID=[]
for c,r in sp:
    ParID.append(O.bodies.append(sphere(c,r)))

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(),
    NewtonIntegrator(),
    PyRunner(command='O.wait()', firstIterRun=200),
    PyRunner(command='deletParticle()', firstIterRun=200),
    PyRunner(command='O.run()', firstIterRun=200),
]

O.run()

Do not use O.wait() and O.run() inside PyRunner command.
While PyRunner command is executed, the simulation is "paused" (contrary to the code in OP).

###
O.engines = [
    ...
    NewtonIntegrator(),
    PyRunner(command='deletParticle()', firstIterRun=200),
]
###

Cheers
Jan

Dear Jan,

Thanks for your reply.
###
O.engines = [
    ...
    NewtonIntegrator(),
    PyRunner(command='deletParticle()', firstIterRun=200),
]
###
This method doesn't work. The simulation will stop .

Best Regards,
GUO

> This method doesn't work. The simulation will stop .

please be more specific [1]

Just in case, a complete MWE below. When I run it, I have no problem.
After 200 steps, I got bunch of "Particle XXXX is deleted." messages is printed in terminal.
In 3D view I see most of particles from middle plane deleted.

###
def deletParticle():
    for bid in ParID:
        x=O.bodies[bid].state.pos[0]
        if (0.2<x<0.3):
            O.bodies.erase(bid)
            ParID.remove(bid)
            print(f'Particle {bid} is deleted.')

sp=pack.SpherePack()
sp.makeCloud([0,0,0],[0.5,0.5,0.5],num=10000)
ParID=[]
for c,r in sp:
    ParID.append(O.bodies.append(sphere(c,r)))

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(),
    NewtonIntegrator(),
    PyRunner(command='deletParticle()', firstIterRun=200),
]

O.run()
###

Cheers
Jan

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

Yes. You are right. Sorry for forgetting to check this code.
Actually, this is a simplified code. Now I know the reason. If I use the periodic boundary, there is an error.

WidthX=0.7 # box width
HeightZ=0.7 # box height, not soil sample's
LengthY=0.25 # pipe length
################## set periodic cell ##################
# O.periodic = True
# O.cell.hSize=Matrix3(WidthX,0,0,0,LengthY,0,0,0,HeightZ)
def deletParticle():
    for bid in ParID:
        x=O.bodies[bid].state.pos[0]
        if (0.2<x<0.3):
            O.bodies.erase(bid)
            ParID.remove(bid)
            print(f'Particle {bid} is deleted.')

sp=pack.SpherePack()
sp.makeCloud([0,0,0],[0.5,0.5,0.5],num=10000)
ParID=[]
for c,r in sp:
    ParID.append(O.bodies.append(sphere(c,r)))

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(),
    NewtonIntegrator(),

    PyRunner(command='deletParticle()', firstIterRun=200),

]

O.run()

> there is an error

please be more specific ([1], point 2). I.e.:
- provide error message
- provide OS and YADE versions
- ...

I got no error, running your script with both commented and uncommented lines:
O.periodic = True
O.cell.hSize=Matrix3(WidthX,0,0,0,LengthY,0,0,0,HeightZ)

Cheers
Jan

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

Dear Jan,

Thank you.

There is no error message. The simulation cannot run anymore.
I find the problem is unstable. Sometimes this code is ok, but sometimes not.

My OS is Ubuntu 20.04. Yade version is 20230205-7089~9ce9332~focal.

Best regards,
GUO, Chang

Do you have the same problems without erasing particles?
If yes, please open a new question for a new problem [1]
Cheers
Jan

Dear Jan,

No. This problem only happens when erasing particles.

For this code below, sometimes there is no bug. But sometimes, it would stop at O.iter=200-203 and I can only force close the terminal.

WidthX=0.7 # box width
HeightZ=0.7 # box height, not soil sample's
LengthY=0.25 # pipe length
################## set periodic cell ##################
# O.periodic = True
# O.cell.hSize=Matrix3(WidthX,0,0,0,LengthY,0,0,0,HeightZ)
def deletParticle():
    for bid in ParID:
        x=O.bodies[bid].state.pos[0]
        if (0.2<x<0.3):
            O.bodies.erase(bid)
            ParID.remove(bid)
            print(f'Particle {bid} is deleted.')

sp=pack.SpherePack()
sp.makeCloud([0,0,0],[0.5,0.5,0.5],num=10000)
ParID=[]
for c,r in sp:
    ParID.append(O.bodies.append(sphere(c,r)))

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(),
    NewtonIntegrator(),

    PyRunner(command='deletParticle()', firstIterRun=200),

]

O.run()

As a suggestion, are you sure you need this ParID list ? Maybe you could do without, for instance with
for b in O.bodies:
  x = b.state.pos
  if ...
     O.bodies.erase(b.id)

which would at least make a simpler script (avoiding to have to wonder how does it work to iterate over an evolving list)

I am using 20230305-7138~9cadf57~jammy1 , Ubuntu 22.04
I have run your code like 10 times, without any problem..
Without being able to reproduce the problem, I cannot help..
Cheers
Jan