Newton / Technical report
Created Modified

Engineering report · XPBD fluids

Newton XPBD Fluid Support

Simulation formulation, rendering, and public API

Newton Physics · CUDA / OpenGL

Upstream
origin/main@5ca5f49c
Feature
eric-heiden/flex-fluid
Hardware
GeForce RTX 4090

Abstract

Newton’s XPBD solver now treats selected particles as fluids using position-based density constraints, bounded cohesion, viscosity, optional vorticity confinement, grouped neighborhoods, and two-way particle–rigid contact. A render-only path fits smoothed anisotropic particles and reconstructs an optically shaded surface in ViewerGL without changing the simulated state. Nine examples cover free surfaces, pumping, multiple fluids and worlds, buoyancy, SDF containers, articulated manipulation, waves, foam, and spray. At approximately 100,000 particles, the measured dam-break simulation reaches a median 90.2 frames per second on an RTX 4090.

1. Scope

The branch adds a fluid particle type to the existing XPBD solver and a corresponding surface representation for ViewerGL.

Fluids use the same particle state and timestep

Particles marked ACTIVE | FLUID participate in a density projection inside SolverXPBD.step(). Cloth, granular particles, rigid contacts, and fluid particles can therefore advance through the same state pair and collision pipeline.

The physical and visual representations are separate

The solver owns particle positions, velocities, density data, and optional diffuse particles. Rendering may smooth positions and fit ellipsoids, but those derived arrays never feed back into the simulation.

Fluid–rigid coupling is two-way

Shape contacts constrain fluid particles while the equal-and-opposite correction is accumulated onto dynamic bodies. The Archimedes screw example uses this coupling to transfer simulated water momentum to a freely rotating wheel.

Nine runnable examples define the supported surface area

The examples exercise containers, pumps, manipulators, mixed fluids, multiple isolated worlds, buoyant objects, wave forcing, and diffuse spray. Each example is registered with an end-state stability check.

2. XPBD fluid simulation

A position-based fluids solve is embedded in each XPBD substep. The implementation uses the same particle arrays, inverse masses, contacts, and state update as Newton’s other XPBD constraints.

Density constraint

For every active fluid particle \(i\), a grouped hash-grid query finds same-world neighbors inside the smoothing radius \(h\). The poly6 kernel estimates density and a unilateral constraint removes compression while allowing a free surface to remain under-dense:

\[ \rho_i = \sum_{j \in \mathcal{N}_i} m_j W_{\mathrm{poly6}}\!\left(\lVert \mathbf{x}_i-\mathbf{x}_j\rVert,h\right), \qquad C_i(\mathbf{x}) = \max\!\left(\frac{\rho_i}{\rho_0}-1,\,0\right). \]

The Jacobi projection computes one multiplier per fluid particle. Here \(w_k\) is inverse mass and \(\varepsilon\) is the relaxation regularizer:

\[ \lambda_i = \frac{-C_i} {\displaystyle \sum_k w_k \left\lVert \nabla_{\mathbf{x}_k} C_i \right\rVert^2+\varepsilon}, \qquad \Delta \mathbf{x}^{\rho}_i = w_i \sum_{j \in \mathcal{N}_i} (\lambda_i+\lambda_j)\frac{m_j}{\rho_0} \nabla W_{\mathrm{spiky}}(\mathbf{x}_i-\mathbf{x}_j,h). \]

fluid_relaxation scales this density correction and a particle-scale bound limits each iteration. A deterministic separation direction resolves coincident samples. The neighbor cap is optional; setting fluid_max_neighbors=0 processes the complete neighborhood.

Substep sequence

  1. Predict. External forces advance particle and body velocities and predict new positions.
  2. Collide. CollisionPipeline generates particle–shape contacts. Texture SDFs provide smooth confinement for mesh containers and pump surfaces.
  3. Project. Each XPBD iteration solves fluid density, particle contacts, rigid contacts, joints, and deformable constraints. Fluid contact corrections can transfer momentum to dynamic bodies.
  4. Update velocity. Corrected positions determine particle velocity; a kernel-weighted neighborhood blend applies viscosity, while optional vorticity confinement restores rotational detail.
  5. Emit diffuse particles. When enabled, low-density crest and trapped-air candidates spawn finite-lifetime foam and spray that advect with the resolved velocity field.

Cohesion and viscosity

Free-surface attraction is a bounded Akinci-style pair bias rather than the negative branch of the density constraint. This keeps isolated particles stable while still forming droplets and strands. Viscosity is applied after projection as a per-substep blend toward the mass- and density-weighted neighborhood velocity:

\[ \mathbf{v}^{\,\prime}_i = \mathbf{v}_i + \nu\left( \frac{\displaystyle\sum_{j \in \mathcal{N}_i} \frac{m_j}{\rho_j}W_{ij}\mathbf{v}_j} {\displaystyle\sum_{j \in \mathcal{N}_i} \frac{m_j}{\rho_j}W_{ij}} - \mathbf{v}_i\right), \qquad 0 \le \nu \le 1. \]

3. Fluid surface rendering

ViewerGL derives a smooth surface from the solved particle state. All smoothing and anisotropy are visual outputs; they do not alter density, contacts, or the next simulation step.

Render-particle fit

SolverXPBD.update_render_particles() computes a weighted neighborhood center and blends it with the simulated position. A regularized covariance fit supplies three clamped ellipsoid axes, allowing sheets and strands to connect without inflating the simulation radius:

\[ \bar{\mathbf{x}}_i = \frac{\sum_j \omega_{ij}\mathbf{x}_j}{\sum_j \omega_{ij}}, \qquad \widetilde{\mathbf{x}}_i = (1-s)\mathbf{x}_i+s\bar{\mathbf{x}}_i, \qquad \boldsymbol{\Sigma}_i = \frac{\sum_j \omega_{ij} (\mathbf{x}_j-\bar{\mathbf{x}}_i) (\mathbf{x}_j-\bar{\mathbf{x}}_i)^{\mathsf T}} {\sum_j \omega_{ij}}. \]
  1. DepthRay-cast anisotropic ellipsoids into linear eye-space depth.
  2. ThicknessAccumulate depth-tested optical thickness from enlarged splats.
  3. SmoothApply a separable depth-aware bilateral filter at framebuffer resolution.
  4. ShadeRecover normals from depth and evaluate refraction, reflection, shadow, and specular light.
  5. DiffuseComposite optional velocity-stretched foam and spray against the fluid depth.

Optical model

Thickness controls Beer–Lambert transmission, while Schlick’s approximation blends the transmitted scene with reflection at grazing angles. The public material parameters expose absorption, index of refraction, normal-incidence reflectance, and specular response:

\[ \mathbf{T}(\tau)=\exp(-\boldsymbol{\sigma}_a\tau), \qquad F(\theta)=F_0+(1-F_0)(1-\cos\theta)^5. \]

Depth, thickness, and bilateral-filter targets use the native framebuffer size. This preserves thin sheets and curved silhouettes in the 1280×720 recordings below.

4. User-level API

Fluid simulation is configured on SolverXPBD; surface generation and optical material settings are explicit rendering calls.

Listing 1. Mark particles as fluid and advance XPBD 
fluid_flags = int(
    newton.ParticleFlags.ACTIVE
    | newton.ParticleFlags.FLUID
)
builder.add_particles(
    pos=positions,
    vel=velocities,
    mass=masses,
    radius=radius,
    flags=fluid_flags,
)

model = builder.finalize()
solver = newton.solvers.SolverXPBD(
    model,
    iterations=2,
    fluid_rest_distance=spacing,
    fluid_rest_density=1000.0,
    fluid_cohesion=0.5,
    fluid_viscosity=0.2,
    fluid_relaxation=0.6,
    fluid_max_neighbors=128,
)

collision = newton.CollisionPipeline(model)
contacts = collision.contacts()

for _ in range(substeps):
    state.clear_forces()
    collision.collide(state, contacts)
    solver.step(state, next_state, None, contacts, dt)
    state, next_state = next_state, state
Listing 2. Build and shade the render surface 
solver.update_render_particles(
    state,
    smoothing=0.6,
    anisotropy_scale=1.0,
)

viewer.log_fluid(
    "/model/fluid",
    solver.render_positions,
    radii=model.particle_max_radius,
    radius_scale=1.8,
    color=(0.113, 0.425, 0.55, 0.8),
    absorption=(0.20, 0.08, 0.03),
    ior=1.333,
    blur_radius_world=0.035,
    anisotropy=solver.render_anisotropy,
    anisotropy_secondary=solver.render_anisotropy_secondary,
    anisotropy_tertiary=solver.render_anisotropy_tertiary,
)
Principal fluid controls on SolverXPBD.
ArgumentRole
fluid_rest_distanceParticle spacing at which a regular lattice is calibrated to rest density.
fluid_smoothing_lengthNeighbor-kernel support radius; defaults to 1.8 times rest distance.
fluid_rest_densityTarget density in kg/m³; may be inferred from particle mass and rest spacing.
fluid_cohesionBounded free-surface attraction in the range [0, 1].
fluid_viscosityPer-substep blend toward neighborhood velocity in the range [0, 1].
fluid_relaxationScale applied to each density-projection correction.
fluid_max_neighborsOptional work bound for over-compressed neighborhoods; zero disables the cap.
max_diffuse_particlesCapacity for visual foam and spray; zero disables diffuse emission.

5. Example simulations

Six runnable XPBD examples are shown at 1280×720 and 30 fps. Every clip contains ten seconds of simulation using the branch’s native ViewerGL fluid surface.

Figure 1. Archimedes screw, 130,001 particles. Fluid–rigid contact carries water upslope through the helical channel; discharged water transfers momentum to a freely rotating paddle wheel before returning to the reservoir.
Figure 2. Dam break, 99,960 particles. A collapsing column impacts two rigid obstacles, testing incompressibility, free-surface cohesion, splashing, and contact response.
Figure 3. Interactive tank, 100,300 particles. Rigid bodies with different densities float, settle, and sink under two-way particle contact.
Figure 4. Three-fluid tank, 94,518 particles. Three independently configured particle groups share one container and interact with buoyant rigid bodies.
Figure 5. Wave pool, 100,491 particles. A kinematic paddle drives sustained waves across a sloped beach while mixed-density bodies test buoyancy and contact stability.
Figure 6. Robot cup transfer, 100,548 particles. A Franka arm carries an SDF-confined fluid through an 8.8-second manipulation trajectory without losing containment.

6. Performance

Measurements use Warp 1.17, CUDA Toolkit 12.9, driver 13.2, and an RTX 4090. Solver trials exclude ViewerGL; capture trials include two 60 Hz simulation steps, rendering, readback, and H.264 encoding per 30 Hz output frame.

Table 1. Warm XPBD dam-break throughput, approximately 100,000 particles, 500 requested frames.
TrialMeasured framesElapsedThroughput
14975.69 s87.4 fps
24975.51 s90.2 fps
34975.46 s91.0 fps
Median4975.51 s90.2 fps
Table 2. End-to-end 1280×720 report capture throughput.
SceneParticlesOutput framesWall timeOutput fps
Archimedes screw130,00130011.70 s25.7
Dam break99,96030010.84 s27.7
Interactive tank100,3003008.72 s34.4
Three-fluid tank94,5183008.65 s34.7
Wave pool100,49130023.69 s12.7
Robot cup transfer100,548300151.57 s2.0

The cup-transfer result is dominated by twelve inverse-kinematics iterations on each 60 Hz simulation frame and should not be interpreted as raw fluid-solver throughput. JSON sidecars beside each video retain the exact capture metadata.

7. Correctness and stability

Focused coverage spans density projection, fluid material terms, coupling, world isolation, surface generation, diffuse particles, and end-to-end examples.

Table 3. Validation performed on the integrated branch.
AreaResultCoverage
XPBD fluid solver61 pass · 1 expected skipCPU and CUDA density, cohesion, viscosity, diffuse particles, SDF confinement, world isolation, capture, reordering, momentum, and non-finite recovery.
Fluid surface API5 passRender-particle fitting, material grouping, lifecycle, visibility, and point fallback.
Fluid examples9 passNine CUDA examples with reduced test particle counts and their built-in physical assertions.
Archimedes screw240 frames passFinite particle and body state, water lift, passive wheel rotation, basin containment, and revolute-bearing integrity.
Report captures6 passSix H.264 simulations at 1280×720, 30 fps, and ten seconds each.

Representative commands

uv run --extra dev -m newton.tests -k test_solver_xpbd_fluid
uv run --extra dev -m newton.tests -k test_viewer_fluid
uv run --extra dev -m newton.tests -k test_viewer_layers
uv run --extra dev -m newton.tests -k test_shape_colors
uv run --extra dev -m newton.tests -k test_viewer_picking
uv run --extra dev -m newton.tests -k TestFluidExamples
uv run -m newton.examples fluid_xpbd_archimedes_screw --viewer gl --headless --num-frames 10
uvx pre-commit run -a

8. Conclusion

The branch adds a position-based fluid model directly to SolverXPBD, including density projection, cohesion, viscosity, vorticity, diffuse particles, multi-world isolation, and two-way rigid coupling. ViewerGL converts the solved particles into a native-resolution anisotropic surface through an explicit public API. Nine examples—including the Archimedes screw pump—exercise the implementation, while the measured dam break sustains a median 90.2 simulation fps at approximately 100,000 particles.