Velocity-Based Monte Carlo Fluids

TL;DR

This work presents a velocity-based Monte Carlo solver for incompressible fluid dynamics that overcomes key limitations of prior vorticity-based Monte Carlo methods. By applying operator splitting to the Navier–Stokes equations and developing pointwise Monte Carlo estimators for advection, external forces, diffusion, and projection, the method integrates velocity-based techniques such as buoyancy modeling and PIC/FLIP within a Monte Carlo framework. A central innovation is reformulating the projection and diffusion steps as integral problems and employing walk-on-boundary methods to handle complex boundaries, enabling accurate pressure gradient and diffusion evaluations in both bounded and unbounded domains. While computationally demanding, the approach yields correct physics in scenarios where vorticity-based solvers falter, and it is compatible with GPU-accelerated ray tracing to accelerate boundary interactions. The work lays a foundation for incorporating variable diffusion, differentiable inverse problems, and fully Monte Carlo free-surface liquids in future developments, offering a flexible, boundary-agnostic path toward realistic fluid animation.

Abstract

We present a velocity-based Monte Carlo fluid solver that overcomes the limitations of its existing vorticity-based counterpart. Because the velocity-based formulation is more commonly used in graphics, our Monte Carlo solver can be readily extended with various techniques from the fluid simulation literature. We derive our method by solving the Navier-Stokes equations via operator splitting and designing a pointwise Monte Carlo estimator for each substep. We reformulate the projection and diffusion steps as integration problems based on the recently introduced walk-on-boundary technique [Sugimoto et al. 2023]. We transform the volume integral arising from the source term of the pressure Poisson equation into a form more amenable to practical numerical evaluation. Our resulting velocity-based formulation allows for the proper simulation of scenes that the prior vorticity-based Monte Carlo method [Rioux-Lavoie and Sugimoto et al. 2022] either simulates incorrectly or cannot support. We demonstrate that our method can easily incorporate advancements drawn from conventional non-Monte Carlo methods by showing how one can straightforwardly add buoyancy effects, divergence control capabilities, and numerical dissipation reduction methods, such as advection-reflection and PIC/FLIP methods.

Figures (6)

• Figure 1: Our method can correctly simulate scenes with different boundary types (resolution $256^3$). From the left, we simulate the colored smoke advected with the velocity field induced by vortex pairs moving with no obstacles (left), in an unbounded domain with a moving square obstacle (middle), and in a domain bounded by an obstacle (right). For each simulation, we show the time evolution from the left to the right, but the time stamps differ for each setup. 1m17s
• Figure 2: Variants of our method with different configurations tested on the Kármán vortex street scene with resolution $128\times256$. Except for (a), all results are generated with our Monte Carlo method. (b) to (d) compare the caching alternatives, (e) to (g) demonstrate simulations with various Reynolds numbers $\textit{Re}$, (h) shows the application of the advection-reflection method to yield reduced numerical dissipation, and (i) to (k) show the application of PIC/FLIP methods. Note that (i) to (k) have no density diffusion because they employ density advection by particles. Given their agreement with the result of the traditional grid-based method (a), we consider all variants of our method to produce qualitatively reasonable results. 1m56s
• Figure 3: Buoyancy and divergence control (resolution $256^2$). We simulate a hot smoke plume rising towards a blue circular obstacle in 2D using the Boussinesq buoyancy model (left). We can also add a velocity sink to the scene, indicated with a red dot, causing the smoke to be sucked into the sink (right). 3m49s
• Figure 4: Buoyancy simulation in 3D. We simulate a smoke plume rising toward a bunny-shaped obstacle (top, resolution $128^3$) and a sphere obstacle (middle, resolution $128^3$) in 3D using the Boussinesq buoyancy model. We can also incorporate a more complicated smoke source shape as well (bottom, resolution $288\times128\times128$). We rendered the images with the principled volume shader in Blender blender2023. 0m15s
• Figure 5: Projection error without solid boundaries. Our projection step projects out the divergent velocity mode from the velocity field, converting the top left field into the one on the bottom left (resolution $256^2$). We visualize the velocity field using arrows and we show its magnitude using the darkness of the red color as well. We measure the root mean squared errors measured against the analytical solution and plot them by altering the number of volume samples $N_V$ with a fixed number of area samples $N_A=10^7$ in the middle and by altering $N_A$ with $N_V=10^7$ on the right. We observe the inverse square root convergence rate in both cases with our formulation with the default sampling strategy (a) as described in \ref{['sec:projection_wo_boundaries']}, while the bias due to the volume term eventually dominates the error on the plot for $N_A$ as we increase $N_A$. We also tested the estimator by (separately) replacing the volume term importance sampling according to the PDF proportional to $1/r^{(d-1)}$ with a uniform sampling (b), and disabling the global velocity shift described in \ref{['eq:pressure_grad_final']} (c). Both of these alternatives either converge slower than the proposed method or fail to converge.