Stirred column: a pitched-blade impeller that actually sweeps

A rotating analytic wall whose geometry actually sweeps, driving 22 000 grains through a Lacey mixing curve — with a no-free-energy guardrail that the first, frictionless version of this page failed by a factor of five.

dem

sdf

csg

moving-walls

granular

mixing

Author

Peclet

Published

August 30, 2026

Runs on a CPU in about ten minutes — 22 000 grains, no fluid.

What you’ll learn

The rotating drum got away with a trick: a barrel is a body of revolution, so you can spin it by giving a static wall a rigid-body velocity field and never moving a single triangle. A stirrer blade cannot be faked that way. Its geometry has to sweep, or there is nothing to push the grains.

This page uses peclet.dem’s two wall controls together, which is the whole point:

set_wall_transform(wall, translation, quat) places an analytic wall rigidly in the world. It moves the geometry — the CSG tree the narrow phase evaluates. Calls are absolute: the placement is composed onto the tree’s authored root transform, not onto the previous frame’s result, so a thousand revolutions accumulate no drift.
set_wall_velocity(wall, lin_vel, ang_vel, center) sets the surface velocity a grain in contact feels.

Neither is inferred from the other. That is deliberate — it is what lets a drum have a velocity without moving, and it is why this example must drive both from the same ang_vel every step. Getting one without the other is a specific, diagnosable bug: transform without velocity gives a blade that shoves grains but exerts no tangential drag; velocity without transform gives the drum trick applied to a shape it does not fit, and the blades stand still while the grains swirl.

You will also see a container and an impeller written as CSG trees rather than voxel grids, and a quantitative mixing measure — the Lacey index — defined explicitly rather than eyeballed.

import importlib.util, os, subprocess, sys
# Set the thread count BEFORE the module is imported: Kokkos reads it at initialisation.
os.environ.setdefault("OMP_NUM_THREADS", "8")
os.environ.setdefault("OMP_PROC_BIND", "false")
_local = os.environ.get("PECLET_LOCAL_BUILD")
if _local:
    for p in _local.split(os.pathsep):
        sys.path.insert(0, p)
elif importlib.util.find_spec("peclet") is None:
    subprocess.run([sys.executable, "-m", "pip", "install", "-q", "peclet", "scikit-image"],
                   check=True)

Requires a recent peclet

set_wall_transform and wall_sdf_at are new, and the CSG authoring layer (peclet.core.geom) is newer than the current PyPI release; on Colab you need a source build of peclet-core and peclet-dem until the next release goes out. The page is frozen, so it renders and reads correctly either way.

import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from skimage import measure

from peclet.core import geom
from peclet import dem as pdem

plt.rcParams.update({"figure.dpi": 130, "font.size": 9, "axes.axisbelow": True,
                     "figure.facecolor": "white", "savefig.bbox": "tight"})
print("threads:", os.environ["OMP_NUM_THREADS"])

threads: 8

Step 1 — The column and the impeller, as CSG trees

Everything is in grain-radius-ish units: the grain radius is RP, gravity points along $-y$, and the column axis is $y$.

The column is an inverted solid: capsule ∩ box. core’s capsule is a $y$-axis primitive with hemispherical caps; intersecting it with a box slices those caps off and leaves a flat-bottomed cylinder. Both leaves are distance-exact and the intersection of two convex bodies keeps that exactness on the inside, which is where the grains live. invert=True turns the solid into a container: the wall SDF becomes positive in the void.

The impeller is a union: a shaft capsule plus four flat blades, each a box rotated by 45° about its own radial axis (the pitch) and then swung to its azimuthal station. No inversion — the wall SDF is positive outside the metal, again where the grains are.

RP     = 0.4       # grain radius
R_COL  = 20.0      # column inner radius
H_COL  = 24.0      # column height (flat bottom at y = 0)
R_IMP  = 16.0      # blade tip radius
Y_IMP  = 3.5       # impeller centre height
PITCH  = np.deg2rad(45.0)
NBLADE = 4
G      = 9.81

def qmul(a, b):
    ax, ay, az, aw = a; bx, by, bz, bw = b
    return (aw*bx + ax*bw + ay*bz - az*by, aw*by - ax*bz + ay*bw + az*bx,
            aw*bz + ax*by - ay*bx + az*bw, aw*bw - ax*bx - ay*by - az*bz)

def qaxis(axis, angle):
    """Quaternion (x, y, z, w) for a rotation of `angle` about `axis`."""
    s = np.sin(0.5 * angle); n = np.asarray(axis, float); n /= np.linalg.norm(n)
    return (float(n[0]*s), float(n[1]*s), float(n[2]*s), float(np.cos(0.5*angle)))

def qrot(q, v):
    x, y, z, w = q; v = np.asarray(v, float)
    qv = np.array([x, y, z]); t = 2 * np.cross(qv, v)
    return v + w * t + np.cross(qv, t)

def column_tree():
    b = geom.SceneBuilder()
    side = b.add_leaf("capsule", [R_COL, H_COL], translation=[0.0, 0.5*H_COL, 0.0])
    slab = b.add_leaf("box", [R_COL*1.5, 0.5*H_COL, R_COL*1.5], translation=[0.0, 0.5*H_COL, 0.0])
    root = b.add_intersection(side, slab)          # flat-ended cylinder
    ni, nr, _, _ = b.encode()
    return np.asarray(ni, np.int32), np.asarray(nr, np.float32), root

def stirrer_tree():
    b = geom.SceneBuilder()
    r_hub, t_b, w_b = 1.2, 0.35, 1.8               # hub radius, blade half-thickness, half-width
    L_b = 0.5 * (R_IMP - r_hub)                    # blade radial half-length
    root = b.add_leaf("capsule", [1.0, 0.5*H_COL], translation=[0.0, 0.5*H_COL, 0.0])   # shaft
    for k in range(NBLADE):
        theta = 2*np.pi*k/NBLADE
        q = qmul(qaxis([0, 1, 0], theta), qaxis([1, 0, 0], PITCH))   # swing, then pitch
        c = qrot(qaxis([0, 1, 0], theta), [r_hub + L_b, 0.0, 0.0])
        blade = b.add_leaf("box", [L_b, t_b, w_b],
                           translation=[float(c[0]), Y_IMP + float(c[1]), float(c[2])],
                           rotation=list(q))
        root = b.add_union(root, blade)
    ni, nr, _, _ = b.encode()
    return np.asarray(ni, np.int32), np.asarray(nr, np.float32), root

wall_sdf_at returns exactly what the narrow phase reads, so it is the honest way to check a placement — and the easiest way to draw the impeller.

Code

def make_sim(capacity=40000):
    s = pdem.Simulation(capacity)
    s.set_domain((-R_COL-4, -4.0, -R_COL-4), (R_COL+4, H_COL+4, R_COL+4))
    s.enable_periodicity(False, False, False)
    s.set_gravity(0.0, -G, 0.0)
    s.set_sphere_shape(1.0)
    s.set_global_scale(RP)          # ghost band is sized off the real radius; without this the
                                    # default globalScale = 1 makes it 10x the grain
    cni, cnr, croot = column_tree()
    sni, snr, sroot = stirrer_tree()
    w_col = s.add_analytic_wall(cni, cnr, croot, True, 0.2, 0.5)    # invert -> container
    w_stir = s.add_analytic_wall(sni, snr, sroot, False, 0.2, 0.5)
    s.set_material_params(0.2, 0.0, 0.5)   # body-body restitution / tangential / friction
    s.set_solver_iterations(4, 4)
    return s, w_col, w_stir

def stirrer_mesh(sim, wall, box=(-R_IMP-1, R_IMP+1), ylim=(0.0, 8.0), n=110):
    xs = np.linspace(box[0], box[1], n)
    ys = np.linspace(ylim[0], ylim[1], max(24, int(n*(ylim[1]-ylim[0])/(box[1]-box[0]))))
    X, Y, Z = np.meshgrid(xs, ys, xs, indexing="ij")
    pts = np.ascontiguousarray(np.column_stack([X.ravel(), Y.ravel(), Z.ravel()]).astype(np.float32))
    fld = np.asarray(sim.wall_sdf_at(wall, pts)).reshape(X.shape)
    v, f, _, _ = measure.marching_cubes(fld, level=0.0,
                                        spacing=(xs[1]-xs[0], ys[1]-ys[0], xs[1]-xs[0]))
    v = v + np.array([xs[0], ys[0], xs[0]])
    return v, f

sim, w_col, w_stir = make_sim()          # the ONE simulation this page uses
figg = plt.figure(figsize=(7.0, 3.2))
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
for panel, ang in enumerate((0.0, np.deg2rad(45.0))):
    sim.set_wall_transform(w_stir, (0.0, 0.0, 0.0), qaxis([0, 1, 0], ang))
    v, f = stirrer_mesh(sim, w_stir)
    ax = figg.add_subplot(1, 2, panel + 1, projection="3d")
    ax.add_collection3d(Poly3DCollection(v[f], facecolor="#8a8f98", edgecolor="none"))
    tt = np.linspace(0, 2*np.pi, 120)
    for yy in (0.0, 8.0):
        ax.plot(R_COL*np.cos(tt), np.full_like(tt, yy), R_COL*np.sin(tt), color="#4c72b0", lw=0.8)
    ax.set_xlim(-R_COL, R_COL); ax.set_zlim(-R_COL, R_COL); ax.set_ylim(0, 2*R_COL)
    ax.set_box_aspect((1, 1, 1)); ax.set_axis_off(); ax.view_init(elev=22, azim=-58)
    ax.set_title("authored pose" if panel == 0 else "after set_wall_transform(45°)", fontsize=8)
plt.show()
sim.set_wall_transform(w_stir, (0.0, 0.0, 0.0), qaxis([0, 1, 0], 0.0))   # back to the authored pose

Figure 1: The impeller, meshed from the wall SDF the narrow phase actually evaluates, inside the column outline. Left: the authored pose. Right: the same wall after `set_wall_transform` swings it 45°, which is what moves the metal.

Step 2 — Fill and settle

A jittered lattice inside the column, with anything starting inside the impeller discarded — checked with wall_sdf_at, not with a hand-written formula that could disagree with the collision geometry.

def fill(sim, w_stir, fill_top=11.0, seed=4):
    rng = np.random.default_rng(seed)
    d = 2.05 * RP
    xs = np.arange(-R_COL + RP*1.2, R_COL - RP*1.2, d)
    ys = np.arange(RP*1.2, fill_top, d)
    X, Y, Z = np.meshgrid(xs, ys, xs, indexing="ij")
    P = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])
    P = P[np.hypot(P[:, 0], P[:, 2]) < R_COL - 1.3*RP]
    P = np.ascontiguousarray((P + rng.normal(0, 0.06*RP, P.shape)).astype(np.float32))
    P = np.ascontiguousarray(P[np.asarray(sim.wall_sdf_at(w_stir, P)) > RP*1.15])
    sim.set_positions(P)
    sim.set_inv_mass(np.full(len(P), 1.0, np.float32))   # dem assigns unit mass regardless of size
    return P

DT = 2e-3
P0 = fill(sim, w_stir)
sim.set_dt(DT)
t0 = time.time()
for _ in range(700):
    sim.step(DT)                 # dem.step() with NO argument advances nothing
Q0 = np.asarray(sim.get_positions()).copy()
print("%d grains settled in %.0f s   bed y = %.2f .. %.2f   max overlap %.4f (grain radius %.2f)"
      % (len(Q0), time.time() - t0, Q0[:, 1].min(), Q0[:, 1].max(), sim.max_overlap(), RP))

22121 grains settled in 32 s   bed y = 0.39 .. 9.99   max overlap 0.0009 (grain radius 0.40)

Set the body–body material, or the bed leaks

dem’s default body–body material is frictionless with zero restitution. A frictionless deep bed behaves like a liquid: it puts full hydrostatic pressure on the container and the position solve cannot hold it, and grains push through the analytic wall. The first version of this page lost grains through the floor (centres 1.7 radii below it) and reached peak speeds four to six times the blade-tip speed. set_material_params(0.2, 0.0, 0.5) fixes it completely — see the guardrail in Step 5, and ISSUES.md.

Step 3 — The Lacey mixing index

Label the bottom half of the settled bed as tracer, then measure how uniformly it spreads. With $p$ the overall tracer fraction and $x_i$ the tracer fraction in bin $i$,

\[ M = \frac{s_0^2 - s^2}{s_0^2 - s_r^2},\qquad s_0^2 = p(1-p),\qquad s_r^2 = \frac{p(1-p)}{\bar n}, \tag{1}\]

where $s^2$ is the variance of $x_i$ across bins, $s_0^2$ is the fully-segregated variance, $s_r^2$ the fully-random (binomial) one, and $\bar n$ the mean grains per bin. $M = 0$ is segregated, $M = 1$ is random (Lacey 1954). The binning is fixed for the whole run — equal-area annuli in $r$ crossed with equal slices in $y$ — so nothing about the measure moves with the grains.

NR, NY = 5, 8
tracer = Q0[:, 1] < np.median(Q0[:, 1])
p_frac = tracer.mean()
r_edges = np.sqrt(np.linspace(0, R_COL**2, NR + 1))                     # equal-area annuli
y_edges = np.linspace(Q0[:, 1].min() - 1e-6, Q0[:, 1].max() + 1e-6, NY + 1)

def lacey(Q, min_per_bin=50):
    r = np.hypot(Q[:, 0], Q[:, 2])
    b = (np.clip(np.digitize(r, r_edges) - 1, 0, NR - 1) * NY
         + np.clip(np.digitize(Q[:, 1], y_edges) - 1, 0, NY - 1))
    n = np.bincount(b, minlength=NR*NY)
    k = np.bincount(b, weights=tracer.astype(float), minlength=NR*NY)
    ok = n >= min_per_bin
    x = k[ok] / n[ok]
    s0 = p_frac * (1 - p_frac)
    return (s0 - x.var(ddof=1)) / (s0 - s0 / n[ok].mean()), int(ok.sum()), float(n[ok].mean())

M0, n_bins, n_bar = lacey(Q0)
print("%d bins, %.0f grains per bin, tracer fraction p = %.4f" % (n_bins, n_bar, p_frac))
print("M at rest = %.4f  (0 = fully segregated by construction)" % M0)

40 bins, 553 grains per bin, tracer fraction p = 0.5000
M at rest = 0.0614  (0 = fully segregated by construction)

Step 4 — Turn the impeller

Both wall controls, every step, from the same OMEGA. The transform is built by integrating the angular velocity into a quaternion; the velocity field is handed the same ang_vel about the same axis.

OMEGA = 1.5                    # rad / time unit
NREV = 5.0
n_steps = int(NREV * 2 * np.pi / OMEGA / DT)
sample_every = 40
frame_every = 140

theta = 0.0
hist_rev, hist_M, hist_vmax = [0.0], [M0], [0.0]
frames = [Q0.copy()]
t0 = time.time()
for i in range(n_steps):
    theta += OMEGA * DT
    sim.set_wall_transform(w_stir, (0.0, 0.0, 0.0), qaxis([0, 1, 0], theta))   # the GEOMETRY
    sim.set_wall_velocity(w_stir, (0.0, 0.0, 0.0), (0.0, OMEGA, 0.0), (0.0, 0.0, 0.0))  # the SURFACE
    sim.step(DT)
    if i % sample_every == 0:
        Q = np.asarray(sim.get_positions()); V = np.asarray(sim.get_velocities())
        hist_rev.append(theta / (2*np.pi))
        hist_M.append(lacey(Q)[0])
        hist_vmax.append(float(np.linalg.norm(V, axis=1).max()))
    if i % frame_every == 0:
        frames.append(np.asarray(sim.get_positions()).copy())
wall_time = time.time() - t0
Qf = np.asarray(sim.get_positions())
sdf_f = np.asarray(sim.wall_sdf_at(w_col, np.ascontiguousarray(Qf.astype(np.float32))))
print("%d steps, %.1f revolutions, %.0f s  (%.1f ms/step, %d grains)"
      % (n_steps, NREV, wall_time, 1000*wall_time/n_steps, len(Qf)))
print("M: %.4f -> %.4f" % (M0, hist_M[-1]))
print("containment: %d grains outside the column, min wall SDF %+.3f (grain radius %.2f)"
      % (int((sdf_f < 0).sum()), sdf_f.min(), RP))

10471 steps, 5.0 revolutions, 474 s  (45.3 ms/step, 22121 grains)
M: 0.0614 -> 0.9177
containment: 0 grains outside the column, min wall SDF +0.379 (grain radius 0.40)

Code

fig, ax = plt.subplots(figsize=(5.6, 3.0))
ax.plot(hist_rev, hist_M, color="#4c72b0", lw=1.4)
ax.axhline(1.0, color="0.7", lw=0.8, ls=":")
ax.set_xlabel("impeller revolutions"); ax.set_ylabel("Lacey index $M$")
ax.set_ylim(0, 1.02); ax.set_xlim(0, NREV); ax.grid(alpha=0.3)
plt.show()

Figure 2: The Lacey index against impeller revolutions. It rises steeply while the blades cut the initial interface, then more slowly; the dips are real — a pitched blade periodically lifts a coherent slug of the lower layer, which briefly *increases* the between-bin variance before it disperses.

Step 5 — The guardrail: no free energy

A moving wall is an energy source, and a wall whose geometry and velocity field disagree is an unphysical one. The cheapest sharp test is that no grain should end up much faster than the fastest point of the machine that is pushing it: the blade tip, at $\omega R_{\text{imp}}$. The right statistic is the median of the per-sample peak speeds — the worst single sample over thousands of steps will always catch a grain being squeezed between a blade and the wall, which is real physics and not an energy leak. Both are reported.

tip = OMEGA * R_IMP
print("blade-tip speed        %.3f" % tip)
print("peak grain speed       %.3f   (ratio %.3f)" % (max(hist_vmax), max(hist_vmax)/tip))
print("median of the peaks    %.3f   (ratio %.3f)"
      % (np.median(hist_vmax[1:]), np.median(hist_vmax[1:])/tip))

blade-tip speed        24.000
peak grain speed       33.582   (ratio 1.399)
median of the peaks    24.241   (ratio 1.010)

Code

fig, ax = plt.subplots(figsize=(5.6, 2.4))
ax.plot(hist_rev[1:], hist_vmax[1:], color="#c44e52", lw=1.0, label="peak grain speed")
ax.axhline(tip, color="#4c72b0", lw=1.2, ls="--", label=r"blade tip $\omega R_{\rm imp}$")
ax.set_xlabel("impeller revolutions"); ax.set_ylabel("speed")
ax.set_xlim(0, NREV); ax.set_ylim(0, 1.6*tip); ax.grid(alpha=0.3)
ax.legend(fontsize=8, frameon=False)
plt.show()

Figure 3: Peak grain speed over the run against the blade-tip speed. Staying at the tip speed is the statement that the wall is doing work on the grains and not manufacturing energy; the earlier frictionless version of this page sat four to six times higher.

Step 6 — The recirculation pattern

Averaging the grain velocities azimuthally onto an $(r, y)$ grid gives the meridional flow. A down-pumping pitched-blade impeller drives material down at the blades, out along the floor, up the wall and back inward across the surface.

Code

Q = np.asarray(sim.get_positions()); V = np.asarray(sim.get_velocities())
r = np.hypot(Q[:, 0], Q[:, 2])
ur = (Q[:, 0]*V[:, 0] + Q[:, 2]*V[:, 2]) / np.maximum(r, 1e-9)
nr_, ny_ = 12, 12
rb = np.linspace(0, R_COL, nr_ + 1); yb = np.linspace(0, Q[:, 1].max(), ny_ + 1)
ir = np.clip(np.digitize(r, rb) - 1, 0, nr_ - 1); iy = np.clip(np.digitize(Q[:, 1], yb) - 1, 0, ny_ - 1)
b = ir*ny_ + iy
cnt = np.bincount(b, minlength=nr_*ny_).astype(float)
UR = np.bincount(b, weights=ur, minlength=nr_*ny_) / np.maximum(cnt, 1)
UY = np.bincount(b, weights=V[:, 1], minlength=nr_*ny_) / np.maximum(cnt, 1)
mask = cnt.reshape(nr_, ny_) >= 30
RC = 0.5*(rb[:-1] + rb[1:]); YC = 0.5*(yb[:-1] + yb[1:])
RG, YG = np.meshgrid(RC, YC, indexing="ij")
fig, ax = plt.subplots(figsize=(5.0, 3.2))
sp = np.hypot(UR.reshape(nr_, ny_), UY.reshape(nr_, ny_))
ax.quiver(RG[mask], YG[mask], UR.reshape(nr_, ny_)[mask], UY.reshape(nr_, ny_)[mask],
          sp[mask], cmap="viridis", scale=60, width=0.005)
ax.axhspan(Y_IMP - 1.8, Y_IMP + 1.8, xmin=1.2/R_COL, xmax=R_IMP/R_COL,
           color="0.85", zorder=0, label="blade band")
ax.set_xlabel("$r$"); ax.set_ylabel("$y$"); ax.set_xlim(0, R_COL); ax.set_ylim(0, Q[:, 1].max())
ax.legend(fontsize=8, frameon=False, loc="upper right")
plt.show()

Figure 4: Azimuthally averaged meridional velocity of the grains. The pattern — downward at the blades, outward along the base, up the column wall and inward across the free surface — is the single-loop circulation reported for four-bladed mixers in the DEM literature. The comparison is QUALITATIVE: the pattern, not the numbers.

The mixing, animated

Code

slab = np.abs(Q0[:, 2]) < 2.0
figm, axm = plt.subplots(figsize=(4.2, 2.6))

def draw(i):
    axm.clear()
    P = frames[i]
    axm.scatter(P[slab & tracer, 0], P[slab & tracer, 1], s=3.0, c="#c44e52", lw=0)
    axm.scatter(P[slab & ~tracer, 0], P[slab & ~tracer, 1], s=3.0, c="#4c72b0", lw=0)
    axm.plot([-R_COL, -R_COL], [0, 14], color="0.4", lw=1.4)
    axm.plot([R_COL, R_COL], [0, 14], color="0.4", lw=1.4)
    axm.plot([-R_COL, R_COL], [0, 0], color="0.4", lw=1.4)
    axm.set_xlim(-R_COL-1, R_COL+1); axm.set_ylim(-0.5, 14)
    axm.set_aspect("equal"); axm.set_xticks([]); axm.set_yticks([])
    axm.set_title("%.2f revolutions" % (i * frame_every * DT * OMEGA / (2*np.pi)), fontsize=8)
    return []

anim = animation.FuncAnimation(figm, draw, frames=len(frames), interval=80, blit=False)
plt.close(figm)
from IPython.display import HTML
HTML(anim.to_jshtml(fps=12))

Figure 5: A vertical slab through the column, grains coloured by their initial layer. The blades cut the interface on the first pass and the pitched pumping folds the two layers together.

Step 7 — Determinism, honestly

peclet.dem’s step is deterministic only single-threaded; the multithreaded contact solve reduces in a nondeterministic order. The showcase run above used 8 threads, so it is a showcase, not a reproducible trajectory. To say what that costs, here is a scaled-down column run three times: twice at one thread (which must agree bit for bit) and once at eight.

Threads cannot be changed after Kokkos initialises, so this check runs the reduced case in subprocesses.

Code

import json, subprocess, tempfile, pathlib, textwrap

CHECK = textwrap.dedent("""
    import json, os, sys, numpy as np
    for _p in os.environ.get("PECLET_LOCAL_BUILD", "").split(os.pathsep):
        if _p:
            sys.path.insert(0, _p)
    from peclet.core import geom
    from peclet import dem as pdem

    RP, R_COL, H_COL, R_IMP, Y_IMP, G = 0.4, 10.0, 14.0, 8.0, 3.0, 9.81
    PITCH, NBLADE, DT, OM = np.deg2rad(45.0), 4, 2e-3, 1.5

    def qmul(a, b):
        ax, ay, az, aw = a; bx, by, bz, bw = b
        return (aw*bx+ax*bw+ay*bz-az*by, aw*by-ax*bz+ay*bw+az*bx,
                aw*bz+ax*by-ay*bx+az*bw, aw*bw-ax*bx-ay*by-az*bz)

    def qaxis(axis, angle):
        s_ = np.sin(0.5*angle); n = np.asarray(axis, float); n /= np.linalg.norm(n)
        return (float(n[0]*s_), float(n[1]*s_), float(n[2]*s_), float(np.cos(0.5*angle)))

    def qrot(q, v):
        x, y, z, w = q; v = np.asarray(v, float); qv = np.array([x, y, z])
        t = 2*np.cross(qv, v); return v + w*t + np.cross(qv, t)

    b = geom.SceneBuilder()
    root = b.add_intersection(
        b.add_leaf("capsule", [R_COL, H_COL], translation=[0.0, 0.5*H_COL, 0.0]),
        b.add_leaf("box", [R_COL*1.5, 0.5*H_COL, R_COL*1.5], translation=[0.0, 0.5*H_COL, 0.0]))
    cni, cnr = b.encode()[0], b.encode()[1]
    b2 = geom.SceneBuilder()
    r_hub, t_b, w_b = 1.2, 0.35, 1.5
    L_b = 0.5*(R_IMP - r_hub)
    sroot = b2.add_leaf("capsule", [1.0, 0.5*H_COL], translation=[0.0, 0.5*H_COL, 0.0])
    for k in range(NBLADE):
        th0 = 2*np.pi*k/NBLADE
        q = qmul(qaxis([0, 1, 0], th0), qaxis([1, 0, 0], PITCH))
        c = qrot(qaxis([0, 1, 0], th0), [r_hub + L_b, 0.0, 0.0])
        sroot = b2.add_union(sroot, b2.add_leaf(
            "box", [L_b, t_b, w_b],
            translation=[float(c[0]), Y_IMP + float(c[1]), float(c[2])], rotation=list(q)))
    sni, snr = b2.encode()[0], b2.encode()[1]

    sim = pdem.Simulation(20000)
    sim.set_domain((-R_COL-4, -4.0, -R_COL-4), (R_COL+4, H_COL+4, R_COL+4))
    sim.enable_periodicity(False, False, False)
    sim.set_gravity(0.0, -G, 0.0); sim.set_sphere_shape(1.0); sim.set_global_scale(RP)
    wc = sim.add_analytic_wall(np.asarray(cni, np.int32), np.asarray(cnr, np.float32),
                               root, True, 0.2, 0.5)
    ws = sim.add_analytic_wall(np.asarray(sni, np.int32), np.asarray(snr, np.float32),
                               sroot, False, 0.2, 0.5)
    sim.set_material_params(0.2, 0.0, 0.5); sim.set_solver_iterations(4, 4)

    rng = np.random.default_rng(7); d = 2.05*RP
    xs = np.arange(-R_COL + RP*1.2, R_COL - RP*1.2, d)
    ys = np.arange(RP*1.2, 6.0, d)
    X, Y, Z = np.meshgrid(xs, ys, xs, indexing="ij")
    P = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])
    P = P[np.hypot(P[:, 0], P[:, 2]) < R_COL - 1.3*RP]
    P = np.ascontiguousarray((P + rng.normal(0, 0.06*RP, P.shape)).astype(np.float32))
    P = np.ascontiguousarray(P[np.asarray(sim.wall_sdf_at(ws, P)) > RP*1.15])
    sim.set_positions(P); sim.set_inv_mass(np.full(len(P), 1.0, np.float32)); sim.set_dt(DT)
    for _ in range(300):
        sim.step(DT)
    Q0 = np.asarray(sim.get_positions()).copy()
    tracer = Q0[:, 1] < np.median(Q0[:, 1]); pf = tracer.mean()
    NR, NY = 4, 5
    re_ = np.sqrt(np.linspace(0, R_COL**2, NR+1))
    ye_ = np.linspace(Q0[:, 1].min()-1e-6, Q0[:, 1].max()+1e-6, NY+1)

    def lac(Q):
        r = np.hypot(Q[:, 0], Q[:, 2])
        bb = (np.clip(np.digitize(r, re_)-1, 0, NR-1)*NY
              + np.clip(np.digitize(Q[:, 1], ye_)-1, 0, NY-1))
        n = np.bincount(bb, minlength=NR*NY)
        k = np.bincount(bb, weights=tracer.astype(float), minlength=NR*NY)
        ok = n >= 20; x = k[ok]/n[ok]; s0 = pf*(1-pf)
        return float((s0 - x.var(ddof=1))/(s0 - s0/n[ok].mean()))

    th = 0.0
    for _ in range(int(0.5*2*np.pi/OM/DT)):
        th += OM*DT
        sim.set_wall_transform(ws, (0.0, 0.0, 0.0), qaxis([0, 1, 0], th))
        sim.set_wall_velocity(ws, (0.0, 0.0, 0.0), (0.0, OM, 0.0), (0.0, 0.0, 0.0))
        sim.step(DT)
    Q = np.asarray(sim.get_positions())
    print("RESULT " + json.dumps({"n": int(len(Q)), "M": lac(Q),
                                  "cksum": float(np.abs(Q, dtype=np.float64).sum())}))
""")

det = []
try:
    runner = pathlib.Path(tempfile.mkdtemp()) / "check.py"
    runner.write_text(CHECK)
    for threads in ("1", "1", "8"):
        env = dict(os.environ, OMP_NUM_THREADS=threads, OMP_PROC_BIND="false")
        out = subprocess.run([sys.executable, str(runner)], env=env, capture_output=True,
                             text=True, timeout=2400)
        line = [l for l in out.stdout.splitlines() if l.startswith("RESULT ")]
        det.append((threads, json.loads(line[-1][7:]) if line else None))
except Exception as exc:                      # never let the check break the page
    print("determinism check skipped:", exc)

for threads, res in det:
    print("  %s thread(s): %s" % (threads, "FAILED" if res is None else
          "n=%d   M=%.6f   position checksum=%.6f" % (res["n"], res["M"], res["cksum"])))
det_ok = len(det) == 3 and all(r is not None for _, r in det)
if det_ok:
    a, b_, c = (r for _, r in det)
    det_bitwise = (a["cksum"] == b_["cksum"])
    det_dM = abs(c["M"] - a["M"])
    print("  two 1-thread runs identical : %s   (checksum difference %.3e)"
          % (det_bitwise, abs(a["cksum"] - b_["cksum"])))
    print("  8-thread Lacey index differs from 1-thread by %.2e (absolute), on M ~ %.3f"
          % (det_dM, a["M"]))

  1 thread(s): n=2585   M=0.406278   position checksum=31386.681234
  1 thread(s): n=2585   M=0.406278   position checksum=31386.681234
  8 thread(s): n=2585   M=0.408983   position checksum=31224.660763
  two 1-thread runs identical : True   (checksum difference 0.000e+00)
  8-thread Lacey index differs from 1-thread by 2.70e-03 (absolute), on M ~ 0.406

Results

claim	measured
grains	22121
bins × grains per bin (fixed binning)	40 x 553
Lacey index at rest (segregated by construction)	0.0614
Lacey index after 5 revolutions	0.9177
blade-tip speed $\omega R_{\rm imp}$	24.000
median per-sample peak grain speed / tip speed	1.010
worst single sample / tip speed	1.399
grains outside the column	0
minimum wall SDF at the end (grain radius 0.40)	+0.379
max grain–grain overlap	0.0152
run cost	474 s, 45.3 ms/step, 8 threads
two 1-thread runs of the reduced case	bitwise identical
8-thread vs 1-thread Lacey index (reduced case)	2.70e-03 absolute

The headline: a genuinely swept pitched-blade impeller takes the bed from a Lacey index of 0.061 to 0.918 in 5 revolutions, with the fastest grain in a typical sample sitting at 1.01× the blade-tip speed (worst single sample over the whole run 1.40×, a grain caught between a blade and the wall) and not one grain outside the container. The recirculation pattern — down at the blades, out along the base, up the wall, inward across the surface — is the single meridional loop reported for four-bladed mixers by (Remy et al. 2009). That comparison is qualitative: the pattern, not the numbers. Our column, grain size, fill level and blade geometry are not theirs, and nothing on this page claims a quantitative match to a published mixer.

Adapt this yourself

Change the pitch. PITCH = 0 gives a flat (radial) blade: the axial pumping disappears and the Lacey curve flattens, because a radial blade mixes in $r$ and $\theta$ but not in $y$.
Change the blade count. NBLADE = 2 or 6; the tree is rebuilt automatically.
Move the impeller up and down as well as round: set_wall_transform takes a translation, so a reciprocating stirrer is one line, and set_wall_velocity’s lin_vel carries the matching surface velocity.
Break the contract on purpose. Comment out the set_wall_velocity call and rerun: the blades still sweep, but there is no tangential drag, and the guardrail plot shows it immediately.
Add the fluid. The same column tree can be handed to peclet.flow (SceneBuilder.encode() + set_scene) for a resolved or unresolved CFD-DEM version.

Reproduce this

PECLET_LOCAL_BUILD=/path/to/suite/dem/build_l4_omp:/path/to/suite/core/python/build_geom \
OMP_NUM_THREADS=8 OMP_PROC_BIND=false \
  quarto render examples/stirred-column/index.qmd --execute

The showcase run is multithreaded on purpose — it is a picture, not a trajectory comparison. Every numeric claim that compares two runs (Step 7) is made at OMP_NUM_THREADS=1, where peclet.dem is deterministic.

References

Lacey, P. M. C. 1954. “Developments in the Theory of Particle Mixing.” Journal of Applied Chemistry 4 (5): 257–68. https://doi.org/10.1002/jctb.5010040504.

Remy, Brenda, Johannes G. Khinast, and Benjamin J. Glasser. 2009. “Discrete Element Simulation of Free Flowing Grains in a Four-Bladed Mixer.” AIChE Journal 55 (8): 2035–48.

--- title: "Stirred column: a pitched-blade impeller that actually sweeps" subtitle: "A rotating analytic wall whose geometry actually sweeps, driving 22 000 grains through a Lacey mixing curve — with a no-free-energy guardrail that the first, frictionless version of this page failed by a factor of five." author: "Peclet" date: "2026-08-30" categories: [dem, sdf, csg, moving-walls, granular, mixing] jupyter: python3 --- [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/computational-chemical-engineering/peclet-examples/blob/main/examples/stirred-column/index.ipynb){target="_blank"}  Runs on a CPU in about ten minutes — 22 000 grains, no fluid. ## What you'll learn The [rotating drum](../rotating-drum/index.qmd) got away with a trick: a barrel is a body of revolution, so you can spin it by giving a *static* wall a rigid-body velocity field and never moving a single triangle. A stirrer blade cannot be faked that way. Its geometry has to sweep, or there is nothing to push the grains. This page uses `peclet.dem`'s two wall controls together, which is the whole point: 1. **`set_wall_transform(wall, translation, quat)`** places an *analytic* wall rigidly in the world. It moves the **geometry** — the CSG tree the narrow phase evaluates. Calls are **absolute**: the placement is composed onto the tree's *authored* root transform, not onto the previous frame's result, so a thousand revolutions accumulate no drift. 2. **`set_wall_velocity(wall, lin_vel, ang_vel, center)`** sets the **surface velocity** a grain in contact feels. Neither is inferred from the other. That is deliberate — it is what lets a drum have a velocity without moving, and it is why this example must drive **both** from the same `ang_vel` every step. Getting one without the other is a specific, diagnosable bug: transform without velocity gives a blade that shoves grains but exerts no tangential drag; velocity without transform gives the drum trick applied to a shape it does not fit, and the blades stand still while the grains swirl. You will also see a **container and an impeller written as CSG trees** rather than voxel grids, and a quantitative mixing measure — the **Lacey index** — defined explicitly rather than eyeballed. ```{python} #| label: bootstrap #| code-summary: "Environment bootstrap (installs peclet from PyPI on Colab/Binder)" import importlib.util, os, subprocess, sys # Set the thread count BEFORE the module is imported: Kokkos reads it at initialisation. os.environ.setdefault("OMP_NUM_THREADS", "8") os.environ.setdefault("OMP_PROC_BIND", "false") _local = os.environ.get("PECLET_LOCAL_BUILD") if _local: for p in _local.split(os.pathsep): sys.path.insert(0, p) elif importlib.util.find_spec("peclet") is None: subprocess.run([sys.executable, "-m", "pip", "install", "-q", "peclet", "scikit-image"], check=True) ``` ::: {.callout-note} ## Requires a recent `peclet` `set_wall_transform` and `wall_sdf_at` are new, and the CSG authoring layer (`peclet.core.geom`) is newer than the current PyPI release; on Colab you need a source build of `peclet-core` and `peclet-dem` until the next release goes out. The page is *frozen*, so it renders and reads correctly either way. ::: ```{python} import time import numpy as np import matplotlib.pyplot as plt from matplotlib import animation from skimage import measure from peclet.core import geom from peclet import dem as pdem plt.rcParams.update({"figure.dpi": 130, "font.size": 9, "axes.axisbelow": True, "figure.facecolor": "white", "savefig.bbox": "tight"}) print("threads:", os.environ["OMP_NUM_THREADS"]) ``` ## Step 1 — The column and the impeller, as CSG trees Everything is in grain-radius-ish units: the grain radius is `RP`, gravity points along $-y$, and the column axis is $y$. The **column** is an *inverted* solid: `capsule ∩ box`. core's `capsule` is a $y$-axis primitive with hemispherical caps; intersecting it with a box slices those caps off and leaves a flat-bottomed cylinder. Both leaves are distance-exact and the intersection of two convex bodies keeps that exactness on the inside, which is where the grains live. `invert=True` turns the solid into a container: the wall SDF becomes **positive in the void**. The **impeller** is a union: a shaft capsule plus four flat blades, each a box rotated by 45° about its own radial axis (the *pitch*) and then swung to its azimuthal station. No inversion — the wall SDF is positive *outside* the metal, again where the grains are. ```{python} #| label: geometry RP = 0.4 # grain radius R_COL = 20.0 # column inner radius H_COL = 24.0 # column height (flat bottom at y = 0) R_IMP = 16.0 # blade tip radius Y_IMP = 3.5 # impeller centre height PITCH = np.deg2rad(45.0) NBLADE = 4 G = 9.81 def qmul(a, b): ax, ay, az, aw = a; bx, by, bz, bw = b return (aw*bx + ax*bw + ay*bz - az*by, aw*by - ax*bz + ay*bw + az*bx, aw*bz + ax*by - ay*bx + az*bw, aw*bw - ax*bx - ay*by - az*bz) def qaxis(axis, angle): """Quaternion (x, y, z, w) for a rotation of `angle` about `axis`.""" s = np.sin(0.5 * angle); n = np.asarray(axis, float); n /= np.linalg.norm(n) return (float(n[0]*s), float(n[1]*s), float(n[2]*s), float(np.cos(0.5*angle))) def qrot(q, v): x, y, z, w = q; v = np.asarray(v, float) qv = np.array([x, y, z]); t = 2 * np.cross(qv, v) return v + w * t + np.cross(qv, t) def column_tree(): b = geom.SceneBuilder() side = b.add_leaf("capsule", [R_COL, H_COL], translation=[0.0, 0.5*H_COL, 0.0]) slab = b.add_leaf("box", [R_COL*1.5, 0.5*H_COL, R_COL*1.5], translation=[0.0, 0.5*H_COL, 0.0]) root = b.add_intersection(side, slab) # flat-ended cylinder ni, nr, _, _ = b.encode() return np.asarray(ni, np.int32), np.asarray(nr, np.float32), root def stirrer_tree(): b = geom.SceneBuilder() r_hub, t_b, w_b = 1.2, 0.35, 1.8 # hub radius, blade half-thickness, half-width L_b = 0.5 * (R_IMP - r_hub) # blade radial half-length root = b.add_leaf("capsule", [1.0, 0.5*H_COL], translation=[0.0, 0.5*H_COL, 0.0]) # shaft for k in range(NBLADE): theta = 2*np.pi*k/NBLADE q = qmul(qaxis([0, 1, 0], theta), qaxis([1, 0, 0], PITCH)) # swing, then pitch c = qrot(qaxis([0, 1, 0], theta), [r_hub + L_b, 0.0, 0.0]) blade = b.add_leaf("box", [L_b, t_b, w_b], translation=[float(c[0]), Y_IMP + float(c[1]), float(c[2])], rotation=list(q)) root = b.add_union(root, blade) ni, nr, _, _ = b.encode() return np.asarray(ni, np.int32), np.asarray(nr, np.float32), root ``` `wall_sdf_at` returns exactly what the narrow phase reads, so it is the honest way to check a placement — and the easiest way to draw the impeller. ```{python} #| label: fig-geometry #| fig-cap: "The impeller, meshed from the wall SDF the narrow phase actually evaluates, inside the column outline. Left: the authored pose. Right: the same wall after `set_wall_transform` swings it 45°, which is what moves the metal." #| code-fold: true def make_sim(capacity=40000): s = pdem.Simulation(capacity) s.set_domain((-R_COL-4, -4.0, -R_COL-4), (R_COL+4, H_COL+4, R_COL+4)) s.enable_periodicity(False, False, False) s.set_gravity(0.0, -G, 0.0) s.set_sphere_shape(1.0) s.set_global_scale(RP) # ghost band is sized off the real radius; without this the # default globalScale = 1 makes it 10x the grain cni, cnr, croot = column_tree() sni, snr, sroot = stirrer_tree() w_col = s.add_analytic_wall(cni, cnr, croot, True, 0.2, 0.5) # invert -> container w_stir = s.add_analytic_wall(sni, snr, sroot, False, 0.2, 0.5) s.set_material_params(0.2, 0.0, 0.5) # body-body restitution / tangential / friction s.set_solver_iterations(4, 4) return s, w_col, w_stir def stirrer_mesh(sim, wall, box=(-R_IMP-1, R_IMP+1), ylim=(0.0, 8.0), n=110): xs = np.linspace(box[0], box[1], n) ys = np.linspace(ylim[0], ylim[1], max(24, int(n*(ylim[1]-ylim[0])/(box[1]-box[0])))) X, Y, Z = np.meshgrid(xs, ys, xs, indexing="ij") pts = np.ascontiguousarray(np.column_stack([X.ravel(), Y.ravel(), Z.ravel()]).astype(np.float32)) fld = np.asarray(sim.wall_sdf_at(wall, pts)).reshape(X.shape) v, f, _, _ = measure.marching_cubes(fld, level=0.0, spacing=(xs[1]-xs[0], ys[1]-ys[0], xs[1]-xs[0])) v = v + np.array([xs[0], ys[0], xs[0]]) return v, f sim, w_col, w_stir = make_sim() # the ONE simulation this page uses figg = plt.figure(figsize=(7.0, 3.2)) from mpl_toolkits.mplot3d.art3d import Poly3DCollection for panel, ang in enumerate((0.0, np.deg2rad(45.0))): sim.set_wall_transform(w_stir, (0.0, 0.0, 0.0), qaxis([0, 1, 0], ang)) v, f = stirrer_mesh(sim, w_stir) ax = figg.add_subplot(1, 2, panel + 1, projection="3d") ax.add_collection3d(Poly3DCollection(v[f], facecolor="#8a8f98", edgecolor="none")) tt = np.linspace(0, 2*np.pi, 120) for yy in (0.0, 8.0): ax.plot(R_COL*np.cos(tt), np.full_like(tt, yy), R_COL*np.sin(tt), color="#4c72b0", lw=0.8) ax.set_xlim(-R_COL, R_COL); ax.set_zlim(-R_COL, R_COL); ax.set_ylim(0, 2*R_COL) ax.set_box_aspect((1, 1, 1)); ax.set_axis_off(); ax.view_init(elev=22, azim=-58) ax.set_title("authored pose" if panel == 0 else "after set_wall_transform(45°)", fontsize=8) plt.show() sim.set_wall_transform(w_stir, (0.0, 0.0, 0.0), qaxis([0, 1, 0], 0.0)) # back to the authored pose ``` ## Step 2 — Fill and settle A jittered lattice inside the column, with anything starting inside the impeller discarded — checked with `wall_sdf_at`, not with a hand-written formula that could disagree with the collision geometry. ```{python} #| label: fill def fill(sim, w_stir, fill_top=11.0, seed=4): rng = np.random.default_rng(seed) d = 2.05 * RP xs = np.arange(-R_COL + RP*1.2, R_COL - RP*1.2, d) ys = np.arange(RP*1.2, fill_top, d) X, Y, Z = np.meshgrid(xs, ys, xs, indexing="ij") P = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()]) P = P[np.hypot(P[:, 0], P[:, 2]) < R_COL - 1.3*RP] P = np.ascontiguousarray((P + rng.normal(0, 0.06*RP, P.shape)).astype(np.float32)) P = np.ascontiguousarray(P[np.asarray(sim.wall_sdf_at(w_stir, P)) > RP*1.15]) sim.set_positions(P) sim.set_inv_mass(np.full(len(P), 1.0, np.float32)) # dem assigns unit mass regardless of size return P DT = 2e-3 P0 = fill(sim, w_stir) sim.set_dt(DT) t0 = time.time() for _ in range(700): sim.step(DT) # dem.step() with NO argument advances nothing Q0 = np.asarray(sim.get_positions()).copy() print("%d grains settled in %.0f s bed y = %.2f .. %.2f max overlap %.4f (grain radius %.2f)" % (len(Q0), time.time() - t0, Q0[:, 1].min(), Q0[:, 1].max(), sim.max_overlap(), RP)) ``` ::: {.callout-important} ## Set the body–body material, or the bed leaks dem's default body–body material is frictionless with zero restitution. A *frictionless* deep bed behaves like a liquid: it puts full hydrostatic pressure on the container and the position solve cannot hold it, and grains push through the analytic wall. The first version of this page lost grains through the floor (centres 1.7 radii *below* it) and reached peak speeds four to six times the blade-tip speed. `set_material_params(0.2, 0.0, 0.5)` fixes it completely — see the guardrail in Step 5, and `ISSUES.md`. ::: ## Step 3 — The Lacey mixing index Label the bottom half of the settled bed as tracer, then measure how uniformly it spreads. With $p$ the overall tracer fraction and $x_i$ the tracer fraction in bin $i$, $$ M = \frac{s_0^2 - s^2}{s_0^2 - s_r^2},\qquad s_0^2 = p(1-p),\qquad s_r^2 = \frac{p(1-p)}{\bar n}, $$ {#eq-lacey} where $s^2$ is the variance of $x_i$ across bins, $s_0^2$ is the fully-segregated variance, $s_r^2$ the fully-random (binomial) one, and $\bar n$ the mean grains per bin. $M = 0$ is segregated, $M = 1$ is random [@lacey1954mixing]. The binning is **fixed** for the whole run — equal-area annuli in $r$ crossed with equal slices in $y$ — so nothing about the measure moves with the grains. ```{python} #| label: lacey NR, NY = 5, 8 tracer = Q0[:, 1] < np.median(Q0[:, 1]) p_frac = tracer.mean() r_edges = np.sqrt(np.linspace(0, R_COL**2, NR + 1)) # equal-area annuli y_edges = np.linspace(Q0[:, 1].min() - 1e-6, Q0[:, 1].max() + 1e-6, NY + 1) def lacey(Q, min_per_bin=50): r = np.hypot(Q[:, 0], Q[:, 2]) b = (np.clip(np.digitize(r, r_edges) - 1, 0, NR - 1) * NY + np.clip(np.digitize(Q[:, 1], y_edges) - 1, 0, NY - 1)) n = np.bincount(b, minlength=NR*NY) k = np.bincount(b, weights=tracer.astype(float), minlength=NR*NY) ok = n >= min_per_bin x = k[ok] / n[ok] s0 = p_frac * (1 - p_frac) return (s0 - x.var(ddof=1)) / (s0 - s0 / n[ok].mean()), int(ok.sum()), float(n[ok].mean()) M0, n_bins, n_bar = lacey(Q0) print("%d bins, %.0f grains per bin, tracer fraction p = %.4f" % (n_bins, n_bar, p_frac)) print("M at rest = %.4f (0 = fully segregated by construction)" % M0) ``` ## Step 4 — Turn the impeller Both wall controls, every step, from the same `OMEGA`. The transform is built by integrating the angular velocity into a quaternion; the velocity field is handed the same `ang_vel` about the same axis. ```{python} #| label: run OMEGA = 1.5 # rad / time unit NREV = 5.0 n_steps = int(NREV * 2 * np.pi / OMEGA / DT) sample_every = 40 frame_every = 140 theta = 0.0 hist_rev, hist_M, hist_vmax = [0.0], [M0], [0.0] frames = [Q0.copy()] t0 = time.time() for i in range(n_steps): theta += OMEGA * DT sim.set_wall_transform(w_stir, (0.0, 0.0, 0.0), qaxis([0, 1, 0], theta)) # the GEOMETRY sim.set_wall_velocity(w_stir, (0.0, 0.0, 0.0), (0.0, OMEGA, 0.0), (0.0, 0.0, 0.0)) # the SURFACE sim.step(DT) if i % sample_every == 0: Q = np.asarray(sim.get_positions()); V = np.asarray(sim.get_velocities()) hist_rev.append(theta / (2*np.pi)) hist_M.append(lacey(Q)[0]) hist_vmax.append(float(np.linalg.norm(V, axis=1).max())) if i % frame_every == 0: frames.append(np.asarray(sim.get_positions()).copy()) wall_time = time.time() - t0 Qf = np.asarray(sim.get_positions()) sdf_f = np.asarray(sim.wall_sdf_at(w_col, np.ascontiguousarray(Qf.astype(np.float32)))) print("%d steps, %.1f revolutions, %.0f s (%.1f ms/step, %d grains)" % (n_steps, NREV, wall_time, 1000*wall_time/n_steps, len(Qf))) print("M: %.4f -> %.4f" % (M0, hist_M[-1])) print("containment: %d grains outside the column, min wall SDF %+.3f (grain radius %.2f)" % (int((sdf_f < 0).sum()), sdf_f.min(), RP)) ``` ```{python} #| label: fig-lacey #| fig-cap: "The Lacey index against impeller revolutions. It rises steeply while the blades cut the initial interface, then more slowly; the dips are real — a pitched blade periodically lifts a coherent slug of the lower layer, which briefly *increases* the between-bin variance before it disperses." #| code-fold: true fig, ax = plt.subplots(figsize=(5.6, 3.0)) ax.plot(hist_rev, hist_M, color="#4c72b0", lw=1.4) ax.axhline(1.0, color="0.7", lw=0.8, ls=":") ax.set_xlabel("impeller revolutions"); ax.set_ylabel("Lacey index $M$") ax.set_ylim(0, 1.02); ax.set_xlim(0, NREV); ax.grid(alpha=0.3) plt.show() ``` ## Step 5 — The guardrail: no free energy A moving wall is an energy source, and a wall whose geometry and velocity field disagree is an *unphysical* one. The cheapest sharp test is that no grain should end up much faster than the fastest point of the machine that is pushing it: the blade tip, at $\omega R_{\text{imp}}$. The right statistic is the **median** of the per-sample peak speeds — the worst single sample over thousands of steps will always catch a grain being squeezed between a blade and the wall, which is real physics and not an energy leak. Both are reported. ```{python} #| label: guardrail tip = OMEGA * R_IMP print("blade-tip speed %.3f" % tip) print("peak grain speed %.3f (ratio %.3f)" % (max(hist_vmax), max(hist_vmax)/tip)) print("median of the peaks %.3f (ratio %.3f)" % (np.median(hist_vmax[1:]), np.median(hist_vmax[1:])/tip)) ``` ```{python} #| label: fig-guardrail #| fig-cap: "Peak grain speed over the run against the blade-tip speed. Staying at the tip speed is the statement that the wall is doing work on the grains and not manufacturing energy; the earlier frictionless version of this page sat four to six times higher." #| code-fold: true fig, ax = plt.subplots(figsize=(5.6, 2.4)) ax.plot(hist_rev[1:], hist_vmax[1:], color="#c44e52", lw=1.0, label="peak grain speed") ax.axhline(tip, color="#4c72b0", lw=1.2, ls="--", label=r"blade tip $\omega R_{\rm imp}$") ax.set_xlabel("impeller revolutions"); ax.set_ylabel("speed") ax.set_xlim(0, NREV); ax.set_ylim(0, 1.6*tip); ax.grid(alpha=0.3) ax.legend(fontsize=8, frameon=False) plt.show() ``` ## Step 6 — The recirculation pattern Averaging the grain velocities azimuthally onto an $(r, y)$ grid gives the meridional flow. A **down-pumping** pitched-blade impeller drives material down at the blades, out along the floor, up the wall and back inward across the surface. ```{python} #| label: fig-recirc #| fig-cap: "Azimuthally averaged meridional velocity of the grains. The pattern — downward at the blades, outward along the base, up the column wall and inward across the free surface — is the single-loop circulation reported for four-bladed mixers in the DEM literature. The comparison is QUALITATIVE: the pattern, not the numbers." #| code-fold: true Q = np.asarray(sim.get_positions()); V = np.asarray(sim.get_velocities()) r = np.hypot(Q[:, 0], Q[:, 2]) ur = (Q[:, 0]*V[:, 0] + Q[:, 2]*V[:, 2]) / np.maximum(r, 1e-9) nr_, ny_ = 12, 12 rb = np.linspace(0, R_COL, nr_ + 1); yb = np.linspace(0, Q[:, 1].max(), ny_ + 1) ir = np.clip(np.digitize(r, rb) - 1, 0, nr_ - 1); iy = np.clip(np.digitize(Q[:, 1], yb) - 1, 0, ny_ - 1) b = ir*ny_ + iy cnt = np.bincount(b, minlength=nr_*ny_).astype(float) UR = np.bincount(b, weights=ur, minlength=nr_*ny_) / np.maximum(cnt, 1) UY = np.bincount(b, weights=V[:, 1], minlength=nr_*ny_) / np.maximum(cnt, 1) mask = cnt.reshape(nr_, ny_) >= 30 RC = 0.5*(rb[:-1] + rb[1:]); YC = 0.5*(yb[:-1] + yb[1:]) RG, YG = np.meshgrid(RC, YC, indexing="ij") fig, ax = plt.subplots(figsize=(5.0, 3.2)) sp = np.hypot(UR.reshape(nr_, ny_), UY.reshape(nr_, ny_)) ax.quiver(RG[mask], YG[mask], UR.reshape(nr_, ny_)[mask], UY.reshape(nr_, ny_)[mask], sp[mask], cmap="viridis", scale=60, width=0.005) ax.axhspan(Y_IMP - 1.8, Y_IMP + 1.8, xmin=1.2/R_COL, xmax=R_IMP/R_COL, color="0.85", zorder=0, label="blade band") ax.set_xlabel("$r$"); ax.set_ylabel("$y$"); ax.set_xlim(0, R_COL); ax.set_ylim(0, Q[:, 1].max()) ax.legend(fontsize=8, frameon=False, loc="upper right") plt.show() ``` ## The mixing, animated ```{python} #| label: fig-movie #| fig-cap: "A vertical slab through the column, grains coloured by their initial layer. The blades cut the interface on the first pass and the pitched pumping folds the two layers together." #| code-fold: true slab = np.abs(Q0[:, 2]) < 2.0 figm, axm = plt.subplots(figsize=(4.2, 2.6)) def draw(i): axm.clear() P = frames[i] axm.scatter(P[slab & tracer, 0], P[slab & tracer, 1], s=3.0, c="#c44e52", lw=0) axm.scatter(P[slab & ~tracer, 0], P[slab & ~tracer, 1], s=3.0, c="#4c72b0", lw=0) axm.plot([-R_COL, -R_COL], [0, 14], color="0.4", lw=1.4) axm.plot([R_COL, R_COL], [0, 14], color="0.4", lw=1.4) axm.plot([-R_COL, R_COL], [0, 0], color="0.4", lw=1.4) axm.set_xlim(-R_COL-1, R_COL+1); axm.set_ylim(-0.5, 14) axm.set_aspect("equal"); axm.set_xticks([]); axm.set_yticks([]) axm.set_title("%.2f revolutions" % (i * frame_every * DT * OMEGA / (2*np.pi)), fontsize=8) return [] anim = animation.FuncAnimation(figm, draw, frames=len(frames), interval=80, blit=False) plt.close(figm) from IPython.display import HTML HTML(anim.to_jshtml(fps=12)) ``` ## Step 7 — Determinism, honestly `peclet.dem`'s step is deterministic only **single-threaded**; the multithreaded contact solve reduces in a nondeterministic order. The showcase run above used `{python} os.environ["OMP_NUM_THREADS"]` threads, so it is a *showcase*, not a reproducible trajectory. To say what that costs, here is a scaled-down column run three times: twice at one thread (which must agree bit for bit) and once at eight. Threads cannot be changed after Kokkos initialises, so this check runs the reduced case in subprocesses. ```{python} #| label: determinism #| code-fold: true import json, subprocess, tempfile, pathlib, textwrap CHECK = textwrap.dedent(""" import json, os, sys, numpy as np for _p in os.environ.get("PECLET_LOCAL_BUILD", "").split(os.pathsep): if _p: sys.path.insert(0, _p) from peclet.core import geom from peclet import dem as pdem RP, R_COL, H_COL, R_IMP, Y_IMP, G = 0.4, 10.0, 14.0, 8.0, 3.0, 9.81 PITCH, NBLADE, DT, OM = np.deg2rad(45.0), 4, 2e-3, 1.5 def qmul(a, b): ax, ay, az, aw = a; bx, by, bz, bw = b return (aw*bx+ax*bw+ay*bz-az*by, aw*by-ax*bz+ay*bw+az*bx, aw*bz+ax*by-ay*bx+az*bw, aw*bw-ax*bx-ay*by-az*bz) def qaxis(axis, angle): s_ = np.sin(0.5*angle); n = np.asarray(axis, float); n /= np.linalg.norm(n) return (float(n[0]*s_), float(n[1]*s_), float(n[2]*s_), float(np.cos(0.5*angle))) def qrot(q, v): x, y, z, w = q; v = np.asarray(v, float); qv = np.array([x, y, z]) t = 2*np.cross(qv, v); return v + w*t + np.cross(qv, t) b = geom.SceneBuilder() root = b.add_intersection( b.add_leaf("capsule", [R_COL, H_COL], translation=[0.0, 0.5*H_COL, 0.0]), b.add_leaf("box", [R_COL*1.5, 0.5*H_COL, R_COL*1.5], translation=[0.0, 0.5*H_COL, 0.0])) cni, cnr = b.encode()[0], b.encode()[1] b2 = geom.SceneBuilder() r_hub, t_b, w_b = 1.2, 0.35, 1.5 L_b = 0.5*(R_IMP - r_hub) sroot = b2.add_leaf("capsule", [1.0, 0.5*H_COL], translation=[0.0, 0.5*H_COL, 0.0]) for k in range(NBLADE): th0 = 2*np.pi*k/NBLADE q = qmul(qaxis([0, 1, 0], th0), qaxis([1, 0, 0], PITCH)) c = qrot(qaxis([0, 1, 0], th0), [r_hub + L_b, 0.0, 0.0]) sroot = b2.add_union(sroot, b2.add_leaf( "box", [L_b, t_b, w_b], translation=[float(c[0]), Y_IMP + float(c[1]), float(c[2])], rotation=list(q))) sni, snr = b2.encode()[0], b2.encode()[1] sim = pdem.Simulation(20000) sim.set_domain((-R_COL-4, -4.0, -R_COL-4), (R_COL+4, H_COL+4, R_COL+4)) sim.enable_periodicity(False, False, False) sim.set_gravity(0.0, -G, 0.0); sim.set_sphere_shape(1.0); sim.set_global_scale(RP) wc = sim.add_analytic_wall(np.asarray(cni, np.int32), np.asarray(cnr, np.float32), root, True, 0.2, 0.5) ws = sim.add_analytic_wall(np.asarray(sni, np.int32), np.asarray(snr, np.float32), sroot, False, 0.2, 0.5) sim.set_material_params(0.2, 0.0, 0.5); sim.set_solver_iterations(4, 4) rng = np.random.default_rng(7); d = 2.05*RP xs = np.arange(-R_COL + RP*1.2, R_COL - RP*1.2, d) ys = np.arange(RP*1.2, 6.0, d) X, Y, Z = np.meshgrid(xs, ys, xs, indexing="ij") P = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()]) P = P[np.hypot(P[:, 0], P[:, 2]) < R_COL - 1.3*RP] P = np.ascontiguousarray((P + rng.normal(0, 0.06*RP, P.shape)).astype(np.float32)) P = np.ascontiguousarray(P[np.asarray(sim.wall_sdf_at(ws, P)) > RP*1.15]) sim.set_positions(P); sim.set_inv_mass(np.full(len(P), 1.0, np.float32)); sim.set_dt(DT) for _ in range(300): sim.step(DT) Q0 = np.asarray(sim.get_positions()).copy() tracer = Q0[:, 1] < np.median(Q0[:, 1]); pf = tracer.mean() NR, NY = 4, 5 re_ = np.sqrt(np.linspace(0, R_COL**2, NR+1)) ye_ = np.linspace(Q0[:, 1].min()-1e-6, Q0[:, 1].max()+1e-6, NY+1) def lac(Q): r = np.hypot(Q[:, 0], Q[:, 2]) bb = (np.clip(np.digitize(r, re_)-1, 0, NR-1)*NY + np.clip(np.digitize(Q[:, 1], ye_)-1, 0, NY-1)) n = np.bincount(bb, minlength=NR*NY) k = np.bincount(bb, weights=tracer.astype(float), minlength=NR*NY) ok = n >= 20; x = k[ok]/n[ok]; s0 = pf*(1-pf) return float((s0 - x.var(ddof=1))/(s0 - s0/n[ok].mean())) th = 0.0 for _ in range(int(0.5*2*np.pi/OM/DT)): th += OM*DT sim.set_wall_transform(ws, (0.0, 0.0, 0.0), qaxis([0, 1, 0], th)) sim.set_wall_velocity(ws, (0.0, 0.0, 0.0), (0.0, OM, 0.0), (0.0, 0.0, 0.0)) sim.step(DT) Q = np.asarray(sim.get_positions()) print("RESULT " + json.dumps({"n": int(len(Q)), "M": lac(Q), "cksum": float(np.abs(Q, dtype=np.float64).sum())})) """) det = [] try: runner = pathlib.Path(tempfile.mkdtemp()) / "check.py" runner.write_text(CHECK) for threads in ("1", "1", "8"): env = dict(os.environ, OMP_NUM_THREADS=threads, OMP_PROC_BIND="false") out = subprocess.run([sys.executable, str(runner)], env=env, capture_output=True, text=True, timeout=2400) line = [l for l in out.stdout.splitlines() if l.startswith("RESULT ")] det.append((threads, json.loads(line[-1][7:]) if line else None)) except Exception as exc: # never let the check break the page print("determinism check skipped:", exc) for threads, res in det: print(" %s thread(s): %s" % (threads, "FAILED" if res is None else "n=%d M=%.6f position checksum=%.6f" % (res["n"], res["M"], res["cksum"]))) det_ok = len(det) == 3 and all(r is not None for _, r in det) if det_ok: a, b_, c = (r for _, r in det) det_bitwise = (a["cksum"] == b_["cksum"]) det_dM = abs(c["M"] - a["M"]) print(" two 1-thread runs identical : %s (checksum difference %.3e)" % (det_bitwise, abs(a["cksum"] - b_["cksum"]))) print(" 8-thread Lacey index differs from 1-thread by %.2e (absolute), on M ~ %.3f" % (det_dM, a["M"])) ``` ## Results | claim | measured | |---|---| | grains | `{python} "%d" % len(Q0)` | | bins × grains per bin (fixed binning) | `{python} "%d x %.0f" % (n_bins, n_bar)` | | Lacey index at rest (segregated by construction) | `{python} "%.4f" % M0` | | **Lacey index after `{python} "%.0f" % NREV` revolutions** | **`{python} "%.4f" % hist_M[-1]`** | | blade-tip speed $\omega R_{\rm imp}$ | `{python} "%.3f" % tip` | | **median per-sample peak grain speed / tip speed** | **`{python} "%.3f" % (np.median(hist_vmax[1:])/tip)`** | | worst single sample / tip speed | `{python} "%.3f" % (max(hist_vmax)/tip)` | | grains outside the column | `{python} "%d" % int((sdf_f < 0).sum())` | | minimum wall SDF at the end (grain radius `{python} "%.2f" % RP`) | `{python} "%+.3f" % sdf_f.min()` | | max grain–grain overlap | `{python} "%.4f" % sim.max_overlap()` | | run cost | `{python} "%.0f s, %.1f ms/step, %s threads" % (wall_time, 1000*wall_time/n_steps, os.environ["OMP_NUM_THREADS"])` | | two 1-thread runs of the reduced case | `{python} ("bitwise identical" if det_ok and det_bitwise else "see Step 7")` | | 8-thread vs 1-thread Lacey index (reduced case) | `{python} (("%.2e absolute" % det_dM) if det_ok else "see Step 7")` | The headline: a genuinely swept pitched-blade impeller takes the bed from a Lacey index of `{python} "%.3f" % M0` to **`{python} "%.3f" % hist_M[-1]`** in `{python} "%.0f" % NREV` revolutions, with the fastest grain in a typical sample sitting at **`{python} "%.2f" % (np.median(hist_vmax[1:])/tip)`× the blade-tip speed** (worst single sample over the whole run `{python} "%.2f" % (max(hist_vmax)/tip)`×, a grain caught between a blade and the wall) and **not one grain** outside the container. The recirculation pattern — down at the blades, out along the base, up the wall, inward across the surface — is the single meridional loop reported for four-bladed mixers by [@remy2009four]. That comparison is **qualitative: the pattern, not the numbers.** Our column, grain size, fill level and blade geometry are not theirs, and nothing on this page claims a quantitative match to a published mixer. ## Adapt this yourself - **Change the pitch.** `PITCH = 0` gives a flat (radial) blade: the axial pumping disappears and the Lacey curve flattens, because a radial blade mixes in $r$ and $\theta$ but not in $y$. - **Change the blade count.** `NBLADE = 2` or `6`; the tree is rebuilt automatically. - **Move the impeller up and down** as well as round: `set_wall_transform` takes a translation, so a reciprocating stirrer is one line, and `set_wall_velocity`'s `lin_vel` carries the matching surface velocity. - **Break the contract on purpose.** Comment out the `set_wall_velocity` call and rerun: the blades still sweep, but there is no tangential drag, and the guardrail plot shows it immediately. - **Add the fluid.** The same column tree can be handed to `peclet.flow` (`SceneBuilder.encode()` + `set_scene`) for a resolved or unresolved CFD-DEM version. ## Reproduce this ```bash PECLET_LOCAL_BUILD=/path/to/suite/dem/build_l4_omp:/path/to/suite/core/python/build_geom \ OMP_NUM_THREADS=8 OMP_PROC_BIND=false \ quarto render examples/stirred-column/index.qmd --execute ``` The showcase run is multithreaded on purpose — it is a picture, not a trajectory comparison. Every *numeric* claim that compares two runs (Step 7) is made at `OMP_NUM_THREADS=1`, where `peclet.dem` is deterministic.