Clump angular momentum not conserved

Hello all,

I have a simulation consisting of two spheres (S1, S2) and a clump made of two spheres (C1, C2). The two spheres, S1 and S2, are given the same initial speed but in the opposite directions along the horizontal. The clump is aligned along the vertical and the two spheres, S1 and S2, are made to hit the spheres in the clump at the same time. Could this be a bug in my code or the NewtonIntegrator?

I expect the angular momentum to be conserved in the interactions, but it is not. I have a similar problem with clumps of pfacet [1]. It was reported in a bug report related to pfacet [2]. I am posting again because this appears to be an error related to clumps in general. I have put the code below.

I am using Yadedaily.

Kind regards,
Rohit K. John

[1] https://answers.launchpad.net/yade/+question/697465
[2] https://gitlab.com/yade-dev/trunk/-/issues/210

# script
# ----------------------------------------------------------------------------
from yade.gridpfacet import *
from yade import *
from yade import geom, utils
from yade import plot
import sys, os

# ---------------------------------------------------------------------------- input parameter
# ----------------------------------------------------- clump
sphere_young = 50e9
sphere_poisson = 0.3
sphere_friction = radians(30)

sphere_offset = 0.3
sphere_radius = 5e-2
sphere_init_vel = 1

sphere_mass = 0.5
sphere_density = sphere_mass / (4/3*pi*sphere_radius**3)
# ---------------------------------------------------------------------------- Materials
sphere_mat = 'sphere_mat'

O.materials.append(
    FrictMat(
        young = sphere_young,
        poisson = sphere_poisson,
        density = sphere_density,
        label = sphere_mat,
        frictionAngle = sphere_friction,
    )
)

# ---------------------------------------------------------------------------- Engines
O.engines = [
                ForceResetter(),

                InsertionSortCollider([
                    Bo1_Sphere_Aabb(),
                ]),

                InteractionLoop(
                    [Ig2_Sphere_Sphere_ScGeom()],
                    [Ip2_FrictMat_FrictMat_FrictPhys()],
                    [Law2_ScGeom_FrictPhys_CundallStrack()],
                )
            ]
# ---------------------------------------------------------------------------- objects
# ----------------------------------------------------- clump
clump_sph_1 = sphere([0, sphere_offset, 0], sphere_radius, material = sphere_mat)
clump_sph_2 = sphere([0,-sphere_offset, 0], sphere_radius, material = sphere_mat)

clump_sph_1_ID = O.bodies.append(clump_sph_1)
clump_sph_2_ID = O.bodies.append(clump_sph_2)

clump_clump_ID = O.bodies.clump([clump_sph_1_ID, clump_sph_2_ID])
clump_mmoi = O.bodies[clump_clump_ID].state.inertia
# ----------------------------------------------------- spheres
sp1 = sphere([ 4*sphere_radius, sphere_offset, 0], sphere_radius, material = sphere_mat)
sp2 = sphere([-4*sphere_radius,-sphere_offset, 0], sphere_radius, material = sphere_mat)

sp1_ID = O.bodies.append(sp1)
sp2_ID = O.bodies.append(sp2)

O.bodies[sp1_ID].state.vel = [-sphere_init_vel,0,0]
O.bodies[sp2_ID].state.vel = [ sphere_init_vel,0,0]

# ---------------------------------------------------------------------------- Additional engines
O.engines += [
    NewtonIntegrator(
        gravity = [0,0,0],
        damping = 0
    ),
    PyRunner(
        command = 'graph()',
        iterPeriod = 1000
    ),
]

# ----------------------------------------------------- Calculating state after collision
m = sphere_mass
v0 = sphere_init_vel
l = sphere_offset
iz = clump_mmoi[2]

calculated_final_sphere_lin_vel = 0.5*(-4.0*iz*l*m*v0/(iz + 2.0*l**2*m) + 2.0*l*m*v0)/(l*m)
calculated_final_clump_ang_vel = 4.0*l*m*v0/(iz + 2.0*l**2*m)

calculated_final_clump_ang_mom = clump_mmoi[2]*calculated_final_clump_ang_vel
calculated_final_sphere_ang_mom = 2*sphere_offset*sphere_mass*calculated_final_sphere_lin_vel

print("clump: ", calculated_final_clump_ang_mom)
print("sphere: ", calculated_final_sphere_ang_mom)
print("total: ", calculated_final_clump_ang_mom + calculated_final_sphere_ang_mom)

# ----------------------------------------------------- plotting
plot.plots = {'t':('clump_Lz', 'calculated_clump_final_Lz'), 't1':('Lz', 'calculated_Lz')}

def graph():
    L_sp1 = get_AngMom(sp1_ID)
    L_sp2 = get_AngMom(sp2_ID)

    angMom_clump = O.bodies[clump_clump_ID].state.angMom
    L_tot = angMom_clump + L_sp1 + L_sp2

    plot.addData(
        t = O.time, t1 = O.time,
        Lx = L_tot[0], Ly = L_tot[1], Lz = L_tot[2],
        clump_Lz = angMom_clump[2],
        calculated = calculated_final_clump_ang_mom + calculated_final_sphere_ang_mom,
        calculated_clump_final_Lz = calculated_final_clump_ang_mom
    )

def get_AngMom(id):
    pos_sp1 = O.bodies[id].state.pos
    vel_sp1 = O.bodies[id].state.vel
    angMom_sp1 = O.bodies[id].state.angMom
    mass_sp1 = O.bodies[id].state.mass
    return angMom_sp1 + mass_sp1 * pos_sp1.cross(vel_sp1)

# ----------------------------------------------------------------------------- simulation controls
O.dt = 2e-7#utils.PWaveTimeStep()
O.saveTmp()
plot.plot()