Kinfung Chu
Yoshifumi Kitamura
Abstract
We present a framework that captures sparse Lagrangian flow information from a volume of real liquid and reconstructs its detailed kinematic information in real time. Our framework can perform flow reconstruction even when the liquid is disturbed by an object of unknown movement and shape. Through a large dataset of liquid moving under external disturbance, an agent is trained using reinforcement learning to reproduce the target flow kinematics with only the captured sparse information as inputs while remaining oblivious to the movement and the shape of the disturbance sources. To ensure that the underlying simulation model faithfully obeys physical reality, we also optimize the viscosity parameters in Smoothed Particle Hydrodynamics (SPH) using classical fluid dynamics knowledge and gradient-based optimization. By quantitatively comparing the reconstruction results against real-world and simulated ground truth, we verified that our reconstruction method is resilient to different agitation patterns.
1 INTRODUCTION
Advances in fluid simulation and GPU parallel computing enable the generation of detailed fluid animation even in real time. Smoothed Particle Hydrodynamics (SPH) is a simulation method that is receiving much attention in the computer graphics community. Although numerous pressure solvers [Bender and Koschier 2017; Ihmsen et al. 2014a] and boundary models [Bender et al. 2019; Akinci et al. 2012] have been proposed that allow robust and stable simulation with complex boundaries, room exists for analysis regarding their physical accuracy. Such analysis is crucial in a reconstruction context because, for any physically based reconstruction framework, an accurate model of the governing physical equations is generally required. In this article, we first approach the flow reconstruction problem by deriving an accurate and fast solution to the relatively simpler problem of fluid simulation in which all the conditions are well defined and given. The more complicated reconstruction problem, which entails such uncertainties as unknown boundary conditions, is then solved with the aid of an accurate simulation model.
After identifying that the viscosity parameters in SPH were not quantitatively anchored to kinematic viscosity, we propose a fast gradient-based method to optimize them. Even though kinematic viscosity \(\nu\) plays an important role in describing the fluid flow’s behavior, as indicated by the definition of the Reynolds number \({{\it Re}}=u L / \nu\), accurate modeling of viscosity in SPH is non-trivial due to the unknown relationship between real-world kinematic viscosity \(\nu\) and the viscosity parameters in the SPH simulation, especially in the presence of both fluid and boundary particles. We iteratively search for the optimal parameters until a reference flow can be replicated in simulations. Upon optimization, the same viscosity parameters work well across a wide range of Reynolds numbers.
With optimized simulation parameters, we proceed to solve the real-time reconstruction problem by sampling the flow data from real-world liquid and applying a guiding force to the simulated liquid. In this framework, fluid simulation and fluid sampling are performed concurrently in real time. The motivation behind this concurrent operation lies in the complementary nature of fluid simulation and fluid tracking. Fluid tracking captures the actual motion of the fluid flow with high accuracy, but in low resolution; fluid simulation predicts fluid flow with arbitrary resolution. In our framework’s implementation, the flow data are sampled using an occlusion-free, electromagnetic 3D-motion-tracking system. Our framework performs reconstruction by converting the sampled flow data in the real world into a guiding force in the simulation using reinforcement learning. Compared to the existing re-simulation methods, the main advantages of our reconstruction framework are speed and resilience. The guiding force is computed in real time using a deep neural network (DNN). Our framework is resilient enough to perform reconstruction even if an arbitrary external force is applied to the real-world liquid. We verify that good reconstruction accuracy can be achieved across different patterns of external disturbance. Using our real-time reconstruction framework, unprecedented interactive fluid applications become possible since the effect of users’ disturbance on a physical liquid is reflected immediately in the virtual domain.
In summary, a real-time flow reconstruction problem entails three challenges: the need for a fast and accurate simulator, the difficulty of measuring the fluid flow with occluding objects in real time, and the limited knowledge of the disturbance source. An accurate simulator is implemented using state-of-the-art SPH methods along with our parameter optimization method in Section 4. Real-time occlusion-free flow measurement is accomplished through an electromagnetic motion-tracking system introduced in Section 5. With the effect of the external disturbance captured at sparse sampling points, we explain the generation of a guiding force from the sampled data for flow reconstruction in Section 6. The following are the contributions of this article:
– | a novel real-time framework that reconstructs fluid flow even in the presence of an unknown disturbance; | ||||
– | the first use of 3D-motion-tracking data and reinforcement learning for real-time flow reconstruction; and | ||||
– | the efficient optimization of SPH viscosity parameters using a gradient-based approach. | ||||
Code and trained models for this article are hosted on GitHub ( https://github.com/kennychufk/alluvion-rl).
2 RELATED WORK
SPH [Monaghan 1992] is a prominent mesh-free approach first used by Stam and Fiume [1995] for fluid simulation in computer graphics. It is often used in real-time fluid simulation due to its simplicity and parallelizability. SPH is a relatively simple fluid model since all the physical quantities of a liquid are represented by particles, each of which has a radius of influence. The largest allowed timestep primarily depends on the efficiency of the pressure solver that maintains the incompressibility condition. Various pressure solvers, such as Predictive-Corrective Incompressible SPH (PCISPH) [Solenthaler and Pajarola 2009], Position Based Fluids (PBF) [Macklin and Müller 2013], and Divergence-Free SPH (DFSPH) [Bender and Koschier 2017], were proposed to enlarge the timesteps while keeping the simulation stable. SPH is also a highly extensible method that supports multi-physics simulation [Lenaerts and Dutré 2009], rigid-fluid coupling [Akinci et al. 2012], and even rigid-rigid interaction [Gissler et al. 2019]. A comprehensive review of SPH methods is available [Koschier et al. 2019].
SPH methods adopt a few parameters that deviate from the definitions in physics, where the viscosity parameter is a notable example. Although viscosity in SPH can be naively implemented using the SPH discretization of the Laplace operator, this naive implementation is sensitive to slight changes in particle configuration [Price 2012]. As a result, numerous explicit [Schechter and Bridson 2012] and implicit [Weiler et al. 2018] viscosity models have been proposed to respectively simulate materials of low and high viscosity. These models employ viscosity parameters that are only weakly associated with the physical definitions of kinematic viscosity or dynamic viscosity. The relationship between the parameter and physical viscosity values is not explicitly defined by the authors of the above methods. SPH’s discretization scheme also requires such parameters as smoothing kernel radius and particle volume. These modeling parameters affect the resolution of the simulation and the computational cost, although they do not correspond to any physical property of fluids.
Parameter optimization for fluid simulation was only recently reported in the literature. Roselli et al. [2018] used a genetic algorithm to obtain parameters for the accurate reproduction of Stokes waves. de Anda-Suárez et al. [2018] compared the performance of different evolutionary metaheuristics for optimizing SPH parameters. Roselli et al. and de Anda-Suárez et al. theoretically computed the ground truth that optimizes the parameters. Conversely, Takahashi and Lin [2019] used the video footage of a viscous liquid to estimate its viscosity through an evolutionary algorithm. Differentiable physics, which is recently receiving research attention, has become an alternative direction for parameter optimization. Du et al. [2020] used gradient-based optimization and backpropagation to efficiently optimize the shape of fluidic devices under Stokes flow. Following this recent trend, we use a gradient-based approach with theoretical ground truth to optimize SPH viscosity parameters.
Fluid flow measurement is an active research field in fluid dynamics for the empirical study of fluid phenomena. Particle Image Velocimetry (PIV) [Lourenco et al. 1989] is a well-established method for the high-resolution measurement of fluid velocity. Tracer particles are suspended in the fluid and a cross section of the fluid is illuminated by a laser sheet. Tracer particles illuminated by the laser sheet are then recorded with a high-speed camera. The particle movements are analyzed using cross-correlation to give a velocity field for the cross section. This basic PIV setup can be extended to support three-dimensional measurement using multiple cameras [Elsinga et al. 2006]. Rainbow PIV [Xiong et al. 2017] simplified the setup so that only one camera is required for three-dimensional tracking by encoding the depth information using color. Instead of the velocity field, Ye et al. [2012] recovered the height field and the water’s surface normal using a camera to capture the distortion of an underwater image pattern. The surface reconstruction technique of Ye et al. [2012] was extended using deep learning by Thapa et al. [2020]. For the three-dimensional Lagrangian measurements of a flow field, Particle Tracking Velocimetry (PTV) [Maas et al. 1993] can be used to trace the movements of individual particles. In this article, we require the sampling of real-world fluid data to evaluate the reconstruction accuracy of our framework. We respectively used IM3D+ [Huang et al. 2022], an electromagnetic motion-tracking system, and PIV to capture sparse 3D Lagrangian flow data and high-resolution 2D velocity fields. With IM3D+, ad-hoc perturbation using occluding objects is allowed since the tracking system supports real-time, occlusion-free tracking of light electromagnetic markers that float on liquids.
Re-simulation with real flow data has been extensively studied in the context of weather forecasting [Ghil and Malanotte-Rizzoli 1991; Kalnay 2002]. In the meteorology community, a methodology called data assimilation incorporates such observation data as atmospheric pressure into numerical simulations to predict the future state and interpolate the spatially sparse observation data. On a smaller scale, measurement-coupled simulation can also accurately simulate complex flows such as Kármán vortex streets and blood flows [Hayase 2015]. A similar methodology can be seen in fluid simulation for computer graphics. The stereoscopic videos of a liquid were used by Wang et al. [2009] to reconstruct its fluid surface. In their work, an incompressible solver was used to physically guide the reconstructed surface. In the aforementioned work of Rainbow PIV [Xiong et al. 2017], 3D vector fields were reconstructed from a series of colored particle images with an optimization method based on physical equations. Gregson et al. [2014] proposed a fluid re-simulation method that increases the resolution of the captured fluid flow data based on Navier-Stokes equations, adding details to the captured data while satisfying the physical constraints of an incompressible flow. Eckert et al. [2018] reconstructed the density map of smoke plumes from just a sequence of images using proximal operator-based optimization. Eckert et al. [2019] also proposed a flow reconstruction algorithm that estimates inflow boundary conditions of smoke plumes and employs correction forces to match observed multi-view flow data. Nagasawa et al. [2019] proposed an algorithm to estimate the motion of a blended liquid using the measured viscosities of its constituents.
Fluid control methods influence the behavior of simulated fluids based on certain reference data. Treuille et al. [2003] proposed to control smoke animations through user-defined keyframes. They used a gradient-based optimization scheme to find the guiding forces required to align the flow of the smoke with the keyframes while respecting the governing physical equations. McNamara et al. [2004] improved the computational efficiency of the optimization process and extended the idea to liquid animations. Shi and Yu [2005] proposed using the principle of control theory to guide a mesh-based fluid simulation toward a target state. The discrepancy between the simulated fluid and the target shape is fed back to the fluid simulation as control forces. Thürey et al. [2006] extended fluid control methods to SPH and Lattice-Boltzmann Method (LBM) simulations. For SPH, the target shapes are modeled as control particles that exert corrective forces on nearby fluid particles. We borrowed this design of corrective forces and used reinforcement learning to determine the suitable force parameters at every simulation step to achieve flow reconstruction.
3 FRAMEWORK OVERVIEW
As illustrated in Figure 1, our reconstruction problem is defined as the reconstruction of liquid motion that is moving under the effect of disturbance with only sparse flow data as input while remaining oblivious to the disturbance source’s shape and motion. For the input of the reconstruction problem, although different flow-sampling methods are possible, to simplify the discussion, flow sampling primarily refers to buoy-based Lagrangian sampling in this article. For brevity, we also refer to the source of the disturbance as an agitator.
Fig. 1. Our framework reconstructs the motion of a liquid moving under the effect of external disturbance. Buoys are floated on the liquid to coarsely measure the flow velocity. Without knowing the details of the disturbance source, our framework uses the trajectories of the buoys alone to reconstruct the fluid flow.
An idealistic approach to the reconstruction problem is to solve for all unknowns, including the liquid’s flow and the agitator’s shape and trajectory, from the kinematic information of the buoys and the governing physical laws. However, with the sparsity of information provided by buoys in practice, such an approach will be futile. Instead, we omit the agitator in the reconstruction result and only focus on the liquid. An obvious limitation of our approach is its inability to model the displaced volume when the agitator is submerged. As a result, we assume that the agitator’s volume is no more than 1% of the total liquid volume.
Compared to an ordinary fluid simulation where all the initial and boundary conditions are well-defined, our reconstruction problem is more challenging because the agitator’s information is hidden. Figure 2 visualizes the different dataflows in ordinary simulations and in our reconstruction framework. In ordinary simulations, the boundary condition is defined by the agitator’s trajectory, the container’s shape, and the geometric and mass properties of the buoys. The trajectories of the buoys and the fluid’s movement can then be uniquely determined using a fluid simulator. For reconstruction, the trajectories of the buoys serve as an input to the reconstruction algorithm. The reconstruction framework attempts to recover the effect of the unknown agitator on the fluid based on the input buoy trajectories.
Fig. 2. Flow chart summarizing interaction among buoys, agitator, and fluid in reality and in reconstruction.
Our framework reconstructs the fluid flow by simultaneously performing fluid simulation and guiding the simulation result based on the buoy trajectory inputs. At each simulation step, a guiding force with specific strength and spatial extent is generated for each buoy with a DNN. The guiding force is input to the simulator to affect the motion of the simulated fluid. With sufficient training data, the DNN is expected to produce appropriate guiding forces that result in a high reconstruction accuracy under different circumstances. In our framework’s implementation, the simulator within the reconstruction framework is implemented using SPH. However, the framework does not depend on the kind of simulation method used as long as the guiding forces can be implemented.
The DNN is trained using reinforcement learning to circumvent the need for differentiable SPH simulations. Although we have access to the ground-truth flow data, training the DNN is difficult with traditional supervised learning combined with backpropagation. Automatic differentiation implementations for fluid simulation, such as Taichi [Hu et al. 2020] and SPNets [Schenck and Fox 2018], have been proposed. However, we tested and confirmed that automatic differentiation is impractical for our use case due to the frequent occurrence of gradient explosions when we simulate beyond 20 simulation steps with thousands of particles. As a result, the gradient of the loss cannot be backpropagated to the network parameters through the SPH simulation process. We instead train the DNN using a reinforcement learning algorithm that supports a continuous action space known as Twin Delayed Deep Deterministic Policy Gradient (TD3) [Fujimoto et al. 2018].
DNN training requires a dataset of fluid moving under the effect of external disturbance. Since it is difficult and time-consuming to collect such a dataset from the real world using fluid measurement tools like PIV, we generate it using SPH.
4 FLUID SIMULATION METHOD
Due to the large variety of available SPH simulation methods, the SPH algorithms employed in this research were first specified before further discussion. We solved the pressure Poisson equation defined in Incompressible SPH [Bøckmann et al. 2012] with the weighted Jacobi method. We used the cubic spline kernel for all the SPH computations except for the guiding force in Section 6. Boundary particles were initialized using the sampling method proposed by Kugelstadt et al. [2021]. To better model boundary pressure, boundary particles with moving least squares pressure extrapolation [Band et al. 2018] were used for solid representation. For the viscous term, to avoid the direct SPH discretization of the Laplace operator, the term is commonly evaluated by combining an SPH derivative with a finite difference [Shahriari et al. 2013], resulting in the following velocity update equation for fluid particle i with fluid neighbors j and solid boundary neighbors b: (1) \(\begin{equation} \begin{split}\Delta \mathbf {v}_{i} & = \Delta t \sum _{j} \nu _{\mathcal {F}}\frac{m}{\rho _{j}}\frac{\mathbf {x}_{ij} \cdot \nabla W\left(\mathbf {x}_{ij}, h\right)}{\Vert \mathbf {x}_{ij}\Vert ^2 + \epsilon ^2 h^2} \mathbf {v}_{ij} \\ & + \Delta t \sum _{b} \nu _{\mathcal {B}}\frac{m}{\rho _{0}}\frac{\mathbf {x}_{ib} \cdot \nabla W\left(\mathbf {x}_{ib}, h\right)}{\Vert \mathbf {x}_{ib}\Vert ^2 + \epsilon ^2 h^2} \mathbf {v}_{ib} , \end{split} \end{equation}\) where \(\mathbf {x}_{ij} =\mathbf {x}_{i} - \mathbf {x}_{j}\), \(\mathbf {v}_{ij} =\mathbf {v}_{i} - \mathbf {v}_{j}\), \(\rho _0\) is the rest density, m is the particle mass, and \(\epsilon\) is a small constant we set as 0.1. This viscosity formulation respectively results in two kinematic viscosity parameters, \(\nu _{\mathcal {F}}\) and \(\nu _{\mathcal {B}}\), for fluid neighbors and boundary neighbors. Surface tension was not modeled as it contributes little to the overall flow in our scenarios.
4.1 Viscosity Parameter Optimization
We aim to find the optimal simulation parameters to faithfully model a liquid with known kinematic viscosity \(\nu\) under a wide range of Reynolds numbers. Searching for the optimal viscosity parameters in simulations is important because it is a major component of the Reynolds number. Although such physical quantities as fluid density and gravity can be specified as is in SPH simulations without raising any issues, directly using the kinematic viscosity value does not lead to the desired viscous behavior.
We propose using well-studied flow patterns with analytic solutions as the ground truth for parameter optimization. We chose the Hagen-Poiseuille flow because it can be easily modeled using SPH. An infinitely long cylindrical pipe is modeled using a periodic boundary condition (Figure 3). In other words, the particles that leave the plane at \(y=H\) are moved back to the simulation domain by subtracting its y-coordinate by 2H.
Fig. 3. Hagen-Poiseuille flow is modeled using a cylinder with its two ends connected through a periodic boundary condition. A constant pressure gradient is applied by accelerating each fluid particle with identical acceleration.
To accurately model the liquid volume using SPH, a water-tight representation of the pipe and a set of initial fluid particle positions resulting in the correct and uniform density distribution are required. With \(n_{p}\) fluid particles each having a particle volume V, pipe radius R can be determined geometrically using the equation \(R=\sqrt {n_{p} V / (2\pi H)}\). To model the pipe that can just contain \(n_{p}\) fluid particles, we first initialize boundary particles for a pipe model with height 2H and a sufficiently large radius. With the pressure solver running, the radius of the pipe model is gradually reduced so that both the fluid and the boundary particles redistribute themselves. The optimal pipe representation is achieved when the density of all the fluid particles is the closest to the rest density \(\rho _0\).
Pressure gradient \(\partial p / \partial y\) across the pipe is maintained by accelerating all the fluid particles with constant acceleration a. As time proceeds, the difference between the axial and peripheral velocities increases, giving the distinctive velocity profile shown in Figure 4. The velocity profiles at different time instances are recorded and compared against the analytic solution of the developing Hagen-Poiseuille flow: (2) \(\begin{equation} \begin{split}v_{y}\left(r, t\right) & = \frac{a}{\nu }\left[\frac{1}{4}\left(R^2-r^2\right) \right. \\ & \left. - 2R^2 \sum _{n=1}^{\infty }\frac{1}{\lambda _{n}^3}\frac{J_0\left(\lambda _{n}\frac{r}{R}\right)}{J_1\left(\lambda _{n}\right)}\exp \left(-\left(\frac{\lambda _{n}}{R}\right)^2 \nu t\right)\right], \end{split} \end{equation}\) where \(J_{\alpha }(x)\) is the Bessel function of the first kind of order \(\alpha\) and \(\lambda _{n}\) are the positive roots of \(J_0(x)=0\).
Fig. 4. Identical optimized parameters work in different Reynold numbers.
By evaluating the velocity error with the mean-squared error, we optimized the two viscosity parameters using gradient-based optimization. The gradient of the error with respect to the parameters was evaluated using forward difference. Although the Hagen-Poiseuille flow is a simple laminar flow, the numerical precision of the estimated gradient still deteriorates significantly when the number of simulation steps is too large. This deterioration happens because forward difference is very sensitive to truncation and rounding errors that are accumulated along many simulation steps. Therefore, it is important to evaluate the error when the flow is still developing instead of when it is close to the terminal state. The half-life of the flow is the time required to reach half of the terminal velocity at the center. With half-life \(t_{1/2}\) of the developing flow estimated as \(R^{2}\ln 2 / (\nu \lambda _0^2)\), we evaluated the velocity error and the error gradient at \(t=0.5t_{1/2}\), \(1.125t_{1/2}\), and \(1.5t_{1/2}\). The viscosity parameters are then updated iteratively using an Adam optimizer with a learning rate of \(2 \times 10^{-4}\).
4.2 Optimization Results
In our experiments, we set \(h=2.5 \,\mathrm{m}\mathrm{m}\) and fixed the number of particles to 9,900, resulting in a total liquid volume of approximately \(15.47 \,\,\mathrm{c}\mathrm{m}^{3}\,\). For the periodic boundary condition, we set \(H = 3h=7.5 \,\mathrm{m}\mathrm{m}\) and geometrically deduced that the pipe is perfectly filled by particles when \(R=18.12 \,\mathrm{m}\mathrm{m}\). The optimization result for a liquid of kinematic viscosity \(\nu =1 \,\mathrm{c}\rm {St}=1 \,\mathrm{m}\mathrm{m}^{2}\,\mathrm{s}^{-1}\), a value that roughly equals the viscosity of water at 20 \(^\circ\)C, is shown in Figure 5. All particles were accelerated at \(a=6.73 \times 10^{-4} \,\mathrm{m}\mathrm{m}\mathrm{/}\,\mathrm{s}^{2}\,\) to give \({{\it Re}}=1\) at the terminal state. After optimizing the two parameters for 1,200 iterations, both viscosity parameters converged to a solution with a mean-squared error of \(4.10 \times 10^{-9} \,\mathrm{m}\mathrm{m}^{2}\,\mathrm{/}\,\mathrm{s}^{2}\,\). On GeForce RTX 3090, performing 1,200 optimization steps takes approximately \(28 \,\mathrm{min}\).
Fig. 5. Optimization process of viscosity parameter values for pure water.
With this optimized set of viscosity parameters, we verify that the simulated velocity profiles remain accurate at higher Reynolds numbers. Figure 4 shows that satisfactory simulation accuracy can still be achieved at \(Re=500\).
By repeating the optimization process for different kinematic viscosities, we obtained the relationship between the viscosity parameters and the actual kinematic viscosity values. Both \(\nu _{\mathcal {F}}\) and \(\nu _{\mathcal {B}}\) are directly proportional to the kinematic viscosity value (Figure 6). Linear regression is performed across the data points to obtain two linear conversion equations for \(\nu _{\mathcal {F}}\) and \(\nu _{\mathcal {B}}\). Liquids of arbitrary viscosities can then be easily simulated using the corresponding viscosity parameters.
Fig. 6. Linear relationship between SPH viscosity parameters and \(\nu\) .
Although the Hagen-Poiseuille flow simulation was already implemented using SPH by Takeda et al. [1994] and later by Shahriari et al. [2013], using an analytic solution to optimize viscosity parameters is novel. Our techniques of preparing a water-tight pipe and ensuring accurate gradient estimation are crucial to close resemblance to the theoretical flow and fast convergence. Furthermore, although there have been recent advancements in SPH pressure solvers and boundary representation in the computer graphics community, it has not yet been shown quantitatively whether these new techniques can model fluids with good physical accuracy. To our knowledge, we are the first to present a precise mapping between kinematic viscosity and its simulation parameters.
5 LAGRANGIAN FLOW TRACKING
To sample the flow data from real-world liquids for real-time reconstruction, a fluid tracker that satisfies the following requirements is necessary:
(1) | support for three-dimensional flow measurement; | ||||
(2) | a temporal resolution comparable to the frame rates of the interactive applications, i.e., on the order of milliseconds; | ||||
(3) | support for real-time output of motion data; and | ||||
(4) | no suffering from occlusion when users interact with the fluid. | ||||
Requirement (4) is the most demanding because, as seen in Section 2, the vast majority of existing measurement methods rely on optical sensors. Moreover, in many existing methods, analyzing raw data is computationally expensive, making them only suitable for offline applications. To our knowledge, since no existing flow measurement system satisfies the above requirements, we turn our attention to motion-tracking systems in general. IM3D+ [Huang et al. 2022], an electromagnetic motion tracker, can be designed to support interactive fluid tracking. IM3D+ operates on the principle of electromagnetic resonance and thus avoids occlusion because water and other weakly diamagnetic materials cause almost no distortion of a magnetic field.
IM3D+ consists of the following aspects: (1) a driving coil that produces an oscillating magnetic field; (2) a 2D array of inductive coil sensors that measure the magnetic flux changes; and (3) multiple inductor-capacitor coils (LC coils), each of which has a distinct resonant frequency, as tracking targets. IM3D+ tracks the position of the LC coils at a refresh rate of \(100 \,\mathrm{Hz}\). The mean position and angular errors are \(1.87 \,\mathrm{m}\mathrm{m}\) and \(0.72 \,\mathrm{^{\circ }}\), respectively. To track the Lagrangian flow data of water, we encapsulated each LC coil in a light cylindrical tube to form a buoy that floats on water (Figure 7). The buoy is made cylindrical so that the LC coil tends to return to an upright orientation even when tilted, ensuring a high level of resonance for accurate tracking. All buoys have almost identical mass and length.
Fig. 7. LC coils are encapsulated in translucent tubes with plastic caps on both ends. Each buoy is \(38.5 \,\mathrm{m}\mathrm{m}\) in length and has a radius of \(3.0 \,\mathrm{m}\mathrm{m}\) . The markings on the tubes are for identification purposes during video recording.
The major challenges of using IM3D+ to capture the flow data lie in the sparsity of the sampling points and the fact that buoys may detach from the surrounding flow when the flow velocity is high. The discrepancy between the buoy velocity and the flow field is related to Stokes number Stk, a quantity defined for suspending particles. It is desirable to have a buoy with \(Stk \ll 1\) so that it closely follows the streamlines. However, buoys generally have a large Stokes number due to their large size. Therefore, note that the buoys are not the exact Lagrangian description of the fluid, especially when the flow abruptly decelerates. Nonetheless, they still loosely represent the surrounding fluid flow and capture the effect of the agitator. The reconstruction framework is also expected to take the discrepancy into account and produce a suitable guiding force even when flow detachment occurs. In our current hardware implementation, at most nine buoys can be deployed simultaneously. As it is likely that improved motion-tracking technologies will support more buoys in the future, our evaluations in Section 7 also consider scenarios with more than nine buoys.
6 RECONSTRUCTION METHOD
6.1 Guiding Force Field
We first show how guiding forces can be incorporated into a typical SPH solver before explaining our design of the guiding force field. Guiding forces exerted on the simulated liquid can be considered an external force in the incompressible Navier-Stokes equation: (3) \(\begin{equation} \frac{\mathrm{D}\mathbf {v}(\mathbf {x}, t)}{\mathrm{D}t}=-\frac{1}{\rho }\nabla p(\mathbf {x}, t) + \nu \nabla ^{2} \mathbf {v}(\mathbf {x}, t) + \mathbf {g} + \mathbf {a}(\mathbf {x},t), \end{equation}\) where \(\mathbf {v}\) is the fluid velocity, \(\mathbf {g}\) is the gravitational acceleration, and \(\mathbf {a}\) is the guiding acceleration. Note that the last three terms on the right-hand side of the equation are non-pressure acceleration. The equation is solved through operator splitting [Ihmsen et al. 2014b] as outlined in Algorithm 1. Within each simulation step, after advancing the velocity with the non-pressure acceleration to give \(\mathbf {v}^\ast\), the interim density \(\rho ^\ast\) can be calculated from the discretized continuity equation \((\rho ^\ast - \rho) / \Delta t=-\rho \nabla \cdot \mathbf {v}^\ast\). Pressure acceleration \(a^p=-\frac{1}{\rho } \nabla p\) that cancels the density error \(\rho _0 - \rho ^\ast\) is required to restore incompressibility. This cancellation can be described by the equation \((\rho _0 - \rho ^\ast) / \Delta t= -\rho \nabla \cdot (\Delta t \mathbf {a}^p)\). The resultant pressure Poisson equation (PPE), \(\rho _0-\rho ^\ast = \Delta t^2 \nabla ^{2}p\), can be solved using the weighted Jacobi method. Any divergence in the guiding acceleration will be removed after solving the PPE.
We base our guiding force design on existing fluid control methods. To approach a reference velocity field \(\mathbf {v}_{\text{t}}(\mathbf {x}, t)\) observed at time t in the real world, its difference from the simulated velocity \(\mathbf {v}(\mathbf {x}, t)\) is fed back to the simulation in the form of compensating acceleration [Shi and Yu 2005]: (4) \(\begin{equation} \mathbf {a}(\mathbf {x}, t)=s \left(\mathbf {v}_{\text{t}}(\mathbf {x}, t) - \mathbf {v}(\mathbf {x}, t) \right), \end{equation}\) where s is a constant gain factor. After an iteration of Algorithm 1, the new velocity error \(\Vert \mathbf {v}_{\text{t}}(\mathbf {x}, t) - \mathbf {v}(\mathbf {x}, t+\Delta t) \Vert\) is expected to be smaller than the original error if s is sufficiently small to avoid over-compensation. The algorithm then proceeds to approach the reference field at \(t+\Delta t\) while the previous reference field \(\mathbf {v}_{\text{t}}(\mathbf {x}, t)\) ceases to influence the simulation.
As the target velocity field is only observed at the positions of buoys, it is necessary to extend the spatial influence of each buoy’s guiding force to drive fluid regions that are not observed in the real world. A popular candidate for this spatial extension is Gaussian blur [Treuille et al. 2003; Thürey et al. 2006]. For buoy k at position \(\mathbf {x}_{k}\) with velocity \(\mathbf {v}_{k}\), the acceleration exerted on fluid particle i is given by (5) \(\begin{equation} \mathbf {a}_{k \rightarrow i} = s_k e^{-\frac{\Vert \mathbf {x}_{k} - \mathbf {x}_{i} \Vert ^2}{{h_k}^2}} \left(\mathbf {v}_{k} - \mathbf {v}_{i} \right), \end{equation}\) where \(s_k\) and \(h_k\) denote the gain and the kernel radius for buoy k, respectively. The total guiding acceleration on particle i is simply the superposition of contributions from all the buoys, i.e., \(\mathbf {a}_i = \sum _{k=1}^{n_b} \mathbf {a}_{k \rightarrow i}\). We use a DNN agent trained with reinforcement learning to decide suitable values for \(s_k\) and \(h_k\).
Due to flow detachment, the observed velocity \(\mathbf {v}_{k}\) may not accurately represent the surrounding fluid velocity. Also, as the actual disturbance originates from the agitator while the buoys only capture its effect, the agent is allowed to shift the origin of the guiding force away from the buoy and adjust the guiding velocity. Position and velocity offsets \(\mathbf {x}^\ast _{k}\) and \(\mathbf {v}^\ast _{k}\) are introduced to allow this freedom. Finally, the agent is expected to intermittently switch off the guiding force when the buoy is hit directly by the agitator or when there are too many buoys nearby. However, through our experiments, we found that the DNN agent rarely outputs near-zero values for \(s_k\). This issue is mitigated by multiplying \(s_k\) with a Heaviside step function \(\theta \left(d_k\right)\), where \(d_k\) is a real-valued output of the agent. The acceleration then becomes (6) \(\begin{equation} \mathbf {a}_{k \rightarrow i} = \theta \left(d_k\right) s_{k} e^{-\frac{\Vert \mathbf {x}_{k} + \mathbf {x}^\ast _{k} - \mathbf {x}_{i} \Vert ^2}{{h_k}^2}} \left(\mathbf {v}_{k} + \mathbf {v}^\ast _{k} - \mathbf {v}_{i} \right) . \end{equation}\)
Unlike optimization-based fluid control methods, control theory-based methods including ours do not guarantee that the target velocity is approached with a minimum error every time [Pan and Manocha 2017]. Nonetheless, our method contains optimization elements as we train the agent to search for good parameter values of the guiding force according to the environment. In terms of generality, Gaussian guiding force alone cannot compensate for all types of discrepancy between target and simulation states. For example, it is impossible to generate small-scale vortices [Thürey et al. 2006]. A vortex flow can only be generated by arranging multiple Gaussian forces in a circular pattern, imposing a limit on the minimum length scale of the vortex. Moreover, it was shown by Shi and Yu [2005] that a potential field needs to be used in conjunction with feedback control forces to reconstruct a rapidly changing target field. Nonetheless, similar Gaussian force designs can still be utilized to produce a wide variety of sophisticated target patterns [Treuille et al. 2003; McNamara et al. 2004]. Gaussian force is adopted in our design for the sake of simplicity.
6.2 Reinforcement Learning
Different parameters of the force fields are generated in a data-driven approach because it would be very difficult to rigorously derive the required force field parameters in different flow conditions using physical laws. We employed TD3 [Fujimoto et al. 2018], a robust actor-critic reinforcement learning algorithm for continuous action space, to train the DNN. In the context of reinforcement learning, the DNN that outputs the force parameters is the actor, also commonly known as the policy network. The critic is another DNN that assesses how favorable it is for the policy network to generate certain force parameters under a specific flow condition. We focused our discussion on the policy network as it entails more engineering choices.
To support an arbitrarily large number of buoys, the policy network is designed to be a homogeneous agent [Gupta et al. 2017] so that every buoy shares the same network parameters. The policy network is therefore invoked \(n_b\) times in every simulation step when there are \(n_b\) buoys in the scene. To make informed decisions about the guiding force, the network accepts an observation vector encoded with a collection of buoy and simulation states. Through extensive experiments, we arrived at the design of the observation vector in Figure 8. Besides the velocity \(\mathbf {v}_{k}\) and the quaternion \(\mathbf {q}_{k}\) of buoy k, four exponential moving averages \(\bar{\mathbf {v}}_k\) of the buoy velocity with smoothing factors \(\alpha _1\) to \(\alpha _4\) are included to represent the velocity history. By comparing \(\mathbf {v}_k\) against the mean speed of the buoys \(\sum _{i=1}^{n_b}\Vert{\mathbf {v}_i}\Vert/n_b\), anomalies such as a buoy being hit by the agitator can be inferred. To avoid lopsided guiding forces, it is crucial to convey information about the buoys adjacent to buoy k. Therefore, the kinematic information of the 11 closest buoy neighbors \(j_1\) to \(j_{11}\) sorted in ascending order of distance are included. The proximity information of the buoy to the boundary of the container is encoded as \(mW_{\mathcal {B}}\left(\mathbf {x}_{k}, h\right)/\rho _{0}=\sum _{b} m W\left(\mathbf {x}_{kb}, h\right)/\rho _{0}\) and \(m\nabla W_{\mathcal {B}}\left(\mathbf {x}_{k}, h\right)/\rho _{0}\). Finally, to inform the network about the effect of the guiding forces, fluid velocity \(\mathbf {v}(\mathbf {x}_k)\) and fluid density \(\rho (\mathbf {x}_k)\) are sampled from the reconstructed liquid at \(\mathbf {x}_k\).
Fig. 8. Inputs and outputs of policy network for buoy k. Golden texts represent information of buoy k. Orange texts represent the global buoy information. When the number of buoys \(n_b\) is less than 12, the invalid buoy neighbor vectors such as \(\mathbf {x}_k-\mathbf {x}_{j_{11}}\) and \(\mathbf {v}_k-\mathbf {v}_{j_{11}}\) are replaced by zero vectors. Blue and purple texts represent fluid and boundary states in the simulation, respectively. Information about the guiding force field is represented in red characters.
6.3 Datasets
An entry of a dataset is essentially sequential frames of flow fields, buoy coordinates, and agitator coordinates sampled at \(100 \,\mathrm{Hz}\). Our datasets consist of synthetic datasets for both training and validation, and real-world empirical sets for additional validation. Although our framework does not dictate whether training data should be acquired from the real world or generated through simulations, the former approach is less practical due to the difficulty of measuring a detailed velocity field, especially when occlusions by the buoys and the agitator occur frequently. Using our SPH simulator with optimized viscosity parameters, we can generate a reasonably large set of reliable training data.
For simplicity in both the discussion and implementation, the same cubical container with a length of \(24 \,\mathrm{c}\mathrm{m}\) is used in all datasets. Nonetheless, the observation and action vectors of the policy network are designed to operate regardless of the container’s shape. The container contains either \(3 \,\mathrm{L}\), \(3.25 \,\mathrm{L}\), or \(3.5 \,\mathrm{L}\) of liquid, amounting to a liquid height of \(5.2 \,\mathrm{c}\mathrm{m}\) to \(6.1 \,\mathrm{c}\mathrm{m}\). The liquid is either pure water (\(\nu =1.03 \,\mathrm{c}\rm {St}\), \(\rho _0=997 \,\mathrm{\mathrm{k}\mathrm{g}}\mathrm{/}\,\mathrm{m}^{3}\,\)) or a glycerol-water mixture (\(\nu =18.94 \,\mathrm{c}\rm {St}\), \(\rho _0=1156 \,\mathrm{\mathrm{k}\mathrm{g}}\mathrm{/}\,\mathrm{m}^{3}\,\)).
While only four to nine buoys can be deployed in the empirical datasets, the number of buoys \(n_b\) varies from 4 to 100 in the synthetic datasets. For synthetic sets, each buoy is modeled as a rigid cylinder using the physical properties described in Appendix A. The small radius of the buoy requires a small particle radius of \(0.98 \,\mathrm{m}\mathrm{m}\), equivalent to a kernel radius of \(h=3.91 \,\mathrm{m}\mathrm{m}\), for accurate modeling of the buoyant force.
Different agitator shapes are required to avoid overfitting. We used the Stanford Bunny and six hand-picked objects from ShapeNetCore [Chang et al. 2015] as the agitator shapes. Their volumes range from \(22 \,\,\mathrm{c}\mathrm{m}^{3}\,\) to \(26.5 \,\,\mathrm{c}\mathrm{m}^{3}\,\). The agitator follows one of the parametric curves in Figure 9.
Fig. 9. Trajectories of agitators used for training and evaluation.
For training, we created a main set for final training and a secondary set for case studies. The agitator’s trajectory in the main set is Stars & trefoiloids (Figure 9(b)) while the secondary set uses the simpler trajectory Diagonal oscillations (Figure 9(a)). In the context of reinforcement learning, a training step finishes when the policy network finishes reconstructing a frame. An episode finishes when all frames in an entry of a dataset are reconstructed. For example, as each entry of the secondary training set spans \(10 \,\mathrm{s}\) and the frame rate is \(100 \,\mathrm{Hz}\), an episode is equivalent to 1,000 steps.
6.4 Eulerian Reward Function
The training process requires a reward function \(R(t)\) that evaluates the closeness between the reconstructed flow and the truth flow at time t. We considered an Eulerian approach and a volumetric approach to the design of \(R(t)\).
The Eulerian approach samples the truth velocity field \(\mathbf {v}_{\text{t}}(\mathbf {x}, t)\) and the reconstructed field \(\mathbf {v}(\mathbf {x}, t)\) from the truth flow and the reconstructed flow using a set of uniformly placed points E within the container. The reward function can then be defined as the negated mean-squared error of the reconstructed field: (7) \(\begin{equation} R_{\mbox{Eulerian}}(t)=-\frac{1}{|E|}\sum _{\mathbf {x}\in E} \Vert \mathbf {v}(\mathbf {x}, t) - \mathbf {v}_{\text{t}}(\mathbf {x}, t)\Vert ^2. \end{equation}\) From the perspective of fluid dynamics, the Eulerian reward conveys a holistic description of the similarity between the two flow fields.
The volumetric approach evaluates the difference in shape between the two liquids. Isosurfaces \(f_{\text{t}}(\mathbf {x}, t)=0\) and \(f(\mathbf {x}, t)=0\) are constructed from the truth liquid surface and the reconstructed surface, respectively. With the sub-level sets \(S(f_{\text{t}}, t)=\left\lbrace \mathbf {x}\mid f_{\text{t}}(\mathbf {x}, t) \le 0 \right\rbrace\) and \(S(f, t)\) denoting the region enclosed by the isosurfaces, the shape discrepancy can be represented by their symmetric difference \(S(f, t)\ominus S(f_{\text{t}}, t)=\left[S(f,t) \setminus S(f_{\text{t}}, t)\right] \cup \left[S(f_{\text{t}}, t) \setminus S(f, t)\right]\). The reward function can then be defined as the negated volume of the symmetric difference: (8) \(\begin{equation} R_{\mbox{volumetric}}(t)=-\iiint _{S(f, t)\ominus S\left(f_{\text{t}}, t\right)}{{\,dV}}. \end{equation}\) In terms of visual perception, the volumetric reward should best capture the perceived similarity between two liquids, or the lack thereof.
The effectiveness of the two reward functions was evaluated using the secondary training set of the trajectory Diagonal oscillations (Figure 9(a)). Volumetric rewards were computed using OpenVDB [Museth 2013]. After training for 400 episodes, a gradual increase in average reward was observed for the Eulerian reward but the average reward remained stagnant throughout the training with volumetric reward. The reconstructed liquids under the volumetric reward were almost stationary, implying that the policy network chose to produce extremely weak force fields to minimize the error in liquid shape. Meanwhile, the Eulerian reward resulted in a policy network that could reconstruct the oscillatory motions in Diagonal oscillations.
Although the volumetric reward assesses similarity at a level closer to the visual perception of liquid motions, it fails to train the network due to its ambiguity. For example, in the case of a vortical flow, the change in the liquid’s shape can be very small. Even for the relatively larger waves produced by the trajectory Diagonal oscillations, the change in the shape is still more subtle than the change in the velocity field. The Eulerian reward signal thus represents the discrepancy between simulation and truth states more holistically. Nonetheless, the volumetric approach is still a metric that aligns better with visual perception, so its variant, height field score, will be used for validation.
Finally, it is worth noting that the Eulerian reward also suffers from an ambiguity issue when the flow speed is small. For example, even when the truth liquid was resting in the container, the policy network trained with the Eulerian reward occasionally produced strong forces working against gravity, resulting in a bulging but still liquid. Although the two liquids are very different, the Eulerian error is actually small in this scenario. This ambiguity can be solved by penalizing the policy network when the mean guiding acceleration on the particles is large but the mean buoy speed is small. The updated reward function becomes (9) \(\begin{equation} R(t)=R_{\mbox{Eulerian}}(t)-\beta \left(\frac{1}{n_b}\sum _{k=1}^{n_b}\Vert {\mathbf {v}_k(t)}\Vert\right)^{-1}\frac{1}{n_p}\sum _{i=1}^{n_p}\Vert {\mathbf {a}_i(t)}\Vert, \end{equation}\) where \(\beta\) is the disproportion penalty coefficient that prevents disproportionate guiding forces in slow flow conditions.
7 RESULTS AND DISCUSSION
For easier interpretation of the reconstruction accuracy, we first define two percentage scores for velocity and shape comparison. The Eulerian score is defined as the advantage of the reconstruction result over a stationary baseline: (10) \(\begin{equation} Q_{\mbox{Eulerian}}(t, w)=1-\frac{\int _{t-\frac{w}{2}}^{t+\frac{w}{2}}\sum \limits _{\mathbf {x}\in E^\ast (\tau)} \frac{\Vert \mathbf {v}(\mathbf {x}, \tau)- \mathbf {v}_{\text{t}}(\mathbf {x}, \tau)\Vert ^2}{|E^\ast (\tau)|}\,d\tau }{\int _{t-\frac{w}{2}}^{t+\frac{w}{2}} \sum \limits _{\mathbf {x}\in E^\ast (\tau)}\frac{\Vert \mathbf {v}_{\text{t}}(\mathbf {x}, \tau)\Vert ^2}{|E^\ast (\tau)|}\,d\tau } , \end{equation}\) where w is a window size for temporal smoothing and \(E^\ast (t)\) is a set of sampling points that are at least \(6 \,\mathrm{c}\mathrm{m}\) away from the surface of the agitator. The sampling points close to the agitator are omitted to eliminate the immediate local influence from the agitator. The score is designed to be a function of t so that the accuracy at any instant is well-defined. The overall Eulerian score for an episode of duration T is \(Q_{\mbox{Eulerian}}(T/2, T)\).
Through experiments, we found that height field comparisons can better highlight shape similarity than volumetric comparisons. The height field of the isosurface \(f(x,y,z,t)=0\) is denoted as \(H_f(x,z,t)\). (11) \(\begin{equation} H_f(x,z,t)=\max _{f(x,y,z,t)=0}y. \end{equation}\) The height field score is defined in a similar fashion: (12) \(\begin{equation} Q_{\mbox{height}}(t, w)=1-\frac{\int _{t-\frac{w}{2}}^{t+\frac{w}{2}}\sum \limits _{(x,z) \in E_{xz}}\left(H_f(x, z, \tau) - H_{f_{\text{t}}}(x, z, \tau)\right)^2 \,d\tau }{\int _{t-\frac{w}{2}}^{t+\frac{w}{2}}\sum \limits _{(x,z) \in E_{xz}}\left(H_{f_0}(x, z) - H_{f_{\text{t}}}(x, z, \tau)\right)^2 \,d\tau }, \end{equation}\) where \(f_0(\mathbf {x})=0\) is the surface of a baseline liquid volume resting in the container. \(E_{xz}\) is a set of uniformly placed points in the xz-plane.
For both scores, a value of zero represents zero advantage over a completely stationary liquid while a value of 100% means a perfect match. Both scores have no lower bounds.
We trained the policy network using the main training set with 15 sequences of the trajectory Stars & trefoiloids (Figure 9(b)). We mostly followed the training algorithm and hyper-parameter choices in the original TD3 implementation [Fujimoto et al. 2018] with some exceptions described in Appendix B. A larger kernel radius of \(h=15.625 \,\mathrm{m}\mathrm{m}\) was used to speed up training. With each episode spanning \(20 \,\mathrm{s}\) (2,000 frames), training the networks for six million steps (3,000 episodes) took about 105 hours on GeForce RTX 3090.
Since we aim at reconstructing fluid flow in real time, a certain period in the simulation must be advanced with less computation time. We found the highest simulation resolution possible by gradually reducing h until this real-time requirement was not satisfied. With nine buoys and \(3.25 \,\mathrm{\mathrm{k}\mathrm{g}}\) of liquid, we measured the computation time required to reconstruct \(2 \,\mathrm{m}\mathrm{s}\) of liquid flow and summarized the results in Table 1. We found that real-time reconstruction was possible when \(h=11 \,\mathrm{m}\mathrm{m}\) since it required \(1.75 \,\mathrm{m}\mathrm{s}\) of computation time to reconstruct \(2 \,\mathrm{m}\mathrm{s}\). Therefore, all evaluations are performed with this choice of h.
Table 1. Average Computation Time Required to Reconstruct \(2 \,\mathrm{m}\mathrm{s}\) of a \(3.25 \,\mathrm{\mathrm{k}\mathrm{g}}\) Liquid Flow under Different Simulation Resolutions
7.1 Validation using Synthetic Data
To test the generalizability of the trained policy, we computed the learning curve with different validation sets. As shown in Figure 10, along with the decreasing trend of the Eulerian error for the training set Stars & trefoiloids, the Eulerian errors for other agitation patterns also decrease in general.
Fig. 10. Learning curves for different agitation patterns. The training set is Stars & trefoiloids. Each curve is normalized so that the maximum error is 1.
The policy network at step number \(5.2 \times 10^{6}\) is selected for further validation because it gives the lowest error for the training set. To obtain a statistical understanding of its performance, a condition with the same trajectory and the same number of buoys was repeated eight times to calculate the mean score and the standard deviation. Across these eight trials, the only variable is the initial positions of the buoys. The performance statistics are summarized in Table 2. The height field scores for Nephroids and Epitrochoids are omitted because the height field fluctuations in the ground truth are insignificant for meaningful evaluations.
| No. of buoys | Diagonal oscillations (Figure 9(a)) | Alternating loops (Figure 9(c)) | Nephroids (Figure 9(d)) | Epitrochoids (Figure 9(e)) | ||
|---|---|---|---|---|---|---|
| Eulerian | Height field | Eulerian | Height field | Eulerian | Eulerian | |
| 4 | 49.1% \(\pm\) 3.7% | 55.6% \(\pm\) 3.6% | 25.0% \(\pm\) 4.2% | 31.2% \(\pm\) 9.8% | 26.7% \(\pm\) 3.9% | 11.4% \(\pm\) 5.1% |
| 8 | 52.9% \(\pm\) 3.1% | 59.1% \(\pm\) 3.6% | 29.2% \(\pm\) 1.8% | 27.0% \(\pm\) 6.5% | 29.9% \(\pm\) 4.8% | 19.3% \(\pm\) 3.6% |
| 22 | 50.8% \(\pm\) 1.9% | 55.7% \(\pm\) 2.4% | 33.1% \(\pm\) 0.5% | 35.7% \(\pm\) 2.7% | 37.7% \(\pm\) 1.6% | 24.6% \(\pm\) 2.2% |
| 66 | 56.9% \(\pm\) 1.2% | 67.7% \(\pm\) 1.1% | 37.3% \(\pm\) 0.7% | 41.5% \(\pm\) 1.0% | 41.7% \(\pm\) 1.9% | 22.7% \(\pm\) 2.7% |
| 100 | 51.3% \(\pm\) 1.3% | 59.9% \(\pm\) 1.1% | 35.4% \(\pm\) 0.7% | 36.0% \(\pm\) 2.1% | 30.2% \(\pm\) 1.9% | 20.7% \(\pm\) 2.5% |
Table 2. Quantitative Performance Validated using Different Agitation Patterns
Effect of the number of buoys on performance: We evaluated the Eulerian scores for Diagonal oscillations with a finer granularity in \(n_b\) and obtained the blue curve in Figure 11. The peak at \(n_b=52\) suggests that the policy network prioritizes the mid-range in \(n_b\) when trained with the main training set, which has a relatively higher variance. Nonetheless, it is possible to better utilize the increased amount of information provided by more buoys. For example, when we trained a policy network against the same trajectory Diagonal oscillations, the Eulerian score plateaued at \(n_b=40\) instead of dropping significantly beyond the mid-range.
Fig. 11. Trend of the score with an increasing number of buoys for Diagonal oscillations (Figure 9(a)). The error bars span one standard deviation.
Flow hindrance by buoys: To assess the intrusiveness of buoys as a flow measurement method, we simulated how the flow speed varies with \(n_b\) when agitated by the same trajectory Diagonal oscillations. The flow speeds at the sampling points \(\mathbf {x}\in E^\ast (t)\) across all frames are calculated and sorted. Different percentiles of the sampled speeds are plotted to give Figure 12. Using linear regression, it is shown that the overall flow speed decreases by \(0.49\%\) to \(0.69\%\) for every 10 buoys introduced to the liquid. With our current hardware limitation of \(n_b=9\), the flow hindrance by the buoys is negligible. However, buoys should not be added indiscriminately due to their slight intrusiveness and diminishing returns.
Fig. 12. Change in overall flow speed under the influence of buoys.
Diagonal oscillations: The reconstruction result with the fourth highest score among the eight trials at \(n_b=8\) is visualized in Figure 13(a). The overall trends of the Eulerian and height field score curves were similar. The initial scores were negative because the liquid in the ground truth was not disturbed until \(t=0.83 \,\mathrm{s}\), making the denominators in Equations (10) and (12) very small initially. At least \(1.0 \,\mathrm{s}\) was required for the establishment of the wave motion in the reconstruction and then the scores turned positive. Both scores started decreasing at \(t=5 \,\mathrm{s}\) but soon recovered as the change in the wave direction was registered by the buoys. A height field accuracy of more than 60% could be maintained for \(5.2 \,\mathrm{s}\).
Fig. 13. Time-averaged ( \(w=0.3 \,\mathrm{s}\) ) Eulerian and height field scores for the fourth best results at \(n_b=8\) . Images in the top row are ground truth while those in the bottom row are reconstructed height fields. Scores were negative in the beginning because the denominators of Equations (10) and (12) were very small initially.
Alternating loops: In the fourth best result at \(n_b=8\) shown in Figure 13(b), the fluctuation in the height field could be approximately reconstructed. Although the scores were lower than Diagonal oscillations in general due to the frequent alternation in the orbiting direction of the agitator, the alternating circular motion could be correctly reconstructed. Both scores dropped approximately \(0.5 \,\mathrm{s}\) after each alternation because it took time for the agitator to accelerate from a momentarily stationary state and reverse the flow. At least \(1.0 \,\mathrm{s}\) was required for the scores to recover to 50%.
Nephroids and Epitrochoids: Example reconstruction results are shown in the supplementary video. For Nephroids, the Eulerian score remained relatively steady and the location of the agitator could even be vaguely identified from the reconstruction animation. Epitrochoids were the most challenging agitation pattern as the truth flow contained mostly local swirls that could not be easily resolved from buoy movements.
Learned strategies: To illustrate the behavior of the policy network, the guiding force vectors \(\mathbf {v}_k + \mathbf {v}_k^\ast\) at the force origins \(\mathbf {x}_k + \mathbf {x}_k^\ast\) are visualized for inspection. Figure 14 shows the guiding forces for Diagonal oscillations at \(n_b=66\). Most vectors aligned with the wave direction when the flow was at its peak velocity. When the wave reached the peak height, vectors close to the wall pointed downward to prevent the liquid from further ascending. The supplementary video shows how the policy network adapted to different numbers of buoys. In particular, the network chose to disable the guiding forces from some buoys at \(n_b=100\) to avoid excessive guiding.
Fig. 14. Visualization of the policy network’s action.
7.2 Validation using Empirical Data
Empirical validation was performed to understand how the tracking noise of the buoys would affect the reconstruction accuracy. Real-world flow data were captured using a PIV system and processed with PIVlab [Thielicke and Stamhuis 2014]. By reproducing the exact trajectory using a robotic manipulator, the 2D Eulerian scores for two trajectories were summarized in Table 3.
| Trajectory | No. of buoys | 2D Eulerian Score |
|---|---|---|
| Diagonal oscillations (Figure 9(a)) | 4 | 32.7% |
| Diagonal oscillations | 5 | 53.1% |
| Diagonal oscillations | 6 | 55.8% |
| Diagonal oscillations | 7 | 50.5% |
| Diagonal oscillations | 8 | 51.2% |
| Diagonal oscillations | 9 | 48.5% |
| Oscillations & loops (Figure 9(f)) | 8 | 47.8% |
| Oscillations & loops | 9 | 38.3% |
Table 3. Quantitative Performance Evaluated Empirically using \(3.25 \,\mathrm{\mathrm{k}\mathrm{g}}\) of the Glycerol-Water Mixture
Diagonal oscillations: The reconstruction result for \(n_b=7\) is shown in Figure 15. Comparing the score curve with the synthetic validation result in Figure 13(a), although the response time was slower in the empirical result, the overall trends of the score curves were comparable. Except for \(n_b=4\), the empirical scores in Table 3 generally fell within the statistical range for synthetic results in Table 2, suggesting the validity of the synthetic validation.
Fig. 15. Variation of time-averaged ( \(w=0.4 \,\mathrm{s}\) ) score over time for Diagonal oscillations (Figure 9(a)) with seven markers.
To show that the empirical score also increases with \(n_b\), we performed multiple trials of reconstructions but with part of the buoys ignored by the policy network. For the empirical data of Diagonal oscillations with nine buoys, as we increased the number of considered buoys \(n_b\) from four to eight, i.e., as the number of ignored buoys decreased from five to one, the score distribution varied as follows: 47.0% \(\pm\) 2.1%, 48.5% \(\pm\) 1.6%, 50.5% \(\pm\) 1.4%, 51.3% \(\pm\) 1.1%, and 51.4% \(\pm\) 1.4%.
Oscillations & loops: The result for \(n_b=8\) is included in the supplementary video. Although the transition from linear oscillation to circular motion could be visually identified from the particle view of the reconstruction animation, the score only recovered slowly during the transition. The slow recovery could be attributed to the 2D nature of the PIV measurement plane, which failed to capture velocity components orthogonal to the plane.
Freehand motion: As it is very likely that a user would test the accuracy of the reconstruction framework by interacting with it with a bare hand, we provided a visual comparison under this freehand condition at the end of the supplementary video. Although a quantitative analysis was not possible due to the frequent occlusion of the laser sheet, the shape of the reconstructed liquid closely resembled the actual liquid most of the time, demonstrating the performance of our framework in practical situations.
8 LIMITATIONS
Our framework struggled to reconstruct the local flow patterns in Epitrochoids due to the scale constraint in creating vortices using Gaussian guiding forces and the limited learning capacity of our framework. A small vortex can only be generated with a collection of closely packed Gaussian forces arranged along the perimeter of the vortex, requiring a high degree of coordination for the policy network. We used a training set with only limited variations in the agitator’s trajectory (Figure 9(b)). Although a more diverse training set is desirable for a generalized agent, the increase in variance impeded the training and thus prevented us from further diversifying the training set. Variance reduction is an active research topic of reinforcement learning. For example, we attempted to reduce the variance by employing a separate critic for each type of agitation pattern, as suggested by Mao et al. [2019], but to no avail. With limited variations in the training set, the learning algorithm focuses on strategies for large-scale patterns instead of the more intricate strategies for local patterns. Nonetheless, our training set could still lead to a policy network that adapts to many unseen agitation patterns.
9 CONCLUSION AND FUTURE WORK
We introduced a novel real-time flow reconstruction framework that provides a high degree of realism by leveraging Lagrangian flow measurement data. The framework employed reinforcement learning with training data generated using SPH. We showed the effectiveness of our framework by quantitatively and visually validating it against both real-world and synthetic ground truth. Our framework can reconstruct a variety of fluid flow features even when the details of external disturbance are unknown.
Although our reinforcement learning algorithm uses a homogeneous agent that shares the same network parameters for each buoy, it is worth studying whether the learning efficiency can be improved by employing multi-agent learning algorithms such as MADDPG [Lowe et al. 2017]. Integrating those algorithms into our framework is non-trivial due to the variability in the number of buoys. However, if such integration is realized, it is likely that the agents can cooperate in a sophisticated manner to reconstruct local flow details.
Our reconstruction framework can be used as a real-time flow visualization tool for educational or scientific purposes. As many flow measurement systems suffer from optical occlusion problems, it is advantageous to use our framework to visualize agitated flows within opaque enclosures. The framework can also be used as a novel interaction interface for games that involve liquid flow. For instance, users can manipulate underwater or floating game objects by pushing or swirling actual water. Users can enjoy the game’s immediate visual feedback while experiencing haptic feedback from tangible water.
APPENDICES
A DATA GENERATION
Buoy properties in Table 4 were used to accurately generate training and validation datasets through SPH simulations.
The second dimension corresponds to the axial direction of the buoy cylinder. The origin is defined as the geometric center of the buoy.
Table 4. Mass and Geometric Properties of Each Buoy
The second dimension corresponds to the axial direction of the buoy cylinder. The origin is defined as the geometric center of the buoy.
B DETAILS IN REINFORCEMENT LEARNING
The hyper-parameters that differ from the original TD3 implementation are summarized in Table 5. We also employ a data enrichment technique that exploits the symmetry of the cubical container for faster training. Each entry from the replay buffer is cloned and transformed eight times with eight unique rotation and mirroring patterns before actual learning.
Table 5. Hyper-parameters used in Reinforcement Learning
Supplemental Material
- . 2012. Versatile rigid-fluid coupling for incompressible SPH. ACM Transactions on Graphics 31, 4 (
July 2012), Article62 , 8 pages.DOI: Google Scholar
Digital Library
- . 2018. MLS pressure boundaries for divergence-free and viscous SPH fluids. Computers & Graphics 76 (2018), 37–46.
DOI: Google ScholarCross Ref
- . 2017. Divergence-free SPH for incompressible and viscous fluids. IEEE Transactions on Visualization and Computer Graphics 23, 3 (
March 2017), 1193–1206.DOI: Google ScholarDigital Library - . 2019. Volume maps: An implicit boundary representation for SPH. In ACM Conference on Motion, Interaction, and Games (MIG ’19).
DOI: Google ScholarDigital Library - . 2012. Incompressible SPH for free surface flows. Computers & Fluids 67 (2012), 138–151.
DOI: Google ScholarCross Ref - . 2015. ShapeNet: An Information-Rich 3D Model Repository.
Technical Report . Stanford University–Princeton University–Toyota Technological Institute at Chicago. https://arxiv.org/abs/1512.03012Google Scholar - . 2018. Optimization of the Parameters of Smoothed Particle Hydrodynamics Method, Using Evolutionary Algorithms. Springer International Publishing, Cham, 153–167.
DOI: Google ScholarCross Ref - . 2020. Functional optimization of fluidic devices with differentiable stokes flow. ACM Transactions on Graphics 39, 6 (
Nov. 2020), Article197 , 15 pages.DOI: Google ScholarDigital Library - . 2018. Coupled fluid density and motion from single views. Computer Graphics Forum, Vol. 37. Blackwell Publishing Ltd, 47–58.
arxiv:1806.06613 . https://arxiv.org/abs/1806.06613v1Google Scholar - . 2019. ScalarFlow: A large-scale volumetric data set of real-world scalar transport flows for computer animation and machine learning. ACM Transactions on Graphics 38, 6 (
Nov. 2019), Article239 , 16 pages.DOI: Google ScholarDigital Library - . 2006. Tomographic particle image velocimetry. Experiments in Fluids 41, 6 (2006), 933–947.Google ScholarCross Ref
- . 2018. Addressing function approximation error in actor-critic methods. In Proceedings of the International Conference on Machine Learning. 1582–1591.Google Scholar
- . 1991. Data assimilation in meteorology and oceanography.
Advances in Geophysics , Vol. 33. Elsevier, 141–266.DOI: Google ScholarCross Ref - . 2019. Interlinked SPH pressure solvers for strong fluid-rigid coupling. ACM Transactions on Graphics 38, 1 (
Jan. 2019),Article 5 , 13 pages.DOI: Google ScholarDigital Library - . 2014. From capture to simulation: Connecting forward and inverse problems in fluids. ACM Transactions on Graphics 33, 4 (
July 2014), Article139 , 11 pages.DOI: Google ScholarDigital Library - . 2017. Cooperative multi-agent control using deep reinforcement learning. Autonomous Agents and Multiagent Systems. Springer International Publishing, 66–83.Google ScholarCross Ref
- . 2015. A review of measurement-integrated simulation of complex real flows. Journal of Flow Control, Measurement & Visualization 3, 2 (
March 2015), 51–66.DOI: Google ScholarCross Ref - . 2020. DiffTaichi: Differentiable programming for physical simulation. In Proceedings of the International Conference on Learning Representations.Google Scholar
- . 2022. Reconstruction of dexterous 3D motion data from a flexible magnetic sensor with deep learning and structure-aware filtering. IEEE Transactions on Visualization and Computer Graphics 28, 6 (2022), 2400–2414.
DOI: Google ScholarDigital Library - . 2014a. Implicit incompressible SPH. IEEE Transactions on Visualization and Computer Graphics 20, 3 (
March 2014), 426–435.DOI: Google ScholarDigital Library - . 2014b. SPH fluids in computer graphics. Eurographics 2014 - State of the Art Reports, and (Eds.). The Eurographics Association.
DOI: Google ScholarCross Ref - . 2002. Atmospheric Modeling, Data Assimilation and Predictability. Cambridge University Press.
DOI: Google ScholarCross Ref - . 2019. Smoothed particle hydrodynamics techniques for the physics based simulation of fluids and solids. Eurographics Tutorial, 2019. Google ScholarCross Ref
- . 2021. Fast corotated elastic SPH solids with implicit zero-energy mode control. Proceedings of the ACM on Computer Graphics and Interactive Techniques 4, 3 (
Sept. 2021), Article33 , 21 pages.DOI: Google ScholarDigital Library - . 2009. Mixing fluids and granular materials. Computer Graphics Forum 28, 2 (2009), 213–218.
DOI: .arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1467-8659.2009.01360.x Google ScholarCross Ref - . 1989. Particle Image Velocimetry. Springer, Berlin, 127–199.
DOI: Google ScholarCross Ref - . 2017. Multi-agent actor-critic for mixed cooperative-competitive environments. Neural Information Processing Systems (NeurIPS) (2017).
DOI: Google ScholarCross Ref - . 1993. Particle tracking velocimetry in three-dimensional flows. Experiments in Fluids 15, 2 (
July 1993), 133–146.DOI: Google ScholarCross Ref - . 2013. Position based fluids. ACM Transactions on Graphics 32, 4 (
July 2013), Article104 , 12 pages.DOI: Google ScholarDigital Library - . 2019. Variance reduction for reinforcement learning in input-driven environments. In Proceedings of the International Conference on Learning Representations.
DOI: Google ScholarCross Ref - . 2004. Fluid control using the adjoint method. ACM Transactions on Graphics 23, 3 (
Aug. 2004), 449–456.DOI: Google ScholarDigital Library - . 1992. Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics 30, 1 (1992), 543–574.
DOI: .arXiv:https://doi.org/10.1146/annurev.aa.30.090192.002551 Google ScholarCross Ref - . 2013. VDB: High-resolution sparse volumes with dynamic topology. ACM Transactions on Graphics 32, 3 (
July 2013), Article27 , 22 pages.DOI: Google ScholarDigital Library - . 2019. Mixing sauces: A viscosity blending model for shear thinning fluids. ACM Transactions on Graphics 38, 4 (
July 2019), Article95 , 17 pages.DOI: Google ScholarDigital Library - . 2017. Efficient solver for spacetime control of smoke. ACM Transactions on Graphics 36, 5 (
July 2017), Article162 , 13 pages.DOI: Google ScholarDigital Library - . 2012. Smoothed particle hydrodynamics and magnetohydrodynamics. Journal of Computational Physics 231, 3 (2012), 759–794.
DOI: .Special Issue: Computational Plasma Physics. Google ScholarDigital Library - . 2018. Ensuring numerical stability of wave propagation by tuning model parameters using genetic algorithms and response surface methods. Environmental Modelling and Software 103, C (
May 2018), 62–73.DOI: Google ScholarDigital Library - . 2012. Ghost SPH for animating water. ACM Transactions on Graphics 31, 4 (
July 2012),Article 61 , 8 pages.DOI: Google ScholarDigital Library - . 2018. SPNets: Differentiable fluid dynamics for deep neural networks. In Proceedings of the 2nd Conference on Robot Learning (Proceedings of Machine Learning Research), , , , and (Eds.), Vol. 87. PMLR, 317–335. https://proceedings.mlr.press/v87/schenck18a.htmlGoogle Scholar
- . 2013. Modeling unsteady flow characteristics using smoothed particle hydrodynamics. Applied Mathematical Modelling 37, 3 (2013), 1431–1450.
DOI: Google ScholarCross Ref - . 2005. Taming liquids for rapidly changing targets. In Proceedings of the 2005 ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA ’05). ACM, New York, NY, 229–236.
DOI: Google ScholarDigital Library - . 2009. Predictive-corrective incompressible SPH. ACM Transactions on Graphics 28, 3 (
July 2009), Article40 , 6 pages.DOI: Google ScholarDigital Library - . 1995. Depicting fire and other gaseous phenomena using diffusion processes. In Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH ’95). ACM, New York, NY, 129–136.
DOI: Google ScholarDigital Library - . 2019. Video-guided real-to-virtual parameter transfer for viscous fluids. ACM Transactions on Graphics 38, 6 (
Nov. 2019), Article237 , 12 pages.DOI: Google ScholarDigital Library - . 1994. Numerical simulation of viscous flow by smoothed particle hydrodynamics. Progress of Theoretical Physics 92, 5 (Nov. 1994), 939–960.
DOI: Google ScholarCross Ref - . 2020. Dynamic fluid surface reconstruction using deep neural network. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR ’20).Google ScholarCross Ref
- . 2014. PIVlab—towards user-friendly, affordable and accurate digital particle image velocimetry in MATLAB. Journal of Open Research Software 2 (Oct. 2014).
DOI: Google ScholarCross Ref - . 2006. Detail-preserving fluid control. In Proceedings of the 2006 ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA ’06). Eurographics Association, 7–12. Google ScholarDigital Library
- . 2003. Keyframe control of smoke simulations. ACM Transactions on Graphics 22, 3 (
July 2003), 716–723.DOI: Google ScholarDigital Library - . 2009. Physically guided liquid surface modeling from videos. In ACM SIGGRAPH 2009 Papers (SIGGRAPH ’09). Association for Computing Machinery, New York, NY, Article
90 , 11 pages.DOI: Google ScholarDigital Library - . 2018. A physically consistent implicit viscosity solver for SPH fluids. Computer Graphics Forum 37, 2 (2018), 145–155.
DOI: .arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.13349 Google ScholarCross Ref - . 2017. Rainbow particle imaging velocimetry for dense 3D fluid velocity imaging. ACM Transactions on Graphics 36, 4 (
July 2017), Article36 , 14 pages.DOI: Google ScholarDigital Library - . 2012. Angular domain reconstruction of dynamic 3D fluid surfaces. In Proceedings of the 2012 IEEE Conference on Computer Vision and Pattern Recognition. 310–317.
DOI: Google ScholarCross Ref
Index Terms
- Real-Time Reconstruction of Fluid Flow under Unknown Disturbance
Recommendations
Divergence-Free SPH Fluid Simulation Using Density Constraint Condition
In this paper, a novel, incompressible fluid simulation framework based on the divergence-free Smoothed Particle Hydrodynamics model is presented. The novel SPH model combines a system of non-linear density constraint conditions and the divergence-free ...
Synthesizing Solid-Induced Turbulence for Particle-Based Fluids
CADGRAPHICS '13: Proceedings of the 2013 International Conference on Computer-Aided Design and Computer GraphicsSimulating the accompanying turbulent details of fluid-solid coupling is still challenging, as numerical dissipation always plagues current fluid solvers. In this paper, we propose a novel particle-based method to simulate turbulent details generated ...
Realistic and stable simulation of turbulent details behind objects in smoothed-particle hydrodynamics fluids
This paper presents a novel realistic and stable turbulence synthesis method to simulate the turbulent details generated behind objects in smoothed particle hydrodynamics SPH fluids. Firstly, by approximating the boundary layer theory on the fly in SPH ...
Comments