VeloxSim-PBD

Position-Based Fluids as a Newton solver, GPU-accelerated via NVIDIA Warp.

SolverPBF is a drop-in newton.solvers.SolverBase subclass that runs the Macklin & Müller 2013 PBF density constraint with the Macklin et al. 2016 XPBD λ-accumulator on top of Newton's particle pipeline.

Install

pip install -e .
# for the visualisation script:
pip install -e ".[viz]"

Requires Python ≥ 3.10, a CUDA-capable GPU, NVIDIA Warp ≥ 1.1, and newton-physics.

Quick start — 3D dam break

# 5 s sim (300 frames) → snapshots_dambreak_3d_newton/frames.npz
python examples/dambreak_3d.py --headless --frames 300

# Standalone interactive HTML viewer (opens in browser)
python examples/view_dambreak.py snapshots_dambreak_3d_newton/frames.npz --html \
    --html-out snapshots_dambreak_3d_newton/dambreak_viewer.html

The viewer is a single self-contained HTML file (three.js loaded from a CDN, mesh data inlined as gzipped + quantised binary). Drag to rotate, scroll to zoom, scrub the slider through frames.

Using SolverPBF in your own scene

import newton
import warp as wp
from veloxsim_pbd.newton_pbf import PBFConfig, SolverPBF

builder = newton.ModelBuilder()
builder.add_ground_plane()
builder.add_particle_grid(
    pos=wp.vec3(0.0, 0.0, 0.05),
    rot=wp.quat_identity(),
    vel=wp.vec3(0.0),
    dim_x=20, dim_y=20, dim_z=40,
    cell_x=0.008, cell_y=0.008, cell_z=0.008,
    mass=0.000256,           # ρ₀ · (2r)³  for r = 4 mm, ρ = 1000
    jitter=0.0,
    radius_mean=0.004,
)
model = builder.finalize()

cfg = PBFConfig(
    particle_radius=0.004,
    kernel_radius=0.016,        # 4 · r
    rest_density=1000.0,
    iterations=10,
    pbf_compliance=1.0e-6,
    s_corr_k=1.0e-4,
    xsph_viscosity=0.08,
    # Wall-ghost density (PhysX-style) — particles within kernel range
    # of these planes get a virtual neighbour density added so the floor
    # layer compacts hydrostatically without fake-attractive artefacts.
    wall_z_min=0.0,
    wall_x_min=-0.4, wall_x_max=0.4,
    wall_y_min=-0.2, wall_y_max=0.2,
)
solver = SolverPBF(model, cfg)

state_0 = model.state()
state_1 = model.state()
contacts = model.contacts()

dt = 1.0 / 60.0
for _ in range(300):
    model.collide(state_0, contacts)
    solver.step(state_0, state_1, None, contacts, dt)
    state_0, state_1 = state_1, state_0

How it works

The solver runs the XPBD loop per substep:

predict (Newton integrator) → [iter]: density Δλ → Δp → apply → … → finalize → XSPH → vorticity

Key design choices:

1. XPBD over PBF — λ accumulates across iterations within a substep, so corrections converge rather than restart. Kills the classic PBF constraint-drift "boiling" at rest.
2. Constraint averaging (Macklin 2014 §4.2) — per-particle Δp is divided by the number of constraints affecting that particle (1 + num_neighbours). This is what makes parallel Jacobi iteration actually stable on GPU; without it you need either tiny step sizes or heavy under-relaxation.4
3. Wall-ghost density (PhysX-style) — a particle within kernel range of a wall gets m·W(d_to_wall) added to its density only (not to the constraint gradient), so the floor layer compacts hydrostatically without the fake-attractive-force artefacts of other approaches.
4. Unilateral constraint (Macklin 2014 eq. 26): C = max(ρ/ρ₀ − 1, 0) — over-density only. Under-density at the free surface gives λ = 0, no pull-together.

Config reference (src/veloxsim_pbd/newton_pbf/config.py)

param	default	purpose
particle_radius	—	resolution (you must set)
kernel_radius	—	SPH kernel support, typically 4·r
rest_density	1000.0	ρ₀ in kg/m³
iterations	4	XPBD iterations per substep
pbf_compliance (α)	1e-6	softer = more compressible
solver_relaxation (ω)	1.0	SOR factor (1…2)
s_corr_k	0.0	artificial pressure against clustering
xsph_viscosity	0.0	bulk velocity smoothing
vorticity_epsilon	0.0	vorticity-confinement strength
max_dp_clamp_h	0.0	per-iteration Δp cap, fraction of kernel_radius (0 = off)
wall_x_min, wall_x_max, wall_y_*, wall_z_*	None	bounds for the wall-ghost density term

Layout

src/veloxsim_pbd/
  surface_kernels.py        GPU anisotropic density field for surfacing
  newton_pbf/
    __init__.py             public API: PBFConfig, SolverPBF
    config.py               PBFConfig dataclass
    solver.py               SolverPBF — newton.solvers.SolverBase subclass
    kernels.py              Warp kernels (density / Δp / vorticity / XSPH)
    _newton_vendored.py     vendored Newton helpers (apply_deltas, contacts)

examples/
  dambreak_3d.py            3D corner-column dam-break demo
  view_dambreak.py          interactive viewer + standalone HTML export

References

• Macklin & Müller 2013. Position Based Fluids.
• Macklin, Müller, Chentanez, Kim 2014. Unified Particle Physics for Real-Time Applications.
• Macklin, Müller, Chentanez 2016. XPBD: position-based simulation of compliant constrained dynamics.
• Fedkiw, Stam, Jensen 2001. Visual Simulation of Smoke (vorticity confinement).
• NVIDIA Flex (PhysX 5 source).
• Newton — the robotics physics framework this integrates with.

License

MIT — see LICENSE.