Velocity-Based Monte Carlo Fluids
Ryusuke Sugimoto Christopher Batty Toshiya Hachisuka



<ccs2012>
   <concept>
       <concept_id>10002950.10003714.10003727.10003729</concept_id>
       <concept_desc>Mathematics of computing Partial differential equations</concept_desc>
       <concept_significance>300</concept_significance>
       </concept>
   <concept>
       <concept_id>10002950.10003714.10003738</concept_id>
       <concept_desc>Mathematics of computing Integral equations</concept_desc>
       <concept_significance>300</concept_significance>
       </concept>
   <concept>
       <concept_id>10010147.10010371.10010352.10010379</concept_id>
       <concept_desc>Computing methodologies Physical simulation</concept_desc>
       <concept_significance>500</concept_significance>
       </concept>
   <concept>
       <concept_id>10010147.10010371.10010372.10010374</concept_id>
       <concept_desc>Computing methodologies Ray tracing</concept_desc>
       <concept_significance>100</concept_significance>
       </concept>
 </ccs2012>


[500]Computing methodologies Physical simulation
[300]Mathematics of computing Partial differential equations
[300]Mathematics of computing Integral equations
[100]Computing methodologies Ray tracing




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 <cit.>.
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 RiouxLavoie2022:Fluids 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.=-1



    [b]0.52
        
    < g r a p h i c s >

        
    < g r a p h i c s >

        
    < g r a p h i c s >

        
    < g r a p h i c s >

        
    < g r a p h i c s >



        

        
    < g r a p h i c s >

        
    < g r a p h i c s >

        
    < g r a p h i c s >

        
    < g r a p h i c s >

        
    < g r a p h i c s >



        

        
    < g r a p h i c s >

        
    < g r a p h i c s >

        
    < g r a p h i c s >

        
    < g r a p h i c s >

        
    < g r a p h i c s >



            
    
    [b]0.47
        
    < g r a p h i c s >

        
    < g r a p h i c s >

        
    < g r a p h i c s >

        

        

        
        
   
    
    Our velocity-based Monte Carlo fluid solver can naturally handle scenes for which the existing vorticity-based Monte Carlo method RiouxLavoie2022:Fluids yields incorrect results. Our solver allows the red and blue smoke densities associated with a pair of vortices to flow between two obstacles (left, top), similar to the conventional non-Monte Carlo velocity-based method <cit.> (left, middle), while the vorticity-based method produces an incorrect result, in which smoke deviates around the outside of the obstacles (left, bottom). Our solver also readily supports commonly available features in traditional velocity-based Eulerian grid fluid solvers, such as buoyancy effects: we simulate a smoke plume rising from a complex-shaped source (right, top) and past a sphere-shaped or bunny-shaped obstacle (right, bottom), where the motion is driven solely by buoyancy forces.
    


Laguerre-Gaussian light induction of orbital currents and Kapitza stabilization in superconducting circuits 
    Alexander Balatsky
    April 29, 2024
============================================================================================================





§ INTRODUCTION

Researchers in computer graphics have recently revisited Monte Carlo PDE solvers <cit.> due to their attractive properties. 
<cit.> showed that geometry processing tasks can reap benefits from Monte Carlo solvers, such as flexibility with respect to geometric boundary representations, robustness to noise in the input geometry, support for pointwise estimation of the solution, and trivial parallelization.
Monte Carlo methods have already been used in light transport simulation for decades since the work of <cit.> due to these same properties.=-1

Looking beyond light transport and geometry processing, RiouxLavoie2022:Fluids developed the first Monte Carlo method to perform smoke simulation based on the Navier-Stokes equations. 
They discussed theoretical and practical aspects of Monte Carlo-based fluid simulation; most advantages in the geometry processing context carry over directly to fluid simulation applications. 
However, their method relies on a vorticity-based formulation, whereas velocity-based formulations are more common in fluid animation.

Since the introduction of velocity-based Eulerian grid fluid simulation techniques to graphics <cit.>, researchers have improved on many aspects of the velocity-based formulation, adding features such as 
buoyancy modeling for smoke <cit.> and velocity divergence control for expansion/contraction and artistic effects <cit.>, as well as 
PIC/FLIP solvers <cit.> and advection-reflection solvers <cit.> for reducing numerical dissipation.
These and many other improvements cannot be straightforwardly incorporated into the Monte Carlo fluid solver of RiouxLavoie2022:Fluids because of its vorticity-based formulation. Furthermore, <cit.> recently showed that existing vorticity-based methods can fail to simulate harmonic velocity fields, leading to incorrect results when the fluid domain is not simply connected (e.g., is disjoint or has holes). Conveniently, however, velocity-based methods can capture these physics without any change.

We therefore introduce a 
Monte Carlo fluid solver that relies on a velocity-based formulation. 
Similarly to the vorticity-based method RiouxLavoie2022:Fluids, our method does not require the boundaries of solid objects to be explicitly discretized; for example, there is no need for a cut-cell method <cit.> or a conforming mesh <cit.> as in traditional solvers, and the only requirement on the geometry is that it support ray intersection queries.
We use the recently introduced walk-on-boundary method <cit.>, which is a Monte Carlo ray tracing solver for PDEs, to perform pressure projection in our solver, and thus, our solver can be trivially accelerated by GPU ray tracing.
Our method can additionally incorporate all of the aforementioned techniques developed for velocity-based simulations: buoyancy, divergence control, PIC/FLIP, and advection-reflection. 

To solve the Navier-Stokes equations, we apply operator splitting <cit.>, and the core of our solver depends on reformulating the projection and diffusion steps as integration problems. In the absence of obstacles, we derive a non-recursive integral for each step and apply Monte Carlo integration. The projection step requires a careful evaluation of the volume integral term arising from the Poisson equation's source term, and we describe in detail how to transform it into a form that is more amenable to numerical evaluation. Domains with boundaries require solving a boundary integral equation, which we treat using the walk-on-boundary method <cit.>. The diffusion equation for the diffusion step is a time-dependent equation, and we apply a time-dependent walk-on-boundary method <cit.> for problems with boundaries; we are the first to apply it to a computer graphics application. After constructing Monte Carlo solvers for each substep, we combine them to design a Navier-Stokes solver with flexible options for how we store the intermediate velocity field in a discretized cache.
To summarize, our contributions include:

    
  * a velocity-based Monte Carlo fluid simulator that uses operator splitting,
    
  * an integral formulation for the projection step with careful Poisson source term handling,
    
  * extension of the projection step to problems with boundaries using the walk-on-boundary method,
    
  * application of the diffusion walk-on-boundary method for the diffusion step, 
    
  * adaptations of more advanced velocity-based solver techniques to our Monte Carlo formulation.





§ METHOD

Our method aims to numerically solve for a velocity field that satisfies the incompressible Navier-Stokes equations with constant density and viscosity:

    ∂/∂ t = - (·)   - 1/ +  + ,
     = 0,

where  is velocity,  is pressure, t is time,  is acceleration due to external forces,  is kinematic viscosity, and  is density.




  
Vorticity-Based Monte Carlo


Let us recap the method by RiouxLavoie2022:Fluids to highlight the challenges in developing a Monte Carlo solver for the velocity-based formulation in <ref>.
Their method uses the vorticity representation of the fluid. For instance, for inviscid fluid in 2D, the vorticity-based formulation simplifies the Navier-Stokes equations to the vorticity transport equation, D/Dt=0, where D/Dt is the material derivative.
We can write the velocity field used to advect the vorticity field as an integral, and combined with a semi-Lagrangian advection gives the equation they used to advance the vorticity field in each time step in the absence of obstacles:

    (, t) ≈( - Δ t ∫_ℝ^2(, ) ×(, t-Δ t) , t-Δ t),

where (, ) is the Biot-Savart kernel. They applied Monte Carlo integration to evaluate the integral.

The vorticity-based formulation offered a straightforward application of the Monte Carlo method. It is nontrivial to design a similar scheme for the velocity-based formulation in <ref> because the Navier-Stokes equations, in their velocity form, couple the velocity and the pressure variables in an intricate way: the pressure serves as a Lagrange multiplier that makes the velocity field incompressible. We address this challenge by adopting operator splitting. 



  
Operator splitting with pointwise estimators The basic formulation of our method follows the standard operator splitting framework <cit.>. We solve <ref> by taking discrete time steps with step size Δ t. For each time step, we decompose the Navier-Stokes equations into four substeps:

    
  * advection: ∂/∂ t = - (·),
    
  * external force integration: ∂/∂ t =,
    
  * diffusion: ∂/∂ t =, and
    
  * projection: ∂/∂ t = - 1/ such that =0.

In contrast to traditional discretization-based solvers, we develop a pointwise (Monte Carlo) estimator for each substep. We design a solver for each substep that can estimate the velocity at any spatial point we are interested in, assuming we can likewise query the velocity estimates from the previous steps at any spatial point.

Basing our approach on this pointwise formulation offers unique advantages.
Our formulation is agnostic to the underlying discretization of the velocity field, which allows for flexible integration of velocity-based techniques used by traditional discretization-based solvers, including the common Stable Fluids-style <cit.> grid-based methods and grid-particle hybrid PIC/FLIP methods, into our Monte Carlo formulation.
While the vorticity-based Monte Carlo method RiouxLavoie2022:Fluids has similar advantages of having flexible discretization options, their method cannot benefit from common velocity-based techniques.
Moreover, we will show in <ref> that our formulation allows us to remove some of the excessive grid interpolation errors that are otherwise introduced in traditional solvers.

Below, we formulate a pointwise estimate of the solution to each substep.
Later, in <ref>, we will describe how we employ these estimators in a few different implementations, some using uniform grid velocity storage and others using hybrid grid-particle velocity storage.
To simplify notation, we use _0 to _4 to indicate the velocity field at each substep within a single time step. The initial velocity field for the current time step _0 is defined as the output velocity field _4 from the last time step. The advection step takes _0 as its input and outputs _1, and so on.



 §.§ Advection

We employ semi-Lagrangian advection, which estimates the new velocity at a point by tracing trajectories backward in time through the velocity field. For example, a forward Euler time discretization for the velocity at point  yields the update rule

    _1() ←_0( -  _0()).

This semi-Lagrangian update satisfies our pointwise evaluation requirement: we can estimate the velocity after the step by evaluating the input velocities at a few points. While the time-discretized equation above is identical to what appears in grid-based fluid solvers, we emphasize that the point where we evaluate the advected velocity, , does not necessarily need to be aligned with discretized locations (e.g., grid points), and when we evaluate the pre-advection velocity at points  and -  _0(), we are not restricted to evaluation by interpolation of some grid data. Thus, the pointwise form offers flexibility that is unavailable with classical solvers.
In practice, rather than forward Euler, we use a higher-order time integration scheme (third-order Runge-Kutta) to improve the accuracy of the traced trajectories, and the adaptations of error-correcting schemes such as MacCormack <cit.> and BFECC <cit.> are also trivial; all of these schemes likewise support pointwise evaluation. 



 §.§ External Force Integration

In the presence of external forces leading to acceleration, such as a buoyancy force, we update the post-advection velocity using a forward Euler discretization:

    _2() ←_1() + Δ t ().

This expression is again a pointwise update of the velocity field, where  is not necessarily at a discretized location.



 §.§ Projection

As the diffusion step is necessary only for viscous fluids and we can simply let _3←_2 for inviscid fluids, we first discuss the projection step. The objective of the projection step is to find the pressure field that projects out the divergent velocity mode of the input field. The post-projection velocity at  is given by

    _4() ←_3() - (),

where  satisfies the Poisson equation

    () = _3(),

and = /. Since only  is required for the update in <ref>, the task is reduced to estimating the gradient of the solution of the Poisson equation; the value of  itself is not needed. 
In each subsection below, we discuss the integral formulation followed by its evaluation with Monte Carlo integration.



  §.§.§ Without solid boundaries

To simplify the problem, we will first consider a case without any solid boundaries (i.e., the domain is unbounded with no obstacles inside it).
Under this assumption, it is well-known that we can write the solution to the Poisson equation using the fundamental solution  of the Laplace operator.  
The fundamental solution  satisfies (, ) + δ( - ) = 0 in the unbounded domain, where δ is the Dirac delta function.
In 2D, (, ) =  -1/2πlog r, and in 3D, (, ) = 1/4π r, where r = ‖ - ‖_2.
We can then write the solution to <ref> as a convolution of the source term of the Poisson equation and the fundamental solution:

    () = -∫_(, ) _·_3(),

where dimension =2, 3. We used the subscript  to indicate that we perform differentiation with respect to , and we will use similar notation throughout this section.
By applying the Laplacian to <ref>, one can confirm that  satisfies <ref>.
To evaluate the gradient of , we take the gradient of <ref>:

    _() = -∫__(, ) _·_3().

This integral is still not suitable for numerical computation for our purposes, as we describe below. We therefore extend beyond the prior work by proposing a further transformation.

Since we typically consider storing information and performing volume integrals over a bounded simulation domain, we replace the (infinite) integral domain  in <ref> with a bounded simulation domain  by assuming that the velocity divergence ·_3 is zero outside the simulation domain.

Attempting to use <ref> in this form for Monte Carlo integration still requires explicit evaluation of velocity divergence inside the domain, whereas we desire a solver that takes only a pointwise-evaluated velocity field as input.
There are a few possible approaches to obtain the required divergence. First, we could apply finite differences to the velocity field by accepting some errors. 
Second, we could differentiate the substep that precedes projection so that it outputs the necessary velocity divergence, in addition to the velocity itself. 
We propose instead a velocity-only design that we believe fits better in our Monte Carlo framework.

To eliminate the dependency on velocity divergence in <ref>, we use  the identity _ = -_ and apply integration by parts:

    _()     = ∫_{_(, )}_·_3()
       
           =   -∫_𝐇(, ) _3() 
       - ∫_{_(, )}()·_3() ,

where  is the outward unit normal. 
The Hessian of the fundamental solution denoted 𝐇, is given by

    𝐇(, )    = 𝐒(, ) - δ()/,    𝐒(, )    =  1/|∂ B| r^+2( - r^2),

where |∂ B| is the surface area of a unit (-1)-sphere, =  -, r=‖‖_2 and  is the identity matrix. While we use the notation above for ease of understanding, technically, it needs to be understood in the sense of the generalized (distributional) derivative, and interested readers can refer to the discussion by <cit.> and <cit.> for further details. Substituting <ref> into <ref>, we get

    _() =    -∫_𝐒(, )  _3()  + 1/_3()   
       - ∫_{_(, )}()·_3()
    
     =     -∫_𝐒(, ) {_3() - _3()}
       - ∫_{_(, )}()·{_3()-_3()}.

Again, technically, the domain of the first integral should exclude an infinitesimal ball around  <cit.>, but we use this simple notation for readability.
In <ref>, to get the final expression, we replaced the original velocity field with a velocity field that is globally shifted by the constant velocity at point , _3(). The computed pressure gradient remains unchanged because the divergence of the shifted velocity field is the same as the original one. This global shift cancels the zeroth order term of the Taylor expansion of 𝐮_3(𝐲) around 𝐱, and the remaining terms of the Taylor expansion multiplied by the function 𝐒 will have 
a singularity of O(1/r^-1) instead of the original O(1/r^) as r→ 0. This lower-order singularity can be handled by an appropriate importance-sampling strategy, as we discuss next.






  
Monte Carlo Estimation
Now that we have an integral representation of the pressure gradient, we can define a Monte Carlo estimator for the pressure gradient. We sample N_V points ^i inside the simulation domain  and N_A points ^j on the simulation domain boundary  using sampling strategies with probability density functions (PDFs) _V and _A, respectively, to get

    _() ≈E_V() + E_A()

where

    E_V()   = -1/N_V∑_i=1^N_V𝐒(, ^i)/_V(^i|){_3(^i) - _3()}
    E_A()   = - 1/N_A∑_j=1^N_A_(, ^j)/_A(^j|) (^j)·{_3(^j)-_3()}.

This formulation is an unbiased estimator for the pressure gradient if the probability that we sample any point with a nonzero contribution is nonzero. If the integrand is non-negative everywhere, importance sampling according to the PDF that is roughly proportional to the integrand can decrease the variance of the estimator. While our integrands may contain negative values, we follow this idea to design our sampling strategy. 
In our implementation, for <ref>, to handle the singular integral, we draw samples such that PDF _V is proportional to 1/r^(-1) inside of the smallest ball around  that fully contains the entire simulation domain, and set zero contributions from the samples outside the simulation domain. This approach is equivalent to evaluating the integral we get by extending the original integral domain to the ball and putting zero integrands in the extended part. We also use antithetic sampling here: in addition to a point ^i sampled this way, we always sample the symmetric point with respect to , 2 - ^i to further reduce variance for smooth velocity fields.
For <ref>, we uniformly sample points on the simulation boundary. Properly sampling the positive and negative contributions separately may further reduce the variance of the estimators <cit.>, but we leave it as future work. 
Equation <ref> only requires the ability to evaluate the velocity at the position of sample points to perform projection, in contrast to traditional solvers that typically require a globally coupled linear system solve.




  §.§.§ With solid boundaries.
 
Next, we consider the case when solid object boundaries are involved (<ref>).
We still solve <ref> for , but with free slip boundary conditions, 

    ∂/∂ = · (_3 - )

on solid boundaries, where  is the unit outward normal from the fluid domain Ω, ∂/∂=· is a normal derivative, and  is the solid velocity. If the simulation domain is bounded, this problem is an interior Neumann problem for the Poisson equation, and otherwise, it is an exterior Neumann problem; we need a Monte Carlo solver that is capable of handling these problems.

The basic walk-on-spheres method <cit.> can only handle Dirichlet problems, and its applicability to unbounded domain problems leaves some concerns about the termination of paths without additional bias.  RiouxLavoie2022:Fluids arbitrarily terminated their walk-on-spheres paths after a few steps in their simulation, leaving some additional bias. <cit.> suggested inverting the unbounded domain into a bounded one using the Kelvin transform so that paths always terminate, but using this approach would require very different treatment for bounded and unbounded domains.

To solve Neumann problems, the walk-on-stars method sawhney2023walk extended the walk-on-spheres method with carefully-designed recursion relationships and tailored sampling techniques based on the extensions by <cit.> and <cit.>.
For a Neumann problem in an unbounded domain, however, the Kelvin transform <cit.> is still required; it transforms a Neumann problem into a Robin problem, to which the application of walk-on-stars has not yet been attempted.
Among others, the walk-on-boundary method <cit.>, introduced recently to graphics by <cit.>, can be applied to Neumann boundary problems in both bounded and unbounded domains under a unified framework. The method builds upon a boundary integral equation analogous to the rendering equation for light transport simulation <cit.>, and uses ray-tracing to solve the boundary integral equation, and we use this method in our solver. 




Following <cit.>, we can write the gradient of the solution to the Poisson equation of <ref> with the boundary condition of <ref> using the single layer potential formulation and define an integral equation for the unknown density function.
For our application, these equations contain the volume integral terms arising from the source term in the Poisson equation of <ref> specific to our problem. 
Similar to <ref>, we transform these volume integral terms as follows.
If the domain  is unbounded, we suggest replacing it with a bounded simulation domain , assuming zero velocity divergence outside of the simulation domain; otherwise, we let =. With this definition, we have  ⊆ for the boundaries, and we distinguish these two sets of boundaries in the following equations. We also remove the explicit dependencies on the velocity divergence by integration by parts similarly to <ref> to get

    _() =
           ∫_{_(, )}[() - ()·{_3() - _3()}]
       -∫_𝐒(, ) {_3() - _3()}
       -∫_\{_(, )}()·{_3() - _3()}

for ∈ and 

    () =
             -∫_ 2(, ) [() - ()·{_3() - _3()}]
       +∫_ 2()𝐒(, ) {_3() - _3()}
       +∫_\ 2(, ) ()·{_3() - _3()}
        + 2()·{_3() - ()}

for ∈.
In both <ref> and <ref>, we have the boundary integral terms that involve the unknown density function  on the solid boundaries , and all the other integrals involve only known quantities.
Note <ref> is a generalization of the strategy described in <ref> as we can recover <ref> by dropping the solid boundary integral term.




  
Monte Carlo Estimation
Based on <ref> and <ref>, we get a biased walk-on-boundary Monte Carlo estimator <cit.>.
Compared to the simplest formulation for the Laplace equation <cit.>, we have some additional non-recursive terms in <ref> and <ref> as a result of the transformations discussed above; we must consider how to sample these added terms. We choose to sample the non-recursive contributions on ∂ using ray intersection sampling for importance sampling and estimate the other non-recursive terms similarly to <ref>. To improve efficiency, we additionally employ boundary value caching, which lets us share the walk-on-boundary subpaths among evaluation points akin to the virtual point light (VPL) method <cit.> in rendering. We describe the details of our particular sampling strategy in the supplemental note. Our Monte Carlo method is compatible with a variety of sampling techniques, so further variance reduction is likely possible. 




 §.§ Diffusion

Simulating viscous fluids requires an additional diffusion step before the projection step. The set of equations we solve here is the constant-coefficient diffusion equation, 

    ∂(, s)/∂ s   = (, s)     for    ∈, s ∈ (0, Δ t),
    (, s)    = _s()    for    ∈, s ∈ (0, Δ t),  and
    (, 0)    = _2()    for    ∈,

and the output of this diffusion step is _3() = (, Δ t). 
The time s here represents the time within each diffusion time step. We use the overbar to indicate that the velocity variable  takes time s as defined here, and differs from the other sections.
The solid velocity gives the no-slip boundary condition, and the output from the preceding advection step gives the initial conditions.

When there are no solid boundaries or the time step size is small enough to ignore boundary effects, we can omit the second equation describing the boundary condition. In such a case, evaluating a convolution of the input velocity field with a Gaussian function would give the exact solution to the diffusion equation. 
Though they apply it to the vorticity field rather than the velocity field, RiouxLavoie2022:Fluids use such a convolution-based approach to model diffusion.
This approach fails when either the time step size is relatively large or the boundary effect is important.
Traditional grid-based solvers similarly cannot rely on a simple Gaussian filter in the presence of large time steps or obstacles, and must instead solve a globally coupled linear system for viscosity <cit.>.=-1

To address this problem, we utilize the walk-on-boundary method for the diffusion equation <cit.> to solve <ref>. Similarly to <ref>, the diffusion walk-on-boundary method is a pointwise estimator and is a natural generalization of the case when there are no boundaries. One difference is that the boundary integral equation and the walk-on-boundary steps are now defined in the space-time domain, and we must sample paths from the point where we want to evaluate the solution with time Δ t towards the initial time within the time step, in the negative time direction.

Notably, for the diffusion walk-on-boundary, the recursion depth does not need to be predefined. By sampling the space-time paths using a PDF proportional to the integral kernel of the integral equation for the diffusion equation, we can terminate recursions when the sampled time is negative.
In our experiments, we observed that the diffusion step is cheaper than the projection step, and we did not attempt any algorithmic improvements, such as the reuse of subpaths. However, there remain many choices for sampling and various efficiency improvement strategies should apply here, too.
The supplemental material summarizes the diffusion walk-on-boundary method for readers' convenience.



§ RESULTS




We explain the results of our practical implementation of the method with discretization structures and other possible extensions.
We implemented our method using CUDA for GPU parallelism and the NVIDIA OptiX ray tracing engine <cit.> via the OWL wrapper library <cit.> to accelerate the ray intersection queries.
For all 2D results, we use consistent parameters as listed in the supplemental note except for <ref>.
The figures have time stamps for the corresponding animations in the supplemental video.=-1



  
Caching of velocity fields
While our fundamental formulation is agnostic to a choice of intermediate storage of velocity fields, or absence of it by the expensive recursive application of estimators to the initial time step (see <ref>), our basic implementation uses a simple uniform grid structure to store the velocities at grid nodes at each time step, and we refer to such a strategy a caching strategy in contrast to the fully time-recursive (cacheless) alternative. 
Since we do not evaluate any finite differences on the stored data, we do not need to use a staggered grid <cit.>.
When we query a velocity from the cache, we bilinearly or trilinearly interpolate the values, and for points outside the cached domain, we take the velocity at the nearest cache point. As with classical grid-based solvers, interpolation- and advection-induced errors sometimes cause us to query values from points inside of solid obstacles; for this issue, we assign velocities inside solid obstacles to take velocities from the solid itself, giving no-slip behavior on the boundary.
Alternatively, one could fill the velocities inside solid obstacles with some form of extrapolation to approximate free slip behavior, for example, but we use the first approach throughout this paper.

While we could cache the velocity field after each substep, the pointwise estimation capability of our approach lets us avoid some of the excessive velocity caching between substeps and reduce the associated errors. For example, we need no caching between the advection and the following projection step for inviscid simulation if we query the advected velocity values directly on the fly whenever the projection estimator needs such a pointwise input velocity. Alternatively, we can also design a method without velocity caching between the projection and the next advection so the velocity field is always advected according to the pointwise divergence-free velocity field.
We implemented these two options in addition to the option where we cache the resulting velocity field both after advection and after projection in <ref> (b) to (d). Our experiments showed consistent results among the three options; further investigations are needed to understand when having fewer caching errors can make critical differences.




  
Boundary conditions <ref> shows that our method can correctly simulate the velocity field due to vortex pairs under different boundary conditions. In <ref>, we further show our method's application to a domain with two disjoint obstacles. In such cases, vorticity-based methods typically use an incorrect boundary condition, and the vorticity-based Monte Carlo method RiouxLavoie2022:Fluids fails to simulate the physics correctly. Applying the modification proposed by <cit.> for vorticity-based simulation to a Monte Carlo method is not straightforward and has not been demonstrated. Hence, ours is currently the only Monte Carlo method that can produce correct simulation results in these scenarios. While the relative performance depends heavily on parameter choices, for the specific simulations in <ref>, ours took 101.0s per time step while the vorticity-based method took 87.6s, with a difference of about 15%.



  
Buoyancy
In addition to the velocity itself, we can additionally simulate the advection and diffusion of the temperature field T and the smoke concentration field s in a manner similar to the velocity field to add buoyancy effects (Figures <ref>, <ref>, and <ref>). We use the Boussinesq approximation, which assumes the fluid density variation to be negligible, so that the buoyancy force can be computed using 

where  is the gravitational constant vector, α and β are positive parameters, and _ambient is the ambient temperature constant <cit.>.
Naively incorporating buoyancy into the vorticity-based Monte Carlo formulation by RiouxLavoie2022:Fluids would require the curl of  <cit.>, which would, in turn, require finite difference approximations with some additional errors and ruin the cache-structure-agnostic nature of pointwise estimators.





  
Divergence control
As our formulation is based on the standard pressure solve framework, we can easily add velocity sinks and sources in our simulation for better artistic control (<ref> (b)) or other expansion/contraction effects, by modifying the pressure solve step similar to the grid-based formulation <cit.>.
We add another term to the right-hand side of <ref>, which leads to an additional integral term in our formulation.
Our implementation supports Dirac delta-type sinks and sources, and we sample the location of sinks and sources directly. It should also generalize straightforwardly to more general volumetric sinks and sources with an appropriate sampling technique. This divergence control technique is yet another advantage of adopting a velocity-based formulation, as it is impossible under vorticity formulations.



  
Viscous fluids
We test our solver in the case of viscous fluids using the von Kármán vortex street scenario (<ref> (e) to (g)). For this simulation, as we increase the viscosity coefficient, or equivalently, as we decrease the Reynolds number, the flow is expected to become less turbulent. We observe that our simulation qualitatively produces the expected type of flow for each Reynolds number configuration we tested. For all the results other than those in <ref> (e) to (g), we disabled the diffusion step. 



  
Advection-reflection solver
In contrast to the case of high viscosity, we may want to introduce less artificial viscosity to the solver. It is known that even without the diffusion step, grid-based solvers suffer from artificial viscosity in the solver because of the dissipation of energy in the advection and projection steps. To address this problem, <cit.> recently introduced an advection-reflection solver, which still uses the same building blocks of advection and projection steps as the standard advection-projection solvers use, but in a more careful way, to reduce artificial viscosity with the same number of projection steps. We adapted its second-order variant <cit.> to our Monte Carlo solver (<ref> (h)), and we observe our method similarly benefits from the use of the advection-reflection method, producing more turbulent animation.



  
Monte Carlo PIC/FLIP
The numerical diffusion effect we have discussed above is partially due to the semi-Lagrangian advection scheme with some grid data interpolation that our method and many traditional methods employ.
To address this problem in the context of classical grid-based solvers, the PIC/FLIP methods <cit.> have been popularized.
The idea of PIC/FLIP is that we handle the advection of velocity using particles by (forward) Lagrangian advection, transfer the velocity from particles to grid, perform projection on the grid, and update particle velocities by transferring grid data back to the particles. PIC and FLIP differ in their last grid-to-particle transfer step: PIC transfers the velocity itself whereas FLIP transfers the velocity update. FLIP is expected to produce more turbulent and noisy results. One can also blend the two transfer methods. Our velocity-based Monte Carlo framework can utilize this idea straightforwardly (<ref> (i) to (k)). In our version, at each step, we first perform the standard particle-to-grid transfer. Then, instead of advecting the particles using the projected grid velocities at the particles' locations via interpolation, we evaluate the velocity (or velocity update) from the grid data exactly at the positions required for particle advection. Exploiting the pointwise estimation capability in this way lets us avoid extra interpolation errors for this step.=-1



  
Performance
We measured the timing of our simulations in <ref> using their GPU implementations against a basic single-thread CPU implementation of traditional grid-based fluids on a workstation with an Intel Xeon Silver 4316 CPU and an NVIDIA RTX A5000 GPU. Stable fluids with cut-cell boundary handling <cit.> (a) runs at 0.064 seconds per time step, while the fastest implementation of our solver that caches the velocity field after each substep (b) runs at 7.8 seconds per time step, two orders of magnitude slower. If we disable the VPL with this setup, the runtime increases to 109 seconds per time step. For the other caching choices, the runtime for the one without caching after projection (c) increases to 32.1 seconds per time step because we perform the projection at all points used in the RK3 advection, and the one without caching after advection (d) runs at 10.1 seconds per time step due to the increased number of evaluation points for the advection. In general, the walk-on-boundary method for handling obstacles requires a large number of paths, and further efficiency improvements would be desirable to improve the method's practicality. The addition of diffusion steps in (e) to (g) resulted in the increase of runtime by 48% (e) to 84% (g) compared to their baseline counterpart implementation (c).



  
Error analysis
We performed a convergence test for the projection step in the absence of obstacles. The estimator is unbounded in this case, and we can observe in <ref> that the error decreases at the inverse square root rate, as we increase the number of samples to estimate the volume term and the area term, respectively. We also observe that using a simpler formulation or a simpler sampling technique makes the estimator converge slower or completely fail to converge. In <ref>, we compared our simulations with different numbers of samples. We observe that using low sample counts can blow up the simulation, and the noise of the velocity estimate is typically larger around the boundary. As <cit.> discuss, the variance of the walk-on-boundary method is particularly large for interior Neumann problems in non-convex domains, and <ref> shows such a scene.



§ CONCLUSION AND DISCUSSION

We have presented a velocity-based Monte Carlo method for fluid simulation. We developed a pointwise estimator for each substep and demonstrated flexible integration of existing velocity-based techniques into our method.

One could also consider using our velocity-based pointwise estimators to design a cacheless scheme without any spatial discretization whatsoever, similar to the vorticity-based strategy described by RiouxLavoie2022:Fluids. 
As RiouxLavoie2022:Fluids pointed out, this kind of fully recursive strategy offers truly pointwise estimation at the cost of large variance and thus has a high computational cost. While of theoretical value, this approach is impractical today because the estimation cost grows exponentially as the simulation proceeds in time, depending on the specific path sampling strategy.


Our method's computational speed is not yet comparable to traditional methods as a relatively large number of samples is required to produce results with low enough error. Moreover, having boundaries in the scene necessitates the use of an advanced Monte Carlo solver for (Neumann) boundary value problems, with corresponding increased variance, leading to higher costs.


In addition to efficiency improvements, our method will benefit from any future improvements to the walk-on-boundary method. Support for spatially varying diffusion equations under the walk-on-boundary method, analogous to what was proposed for the walk-on-spheres method Sawhney:2022:GFMC, would allow us to support variable fluid densities in our framework.
Applying differentiable rendering techniques similar to <cit.> would enable the use of our solver for inverse problems.
Supporting double-sided boundary conditions by extending the boundary integral formulation would also be a useful extension.

We designed our solver with a velocity-only formulation, but it may be beneficial to revisit this choice as (grid-based) methods utilizing velocity gradient information <cit.> have recently shown promising improvements. Other ideas in the velocity-based fluid literature, such as the method of characteristic mapping <cit.> to lessen advection errors or energy-preserving integrators <cit.> to reduce splitting errors, are also worth investigating in our Monte Carlo context.

Lastly, we believe our method paves the way for a full Monte Carlo simulation of liquids with free surfaces, and we look forward to the continued advancement of Monte Carlo techniques for physical simulations more broadly.


This research was partially funded by NSERC Discovery Grants (RGPIN-2021-02524 & RGPIN-2020-03918), CFI-JELF (Grant 40132), and a grant from Autodesk. 
Ryusuke Sugimoto was partially funded by David R. Cheriton Graduate Scholarship. 
This research was enabled in part by support provided by SHARCNET and the Digital Research Alliance of Canada.
We thank Joel Wretborn, Brooke Dolny, Rikin Gurditta, and Clara Kim for proofreading.
We would like to thank the anonymous reviewers for their constructive evaluations and feedback.











ACM-Reference-Format
















 









§ WALK-ON-BOUNDARY FOR PROJECTION

For <ref> and <ref>, we defined the normal directions as pointing outward from the fluid domain to describe both interior and exterior problems with a single pair of equations, whereas the previous work  <cit.> had two separate pairs of equations for interior and exterior problems. Hence, our equations for exterior problems differ from <cit.> and <cit.> in their signs.



  
Monte Carlo Estimation
Based on <ref> and <ref>, we use a walk-on-boundary Monte Carlo estimator <cit.>.
Note that <ref> contains the unknown density function  itself on the left-hand side of the equation and in one of the integrands on the right-hand side, making it a recursive integral equation similar to the rendering equation <cit.> used in path tracing for light transport simulation rather than a simple non-recursive one as in <ref>. Following <cit.>, we truncate the recursion after M steps and multiply the contribution of the longest path by 0.5 to design a biased estimator.

We can consider various sampling strategies to estimate the solution to the truncated recursive integral equation with Monte Carlo methods. Following <cit.>, we use a forward estimator to solve this in our implementation: to sample a path consisting of multiple boundary points, we randomly sample a boundary point first, take a random walk on randomly sampled boundary points, and finally connect them to each individual evaluation point. 
In the below, for _U(_0), we use uniform sampling to sample a point _0 on the solid boundary . Then, for _R(_m+1|_m), we uniformly randomly sample a direction from a unit hemisphere for a line that goes through the point _m and uniformly randomly sample one of the intersection points between the line and the solid boundary  to sample the next point _m+1. Recursively applying this sampling strategy lets us sample a series of points _1, _2, ..., _M+1 to form a path, given previous points in the path. We use this strategy because the PDF associated with this specific forward sampling strategy is known to be proportional to the integral kernel  we have in <ref>.

Compared to the simplest formulation for the Laplace equation <cit.>, we have some additional non-recursive terms in <ref> and <ref> as a result of the transformations we discussed in the main paper, and we need to consider how to sample such terms, too. We choose to sample the non-recursive contributions on ∂ using the forward estimator as well for the purpose of importance sampling while estimating the other terms similarly to <ref>. To do so, we first define two estimators for length m contributions, μ^1_m() and μ^2_m(), recursively, based on <ref>: we define the base case for length 1 contributions as

    μ^1_1(_1) = 2(_1, _0)/_R(_1|_0) (_0)·{_3(_0)-_3(_1)},


    μ^2_1(_0) = 2(_0) ·{ -E_V(_0) - E_A(_0)+ _3(_0) - (_0)},

where μ^1_1(_1) corresponds to the non-recursive contribution from the first integral in <ref> and μ^1_2(_0)corresponds to the rest of the non-recursive contributions in <ref>.
Then, we define the contributions from longer paths as

    μ^1_m+1(_m+1)   = -2(_m+1, _m)/_R(_m+1|_m)μ^1_m(_m),
    μ^2_m+1(_m)   = -2(_m, _m-1)/_R(_m|_m-1)μ^2_m(_m-1).

Note in particular, when we use the above-mentioned uniform line intersection sampling for _R(|), 

    2(, )/_R(|) = κ(, ) ·sgn(() ·( - )),

where κ(, ) is the number of intersection points that the line that goes through the two points has, excluding point , and sgn is the sign function.
Then, based on <ref>, we can estimate the pressure gradient by taking N_P sample paths in addition to the terms in <ref>:

    _()
    
    =    E_V() + E_A()
       +  1/N_P∑_k=1^N_P[  -_(, _0^k)/2_U(_0^k)(_0^k)·{_3(_0^k) - _3()}. 
       + _(, _M+1^k)/2_U(_0^k)μ^1_M(_M+1^k) +  _(, _M^k)/2_U(_0^k)μ^2_M(_M^k)
       + .∑_m=1^M-1_(, _m+1^k)/_U(_0^k)μ^1_m(_m+1^k + _(, _m^k)/_U(_0^k)μ^2_m(_m^k)].

On the right-hand side, the first line estimates the last two integrals in <ref>, the second line estimates the non-recursive terms in the first integral, the third line estimates the longest path contributions, and the last line estimates the shorter path contributions. Note that while we do not explicitly indicate so, the contributions from paths of any length μ^1_m(_m+1^k) and μ^2_m(_m^k) implicitly depend on the sampled boundary point _0^k due to the sampling strategy we employ.

To use the estimator <ref>, the most naive approach would be to generate sample paths for each individual evaluation point, which has a high computational cost. Instead, we use a boundary value caching strategy similar to the virtual point light method in rendering <cit.>. We first compute the contributions of N_P subpaths up to the point before we connect them to the evaluation points and cache them at M+1 boundary points per path. We then connect all of them to all evaluation points. This significantly increases the effective number of sample paths per evaluation point without increasing the computational cost too much while introducing a correlation of estimates between evaluation points. Even this correlation can be preferable for our application because it guarantees that the contributions from these paths changes smoothly across evaluation points. 

<cit.> apply a similar caching technique to their Neumann problem walk-on-boundary solver based on the same single-layer boundary integral formulation as well, but they use a backward estimator, which works only with a less efficient resampled importance sampling strategy. In contrast, our method can utilize a forward estimator with more efficient line intersection importance sampling. This is at the cost of increasing storage per cache point by a small constant multiplication factor and increasing computation per evaluation point when we sum up the contributions from all cache points.

Similar to <ref>, for the estimation of volume integrals E_V() in <ref> and <ref>, we use the importance sampling strategy so PDF _V is proportional to 1/r^(d-1), and let 𝐒(, ) = 0 for all  outside the simulation domain. Additionally, we need to set 𝐒(, ) = 0 for all  inside of solid obstacles, too. We perform this insidedness test by casting a ray from point  in a random direction and taking a sum of its intersection signs. This strategy can be considered a Monte Carlo estimator for the generalized winding number <cit.>. We use uniform sampling for E_A().


For all 2D examples in the paper except <ref>, we use N_V = N_A = 5·10^5 for direct contributions to each evaluation point. For indirect (path) contributions, we use N_P = 5·10^5 paths with path length 4. For each path, the volume and area term sample counts for the pseudo boundary are N_V=N_A=10.





§ WALK-ON-BOUNDARY FOR DIFFUSION

To get scalar diffusion equations from <ref>, we first absorb the viscosity coefficient ν into the variable to define (, s)=(, s/) and redefine the range of time s accordingly:

    ∂(, s)/∂ s   = (, s)     for    ∈, s ∈ (0, Δ t)
    (, s)    = (, s) = (, s/ν)    for    ∈, s ∈ (0, Δ t), and
    (, 0)    = (, 0) =  (, 0)    for    ∈,

where the asterisk in  indicates that it is a given function.  
This is still a vector-valued equation, but because we only need to deal with simple Dirichlet boundary conditions, we can solve the vector diffusion equation <ref> component-wise. Thus, we describe the scalar diffusion equation solver for one scalar component of 𝐰, w below.

We will summarize the diffusion walk-on-boundary solver from the book by <cit.> here, focusing on its use in our context. For more general cases not described here, such as Neumann problems and nonhomogeneous problems, readers should refer to the book. The idea of the diffusion walk-on-boundary method is very similar to the one for the Poisson equation: we convert the partial differential equation into an integral equation and solve it using a ray-tracing-style solver. However, with the additional time dependency, we need to additionally consider the time variation in the diffusion equation and take walks in the space-time domain.
Note that this solver does not discretize the time domain within each time step, unlike most traditional methods that discretize the time with a finite difference.

First, we define the fundamental solution for the diffusion equation <ref>, also known as the heat kernel,

    (, s; , τ) = (t') (4π t')^-d/2e^-r^2/4t',

where t' = s-τ and (·) is the Heaviside step function.


Using the fundamental solution, we can write the solution to the diffusion equation in the form of the double layer potential and an additional initial condition term for ∈ and s∈ [0, νΔ t]:

    (, s) =    -∫_0^s ∫_(, s; , τ) (, τ) 
       +∫_(, s; , 0 ) (, 0) ,

where (·, ·) is an unknown density function defined for ∈ and s∈ [0, νΔ t].
We can also derive a recursive boundary integral equation for (·, ·) by taking the limit → in <ref>:

    (, s) =   ∫_0^s ∫_ 2(, s; , τ) (, τ) 
       + 2 (, s) - 2∫_(, s; , 0 ) (, 0) .




  
Monte Carlo Estimation We design a Monte Carlo estimator based on <ref> and <ref>.
First, we define an estimator for the initial condition term with N_I samples using PDF _I as

    E_I(, s)= 1/N_I∑_i=1^N_I𝐙(, s;^i, 0)/_I(^i|, s)(^i, 0).

Using this estimator, we can define the estimator for <ref> with N_D path samples using PDF _D as

    (, s) = E_I(, s) +  1/N_D∑_j=1^N_D  -(, s; ^j, τ^j)/_D(^j, τ^j|, s)(^j, τ^j)

and the one for <ref> with 1 sample with PDF _E as

    (, s) =  2(, s; , τ)/_E(, τ|, s)(, τ)   + 2 (, s) -2E_I(, s),

which is a recursive estimator: (, τ) included on the right-hand side should be estimated recursively.
We use a backward estimator and start generating paths for this recursive estimator from each evaluation point (, νΔ t) toward time 0. Note that _D and _E only need to sample points with time τ < s because the integral kernel  is zero for τ≥ s. Intuitively, this is because the solution to the heat equation depends only on previous times.
Using this property, we sample points in the space-time domain so that times for the sampled sequence of points strictly decrease.
We terminate the recursion when the sampled time τ is negative because the original integral domain does not contain the negative part. While the walk-on-boundary method for the Poisson equation is biased, the walk-on-boundary method for the diffusion equation gives an unbiased solution estimate. 

As for the specific sampling strategy, for the initial condition sampling _I(^i|, s), we sample the points by _i← + √(t γ_d/2)ω⃗, where γ_d/2 is a sample drawn from the Gamma distribution with shape parameter d/2 and scale parameter 1, and ω⃗ is a uniformly random direction sampled on a unit sphere. With this sampling strategy, we get

    𝐙(, s;^i, 0)/_I(^i|, s) = 1.
Similar to <ref>, we set 𝐙(, s;^i, 0) = 0 for all  inside of solid obstacles in <ref>.

For _D(, τ|, s) and _E(, τ|, s), we sample  using the uniform line intersection sampling as in <ref> and sample time τ, given t, by τ← s - ‖ - ‖/4γ_d/2 to get

    (, s; , τ)/_D(, τ|, s)   = -κ(, ) ·sgn(() ·( - )),   and
    
     2(, s; , τ)/_E(, τ|, s)   = -κ(, ) ·sgn(() ·( - )).

In the above, the left-hand side scaling factors of the two expressions differ by 2 because the PDF _D and _E differ even though we use the same strategy: point  for _D lies within the domain, whereas point  for _E lies on a boundary. 

To draw samples from the Gamma distribution, we use  for 2D, where α is a uniformly random sample in the range (0, 1), and γ_3/2←γ_1 + ξ^2/2 for 3D, where ξ is a standard normal sample.


For problems without boundaries, we drop the second term in <ref>, and the remaining first term estimates the convolution of the input field and a Gaussian function.

In our implementation, we use N_I=5·10^5 initial conditions samples at each evaluation point, with N_D = 5·10^5 paths with N_I = 10 initial condition samples per bounce. Unlike the walk-on-boundary method for the projection step, we do not share the subpaths among evaluation points as we found that the diffusion walk-on-boundary is typically cheaper, and there was not much need to improve its efficiency relative to the projection.
